>>> import numpy as np
>>> from scipy import stats
>>> a = np.r_[1., 2., np.nan, 4., 5.]
>>> stats.nanmean(a)
2.9999999999999996
>>> np.nansum(a)/np.sum(~np.isnan(a))
3.0
浮動小数点表現の制限を認識しています。不器用な表現が「より良い」結果をもたらすように見える理由に興味があります。
>>> import numpy as np
>>> from scipy import stats
>>> a = np.r_[1., 2., np.nan, 4., 5.]
>>> stats.nanmean(a)
2.9999999999999996
>>> np.nansum(a)/np.sum(~np.isnan(a))
3.0
浮動小数点表現の制限を認識しています。不器用な表現が「より良い」結果をもたらすように見える理由に興味があります。
まず第一に、scipy.nanmean()
これは私たちが何と比較しているのかを知るためです:
def nanmean(x, axis=0):
x, axis = _chk_asarray(x,axis)
x = x.copy()
Norig = x.shape[axis]
factor = 1.0-np.sum(np.isnan(x),axis)*1.0/Norig
x[np.isnan(x)] = 0
return np.mean(x,axis)/factor
数学的には、2つの方法は同等です。数値的には、それらは異なります。
あなたの方法は単一の部門を含みます、そしてそれはそう起こります:
1. + 2. + 4. + 5.
として表すことができます。float
と4.
)は2の累乗です。これは、除算の結果が正確であることを意味し3.
ます。
stats.nanmean()
最初にの平均を計算し、[1., 2., 0., 4., 5.]
次にそれを考慮して調整する必要がありますNaNs
。たまたま、この平均(2.4
)を正確にaとして表すことはできないfloat
ため、この時点からの計算は不正確になります。
あまり考えていませんが、役割が逆になりstats.nanmean()
、他の方法よりも正確な結果が得られる例を構築できるかもしれません。
私を驚かせるのは、それstats.nanmean()
が単に次のようなことをしないということです。
In [6]: np.mean(np.ma.MaskedArray(a, np.isnan(a)))
Out[6]: 3.0
これは、現在のアプローチよりも優れたアプローチのように思えます。
答えは次のコードにありstats.nanmean
ます:
x, axis = _chk_asarray(x,axis)
x = x.copy()
Norig = x.shape[axis]
factor = 1.0-np.sum(np.isnan(x),axis)*1.0/Norig
x[np.isnan(x)] = 0
return np.mean(x,axis)/factor
1.0 - np.sum
合計の減算と関係があると思います。
@eumiroが述べたように、stats.nanmeanは、単純なワンライナーの方法とは異なる迂回法で平均を計算しました。
同じ参照コードから、
np.sum(np.isnan(x),axis)
numpy.int32
*を掛けると1.0
、結果が整数で結果に差が生じた場合に得られるものとは対照的に、浮動小数点近似が得られます。
>>> numpy.int32(1)*1.0/5
0.20000000000000001
>>> int(numpy.int32(1))*1.0/5
0.2
>>> type(np.sum(np.isnan(x),axis))
<type 'numpy.int32'>