>>> 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
在数学上,这两种方法是等价的。在数值上,它们是不同的。
您的方法涉及单个部门,并且碰巧:
1. + 2. + 4. + 5.
) 可以精确地表示为 a 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'>