错误的numpy平均值？

我通常使用巨大的模拟。有时候，我需要计算一组粒子的质心。我注意到在很多情况下，numpy.mean（）返回的平均值是错误的。我可以看出，这是由于累加器饱和所致。为了避免这个问题，我可以将所有粒子中的所有粒子进行总和分解，但这是不舒服的。任何人都有和想法如何以优雅的方式解决这个问题？错误的numpy平均值？

只是为了piking了你的好奇心，下面的例子产生类似于我在模拟观察的东西：

import numpy as np 
a = np.ones((1024,1024), dtype=np.float32)*30504.00005

如果您检查最大值和最小值，您可以：

a.max() 
30504.0 
a.min() 
30504.0

然而，平均值为：

a.mean() 
30687.236328125

你可以弄清楚，什么是错的这里。使用dtype = np.float64时不会发生这种情况，所以应该很好地解决单精度问题。

来源

2013-07-04 Alejandro

如果这些答案中的任何一个解决了您的问题，您应该接受它。 – tacaswell

这不是一个NumPy问题，它是一个浮点问题。同样发生在C：

float acc = 0; 
for (int i = 0; i < 1024*1024; i++) { 
    acc += 30504.00005f; 
} 
acc /= (1024*1024); 
printf("%f\n", acc); // 30687.304688

（Live demo）

的问题是，浮点具有有限的精度;随着累加器值相对于添加到其中的元素增加，相对精度下降。

一个解决方案是通过构造一个加法器树来限制相对增长。这里的C中的例子（我的Python是不够好...）：

float sum(float *p, int n) { 
    if (n == 1) return *p; 
    for (int i = 0; i < n/2; i++) { 
     p[i] += p[i+n/2]; 
    } 
    return sum(p, n/2); 
} 

float x[1024*1024]; 
for (int i = 0; i < 1024*1024; i++) { 
    x[i] = 30504.00005f; 
} 

float acc = sum(x, 1024*1024); 

acc /= (1024*1024); 
printf("%f\n", acc); // 30504.000000

（Live demo）

来源

2013-07-04 06:24:23

谢谢奥利，我知道这不是numpy的问题。我认为有一个函数可以自己分割累加器以避免这个问题（在numpy中实现）应该很有趣。 – Alejandro

@Alejandro：查看更新后的答案。 –

谢谢奥利，我喜欢你的方法。这是非常有用的 – Alejandro

你可以调用np.mean与dtype关键字参数，指定累加器的类型（其默认与浮点数组的数组类型相同）。

所以调用a.mean(dtype=np.float64)将解决你的玩具的例子，也许你的问题与更大的数组。

来源

2013-07-04 06:29:20 Jaime

是的，它是在问题中说明的。正如你所说，np.float64解决了这个问题。但是，在不改变dtype的情况下手工计算平均值时可以解决问题。如果你采用少量的数据子集并计算部分求和，即使采用单精度，你也可以得到更好的结果 – Alejandro

正确的做法是使用（Welford的方法）[http://stackoverflow.com/questions/895929/how -do-i -definition-the-standard-deviation-stddev-of-a-set-of-values/897463＃897463]或类似的变体，但没有类似的东西在numpy中实现。让你的'np.float64'数组最好的事情是告诉'np.mean'使用'dtype'关键字使用'np.float64'累加器。 – Jaime

快速和肮脏的答案

assert a.ndim == 2 
a.mean(axis=-1).mean()

这给了预期的结果为1024 * 1024矩阵，当然，这不会是更大的阵列真的......

如果计算将平均不是你的代码中的瓶颈我会在python中实现自己的特别算法：但是细节取决于你的数据结构。

如果计算均值是一个瓶颈，那么一些专门的（并行）还原算法可以解决这个问题。

编辑

这种方法可能看起来很可笑，但将肯定缓解这个问题，是几乎一样有效.mean()本身。

In [65]: a = np.ones((1024,1024), dtype=np.float32)*30504.00005 

In [66]: a.mean() 
Out[66]: 30687.236328125 

In [67]: a.mean(axis=-1).mean() 
Out[67]: 30504.0 

In [68]: %timeit a.mean() 
1000 loops, best of 3: 894 us per loop 

In [69]: %timeit a.mean(axis=-1).mean() 
1000 loops, best of 3: 906 us per loop

给人一种更明智的办法需要对数据结构的一些更多的信息，它的大小和目标architeture。

来源

2013-07-04 06:51:12

可以部分地通过纠正这种内置math.fsum，跟踪下来的部分和（该文档包含一个链接到AS配方原型）：

>>> fsum(a.ravel())/(1024*1024) 
30504.0

据我所知，numpy没有模拟量。

来源

2013-07-04 08:14:44

+1表示精度，但在我的机器上比'a.mean（）'或'a.mean（axis = -1）.mean（）'慢100多倍。 –

确定它是纯Python。即使这种事情变得越来越模糊，与仅仅总结事情相比，仍然有相当多的工作要做。但问题当然是这样做是否会在你的真实代码中造成瓶颈 - 你在原文中提到'有时':-)。 –

'math.fsum'在C中实现，AS配方只是一个参考。 AS python的代码可能会慢几千倍......因为OP说的是“巨大”的问题，我认为虽然速度是一个问题，但在这里我是孤身一人。在交易速度和小内存占用的准确性方面没有任何错误...... –

错误的numpy平均值？

回答

相关问题