numpy.correlateを使用して自己相関を行うにはどうすればよいですか？

Question

数字のセットの自己相関を行う必要があります。これは、理解しているように、数字とセット自体の相関にすぎません。

numpyの相関関数を使用して試してみましたが、結果は信じられません。ほとんどの場合、最初の数が最大ではないベクトルを与えるため、 。

つまり、この質問は実際には2つの質問です。

1. numpy.correlate は正確に何をしていますか？
2. 自己相関を行うためにどのように（または他の何か）使用できますか？

Answer 1

最初の質問に答えるために、numpy.correlate(a, v, mode)はaの畳み込みをvの逆で実行し、指定されたモードでクリップされた結果を提供します。 畳み込みの定義、C（t）= <！>＃8721; _{-<！>＃8734; <！> lt; i <！> lt; <！>＃8734;} a _i v _{t + i} ここで、-<！>＃8734; <！> lt; t <！> lt; <！>＃8734;、-<！>＃8734;からの結果を許可します。 <！>＃8734;ですが、明らかに無限に長い配列を格納することはできません。そのため、クリップする必要があります。それがモードの出番です。3つの異なるモードがあります。完全、同じ、<！>アンプ、有効：

• <！> quot; full <！> quot;モードは、tとnumpy.convolveの両方に重複があるすべてのnumpy.correlateの結果を返します。
• <！> quot; same <！> quot; modeは、最短ベクトルと同じ長さの結果を返します（xまたは<=>）。
• <！> quot; valid <！> quot;モードは、<=>と<=>が完全に重なり合っている場合にのみ結果を返します。 <=>のドキュメントで詳細を確認できます。モード。

2番目の質問については、<=> が自己相関を与えていると思いますが、もう少しだけ与えています。自己相関は、特定の時間差で信号または関数がそれ自体にどれだけ似ているかを見つけるために使用されます。時間差が0の場合、信号はそれ自体と同一であるため、自己相関は最高になるはずです。したがって、自己相関結果配列の最初の要素が最大になると予想しました。ただし、相関は時間差0で開始されません。負の時間差で開始され、0に近くなり、その後正になります。つまり、次のことを期待していました：

autocorrelation（a）= <！>＃8721; _{-<！>＃8734; <！> lt; i <！> lt; <！>＃8734;} a _i v _{t + i} ここで、0 <！> lt; = t <！> lt; <！>＃8734;

しかし、得られたのは：

autocorrelation（a）= <！>＃8721; _{-<！>＃8734; <！> lt; i <！> lt; <！>＃8734;} a _i v _{t + i} ここで、-<！>＃8734; <！> lt; t <！> lt; <！>＃8734;

行う必要があるのは、相関結果の最後の半分を取得することです。これは、探している自己相関であるべきです。それを行う簡単なpython関数は次のようになります。

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

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

Answer 2

自己相関には、統計と畳み込みの2つのバージョンがあります。少し詳細を除いて、両方とも同じことを行います。統計バージョンは、間隔[-1,1]にあるように正規化されます。統計を行う方法の例を次に示します。

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

Answer 3

代わりに numpy.corrcoef 関数を使用しますの統計的相関を計算する numpy.correlate tのラグ：

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

Answer 4

同じ問題に遭遇したばかりなので、数行のコードを共有したいと思います。実際、今までに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 ]、次のように正しい結果が得られるようです：

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="same"ではなくmode="full"を使用したことに注意してください。

Answer 5

このトピックに混乱を招くものが2つあると思います：

1. 統計的v.s.信号処理の定義：他の人が指摘したように、統計では、自己相関を[-1,1]に正規化します。
2. 部分的な対非部分的な平均/分散：時系列がラグ<！> gt; 0でシフトするとき、それらのオーバーラップサイズは常に<！> lt;元の長さ。元の平均と標準（非部分）を使用しますか、それとも常に変化するオーバーラップ（部分）を使用して新しい平均と標準を計算しますか？ （おそらくこれには正式な用語がありますが、ここでは<！> quot; partial <！> quot;を使用します）。

部分的なv.sで、1d配列の自己相関を計算する5つの関数を作成しました。非部分的な区別。統計からの式を使用するものもあれば、信号処理の意味で相関する使用もあり、これはFFTを介して行うこともできます。ただし、すべての結果は統計の定義における自己相関であるため、それらが互いにどのようにリンクされているかを示しています。以下のコード：

import numpy
import matplotlib.pyplot as plt

def autocorr1(x,lags):
    '''numpy.corrcoef, partial'''

    corr=[1. if l==0 else numpy.corrcoef(x[l:],x[:-l])[0][1] for l in lags]
    return numpy.array(corr)

def autocorr2(x,lags):
    '''manualy compute, non partial'''

    mean=numpy.mean(x)
    var=numpy.var(x)
    xp=x-mean
    corr=[1. if l==0 else numpy.sum(xp[l:]*xp[:-l])/len(x)/var for l in lags]

    return numpy.array(corr)

def autocorr3(x,lags):
    '''fft, pad 0s, non partial'''

    n=len(x)
    # pad 0s to 2n-1
    ext_size=2*n-1
    # nearest power of 2
    fsize=2**numpy.ceil(numpy.log2(ext_size)).astype('int')

    xp=x-numpy.mean(x)
    var=numpy.var(x)

    # do fft and ifft
    cf=numpy.fft.fft(xp,fsize)
    sf=cf.conjugate()*cf
    corr=numpy.fft.ifft(sf).real
    corr=corr/var/n

    return corr[:len(lags)]

def autocorr4(x,lags):
    '''fft, don't pad 0s, non partial'''
    mean=x.mean()
    var=numpy.var(x)
    xp=x-mean

    cf=numpy.fft.fft(xp)
    sf=cf.conjugate()*cf
    corr=numpy.fft.ifft(sf).real/var/len(x)

    return corr[:len(lags)]

def autocorr5(x,lags):
    '''numpy.correlate, non partial'''
    mean=x.mean()
    var=numpy.var(x)
    xp=x-mean
    corr=numpy.correlate(xp,xp,'full')[len(x)-1:]/var/len(x)

    return corr[:len(lags)]


if __name__=='__main__':

    y=[28,28,26,19,16,24,26,24,24,29,29,27,31,26,38,23,13,14,28,19,19,\
            17,22,2,4,5,7,8,14,14,23]
    y=numpy.array(y).astype('float')

    lags=range(15)
    fig,ax=plt.subplots()

    for funcii, labelii in zip([autocorr1, autocorr2, autocorr3, autocorr4,
        autocorr5], ['np.corrcoef, partial', 'manual, non-partial',
            'fft, pad 0s, non-partial', 'fft, no padding, non-partial',
            'np.correlate, non-partial']):

        cii=funcii(y,lags)
        print(labelii)
        print(cii)
        ax.plot(lags,cii,label=labelii)

    ax.set_xlabel('lag')
    ax.set_ylabel('correlation coefficient')
    ax.legend()
    plt.show()

出力図は次のとおりです。

5行すべてが表示されないのは、3行が重複しているためです（紫色で）。オーバーラップはすべて非部分的な自己相関です。これは、信号処理メソッド（np.correlate、FFT）からの計算では、オーバーラップごとに異なる平均/ stdが計算されないためです。

また、FFTを行う前に時系列に0をパディングしなかったため、循環FFTであるため、fft, no padding, non-partial（赤い線）の結果が異なることに注意してください。理由を詳しく説明することはできません。それは他の場所から学んだことです。

Answer 6

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

自己相関の数学的特性のみに基づいて信号の自己相関を計算できる数行のコードを皆さんと共有したいと思いました。つまり、自己相関は次の方法で計算できます。

1. 信号から平均値を減算し、不偏信号を取得します

2. 不偏信号のフーリエ変換を計算

3. 不偏信号のフーリエ変換の各値の二乗ノルムを取得して、信号のパワースペクトル密度を計算します

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)

Answer 7

このような自己相関にはtalib.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)

Answer 8

フーリエ変換と畳み込み定理の使用

時間の複雑さは N * log（N）です

def autocorr1(x):
    r2=np.fft.ifft(np.abs(np.fft.fft(x))**2).real
    return r2[:len(x)//2]

これは正規化された偏りのないバージョンで、 N * log（N）

でもあります

def autocorr2(x):
    r2=np.fft.ifft(np.abs(np.fft.fft(x))**2).real
    c=(r2/x.shape-np.mean(x)**2)/np.std(x)**2
    return c[:len(x)//2]

A。Levyが提供する方法は機能しますが、PCでテストしましたが、時間の複雑さは N * N

のようです

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

Answer 9

OPの質問に対する本当の答えは、Numpy.correlateドキュメントからのこの抜粋に簡潔に含まれていると思います：

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

これは、「mode」定義なしで使用すると、Numpy.correlate関数は、2つの入力引数に同じベクトルが与えられた場合（つまり、自己相関の実行に使用された場合）、スカラーを返すことを意味します。

Answer 10

パンダのないシンプルなソリューション：

import numpy as np

def auto_corrcoef(x):
   return np.corrcoef(x[1:-1], x[2:])[0,1]

Answer 11

pandas datetime Seriesのリターンを指定して、統計的な自己相関をプロットします。

import matplotlib.pyplot as plt

def plot_autocorr(returns, lags):
    autocorrelation = []
    for lag in range(lags+1):
        corr_lag = returns.corr(returns.shift(-lag)) 
        autocorrelation.append(corr_lag)
    plt.plot(range(lags+1), autocorrelation, '--o')
    plt.xticks(range(lags+1))
    return np.array(autocorrelation)

Answer 12

私は計算生物学者であり、確率過程の時系列のカップル間の自動/相互相関を計算しなければならなかったとき、np.correlateが必要な仕事をしていないことに気付きました。

実際、<=>に欠けていると思われるのは、距離<！>＃120591;。 可能なすべての時点の平均 です。

必要なことを行う関数を定義する方法は次のとおりです。

def autocross(x, y):
    c = np.correlate(x, y, "same")
    v = [c[i]/( len(x)-abs( i - (len(x)/2)  ) ) for i in range(len(c))]
    return v

これまでの回答のいずれも、この自動/相互相関のインスタンスをカバーしていないようです。

Answer 13

numpy.correlateの代替は statsmodelsで入手できます。 tsa.stattools.acf（）。これにより、OPで説明されているような自己相関関数が連続的に減少します。実装は非常に簡単です：

from statsmodels.tsa import stattools
# x = 1-D array
# Yield normalized autocorrelation function of number lags
autocorr = stattools.acf( x )

# Get autocorrelation coefficient at lag = 1
autocorr_coeff = autocorr[1]

デフォルトの動作は40 nlagsで停止しますが、これは特定のアプリケーションのnlag=オプションで調整できます。ページの一番下に、関数の背後にある統計情報の引用があります。