numpy平均値が間違っていますか？

Question

私は通常、巨大なシミュレーションを行います。場合によっては、粒子の集合の中心を計算する必要があります。私は多くの状況で、numpy.mean（）が返す平均値が間違っていることに注意しました。私はそれがアキュムレータの飽和に起因していることを理解することができます。この問題を避けるために、粒子の小さな集合の中のすべての粒子に対して総和を分割できますが、それは不快です。誰もがこの問題を優雅なやり方で解決する方法を持っていますか？あなたは最大値と最小値をチェックすると、あなたが得る

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

：

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

はちょうどあなたの好奇心をpikingために、次の例では、私は私のシミュレーションで観察するものと同様のものを生み出します

ただし、平均値は

a.mean() 
30687.236328125 

です。何かが間違っていることがわかりますここに。これは、dtype = np.float64を使用しているときには起こらないので、問題を単精度で解決するのがいいはずです。

Answer 1

これは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）

問題は、浮動小数点精度が限られているということです。アキュムレータ値が追加される要素に対して相対的に増加すると、相対精度が低下します。

1つの解決策は、加算器ツリーを構築することによって、相対的な成長を制限することです。ここでの例では、（私のPythonが...十分ではありません）Cにあります：

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）

Answer 2

あなたは（アキュムレータのタイプを指定dtypeキーワード引数、とnp.meanを呼び出すことができます浮動小数点型配列の配列と同じ型にデフォルト設定されています）。

だから、a.mean(dtype=np.float64)を呼び出すと、あなたのおもちゃの例と、おそらく大きな配列の問題が解決されます。

Answer 3

迅速かつ汚い答え

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 

もっと賢明な回答をするには、データ構造、サイズ、ターゲットアーキテクチャに関する情報が必要です。

Answer 4

あなたは、部分的に使用することによってこの問題を解決することができます内蔵の部分和をダウン追跡math.fsum、（ドキュメントは、ASレシピのプロトタイプへのリンクを含んで）：私の知る限り承知しているよう

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

を、numpyにはアナログがありません。