2009-03-13 11 views
69

数字の集合の自己相関を行う必要があります。私はそれがセットと自分自身との相関関係であることを理解しています。numpy.correlateを使用して自己相関を行うにはどうすればよいですか?

私はnumpyの相関関数を使って試してみましたが、ほとんどの場合、最初の数字がでなく、であるべきであるので、結果は信じられません。

ので、この質問は、実際には2つの質問されていますnumpy.correlateは正確に何を

  1. をしているのですか?
  2. 自己相関を実行するにはどうすればよいですか?
+0

も参照してください:以下をさhttp://stackoverflow.com/questions/12269834/is-there-any-numpy-autocorrellation-function-with-standardized-output正規化自己相関の詳細については。 – amcnabb

答えて

79

numpy.correlate(a, v, mode)vの逆でaの畳み込みを実行し、指定されたモードによって切り取られた結果を与えている、あなたの最初の質問に答えるために。 definition of convolution、C(T)=Σ-∞<私<∞ V トン+ I-∞<トン<∞は、-∞から∞に結果することができますが、あなたは明らかにできません無限に長い配列を格納します。だから、クリッピングする必要がある、とモードの出番つまり3種類のモードがあります。フル、同じ、&有効:

  • 「フル」モードでは、avの両方を持っているすべてのtの結果を返します。いくつかの重複。
  • "同じ"モードは、最短ベクトル(aまたはv)と同じ長さの結果を返します。
  • "有効"モードは、avが完全に重複している場合にのみ結果を返します。 numpy.convolvedocumentationは、モードの詳細を示します。あなたの2番目の質問については

、私はnumpy.correlateそれだけでだけでなく、あなたにもう少しを与えている、あなたの自己相関を与えることだと思います。自己相関は、信号または関数が、ある時間差でどれほど類似しているかを見つけるために使用されます。時差が0の場合、自己相関は信号がそれ自身と同一であるため最高でなければならないため、自己相関結果配列の最初の要素が最大になると予想しました。しかし、相関は0の時間差で開始していません。負の時間差で開始し、0になり、次に正になります。それはあなたが期待していた、次のとおりです。

自己相関(A)=Σ-∞< I∞ < V トン+ I = T <∞

しかし、あなたが得ましたた:

自己相関(A)=Σ-∞<私<∞ V t + iここで、-∞< t <∞

相関結果の最後の半分を取る必要があり、それはあなたが探している自己相関でなければなりません。それを行うための単純なPythonの関数は次のようになります。

def autocorr(x): 
    result = numpy.correlate(x, x, mode='full') 
    return result[result.size/2:] 

あなたは、もちろん、xが実際に1次元配列であることを確認するために、エラーチェックが必要になります。また、この説明はおそらく数学的に厳密なものではありません。私は畳み込みの定義がそれらを使用するので無限大を投げていますが、必ずしも自己相関には適用されません。したがって、この説明の理論的な部分はやや怪しいかもしれませんが、実用的な結果が役立つことを願っています。 Thesepagesは、自己表記法や重い概念を気にする必要がなければ、非常に役立ちます。

+4

現在のnumpyのビルドでは、モード 'same'を指定して、A. Levyが提案したものを正確に達成することができます。関数の本体は 'return numpy.correlate(x、x、mode = 'same')' –

+10

@DavidZwickerを読み取ることができますが、結果は異なります。 'np.correlate(x、x、mode = 'full')[len(x)// 2:]!= np.correlate(x、x、mode = 'same')'。たとえば、 'x = [1,2,3,1,2]; np.correlate(x、x、mode = 'full'); '{' >>>配列([2,5,11,13,19,13,11,5,2]) '}' np.correlate( x、x、mode = 'same'); '{' >>> array([11、13、19、13、11])}正しいものは 'np.correlate(x、x、mode = 'full')[len(x)-1:];' {'>>> array([19,13,11,5,2])です。 '}参照**最初のアイテム**は**最大**です**。 – Developer

+7

この回答は非正規化自己相関を与えることに注意してください。 – amcnabb

12

自動相関には統計と畳み込みという2つのバージョンがあります。少し詳細を除いて、どちらも同じです:前者は区間[-1,1]上にあるように正規化されています。ここでは、統計的なものを行う方法の例です:

def acf(x, length=20): 
    return numpy.array([1]+[numpy.corrcoef(x[:-i], x[i:]) \ 
     for i in range(1, length)]) 
+5

'' corrcoef''の戻り値が2x2の行列であるため、 '' numpy.corrcoef [x:-i]、x [i:])[0,1] ''を第2行に入れます。 – luispedro

+0

どういう違いがありますか統計的畳み込み自己相関と畳み込み自己相関の間にあるか? –

9

私はちょうど同じ問題に遭遇したとして、私はあなたと数行のコードを共有したいと思います。実際には、現在のところ、stackoverflowの自己相関に関するいくつかのかなり類似した投稿があります。あなたはa(x, L) = sum(k=0,N-L-1)((xk-xbar)*(x(k+L)-xbar))/sum(k=0,N-1)((xk-xbar)**2)として自己相関を定義した場合、[これはIDLのa_correlate機能に与えられた定義であり、それは私が質問#12269834の解答2に見るものと一致]、次は正しい結果を与えるようだ:

import numpy as np 
import matplotlib.pyplot as plt 

# generate some data 
x = np.arange(0.,6.12,0.01) 
y = np.sin(x) 
# y = np.random.uniform(size=300) 
yunbiased = y-np.mean(y) 
ynorm = np.sum(yunbiased**2) 
acor = np.correlate(yunbiased, yunbiased, "same")/ynorm 
# use only second half 
acor = acor[len(acor)/2:] 

plt.plot(acor) 
plt.show() 

ように私はこれをsin曲線と均一なランダム分布でテストしていますが、どちらの結果も期待通りです。他のものと同じようにmode="full"の代わりにmode="same"を使用しました。私はOPの質問に本当の答えは簡潔Numpy.correlateのドキュメントからこの抜粋に含まれていると思います

def autocorr(x, t=1): 
    numpy.corrcoef(numpy.array([x[0:len(x)-t], x[t:len(x)]])) 
13

mode : {'valid', 'same', 'full'}, optional 
    Refer to the `convolve` docstring. Note that the default 
    is `valid`, unlike `convolve`, which uses `full`. 

その2つの入力引数(IE用と同じベクトルを与えられたとき、これは、ノー「モード」の定義で使用する場合、Numpy.correlate関数は、スカラーを返します、ということを意味 - とき自己相関を実行するために使用される)。

+0

「相関係数」は、信号処理で使用される自己相関ではなく、統計に使用された自己相関ではありませんか? https://en.wikipedia.org/wiki/Autocorrelation#Signal_processing –

+0

@DanielPendergast私は信号処理に精通していません。 numpy docsから: "Pearsonのproduct-moment相関係数を返します。"それは信号処理版ですか? –

0

:トンの遅れのための統計的相関を計算する代わりにnumpy.correlatenumpy.corrcoef機能を使用して

1

私はタリブを使用します。このような自己相関のためのCORRELは、私はあなたが他のパッケージと同じことを行うことができ疑う:

def autocorrelate(x, period): 

    # x is a deep indicator array 
    # period of sample and slices of comparison 

    # oldest data (period of input array) may be nan; remove it 
    x = x[-np.count_nonzero(~np.isnan(x)):] 
    # subtract mean to normalize indicator 
    x -= np.mean(x) 
    # isolate the recent sample to be autocorrelated 
    sample = x[-period:] 
    # create slices of indicator data 
    correls = [] 
    for n in range((len(x)-1), period, -1): 
     alpha = period + n 
     slices = (x[-alpha:])[:period] 
     # compare each slice to the recent sample 
     correls.append(ta.CORREL(slices, sample, period)[-1]) 
    # fill in zeros for sample overlap period of recent correlations  
    for n in range(period,0,-1): 
     correls.append(0) 
    # oldest data (autocorrelation period) will be nan; remove it 
    correls = np.array(correls[-np.count_nonzero(~np.isnan(correls)):])  

    return correls 

# CORRELATION OF BEST FIT 
# the highest value correlation  
max_value = np.max(correls) 
# index of the best correlation 
max_index = np.argmax(correls) 
5

あなたの質問1は、すでに広く、ここでいくつかの優れた答えで議論されてきました。

私は自己相関の数学的性質に基づいて信号の自己相関を計算することを可能にする数行のコードをあなたと共有することを考えました。すなわち、自己相関は次のように計算することができるされている。

  1. はフーリエ公平信号

  2. 計算の変換を計算公平信号

  3. の信号から平均値を減算し得不偏信号のフーリエ変換の各値の2乗ノルムをとることにより、信号のパワースペクトル密度を計算する。

  4. 逆フーリエ変換パワースペクトル密度

  5. はフーリエ公平信号の二乗和によりパワースペクトル密度の変換、そして得られたベクター

これを行うためのコードの半分だけを取る逆正規化

def autocorrelation (x) : 
    """ 
    Compute the autocorrelation of the signal, based on the properties of the 
    power spectral density of the signal. 
    """ 
    xp = x-np.mean(x) 
    f = np.fft.fft(xp) 
    p = np.array([np.real(v)**2+np.imag(v)**2 for v in f]) 
    pi = np.fft.ifft(p) 
    return np.real(pi)[:x.size/2]/np.sum(xp**2) 
+0

これには何か問題がありますか?他の自己相関関数と結果を一致させることはできません。関数は類似しているように見えますが、やや渋滞しているようです。 – pindakaas

+0

@pindakaasあなたはもっと具体的になりますか?あなたが見つけた矛盾と、どの機能に関する情報を提供してください。 – Ruggero

+0

'p = np.abs(f)'を使わないのはなぜですか? – dylnan

関連する問題