2019年12月5日木曜日

rEDM を用いた Empirical Dynamic Modeling: (1) 動態の再構成

20191205_Blogger0004

時系列データを用いた動態の再構成

1. 時系列データ

時系列データはある系のある変数を経時観察・経時測定して得られるデータです。例えば、毎日自分の (=興味ある系) 体重 (= 興味ある変数) を測定し続ければそれは立派な時系列データとなります。

時系列データを図示する際には、基本的には横軸に時間、縦軸に変数をプロットする場合がほとんどです。例えば、以下の図は R に標準で用意されているデータセット Nile (ナイル川の流量変化の時系列) を時系列データとして表示したものです。横軸に時間がとってあり、縦軸に変数がプロットされていますね。

enter image description here

2. 「時系列データから系の動態を描く」 という概念

次回以降に解説予定の Empirical Dynamic Modeling という時系列解析では 「時系列データから系の動態を描く」 という概念が重要になります。そこで、この概念について 「バタフライアトラクター」 (もしくはローレンツアトラクターとも呼ばれる) という有名な 3 変数の微分方程式を例に説明します (Lorenz 1963)。バタフライアトラクターとは以下の微分方程式において、p=10p = 10, r=28r = 28, b=8/3b = 8/3 としたときに現れる動態のことを言います。

dxdt=pX+pYdydt=XZ+rXYdzdt=XYbZ\frac{dx}{dt} = -pX + pY\\ \frac{dy}{dt} = -XZ + rX - Y\\ \frac{dz}{dt} = XY - bZ

この方程式において、p=10p = 10, r=28r = 28, b=8/3b = 8/3 としたときに生成される時系列データは以下に図示するとおりです。
enter image description here

さて、上の 3 つの時系列データはそれだけで動態の特徴を表しています。しかし、この動態の表し方はこれだけでしょうか?見方を変えるだけで同じデータを全く違う方法で表すことが可能です。例えば、XX, YY, ZZ を 3 次元空間の xx軸・yy軸・zz軸にプロットすると以下のようになります (XXZZをプロットして色でYYを表しています)。

enter image description here

この蝶 (バタフライ) のような形 (= バタフライアトラクター) は、時間情報以外は元の 3 つの時系列プロットで表したものと同じ情報を含んでいます。どの時間・どの場所に点がくるか、という情報はここには表示されておらず、例えばアニメーションで点の動きをこのプロット内に示せば、元の情報を余すことなく表現したことになります。

さて、今回は 3 つの時系列データを元にして 3 次元空間に一つの動態を描きましたが、見方を変えると 「3 つの時系列データは 3 次元空間の中にある 1 つの動態から得られたもの」 と見ることもできるでしょう。つまり、バタフライアトラクター上の点の動きを xx 軸に射影すれば XX の時系列データが得られ、yy, zz 軸に射影すればそれぞれ YY, ZZ の時系列データが得られます。時系列データは (実際に自分たちが観測できるかどうかは別として) 背後にある 「真の動態」 を何らかの軸に射影して得られたもの、と捉えることができます。例えば、人間の体重は背後にある様々な要因 (食べた物・運動量・ストレスなどなど) が複雑に絡み合った (おそらく観察できない) 動態を体重計を用いることで、その動態を一つのある軸に射影した結果得られたデータ、と捉えることもできます。

以上の解説は YouTube にある以下のアニメーションをみることでよりよく理解できるかもしれません (ただし、英語です)。
Introduction to Empirical Dynamic Modeling

3. 時間遅れ時系列による状態空間の再構成

自分が興味ある系の動態を 「きちん」 と描くことができれば非常に有用です (「きちんと」 とは何なのかはあとで考えることにします)。例えば、過去の系の挙動を参照して近い将来何が起こりそうかを予測したり、どのくらいの外力を加えると系の状態が元に戻らなくなるかを予測したりすることもできるかもしれません。

しかし、興味ある系が野外の系であった場合、その動態をきちんと描くことは容易ではありません。なぜなら、これまでに説明した方法で動態をきちんと描くためにはバタフライアトラクターの例のように系の動態に関わっている全ての変数の時系列データを得る必要があるからです。野外の系の場合は、その系の動態に関わっている変数がいくつあるのか事前には分かりません。2 個かもしれないし、10 個かもしれないし、はたまた 100 個かもしれません。また、仮に何個の変数を測定すればいいかが分かったとしても、それらが全て技術的に測定可能かどうかは別問題です。さらに技術的には測定可能でも、労力的・金銭的に測定が不可能ということもありそうです。このようなことを考えると系の動態に関連する変数を全て測定し、きちんと興味ある系の動態を描くことはほとんど不可能のように思えてきます。

このような状況下で系の状態をきちんと描く (= 再構成と呼ぶことにします) 手助けとなるある数学定理が存在します。Takens の埋め込み定理 (Takens 1981) と呼ばれる定理は、時間遅れ座標系を取ることで 1 変数からでも系の動態を再構成できることを示されています (ちなみに 「埋め込み定理」 は Takens のもの以外にもたくさんあります)。

具体的にどういうことか見てみましょう。

まず、X(t)X(t) を元に時間遅れ座標系を作ってみます。ここでは X(t1)X(t-1)X(t2)X(t-2) を作ってみましょう。その後、3 次元空間に元の {X(t),Y(t),Z(t)}\{X(t), Y(t), Z(t)\} の代わりに {X(t),X(t1),X(t2)}\{X(t), X(t-1), X(t-2)\} をプロットしてみます。すると、下の左のような動態が現れます。同様のことを Y(t)Y(t) についても行うと、今度は下の右の図のような動態が現れます。

enter image description here

XX, YY, ZZ を 3 次元空間にプロットした図と見比べてみてください。時間遅れ時系列のプロットにより描かれる形が元の形と非常によく似ていることが分かると思います。

すなわち、たとえ興味ある系の中から 1 変数しかモニタリングできなかったとしても、時間遅れ埋め込み (Time-delay embedding) の方法により系の動態を再構成できる可能性がある、ということです。

4. なぜ再構成できるのか?

時間遅れを取ってそれをプロットするだけで系の動態が再構成できるというのは (特に実証研究者にとっては) 驚異的です。なぜこんなことが可能なのでしょうか?定理の証明を理解するには微分幾何学の知識が必要で、(自分も含めて) 実証研究を行う生態学者には少しハードルが高いです。そこで Munch et al. (2019) で紹介されている手順で定理を解釈してみます。

まず、測定可能な変数 yy と測定不可能な変数 zz を考え、その時間変化は以下の式で表されるとします。
yt+1=F(yt,zt)zt+1=G(yt,zt)y_{t+1} = F(y_t, z_t)\\ z_{t+1} = G(y_t, z_t)
つまり、yt+1y_{t+1}yty_tztz_t が何らかの関数 FF で変換されて生成され、zt+1z_{t+1}yty_tztz_t が何らかの関数 GG で変換されて生成される状況です。この状況で Takens の埋め込み定理は yy だけで動態を再構成できる、と言っているわけです。

さて、zt+1=G(yt,zt)z_{t+1} = G(y_t, z_t) なので、zt=G(yt1,zt1)z_{t} = G(y_{t-1}, z_{t-1}) です。これを yt+1y_{t+1} を表す式に入れると以下のようになります。
yt+1=F(yt,zt)=F(yt,G(yt1,zt1))y_{t+1} = F(y_t, z_t) = F(y_t, G(y_{t-1}, z_{t-1}))

この時点で、もし運が良ければ yt+1=F(yt,zt)y_{t+1} = F(y_t, z_t) を利用して zt1z_{t-1} を別の表し方で書けるかもしれません。つまり、時間を一つ戻して yt=F(yt1,zt1)y_t = F(y_{t-1}, z_{t-1}) とし、これは yty_t, yt1y_{t-1}, zt1z_{t-1} からなる方程式なので、zt1z_{t-1} について解くことができるかもしれません。解けたときに zt1=Φ(yt,yt1)z_{t-1} = \Phi(y_t, y_{t-1}) と表せるとしましょう。すると、yt+1y_{t+1} は以下のように表されます。
yt+1=F(yt,zt)=F(yt,G(yt1,Φ(yt,yt1)))y_{t+1} = F(y_t, z_t) = F(y_t, G(y_{t-1}, \Phi(y_t, y_{t-1})))

どうでしょうか?yt+1y_{t+1}yty_tyt1y_{t-1} という時間遅れで表現され、zz が消えました。

もし zt1=Φ(yt,yt1)z_{t-1} = \Phi(y_t, y_{t-1}) と解けない場合は zt1=G(yt2,zt2)z_{t-1} = G(y_{t-2}, z_{t-2})yt+1=F(yt,G(yt1,G(yt2,zt2)))y_{t+1} = F(y_t, G(y_{t-1}, G(y_{t-2}, z_{t-2}))) として、zt2z_{t-2} をうまく yty_t, yt1y_{t-1}, yt2y_{t-2} で表すことができないか考えることになります。この作業を繰り返すことで yyzz で表されていたものを yy の時間遅れにより表すことができます。以上について、詳しくは Munch et al. (2019) の 3-4 ページ目を参照して下さい。

さて、次回のポストでは動態の再構成を利用してどのように時系列データの解析を行っていくか、rEDM という R のパッケージの使い方とともに解説していく予定です。

参考文献

  • Lorenz, E. N. (1963) Deterministic Nonperiodic Flow, Journal of Atmospheric Sciences, Vol.20, pp.130-141
  • Takens, F. (1981) Detecting strange attractors in turbulence. In D. Rand & L.-S. Young (Eds.), Lecture Notes in Mathematics (Vol. 898, pp. 366–381).
  • Munch, S. B., Brias, A., Sugihara, G., Rogers, T. L. (2019) Frequently asked questions about nonlinear dynamics and empirical dynamic modelling. ICES Journal of Marine Science, doi:10.1093/icesjms/fsz209

Written with StackEdit.

0 件のコメント:

コメントを投稿

Empirical Dynamic Modeling with rEDM: (2) Near-future forecasting (Simplex projection)

20231015_Blogger0021 *This is an English version of my previous post (first translated by ChatGPT, checked and edited ...