読者です 読者をやめる 読者になる 読者になる

統計学入門読み始め

P57 図3.20の自己相関係数計算

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt

data = np.array([4.96,11.15,14.37,9.25,1.01,-0.38,7.49,16.83,11.21,3.11,
            3.03,8.7,16.47,14.29,1.89,-7.99,-5.91,6.58,12.30,14.35,
            4.65,-1.31,-0.71,8.65,16.97,16.11,4.64,-1.01,0.78,8.20,
            12.95,9.55,8.61,4.88,-1.45,0.99,3.35,2.71,4.68,2.71,
            7.11,13.22,16.01,14.78,3.08,-7.12,-0.53,11.28,17.80,13.18,
            5.51,-2.82,-7.13,1.91,15.28,16.08,9.64,2.60,1.16,3.51])
#自己相関係数
def serial_corr(data,shift):
    #データ個数を算出
    n = len(data)
    #平均を算出
    m = data.mean()
    #分母を算出
    bunbo = sum(pow(data-m,2))/n
    #分子を算出
    bunsi = 0
    for i in range(n-shift):
        bunsi += (data[i]-m)*(data[i+shift]-m) 
    return (bunsi/(n-shift))/bunbo

plt.plot(data)
plt.show()
print(serial_corr(data,1))
print(serial_corr(data,2))
print(serial_corr(data,3))
print(serial_corr(data,4))
print(serial_corr(data,5))
print(serial_corr(data,6))
print(serial_corr(data,7))
print(serial_corr(data,8))
print(serial_corr(data,9))
print(serial_corr(data,10))
korero = [serial_corr(data,s) for s in range(1,11)]
plt.plot(korero)
plt.show()

<結果>

0.53209073795
-0.363539399131
-0.816080494564
-0.468670102984
0.26227969162
0.627346907724
0.328015213105
-0.255746034025
-0.545368160167
-0.259160467086

f:id:bitop:20160529154949p:plain
本とインデックスがずれている
f:id:bitop:20160529155109p:plain