为什么scipy.stats.nanmean会给出numpy.nansum不同的结果？

Question

>>> 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

我知道浮点表示的局限性.只是好奇为什么更笨拙的表达似乎给出"更好"的结果.

Answer 1

首先,这是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

从数学上讲,这两种方法是等价的.在数字上,它们是不同的.

你的方法涉及一个单独的划分,它发生了:

• numerator(1. + 2. + 4. + 5.)可以完全表示为a float; 和
• 分母(4.)是2的幂.

这意味着划分的结果是准确的3..

stats.nanmean()首先计算平均值[1., 2., 0., 4., 5.],然后调整它来计算NaNs.碰巧的是,这个mean(2.4)不能完全表示为a float,所以从这一点来说计算是不精确的.

我没有多想过,但是有可能构建一个角色可以反转的例子,并且stats.nanmean()可以提供比其他方法更准确的结果.

让我感到惊讶的是,stats.nanmean()不仅仅是这样做:

In [6]: np.mean(np.ma.MaskedArray(a, np.isnan(a)))
Out[6]: 3.0

在我看来,这对于它目前所做的事情来说是一种更好的方法.