スポンサードリンク
これまでに、カオス時系列データに関する勉強として、レスラー方程式やローレンツ方程式を扱ってきましたが、その中で得られた時系列データのデータの振れ幅（x(t)-x(t-1)）について何となく気になりました。
これまでのカオス時系列データに関する記事は以下です。
＞＞＞ローレンツ方程式をRunge-Kutta法（python）で計算
＞＞＞Euler法とRunge-Kutta法のカオス時系列データの比較
今回は得られた時系列データが理論にあっているか（べき乗則）の確認として、べき分布を表示してみました。
念のため、今回計算に用いた時系列データは以下のグラフです。
ついでに時系列データはローレンツ方程式をRunge-Kutta法（python）で計算の記事にあるコードを用いれば、簡単に計算できます。
べき乗則とは？
統計モデルの一つ。もっとも一般的な定義は
\begin{eqnarray}
f(x)=ax^{k}+o(x^{k})
\end{eqnarray}
oはランダウの記号、kはスケーリング指数です。
両対数グラフにすると、線形になるのが特徴。
ついでに片対数グラフで線形になるのは、指数関数です。
そしてこのべき乗則は自然現象や社会現象でもみられます。有名なところで言えば、株式市場の株価や地震の発生回数とその規模などもこれに従います。
なので、例えば株式市場で言えば、ちっちゃな下落や上昇はたくさんあるが歴史的な株価の暴落は数える程しか起きないと言うのがそうであります。
地震で言えば、震度１や２の地震は頻繁に発生するが、震度６や７の地震はか数えるほどしか起きないということです。
要は、事象の発生回数とその規模の相関関係を示す法則と言えます。
計算方法
結果を表示する前に今回の計算の流れ（方法）を述べます。計算に用いたコードはこの記事の最後に載せます。
１、Δt間での変化量ごとにグルーピングする（x(t)-x(t-1)の大きさごとにいくつかに分類）。地震の例で言えば、これが震度１、震度２、震度３のような大きさの定義になります。
２、Δtの分け方については、基本は１ステップとし、比較のためにもう少し変化の勾配がきつくなるようにΔtを設定した。
３、全ての時系列に対して行い、規模ごとに分類し配列にぶち込む。
４、Δx/Δt=1,1.1の時を最小の規模（地震で言えば震度１）とし、その1.5倍を次の規模（震度２）、２倍（震度３）、３倍と増やしていき、最小規模の６倍の規模までカウントするようにします。
５、なお今回の計算では、Δx/Δt、事象規模の間隔（1.5倍、２倍…等の分け方）については直感的に決めました。もしかしたら最適な決め方があるのかもしれません。
上の事象規模の分類の仕方ですと、最小規模の1.1〜1.4倍の事象も1.5倍の事象として分類していることになります。
なので精度を挙げるなら、ここの間隔（地震で言えば震度の分類）をより細かくするべきでしょう。
スポンサードリンク
べき分布の結果
以下に計算結果を示します。
Δt＝1と1.1の時のデータです。もう少し線形なグラフが得られるかなと予測したのですが、全然違いました。ガウス分布の両対数グラフに近いかもしれません。
上でも述べたように、発生規模の分類の仕方に問題があるかもしれません。
あとローレンツ力学系はカオス系なので、初期値鋭敏性を持っています。
なので違う初期条件や計算時間をもっと長くしたりするなどしたら、結果はまた変わってくるかもしれません。
ついでに片対数（縦軸対数）グラフです。
こっちの方が線形ですね笑
ついでに、両軸とも対数にしていないグラフは以下です。
スポンサードリンク
計算に用いたコード
先に謝ります。ひどいです。いや本当にひどいコード。
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 |
# -*- coding: utf-8 -*- import numpy as np import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D # lorenz equation def lorenz(x, y, z, p=10, r=28, b=8/3): x_dot = -p*x +p*y y_dot = -x*z +r*x-y z_dot = x*y -b*z return np.array([x_dot, y_dot, z_dot]) t=0 ##### 初期時間###### dt = 0.05 ##### 微分間隔###### stepCnt =50000 ######時間ステップ## # Need one more for the initial values xs = np.empty((stepCnt + 1,)) ys = np.empty((stepCnt + 1,)) zs = np.empty((stepCnt + 1,)) # Setting initial values xs[0], ys[0], zs[0] = (0.1, 0.5, 1.2) k0=[0,0,0] k1=[0,0,0] k2=[0,0,0] k3=[0,0,0] listData=[0] width=0.1 d=[] d2=[] d3=[] d4=[] d5=[] d6=[] d7=[] d8=[] d9=[] d10=[] d11=[] d12=[] d02=[] d22=[] d32=[] d42=[] d52=[] d62=[] d72=[] d82=[] d92=[] d102=[] d112=[] d122=[] delta=1 delta2=1.1 for i in range(stepCnt): x,y,z=xs[i],ys[i],zs[i] ########ルンゲクッタ#################### k0 = dt * lorenz(x,y,z) k1 = dt * lorenz(x+k0[0]/2., y+k0[1]/2., z+k0[2]/2.) k2 = dt * lorenz(x+k1[0]/2., y+k1[1]/2., z+k1[2]/2.) k3 = dt * lorenz(x+k2[0], y+k2[1], z+k2[2]) dx = (k0[0]+2.0*k1[0]+2.0*k2[0]+k3[0])/6.0 dy = (k0[1]+2.0*k1[1]+2.0*k2[1]+k3[1])/6.0 dz = (k0[2]+2.0*k1[2]+2.0*k2[2]+k3[2])/6.0 xs[i+1] = xs[i] + dx ys[i+1] = ys[i] + dy zs[i+1] = zs[i] + dz ######## listData.append(i) ##########xの変化規模の分類####################### if delta > abs(xs[i]-xs[i-1]) >= 0: d.append(i) elif 1.5*delta > abs(xs[i]-xs[i-1]) >= delta: d2.append(i) elif 2*delta > abs(xs[i]-xs[i-1]) >= 1.5*delta: d3.append(i) elif 2.5*delta > abs(xs[i]-xs[i-1]) >= 2*delta: d4.append(i) elif 3*delta > abs(xs[i]-xs[i-1]) >= 2.5*delta: d5.append(i) elif 3.5*delta > abs(xs[i]-xs[i-1]) >= 3*delta: d6.append(i) elif 4*delta >abs(xs[i]-xs[i-1]) >= 3.5*delta: d7.append(i) elif 4.5*delta >abs(xs[i]-xs[i-1]) >= 4*delta: d8.append(i) elif 5*delta >abs(xs[i]-xs[i-1]) >= 4.5*delta: d9.append(i) elif 5.5*delta >abs(xs[i]-xs[i-1]) >= 5*delta: d10.append(i) elif 6*delta >abs(xs[i]-xs[i-1]) >= 5.5*delta: d11.append(i) elif abs(xs[i]-xs[i-1]) >= 6*delta: d12.append(i) if delta2 > abs(xs[i]-xs[i-1]) >= 0: d02.append(i) elif 1.5*delta2 > abs(xs[i]-xs[i-1]) >= delta2: d22.append(i) elif 2*delta2 > abs(xs[i]-xs[i-1]) >= 1.5*delta2: d32.append(i) elif 2.5*delta2 > abs(xs[i]-xs[i-1]) >= 2*delta2: d42.append(i) elif 3*delta2 > abs(xs[i]-xs[i-1]) >= 2.5*delta2: d52.append(i) elif 3.5*delta2 > abs(xs[i]-xs[i-1]) >= 3*delta2: d62.append(i) elif 4*delta2 >abs(xs[i]-xs[i-1]) >= 3.5*delta2: d72.append(i) elif 4.5*delta2 >abs(xs[i]-xs[i-1]) >= 4*delta2: d82.append(i) elif 5*delta2 >abs(xs[i]-xs[i-1]) >= 4.5*delta2: d92.append(i) elif 5.5*delta2 >abs(xs[i]-xs[i-1]) >= 5*delta2: d102.append(i) elif 6*delta2 >abs(xs[i]-xs[i-1]) >= 5.5*delta2: d112.append(i) elif abs(xs[i]-xs[i-1]) >= 6*delta2: d122.append(i) ##################グラフにするためのデータの用意（配列にぶち込む）############################ lenlist=[np.log10(len(d)),np.log10(len(d2)),np.log10(len(d3)),np.log10(len(d4)) ,np.log10(len(d5)),np.log10(len(d6)),np.log10(len(d7)),np.log10(len(d8)) ,np.log10(len(d9)),np.log10(len(d10)),np.log10(len(d11)),np.log10(len(d12))] slope=[np.log10(delta),np.log10(1.5*delta),np.log10(2*delta),np.log10(2.5*delta),np.log10(3*delta),np.log10(3.5*delta) ,np.log10(4*delta),np.log10(4.5*delta),np.log10(5*delta),np.log10(5.5*delta),np.log10(6*delta),np.log10(6.5*delta)] lenlist2=[np.log10(len(d02)),np.log10(len(d22)),np.log10(len(d32)),np.log10(len(d42)) ,np.log10(len(d52)),np.log10(len(d62)),np.log10(len(d72)),np.log10(len(d82)) ,np.log10(len(d92)),np.log10(len(d102)),np.log10(len(d112)),np.log10(len(d122))] slope2=[np.log10(delta2),np.log10(1.5*delta2),np.log10(2*delta2),np.log10(2.5*delta2),np.log10(3*delta2),np.log10(3.5*delta2) ,np.log10(4*delta2),np.log10(4.5*delta2),np.log10(5*delta2),np.log10(5.5*delta2),np.log10(6*delta2),np.log10(6.5*delta2)] lenlist3=[len(d),len(d2),len(d3),len(d4),len(d5),len(d6) ,len(d7),len(d8),len(d9),len(d10),len(d11),len(d12)] slope3=[delta,1.5*delta,2*delta,2.5*delta,3*delta,3.5*delta ,4*delta,4.5*delta,5*delta,5.5*delta,6*delta,6.5*delta] lenlist4=[len(d02),len(d22),len(d32),len(d42) ,len(d52),len(d62),len(d72),len(d82) ,len(d92),len(d102),len(d112),len(d122)] slope4=[delta2,1.5*delta2,2*delta2,2.5*delta2,3*delta2,3.5*delta2 ,4*delta2,4.5*delta2,5*delta2,5.5*delta2,6*delta2,6.5*delta2] ###########all-time-datasiries(graph-setting)################## ''' plt.plot(slope, lenlist3, c='b', label="$\Delta$t=1",marker="o") plt.plot(slope2, lenlist4, c='r', label="$\Delta$t=1.1",marker="o") ax2 = plt.gca() ax2.set_ylabel("number") ax2.set_xlabel("$\Delta$x-scale") ax2.set_title("power law") plt.legend() plt.show() ''' plt.plot(slope, lenlist, c='b', label="$\Delta$t=1",marker="o") plt.plot(slope2, lenlist2, c='r', label="$\Delta$t=1.1",marker="o") ax1 = plt.gca() ax1.set_ylabel("number") ax1.set_xlabel("$\Delta$x-scale") ax1.set_title("power law") plt.legend() plt.show() |
配列の作成、Δtの場合分けをもっとうまくやればこんなに長くなりません。
奮闘しましたが、書けなかったのでなんとも原始的なやり方になってしまいました。
事象の規模の分類のところもどんだけelif使うんだよって感じですが笑
もう少しpython 勉強します笑。
べき分布の計算をする際に役立ててもらえれば幸いです。
