錯誤的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平均值？

回答

相關問題