2017-01-12 2 views
-1

私は非ベクトル化と半ベクトル化の方法でstdを計算しようとしています。 ベクトル化されていないバージョンのコードは正常に機能しますが、半ベクトル化されたバージョンも同様に機能しますが、生成される結果は同じではありません。
これはバージョン1:このスニペットで同じ標準偏差が生成されないのはなぜですか?

import math 
#unvectorized version --really slow! 
def calc_std_classic(a): 
    batch = a.shape[0] 
    channel = a.shape[1] 
    width = a.shape[2] 
    height = a.shape[3] 
    mean = calc_mean_classic2(data_train) 

    sum = np.zeros((channel)) 
    for i in range(batch): 
     for j in range(channel): 
      for w in range(width): 
       for h in range(height): 
        sum[j] += (abs(a[i,j,w,h] - mean[j])**2) 

    var = (sum/(width*height*batch)) 
    return [(math.sqrt(x)) for x in var ] 

半ベクトル:

def calc_std_classic2(a): 
    batch = a.shape[0] 
    channel = a.shape[1] 
    width = a.shape[2] 
    height = a.shape[3] 
    mean = calc_mean_classic2(data_train) 

    sum = np.zeros((channel)) 
    for i in range(batch): 
     for j in range(channel): 
      sum[j] += np.sum(abs(a[i,j,:,:] - mean[j])**2) 

    var = (sum/(width*height*batch)) 
    return [(math.sqrt(x)) for x in var ] 

、これはその必要に応じて平均算出する方法である:

def calc_mean_classic2(a): 
    #sum all elements in each channel and divide by the number of elements 
    batch = a.shape[0] 
    channel = a.shape[1] 
    width = a.shape[2] 
    height = a.shape[3] 

    sum = np.zeros((channel)) 
    for i in range(batch): 
     for j in range(channel): 
      sum[j] += np.sum(a[i,j,:,:]) 

    return (sum/(width*height*batch)) 

ニシキヘビを使用して生成された出力numpy.std()であり、2つの方法は次のとおりです。

std = np.std(data_train, axis=(0,2,3)) 
std2 = calc_std_classic(data_train) 
std3 = calc_std_classic2(data_train) 

が発生:

std = [ 62.99321928 62.08870764 66.70489964] 
std2 = [62.99321927774396, 62.08870764038716, 66.7048996406483] 
std3 = [62.99321927813685, 62.088707640014405, 66.70489964063101] 

あなたが見ることができるように、すべての3つは、8桁まで同じ結果を生成します。 3番目の方法では残りの数字が異なります。

私はここで間違っていますか?

+0

おそらくnumpyはより安定した合計アルゴリズムを使用します。 –

+0

@WillemVanOnsem:それはすべてですか? 私が行っている手順は正しいですか? – Breeze

+0

投票が遅れましたか?なぜ誰かがこの質問をd​​ownvoteだろうか?どうしたの ? – Breeze

答えて

1

浮動小数点演算エラー伝播には、多くの優れたリソースがあります。しかし、直ちに問題になるのは、numpy.ndarrayの表示が、python floatとは異なる精度で浮動することです。あなたの明示的な場合には

>>> print(np.std(arr, ....)) 
[ 0.28921072 0.2898092 0.28961785] 
>>> print(np.std(arr, ....).tolist()) 
[0.28921072085015914, 0.28980920065339233, 0.28961784922639483] 

:だから、あなたは、同一のデータ構造(例えばlist秒)に変換する必要がありますあなたの結果を比較すること

1はナイーブ使用しているためcalc_std_classiccalc_std_classic2間の差があります合計はa1+a2+....+an、他方はnp.sumです。 np.sumはナイーブな集計になる可能性がありますが、私が知る限りpairwise summationを使用しています。より高い精度が必要な場合は、Kahan summationを実装するか、python-builtin statistics._sumを使用します。

numpyがどのアルゴリズムを使用しているかわからないため、np.stdとあなたの亜種の違いは説明が難しいです。全体がarticle about "Algorithms for calculating variance" on wikipediaです。特に、item - meanの減算のために、単純なインプリメンテーションでアンダーフロー/オーバーフローの問題が発生する可能性があることに注意してください。


一般的なアドバイス:あなたは、その後statisticsを使用し、最高の精度でそれをしたい場合は、それは速く、その後、numpyを使用したい場合は

。 NumPyはほとんどパフォーマンス面に焦点を当てているので、が最も正確なアルゴリズムを実装していない可能性があります。おそらく正確でも高速でもないので、アルゴリズムの研究をしないで純粋な実装を避けてください。

+0

ありがとう、本当にありがとうございました:) – Breeze