スポンサードリンク
多変量時系列解析の勉強その2
前回の記事で「アトラクター」という単語が出てきたので、今回はその用語を中心に書いていこうと思う。
>>>前回の記事「多変量時系列解析に関する論文」。
とりあえずググると、アトラクターとは、「ある力学系がそこに向かって時間発展をする集合のことである。カオスな力学系に対してアトラクターを描写することは、現在においてもカオス理論における一つの研究課題である。」とか「実世界における力学系は散逸的であることが多いであろう。すなわち、もし力学系に運動の駆動力が無ければ、運動は停止するものと考えられる(そのような散逸は、様々な原因による内部摩擦や熱力学的損失、物質の損失などにより生じうる)。散逸と駆動力が組み合わさることにより、初期の摂動を鎮め、その力学系の振る舞いを典型的なものへと落ち着かせる傾向にある。そのような典型的な振る舞いに対応している力学系からなる位相空間の一部分はattracting section または attractee と呼ばれる。」とかウィキペディアには書いてある。
しかしこれだけではもちろん意味不なので、論文を読んでみた。
>>>読んだ論文はこちら。
今回はこちらの内容を元に、勉強したことを書いていこうと思います。
埋め込みとアトラクタ表示
例えばある3変数関数f(x,y,z)があるとする。しかし得られたデータは状態変数xの時系列データ(x(t))であるとする。そしてここから遅れ時間\(/delta\)なる変数を新たに導入し、新たにn次元ベクトル\({\bf V}= (x(t),x(t+\delta),…,x(t+(n-1)\delta))(t = 0, 1, · · · , N − 1 − (n − 1)\delta)\)を生成する。
以下では、例として、連続力学系のRossler方程式を取り上げる。
\[
\left\{
\begin{array}{ll}
\frac{dx}{dt}= -y-z\\
\frac{dy}{dt}= x+ay\\
\frac{dz}{dt}= b+z(x-c)
\end{array}
\right.
\]
なお今回は、a = 0.398, b = 2, c = 4 とし,ルンゲクッタ法を用 いて数値計算を行なった.
まず、全ての状態変数(x、y、z)を用いて軌道を描いたアトラクターが以下である。
図1 全ての状態変数を用いて描いたアトラクター
次に、3つの状態変数のうちxのみが観測可能であると仮定し、遅れ時間\(\delta\)を用いて。アトラクタを再構成したものが以下である。なお遅れ時間\(\delta=1,10,500\)を用いて時系列データを作成した。
図2 遅れ時間[latex]\delta=1[\mathrm{s}][/latex]の再構成アトラクタ
図3 遅れ時間[latex]\delta=10[\mathrm{s}][/latex]の再構成アトラクタ
図4 遅れ時間[latex]\delta=500[\mathrm{s}][/latex]の再構成アトラクタ
図からわかる通り、図3のアトラクタの位相は全ての状態変数が観測できた時に描けるアトラクタ(図1)と位相構造が似ている。
そして今回の計算ではTakensの埋め込み定理は満たされていないが、それでも元の力学系の位相的構造が保存されていることがわかる。
本来、アトラクタの再構成を適切におこなうには、時間遅れ\(\delta\)は時系列データの「自己相関関数」や「パワースペクトラム」を用いるらしいが、今回のような状態変数が少なく、データを視覚的に捉えやすい場合は、必ずしもそうする必要がないことがわかる。
スポンサードリンク
Runge-Kutta法によるR\(\ddot{\mathrm{o}}\)ssler方程式の計算
以下に、今回の計算に用いたプログラムを載せる。
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
|
# -*- coding: utf-8 -*- import numpy as np import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D import csv # rossler equationの関数作成 def lorenz(x, y, z, a=0.398, b=2, c=4): x_dot = -y -z y_dot = x + a*y z_dot = b + z*(x-c) return np.array([x_dot, y_dot, z_dot]) #初期値指定 t=0 dt = 0.05 stepCnt = 7000 xs = np.empty((stepCnt + 1,)) ys = np.empty((stepCnt + 1,)) zs = np.empty((stepCnt + 1,)) 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] # 書き込みファイルの指定 f = open('xdata.csv', 'w') csvWriter = csv.writer(f) #ルンゲ-クッタ法 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 = [] #listの初期化 listData.append(x) csvWriter.writerow(listData) #1行書き込み #グラフ描画設定 fig = plt.figure() ax = fig.gca(projection='3d') ax.plot(xs, ys, zs) ax.set_xlabel("X Axis") ax.set_ylabel("Y Axis") ax.set_zlabel("Z Axis") ax.set_title("Lorenz Attractor") f.close() plt.show() |
なお、上のプログラムは全ての状態変数が観測できた時の、アトラクタを描画する際のプログラムである。
本来なら、このプログラムに時間遅れを加味した時系列データをいい感じに元のデータにappendでもしたり、二次元配列のデータの任意の一部のデータを除去したりできればよかったです。ですがそこまでの実力がなかったので(一応数時間奮闘はしました笑)、時間遅れの再構成のアトラクタは別にプログラムを書いた(とは言ってもデータ処理の部分のみ)。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22
|
import numpy as np import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D import csv import pandas as pd f = open('xdata.csv', 'rb') f = open('xdata-2.csv', 'rb') f = open('xdata-3.csv', 'rb') fig = plt.figure() ax = fig.gca(projection='3d') x = np.loadtxt("xdata.csv", unpack=True) y = np.loadtxt("xdata-2.csv", unpack=True) z = np.loadtxt("xdata-3.csv", unpack=True) ax.plot(x, y, z) ax.set_xlabel("X") ax.set_ylabel("Y") ax.set_zlabel("Z") ax.set_title("Takens embedding '$\delta=1$'") f.close() plt.show() |
xdata.csvは観測できる時系列データです。そして、xdata-2とxdata-3は時間遅れの新規時系列データx(t+\(\delta\))とx(t+2\(\delta\))ということになります。
これらのデータはとても原始的に用意しました笑。xdataの時系列データの時間遅れの分だけ上の行から取り除いたデータを、xdata2にコピペし、さらに、そのxdata2の時系列データから、時間遅れの分だけ上の行から取り除いたデータを、xdata3にコピペする。
最後に時間遅れの分、xdata2,xdata3は、要素数(ファイルの行数)が少ないので、一番データ行数の少ない、xdata3と行数が合うようにデータを消して、3つのファイルの行数を揃えます。
こんな原始的にやらなくても、pythonならこれくらいのこと瞬殺できそうですが、ここはまた今度にします(pythonの無駄遣い笑)。
最後に
今回は「多変量時系列解析」に関してしっかりと理解するため、論文のような計算結果が得られるか確認しました。ルンゲクッタとかはググったら瞬殺できたのでよかったのですが、pythonに慣れていなく、データ処理の所に時間がかかってしまいました。
今後はpythonのデータ処理と、ローレンツ方程式ではどんな計算結果が得られるのか勉強したいと思います。