mean and nanmean give different results on masked arrays

[clintrowe]

I'm working on code to compute some statistics and created a masked array. Even though the current dataset contains no NaNs, future data may. I found that, after masking the array, I obtained differing results with mean than with nanmean. Here's a snippet using a 4x1096 array:

print np.count_nonzero(np.isnan(stats),axis=1)
print np.mean(stats,1)
print np.nanmean(stats,1)
stats = np.ma.masked_invalid(stats)
print stats.size - np.count_nonzero(stats,axis=1)
print np.mean(stats,1)
print np.nanmean(stats,1)

with output:

0
[ 27.42351151 26.58027077 26.62998009 26.5291214 ]
[ 27.42351151 26.58027077 26.62998009 26.5291214 ]
0
[27.423510920392335 26.58027058622263 26.62998082515967 26.52912223426095]
[27.423511505126953 26.580270767211914 26.629980087280273 26.52912139892578]

Python 2.7.13 |Anaconda custom (64-bit)| (default, Dec 20 2016, 23:09:15)
numpy.version.version
'1.11.3'

[ahaldane]

I can reproduce this, and it is a precision issue. Your stats variable is likely of f4 type, eg I reproduce your result with

>>> from numpy.random import rand
>>> stats = rand(4,1096).astype('f4')

It doesn't happen if you use 'f8'.

Here's what is going on: Normally, mean and nanmean do not upcast floating point types to f8 (read their docstrings) meaning they compute the mean of your data with limited f4 precision. That's what is happening in your 1st, 2nd, and 4th computations. (the 4th line just has more digits printed for some reason, which appears to be fixed in master).

However, MaskedArrays happen to have their own totally separate implementation of mean (but non nanmean), which evidently does upcast to f8 and therefore does a more accurate computation in your 3rd line. (So in fact your third line is the most "correct" value).

Probably we want to stop MaskedArray.mean from upcasting. I think I have an updated implementation of MaskedArray.mean lying around somewhere which behaves more like np.mean, probably I should try submitting it again.

For now if you want to avoid these precision issues you can cast all your data to f8.

[clintrowe]

Thanks for the great explanation. Indeed, my stats array was 'f4' and, as you said, the differences disappear if I use 'f8'. Either set of results is more than precise enough for my needs, it was only the differences that arose with the two implementation following masking that I found a little disconcerting.