カオス時系列データのべき乗則の確認(python)



スポンサードリンク

これまでに、カオス時系列データに関する勉強として、レスラー方程式やローレンツ方程式を扱ってきましたが、その中で得られた時系列データのデータの振れ幅(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はスケーリング指数です。

両対数グラフにすると、線形になるのが特徴。

ついでに片対数グラフで線形になるのは、指数関数です。

そしてこのべき乗則は自然現象や社会現象でもみられます。有名なところで言えば、株式市場の株価や地震の発生回数とその規模などもこれに従います。

なので、例えば株式市場で言えば、ちっちゃな下落や上昇はたくさんあるが歴史的な株価の暴落は数える程しか起きないと言うのがそうであります。

地震で言えば、震度1や2の地震は頻繁に発生するが、震度6や7の地震はか数えるほどしか起きないということです。

要は、事象の発生回数とその規模の相関関係を示す法則と言えます。

計算方法

結果を表示する前に今回の計算の流れ(方法)を述べます。計算に用いたコードはこの記事の最後に載せます。

1、Δt間での変化量ごとにグルーピングする(x(t)-x(t-1)の大きさごとにいくつかに分類)。地震の例で言えば、これが震度1、震度2、震度3のような大きさの定義になります。

2、Δtの分け方については、基本は1ステップとし、比較のためにもう少し変化の勾配がきつくなるようにΔtを設定した。

3、全ての時系列に対して行い、規模ごとに分類し配列にぶち込む。

4、Δx/Δt=1,1.1の時を最小の規模(地震で言えば震度1)とし、その1.5倍を次の規模(震度2)、2倍(震度3)、3倍と増やしていき、最小規模の6倍の規模までカウントするようにします。

5、なお今回の計算では、Δx/Δt、事象規模の間隔(1.5倍、2倍…等の分け方)については直感的に決めました。もしかしたら最適な決め方があるのかもしれません。

上の事象規模の分類の仕方ですと、最小規模の1.1〜1.4倍の事象も1.5倍の事象として分類していることになります。

なので精度を挙げるなら、ここの間隔(地震で言えば震度の分類)をより細かくするべきでしょう。


スポンサードリンク

べき分布の結果

以下に計算結果を示します。


Δt=1と1.1の時のデータです。もう少し線形なグラフが得られるかなと予測したのですが、全然違いました。ガウス分布の両対数グラフに近いかもしれません。

上でも述べたように、発生規模の分類の仕方に問題があるかもしれません。

あとローレンツ力学系はカオス系なので、初期値鋭敏性を持っています。

なので違う初期条件や計算時間をもっと長くしたりするなどしたら、結果はまた変わってくるかもしれません。

ついでに片対数(縦軸対数)グラフです。

こっちの方が線形ですね笑

ついでに、両軸とも対数にしていないグラフは以下です。




スポンサードリンク

計算に用いたコード

先に謝ります。ひどいです。いや本当にひどいコード。

配列の作成、Δtの場合分けをもっとうまくやればこんなに長くなりません。

奮闘しましたが、書けなかったのでなんとも原始的なやり方になってしまいました。

事象の規模の分類のところもどんだけelif使うんだよって感じですが笑

もう少しpython 勉強します笑。

べき分布の計算をする際に役立ててもらえれば幸いです。

記事が役に立ったらシェア!

Euler法とRunge-Kutta法のカオス時系列データの比較



スポンサードリンク

モチベーション

先日、非線形力学系ローレンツ方程式のアトラクタ描画をルンゲクッタ法で計算しました。x,y,zに関する時系列データも簡単に計算できたので、描画してみたところ以下のようなグラフになりました。以下はx-tグラフです。

上記のグラフがカオス的な振る舞いを示しているというのは前回の記事に示した通りですが、これがルンゲクッタではなく、違う近似解法の場合どういった挙動を示すのか気になりました。

ネットを見る限りオイラー法とルンゲクッタ法の解の誤差について、示している記事もいくつかあったのですが、1変数の一階微分方程式や単純なバネのモデルのような1変数の二階微分方程式しか見つかりませんでした。

なので、今回は3変数関数のローレンツ方程式を用いて、数値解法が異なることによる、誤差の程度を検証したいと思います(定量評価はしてません)。

Eular法とは

微分方程式の近似解法の中でも比較的計算が簡単なものですが、その分、とても粗い近似解法です。

研究等ではまず使われないですが、数値解放の概念を理解する手助けにはなります。以下のような方程式を考えます。
\begin{eqnarray}
\frac{dy}{dt}=\dot{y}=f(t,y(t))
\end{eqnarray}

次に微分の定義を用いて、
\begin{eqnarray}
\dot{y}=\lim_{\Delta{t} \to 0} \frac{y(t+\Delta t)-y(t)}{\Delta t}
\end{eqnarray}
0<t<<1の時
\begin{eqnarray}
\dot{y}=\frac{y(t+\Delta t)-y(t)}{\Delta t}
\end{eqnarray}
と表せる。これを最初の式を用いて式変形していくと、
\begin{eqnarray}
\dot{y}\Delta t = y(t+\Delta t)-y(t) \\
y(t+\Delta t) = y(t)+\Delta tf(t,y(t))
\end{eqnarray}
となり、t秒の情報(y(t)、f(t,y(t))を用いて、Δt秒後のy(t+Δt)を求めるのがEular法です。

計算ではこのΔtを小さくするほど、微分に近づくことになるので精度はよくなりますが、その分、計算時間が長くなります。

Runge-Kutta法とは

ルンゲクッタ法は、計算に少々時間がかかりますが、その精度の高さから、微分方程式の数値解析には広く用いられています。

今回計算に用いているのは、4次のルンゲクッタ法です。4次のルンゲクッタ法には公式が存在するのですが、その公式の導出をしている人はどうやらいないみたいです。

公式の導出にはオイラー法とホイン法(2次のルンゲクッタ法)を用いなければいけないみたいですが、計算量がえげつないみたいで、そんなものを書こうとするマゾな人はいないらしいです。

簡単に言うとテイラー展開の4次までの項の係数と一致するように計算するらしいです。まあいろんな人が書いているので、ここで詳しく書く必要もないでしょう。

どうしても知りたい人はこのサイトを参考にしてください笑。→http://eman-physics.net/math/differential21.html



スポンサードリンク

計算結果比較

以下に、オイラー法とルンゲクッタ法で計算したx-tグラフを示します。

最初の方(t=0から500)くらいでは、二つの数値解法にはあまり差がありませんが、徐々に差は大きくなり、やがて全く異なる挙動になります。

試しにt=0から500の範囲のグラフを載せます

解に大きなズレはなくこの程度の範囲ならば、オイラー法でも良さそうです。

しかし実際の科学計算ではオイラー法はまず用いられません。研究で使用されている方程式はより高次で、高変数のものが多く、計算範囲ももっと広い場合がほとんどです。

なので、この程度の方程式かつこの程度の時間範囲しか解の正確性を保てないオイラー法はまず使われません。

オイラー法は比較的計算も楽で、簡単に理解できるので、数値解法の玄関的な役割として、教科書等に書かれているだけです。

ついでに、1変数関数であるならば、ある程度計算時間が長くてもオイラー法とルンゲクッタ法の解に差はあまり見られません。 以下は一次元の減衰系振動モデルのx-tグラフです。

一次元減衰振動モデル

 

若干紫っぽくなってるのは、赤と青が重なっているためです。一次元の簡単な運動方程式程度であれば、二つの解放に誤差が出ないことが観察できます。

今回用いたプログラム

  • ローレンツ方程式

 

  • 一次元減衰振動モデル

ルンゲクッタ法とオイラー法で別々にリストを作って、別々にfor文を回しているので、すごい効率の悪いプログラムです。

if文等用いて変数を減らし、一つのfor文ですむように書き直したいと思います。


スポンサードリンク

最後に

今回は3変数関数におけるオイラー法とルンゲクッタ法の数値計算の誤差について述べました。

数値解法は他にもホイン法、ベルレ法、かえる跳び法、修正オイラー法等あるので、それらを用いて計算し、解の比較をしていきたいと思います。

ついでにルンゲクッタ法よりも僕は、計算回数と長期正確性を加味して、速度ベルレ法が一番おすすめです笑。

記事が役に立ったらシェア!

ローレンツ方程式をRunge-Kutta法(python)で計算



スポンサードリンク

前回は、アトラクターの勉強ということで、非線形力学系のレスラー方程式のアトラクター再構成の数値実験を行いました。

>>>前回の記事「Rossler方程式を Runge-Kutta法を用いて計算してみた(言語はpython)。」

>>>前々回の記事「多変量時系列解析に関する論文」

今回はローレンツ方程式について実施した数値実験の結果と、時間遅れの大きさについて論じて行きます。

時間遅れ座標系を用いた、アトラクタの再構成

まず、初めにローレンツ方程式を示します。
\[
\left\{
\begin{array}{ll}
\frac{dx}{dt}= -px+py\\
\frac{dy}{dt}= -xz+rx-y\\
\frac{dz}{dt}= xy-bz
\end{array}
\right.
\]
x,y,z 全ての変数の時系列データが観測できた場合のアトラクタを以下に示します。

図1 (x,y,z)の時系列データから得られるアトラクター

アトラクタ再構成のやり方には

  • 微分座標系(V(t)=(\(x(t),\dot x(t),\ddot y(t)\))を用いる
  • 時間遅れ座標系(V(t)=(\(x(t),x(t+\delta),x(t+2\delta)\)))を用いる

方法があるが、ノイズがある場合には、微分座標系はノイズを拡大させる原因にもなるので、時間遅れ座標系を用いる場合が多いです。

以下に観測できるデータがx(t)、時間遅れの大きさ\(\delta=1,10,100,500\)の場合での再構成したアトラクタを示します。

 

図2 時間遅れの大きさδ=1
図3 時間遅れの大きさδ=10
図4 時間遅れの大きさδ=100

 

 

図5 時間遅れの大きさδ=500

 

 

 

 

 

 

 

 

 

 

そして力学系のアトラクタは平衡点、リミットサイクル、kトーラス、ストレンジアトラクタ、ランダムに分類されます。

状態空間(位相構造)を観察する限りでは、ローレンツ方程式のアトラクタはストレンジアトラクタというものに分類されます。

このアトラクターが「ストレンジアトラクターになっている」というのが大切で、これが観察できれば、この非線形力学系は「カオス的な振る舞いをする」と言えます。

ターケンスの埋め込み定理

しかし、上記のような時間遅れ座標系を用いてアトラクタを再構成するには、どうやら以下のような定理が成り立たないといけないみたいです。

「d 次元のコンパクトな多様体 M と C 1 級微分同相写像 f : M → M,C 1 級関数 g : M → R 1 が与えられたとき, m = 2 d + 1 であ れば,次式の写像 V : M → R m は,埋め込みである」

だめだ、わけワカンねぇ。本当に同じ言語なのかと疑うほどであるが、m=2d+1が成り立てば良いらしい。

ここでmとはアトラクタの次元のことで、ここではm=3になります。

問題はdです。この「ターケンスの埋め込み定理」が成り立つには、d=1であれば良いのですが、僕はこのdが何を意味するのかわかりません。可観測な状態変数の数であるならば、今回は可観測なデータがx(t) のみを仮定しているので、成り立ちますがどうなのでしょう。

ここら辺、追求したいですが、体力が切れてきたので、ここはまた今度にします。

しかし前回読んだ論文では、ターケンスの埋め込み定理はあくまで十分条件らしいので、必ず成り立っていないといけないみたいではないようです。


スポンサードリンク

時間遅れの大きさの最適化

時間遅れの大きさの最適化には以下の二つの方法があります。

  • 線形自己相関関数が最初に0になる点を求める
  • 平均相互情報量が最初に0なる点を求める

以下では平均相互情報量の方法について述べて行きます。まずは定義から

\begin{eqnarray}
\hat{I}_{AB}(T)=\sum_{a_{i},b_{k}} P_{AB}(a_{i},b_{k})I_{AB}(a_{i},b_{k})
\end{eqnarray}
A={x(n)},B={x(n+T)}とすれば、
\begin{eqnarray}
\hat{I(T)}=\sum_{n=1}^N P(x(n),x(n+T))\log2(\frac{P(x(n),x(n+T))}{P(x(n))P(x(n+T))})
\end{eqnarray}
となります。そしてこの\(\hat{I(T)}\) の最初(シグマなので)の最小値をとる時のTを「最適な遅れ時間」として選択します。

そして、この遅れ時間の最適化の計算はプログラムを書いて計算し、また次回以降載せます。

計算するとδ=18くらいになるみたいです。上の再構成アトラクタ(図2〜図5)もδ=10が一番元の位相構造に近いので、適切な手法と言えそうです。

計算に用いたプログラム

以下に計算に用いたプログラムを載せます。

今回は、遅れ時間の時系列データもプログラムできました。しかし、遅れ時間のデータをfor文で別々に書いているので、ここを一つにまとめられたらよかったのです。

それと、遅れ時間も毎回、δの部分を書き換えて計算しているので、δ の配列を作り、配列の要素ごとに計算するというのも、やってみたいです。

アトラクタの軌跡のアニメーションも作ろうと思って計算したのですが、PCが低スペックでフリーズしたこともあり、今回はやめておきました。

そこらへんも、次回以降書いて行きたいと思います。



スポンサードリンク

記事が役に立ったらシェア!

多変量時系列解析に関する論文

最近、人工知能とかデータサイエンスとか、Pythonに興味があるので、勉強を開始しました。とりまテキトーに論文を漁る所から。
辿り着いた論文は”Analyzing multiple nonlinear time series with extended Granger causality” PDFはこちら



スポンサードリンク

論文の内容の前に

その前に日本語で書かれても全く理解出来ないです。何せ専門用語しらないのですから。
今回は私がこの論文を読むにあたり、調べた用語などのアウトプットとアブストラクトの要約を述べようと思います。
  1. 多変量解析
    多数の変量を同時に解析すること。多変量解析があるならば、一変量解析、二変量解析も当然あります。
    データ解析は数値計算をするだけでなく、データの可視化をすることで、データに対する
    理解をする手助けをすることも含めます。
    例えば、一変量解析ならば、ヒストグラムなど。二変量解析ならば、散布図、三変量ならば、3次元の散布図などです。
    人間は4次元以上のものは、グラフのようなデータの可視化ができないので、多変量解析はまずは、次元を落としてやるところから始めます。または、方程式を可視化する方法もあるんだとか。
  2. グレンジャー因果
    ある変数に他の変数を説明変数として加えたときに、平均二乗誤差に変化が出る場合、グレンジャー因果があるという。
    この論文ではこの手法を用いています。
    このグレンジャ因果は線形確率系に用いることで、多変量の関係を見つける手法らしいです。
  3. ターケンスの埋め込み定理
    1つの状態変数から多数の状態変数を復元する際に用いられる。この定理によると、復元したデータの位相空間は復元する前の位相空間(アトラクター)と微分的に同相となります。
  4. 相空間
    一般的な力学系では、運動を記述する際、位置と時間だけのグラフでは、完全に特徴を記述できず、N次元の空間を想定する。この空間のことを相空間あるいは、状態空間という。

アブストラクトの内容

多変量解析において、多変数の相関を特定することは重要である。そりゃ、そうだ!運動量が位置(例えばx)によって変化するのに、xの関数で表そうとしないなんて馬鹿な話はありませんね。この論文では、線形確率系において、使われたグレンジャー因果を非線形系に拡張し、その結果を述べるそうです。
特に3変量やそれ以上の時系列データに対して、変数の相関が見られるかもしくは、他の工程により分解されるかどうかを決定することが可能となる、方法を提案するらしい。
なるほど、非線形系に対してグレンジャー因果を応用したときに、変数の相関を確認することができるかどうかを明らかにしたいそうです。
まずは、グレンジャー因果と多変量解析の勉強をしていくところですかね!

 



スポンサードリンク

記事が役に立ったらシェア!