returning a masked nan would make more sense to me, but should we change the behaviour instead of keeping the old?

Agreed that a masked nan might be more logical, but since np.ma.mean also just returns nan, it is probably best to keep the behaviour as is.

Not a huge fan of this special-casing of index building

Can we just do something like:

indexer = [slice(None)] * asorted.ndim
if asorted.shape[axis] == 0:
    indexer[axis] = []
elif asorted.shape[axis] % 2 == 1:
    indexer[axis] = [h]
else:
    indexer[axis] = [h-1, h]

# ....

return np.ma.mean(asorted[indexer], axis=axis)

Edit: no, that doesn't work since the counts are different on each dimension

I think I meant:

indexer = np.ix_(*(np.arange(x) for x in asorted.shape))
if asorted.shape[axis] == 0:
    indexer[axis] = # empty of the right shape
elif:
    indexer[axis] = np.stack([h-1, h], axis=axis)

What does this return for a fully-masked matrix? nan or masked?

hm right you can index with mixed slices and indices, might be nicer

it returns unmasked nan array of the input shape % the reduced axis

hm right you can index with mixed slices and indices, might be nicer

Doesn't work when indexing both h and h-1 and the same time. But I'd replace the empty slice with a empty array, for consistency. So remove all the slices, rather than all the arrays

not really nicer due to the way the mean works at the end, which is ugly but I don't really want to change it for a bugfix

Fair enough, I might do some maintenance on this another time

@juliantaylor: re index-building, #8708