Numpy float 均值计算精度

Question

我碰巧有一个 numpy 数组：

a.dtype, a.shape
#(dtype('float64'), (32769,))

值为：

a[0]
#3.699822718929953
all(a == a[0])
True

但是：

a.mean()
3.6998227189299517

平均值偏离第 15 位和第 16 位数字。

任何人都可以说明这种差异是如何累积超过 30K 均值的，是否有办法避免它？

以防万一，我的 OS 是 64 位。

Answer 1

这是因为 float64 类型无法以正确的精度存储浮点数的总和。为了解决这个问题，您需要使用更大的数据类型 course^*。 Numpy 有一个 longdouble dtype，你可以在这种情况下使用它：

In [23]: np.mean(a, dtype=np.longdouble)                                                                                                                                                                    
Out[23]: 3.6998227189299530693

另外，注意：

In [25]: print(np.longdouble.__doc__)                                                                                                                                                                       
Extended-precision floating-point number type, compatible with C
    ``long double`` but not necessarily with IEEE 754 quadruple-precision.
    Character code: ``'g'``.
    Canonical name: ``np.longdouble``.
    Alias: ``np.longfloat``.
    Alias *on this platform*: ``np.float128``: 128-bit extended-precision floating-point number type.

_{*阅读评论了解更多详情。}

Answer 2

平均值（根据定义）：

a.sum()/a.size

不幸的是，将所有这些值相加并相除会累积浮点错误。它们的大小通常约为：

np.finfo(np.float).eps
Out[]: 2.220446049250313e-16

是的，e-16，关于你从哪里得到它们。您可以通过使用更高准确度的 float 来减小误差，例如 float128（如果您的系统支持），但只要您对大量 float 求和，它们就会不断累积一起。如果您真的想要身份，则必须对其进行硬编码：

def mean_(arr):
    if np.all(arr == arr[0]):
        return arr[0]
    else:
        return arr.mean()

实际上，您永远不想在浮点数之间使用 ==。通常在 numpy 中我们使用 np.isclose 或 np.allclose 来比较浮点数正是出于这个原因。有方法使用其他包和利用神秘的机器级方法计算数字来获得（接近）完全相等，但很少值得性能和清晰度受到影响。

Answer 3

这是最大误差范围的粗略近似值。这不代表平均误差，可以通过更多分析加以改进。

考虑使用四舍五入到偶数的浮点算法计算和：

sum = 0;
for (i = 0; i < n; ++n)
    sum += a[i];

其中每个 a[i] 在 [0, m] 中。

令ULP(x)表示浮点数x中的最小精度单位。（例如，在具有 53 位尾数的 IEEE-754 二进制 64 格式中，如果不大于 |x| 的 2 的最大幂为 2^p，然后ULP(x) = 2^p−52。使用最近舍入法，结果为 x 的任何运算中的最大误差为 ½ULP(x).

如果忽略舍入误差，i次迭代后sum的最大值为i•m。因此，迭代 i 中加法的误差界限为 ½ULP(i•m) . （实际上 i=1 为零，因为这种情况加起来为零，没有错误，但我们忽略了这个近似值。）那么所有加法的边界总和是i 从 1 到 n 的 ½ULP(i•m) 的总和。这大约是 ½•n•(n+1)/2•ULP(m) = ¼•n•(n+1)•ULP(m)。（这是一个近似值，因为它将 i 移动到 ULP 函数之外，但 ULP 是一个不连续的函数。它是“近似线性的”，但有跳跃。由于跳跃是由以下因素决定的二，近似值最多可以偏离两倍。）

因此，对于 32,769 个元素，我们可以说总舍入误差最多为 ¼•32,769•32,770•ULP(m)，大约为 2.7•10⁸乘以最大元素值的ULP。 ULP是2⁻⁵²乘以不小于m的2的最大幂，所以大约是2.7•10⁸•2⁻⁵² = 6•10⁻⁸次m.

当然，32,768 个总和（不是 32,769，因为第一个必然没有错误）在同一个方向上偶然出现的可能性微乎其微，但我猜想有人可能会设计出一系列接近该值的值.

一个实验

这是（蓝色）10,000 个求和数组样本的平均误差图表，大小为 100 到 32,800 x 100s，元素是从 [0, 1) 上的均匀分布中随机抽取的。通过比较使用 float (IEEE-754 binary32) 计算的总和与使用 double (IEEE-754 binary64) 计算的总和来计算误差。（样本都是 2⁻²⁴ 的倍数，并且 double 具有足够的精度，因此总和高达 2²⁹ 这样值是准确的。）

绿线是cn√n和c 设置为匹配蓝线的最后一个点。我们看到它长期跟踪蓝线。在平均和超过 2 的幂的点，平均误差会暂时增加得更快。在这些点上，总和进入了一个新的 binade，并且由于 ULP 的增加，进一步的添加具有更高的平均误差。在 binade 的过程中，这个固定的 ULP 相对于 n 减少，使蓝线回到绿线。

Numpy float 均值计算精度

Numpy float mean calculation precision

python

floating-point

precision

numpy

一个实验