python3で自己相関関数

久々に読んだら，規格化数列がまどろっこしかったので修正し実行速度を再測定．
マシンが変わったので実行速度も変化した．

---

久々に読んだら関数名が意味不明だったので比較した関数の説明を追記．

---

なんか，思った通りの結果が出なかったので計算内容を確認したのと実行速度を測定したのと．

比較した関数

• corr1
  • hadacchi の自作関数
  • numpyのベクトル演算でlagだけずらした時系列間の積をとり，その積のベクトルをmeanメソッドで平均して自己相関を計算する関数
  • 他の3関数と比較して周期関数のピークが減衰せずに出てくる（正弦波の自己相関関数が減衰するとか明らかにおかしい）というメリットはありつつも， デメリットとしてnlagsをあまりに大きい値にしてしまうと，相関を計算するサンプルが少なすぎて相関係数の値が不安定（時系列をサンプルする毎に大きく値が変化してしまう）になる点がある
• corr3
  • corr2は作ったけど意味なかったので消して欠番
  • numpyのベクトル演算までは一緒だがsumメソッドで内積を計算した後，原系列長で除算して平均して自己相関を計算する関数
    • lagだけ時系列をシフトするとリテラル数が減るのに除算の母数が原系列長と一定なので変だと思うんだけど，この関数と比較することで，実はscipyとかstatsmodelsなんかの自己相関関数はこの実装になっていることに気付いた
• statsmodels.api.tsa.acf
  • statsmodelsモジュールにある自己相関関数を計算してくれる関数
  • いかなるlagの相関係数の計算であっても原系列長で除算していることが，corr3との比較で分かった
• scipy.signal.correlate
  • scipyモジュールにある自己相関関数を計算してくれる関数…かと思いきや内積しか計算してくれないから，ユーザが返り値を原系列長とかで除算しないと使えない

検算

コード

import scipy.signal as sig
import numpy as np
import statsmodels.api as sm
import matplotlib.pyplot as plt
import pandas
from IPython.core.display import display

x = np.arange(0,10000,0.1)
y = np.sin(x)
ny = (y-y.mean())/y.std()
dfy = pandas.DataFrame({'y':y})

def corr1(ts1, ts2, nlags=40):
    ts1, ts2 = np.asarray(ts1), np.asarray(ts2)
    nlags += 1
    N      = min(len(ts1),len(ts2))
    nts1   = ((ts1-ts1.mean())/ts1.std())[:N]
    nts2   = ((ts2-ts2.mean())/ts2.std())[:N]
    nlags  = min(N,nlags)
    result = [(nts1*nts2).mean()]
    for lag in range(1,nlags):
        result.append((nts1[lag:]*nts2[:-lag]).mean())
    return np.asarray(result)

def corr3(ts1, ts2, nlags=40):
    ts1, ts2 = np.asarray(ts1), np.asarray(ts2)
    nlags += 1
    N      = min(len(ts1),len(ts2))
    nts1   = ((ts1-ts1.mean())/ts1.std())[:N]
    nts2   = ((ts2-ts2.mean())/ts2.std())[:N]
    nlags  = min(N,nlags)
    result = [(nts1*nts2).sum()/N]
    for lag in range(1,nlags):
        result.append((nts1[lag:]*nts2[:-lag]).sum()/N)
    return np.asarray(result)

acf = sm.tsa.acf(y, nlags=5000)
sigcorr_buf = sig.correlate(ny,ny,mode='same')
sigcorr = sigcorr_buf[len(sigcorr_buf)//2:]  # 負のlagは消す
sigcorr /= len(sigcorr_buf)  # 系列長で規格化

myacf1 = corr1(y,y,nlags=5000)
myacf3 = corr3(y,y,nlags=5000)

display(pandas.DataFrame({ 'acf':acf
                         , 'corr1':myacf1
                         , 'corr3':myacf3
                         , 'sigcorr/N':sigcorr[:5001]
                         }))

result

まとめると，

• statsmodels.api.tsa.acfはlagが1以上の場合も，0の場合と同じ系列長を母数として平均値を計算している
• scipy.signal.correlateはそもそも原信号の規格化すらしてないので，全部手でやる必要がある
• scipy.signal.correlateは負の方向のlagにも計算をするため前半の結果は捨てても良い
  • 実は自己相関もFFTも実数系列を対象に処理すると正負方向に線対称なので，ホントは計算量を半分くらいに削減できるはずである
  • バタフライ演算から直さないといけないので，詳しくは numerical recipes を見る
• 平均をちゃんと取るcorr1（mean）と比べて絶対値で少し小さくなるのは，lagが大きくなると，例えばlag=4000では コンボリュージョンの計算に 10000-4000=6000 個の値しか使えてないのに，平均値の計算をlags=10000で計算しているため
  • なので，正弦波の自己相関を計算した場合，画像のように前者は自己相関のピーク値が減衰していくが，corr1は減衰しない

idx	acf	corr3	corr1	sigcorr/N
0	1.000000	1.000000	1.000000	1.000000
1	0.996869	0.996869	0.996879	0.996869
2	0.987520	0.987520	0.987540	0.987520
3	0.972032	0.972032	0.972061	0.972032
4	0.950536	0.950536	0.950574	0.950536
5	0.923214	0.923214	0.923260	0.923214
6	0.890298	0.890298	0.890352	0.890298
7	0.852067	0.852067	0.852126	0.852067
8	0.808842	0.808842	0.808906	0.808842
9	0.760986	0.760986	0.761055	0.760986
10	0.708902	0.708902	0.708973	0.708902
11	0.653023	0.653023	0.653095	0.653023
12	0.593815	0.593815	0.593886	0.593815
13	0.531767	0.531767	0.531836	0.531767
14	0.467389	0.467389	0.467454	0.467389
15	0.401208	0.401208	0.401268	0.401208
16	0.333761	0.333761	0.333814	0.333761
17	0.265590	0.265590	0.265635	0.265590
18	0.197240	0.197240	0.197276	0.197240
19	0.129250	0.129250	0.129275	0.129250
20	0.062149	0.062149	0.062162	0.062149
21	-0.003546	-0.003546	-0.003547	-0.003546
22	-0.067341	-0.067341	-0.067356	-0.067341
23	-0.128762	-0.128762	-0.128791	-0.128762
24	-0.187364	-0.187364	-0.187409	-0.187364
25	-0.242735	-0.242735	-0.242796	-0.242735
26	-0.294497	-0.294497	-0.294574	-0.294497
27	-0.342314	-0.342314	-0.342406	-0.342314
28	-0.385888	-0.385888	-0.385996	-0.385888
29	-0.424967	-0.424967	-0.425091	-0.424967
...	...	...	...	...
4971	-0.089497	-0.089497	-0.094179	-0.089497
4972	-0.113837	-0.113837	-0.119793	-0.113837
4973	-0.140299	-0.140299	-0.147641	-0.140299
4974	-0.168542	-0.168542	-0.177364	-0.168542
4975	-0.198200	-0.198200	-0.208576	-0.198200
4976	-0.228886	-0.228886	-0.240871	-0.228886
4977	-0.260194	-0.260194	-0.273822	-0.260194
4978	-0.291705	-0.291705	-0.306987	-0.291705
4979	-0.322991	-0.322991	-0.339915	-0.322991
4980	-0.353620	-0.353620	-0.372153	-0.353620
4981	-0.383160	-0.383160	-0.403245	-0.383160
4982	-0.411182	-0.411182	-0.432741	-0.411182
4983	-0.437269	-0.437269	-0.460201	-0.437269
4984	-0.461016	-0.461016	-0.485199	-0.461016
4985	-0.482038	-0.482038	-0.507329	-0.482038
4986	-0.499972	-0.499972	-0.526208	-0.499972
4987	-0.514479	-0.514479	-0.541483	-0.514479
4988	-0.525255	-0.525255	-0.552830	-0.525255
4989	-0.532025	-0.532025	-0.559962	-0.532025
4990	-0.534555	-0.534555	-0.562631	-0.534555
4991	-0.532649	-0.532649	-0.560630	-0.532649
4992	-0.526151	-0.526151	-0.553797	-0.526151
4993	-0.514953	-0.514953	-0.542016	-0.514953
4994	-0.498989	-0.498989	-0.525219	-0.498989
4995	-0.478243	-0.478243	-0.503387	-0.478243
4996	-0.452743	-0.452743	-0.476551	-0.452743
4997	-0.422566	-0.422566	-0.444793	-0.422566
4998	-0.387838	-0.387838	-0.408242	-0.387838
4999	-0.348727	-0.348727	-0.367077	-0.348727
5000	-0.305450	-0.305450	-0.321526	-0.305450

プロット

速度測定

先に信号を100系列作っておいて，自己相関を計算する．

コード

import copy
import scipy.signal as sig
import numpy as np
import statsmodels.api as sm
import time
import matplotlib.pyplot as plt
import pandas

# 上のソースと同じ
def corr1(ts1, ts2, nlags=40):
    ts1, ts2 = np.asarray(ts1), np.asarray(ts2)
    nlags += 1
    N      = min(len(ts1),len(ts2))
    nts1   = ((ts1-ts1.mean())/ts1.std())[:N]
    nts2   = ((ts2-ts2.mean())/ts2.std())[:N]
    nlags  = min(N,nlags)
    result = [(nts1*nts2).mean()]
    for lag in range(1,nlags):
        result.append((nts1[lag:]*nts2[:-lag]).mean())
    return np.asarray(result)

def corr3(ts1, ts2, nlags=40):
    ts1, ts2 = np.asarray(ts1), np.asarray(ts2)
    nlags += 1
    N      = min(len(ts1),len(ts2))
    nts1   = ((ts1-ts1.mean())/ts1.std())[:N]
    nts2   = ((ts2-ts2.mean())/ts2.std())[:N]
    nlags  = min(N,nlags)
    result = [(nts1*nts2).sum()/N]
    for lag in range(1,nlags):
        result.append((nts1[lag:]*nts2[:-lag]).sum()/N)
    return np.asarray(result)

N=12000
nlags = 5000
Xs    = [np.random.normal(size=N) for i in range(100)]

# corr1
%%timeit
for x in Xs:
    cor = corr1(x, x, nlags=nlags)

# 代入
m1corr = copy.copy(cor)

# corr3
%%timeit
for x in Xs:
    cor = corr3(x, x, nlags=nlags)

# 代入
m3corr = copy.copy(cor)

# sm.tsa.acf
%%timeit
for x in Xs:
    cor = sm.tsa.acf(x, nlags=nlags)

# 代入
acfcorr = copy.copy(cor)

# scipy.signal.correlate を acf 風に
%%timeit
for x in Xs:
    nx = (x-x.mean())/x.std()
    sigcorr_buf = sig.correlate(nx,nx,mode='same')
    sigcorr = sigcorr_buf[len(sigcorr_buf)//2:]  # 正のlagのみ
    sigcorr /= len(sigcorr_buf)  # 定数で規格化

# 代入
sigcorr_acf = copy.copy(sigcorr[:nlags+1])

# scipy.signal.correlate を corr1 風に
%%timeit
for x in Xs:
    nx = (x-x.mean())/x.std()
    sigcorr_buf = sig.correlate(nx,nx,mode='same')
    sigcorr = sigcorr_buf[len(x)//2:]
    sigcorr /= np.arange(len(x),len(x)//2,-1)  # 計算に用いた系列長で規格化

# 代入
sigcorr_3 = copy.copy(sigcorr[:nlags+1])

# 式が間違えてないか確認
results = pandas.DataFrame({'c1':m1corr, 'c3':m3corr, 'acf':acfcorr, 'sigcorr_acf':sigcorr_acf, 'sigcorr_3':sigcorr_3})
results[4000:].plot(figsize=[16,9],style=['o','^','v','+','x',''],alpha=0.8,ylim=[0,0.04])  # 差が広がるlagが大きい領域のみplot
plt.show()

result

速度の出力(100系列分の合計時間）
scipy.signal.correlation くっそ早いよ

method	result (1回あたりの実行時間は100で除したもの)
corr1	16.3 s ± 137 ms per loop
corr3	12.1 s ± 38.6 ms per loop
acf	6.96 s ± 175 ms per loop
sigcorr_acf	218 ms ± 101 µs per loop
sigcorr_3	220 ms ± 117 µs per loop

$lag>4000$ の領域の相関関数を比較して，特にsigcorrの関数外操作に間違えがないかチェック．
確かに，corr1とacfとsigcorr_arcは等しく，corr3とsigcorr_3も等しい．