数理解析研究所講究録
第 1919 巻 2014 年 38-60

38

時系列データに対する確率微分方程式モデルの統計的係数決定公式と軌道の予測可能性について
中野直人 (Naoto Nakano)*

Advanced Institute for Materials Research, Tohoku University稲津將 (Masaru Inatsu)

Graduate School of Science, Hokkaido University
楠岡誠一郎 (Seiichiro Kusuoka)

Graduate School of Science, Tohoku University
齊木吉隆 (Yoshitaka Saiki)

Graduate School of Commerce and Management, Hitotsubashi University

1序状態が時々刻々変化していく現象に対して，現在の状態から未来の状態をどれだけ予測できるかはその現象を表すモデルと現在の状態を表すデータの正確さにかかつている．とはいえ，実際はいわゆるモデル化誤差や観測誤差が発生するため，モデルとデータを用いた解析では何らかの不確定な要素が常に入り込むものである．
そのためそのようなデータ解析では，理想的な 「正確さ」 を追求できるとは限らないため，確率論的枠組みが用いられることが多い．たとえば経済学 [6], 気象学[1, 5], 生物モデル [13], 材料科学 [14] などの幅広い研究分野で確率論的モデルによる研究が進められている．確率論的モデルと一口に言ってもアプローチの仕方は多岐にわたり，例えば回帰分析のような予測モデルや隠れマルコフモデルなど様々な手法が開発されているが，ここでは確率論的枠組みの中でも時間連続的に変化するデータに対して用いられる確率微分方程式モデルに焦点を当てることにする$*e$-mail: n-nakano@wpi-aimr.tohoku.ac.jp

39[3, 4, 7, 9, 10].

確率微分方程式によるモデリングは，時系列データが相空間内の状態の変化 (軌道 として表されるものと解釈し，状態の変化に対して決定論的要因と確率論的要$)$

因が両存すると解釈される現象のデータに対して用いられる．しばしば，スケー

ルの大きな変動 (遅い変数) が決定論的寄与を，微細なゆらぎ (速い変数) が確率論的寄与を表すようなモデルが多く，Langevin 方程式がその最もシンプルでかつ重要な例である．確率微分方程式にょる解析というと，与えられた方程式を数学的に解析したり数値シミュレーションしたりすることでその現象の理解を深めるモデル
主体のアプローチもあるが，対象の時系列データが何かの確率微分方程式に従って時間発展したものであるとみなし，その係数をデータから推定するという逆問題的アプローチも存在する．すなわち．長さ $K$ の $N$ 次元の時系列データ$\{X_{t}=(X_{t}^{1}, X_{t}^{2}, \ldots, X_{t}^{N})^{T}\in \mathbb{R}^{N}|t\in \mathbb{T}\},$

は相空間$\mathbb{R}^{N}$$\mathbb{T}=\{t=j\triangle t|j=1, 2, .

.

.

, K\}$

上のマルコフ型確率微分方程式(1.1)$dX_{t}=A(X_{t})dt+\mathbb{S}(X_{t})dW_{t}$の解軌道からと$\mathbb{S}$$\triangle t$

おきにサンプルされた点であると仮定し，そのデータから係数 $A$を決定して解の挙動や統計的性質を調べるというものである．ここで， $x=$$(x^{1}, x^{2}, \ldots, x^{N})^{T}$

は相空間の位置とし， $A(x)=(A^{1}(x), A^{2}(x), \ldots, A^{N}(x))^{T}$ はドリフトベクトル， $W_{t}=(W_{t^{1}}, W_{t}^{2}, \ldots, W_{t}^{M})^{T}$ は独立な $M$ 個の Winner 過程を成分に持つ確率過程であり，
$1$

$T$

は転置を表す． $\mathbb{S}(x)=(S^{ij}(x))(i=1,2,$$\ldots,$

$N,$

$j=$

, 2, . . . , $M)$ は乗法的ノイズの係数行列であり，系にかかるランダム性の強さが状態に依存して変化しうる場合を考慮している．ここでいう系のランダム性とは，データの持つ測定誤差や Langevin 方程式のように確率論的ノイズとして入るランダム

外力だけを指すのではない．高次元力学を低次元相空間に射影したときは，縮約された高次の力学はその低次元空間での力学として表現できなくなるため，確率論的に扱われることになる．このようなシステム外の影響も含めたゆらぎを表す項として乗法的ノイズ項を採用するのである．
ここで，(1. 1) から導かれる Fokker-Planck 方程式は， $p(x, t)$ を (1.1) の解の遷移確率密度函数とすると，以下のように表される．
$\frac{\partial p}{\partial t}(x, t)=-\sum_{i=1}^{N}\frac{\partial}{\partial x^{i}}(A^{\dot{\iota}}(x)p(x, t))+\sum_{i,j=1}^{N}\frac{\partial^{2}}{\partial x^{i}\partial x^{j}}(B^{ij}(x)p(x, t))$

.

(1.2)

40
ここで

$\mathbb{B}$(x) $=(B^{ij}(x))$ は拡散行列とよばれ，である．確率論的には$\mathbb{S}(x)$ではなく$\mathbb{B}(x)$$\mathbb{B}=\frac{1}{2}\mathbb{S}\mathbb{S}^{T}$で表される非負対称行列が (1.1) の解の拡散をつかさどる行列であることが知られている．この Fokker-Planck 方程式 (1.2) の係数は，確率微分方
程式 (1.1) の解と以下のような関係があることが知られている :$A^{i}(x)= \lim_{tarrow 0}\frac{E_{x}[X_{t}^{i}-x^{i}]}{t},$ただし，E、は，$x$$B^{\dot{\iota}j}(x)= \lim_{tarrow 0}\frac{E_{x}[(X_{t}^{i}-x^{i})(X_{t}^{j}-x^{j})]}{2t}$

.

(1.3)

を初期値とする (1.1) の解が導く確率測度に対する期待値である．

(1.3) は，相空間の位置$x$を出発する (1.1) の解からその位置における (1.2) の係数の値を決定することができることを示している．

データ解析の観点では，この解と係数の関係式を利用して時系列データから統計
的に $A(x)$ と$\mathbb{B}(x)$の近似値を計算しようという研究手法が提案されてきた [1, 11].

すなわち，相空間をセルに分割し，各セルに滞在する時系列データを用いて$\{\begin{array}{l}\overline{A}^{i}(\overline{x}_{m})=\frac{\langle X_{t+\Delta t}^{i}-X_{t}^{i}\rangle}{\Delta t},\overline{B}^{ij}(\overline{x}_{m})=\frac{\langle(X_{t+\Delta t}^{i}-X_{t}^{i})(X_{t+\Delta t}^{j}-X_{t}^{j})\rangle}{2\Delta t}\end{array}$と計算するのである (Figure 1) . ただし，セル，$\overline{x}_{m}$(1.4)

は相空間を分割したときの $m$ 番目の$\overline{A}(\overline{x}_{m})=(A^{1}(\overline{x}_{m}),\overline{A}^{2}(\overline{x}_{m}), \ldots,\overline{A}^{N}(\overline{x}_{m}))^{T}$と$\overline{\mathbb{B}}(\overline{x}_{m})=(\overline{B}^{ij}(\overline{x}_{rn}))$はそれぞにおけるドリフトベクトルと拡散行列の計算値である．変数の上に付くバー C) は離散化されたセルやその上で定義される変数に対して用いることにする．
れセル$\overline{x}_{m}$また，は各セルに属する時系列データに対するアンサンブル平均を意味し，例
えばドリフトベクトルの計算では，セル$\overline{x}_{m}$の滞在時刻を$\mathbb{T}_{m}=\{t\in \mathbb{T}|X_{t}\in\overline{x}_{m}\}$とおき， $\{t\in \mathbb{T}_{m}|t<K\triangle t\}$ が空でない限り$\overline{A}^{i}(\overline{x}_{m})=\frac{\langle X_{t+\triangle t}^{i}-X_{t}^{i}\rangle}{\triangle t}=\frac{1}{\#\mathbb{T}_{m}-1}\sum_{t\in T_{m},t<K\Deltat}\frac{X_{t+\triangle t}^{i}-X_{t}^{i}}{\triangle t}$とおこなう．ここでは時系列データにエルゴード性を仮定し，データ長が十分あり
且つその軌道がセル$\overline{x}_{m}$を十分多く通過するのであれば，そのセルに含まれるデー

タに対するアンサンブルを計算することで， $(1.3)$ での期待値の計算の代わりにすることができるとしている．この手法を用いれば時系列データだけから確率微分方
程式モデルを構成することができるわけである．
一方で，(1.4) は (1.3) から自然に導かれる時空離散化を用いており，統計的係数決定公式として理に適うように思えるのだが，実はドリフトベクトルのノルムが41大きい領域で拡散行列の計算誤差が大きくなる傾向となることが今回ゎかった．特に，拡散行列のノルムは系の本質的なノイズの強さを表す量であるが，(1.4) から計算される量ではノルムが必ず過大評価となってしまう．これにより，データから評価した相空間の各位置におけるノイズの強さが実際の相空間における時系列の軌道の振る舞いとは直感的に異なってしまう例があり (\S 3 参照), 実用上好ましくない．また，ドリフトベクトルの計算時にもセルのサイズとデータサンプル幅に制限があることもわかっており (\S 2.1 参照), 時系列データから確率微分方程式の係数を算出するときには幾つか払うべき注意点があることがゎかる．本研究では，それらの注意点を数学的に示し，確率微分方程式モデリングにおける新しい係数決定公式を提示することが目的である．

Figure 1: 相空間分割セルの概念図

本稿の構成は以下の通りである．まず\S 2.1 ではドリフトベクトルの計算手法に対する注意点を例を挙げながら説明し，その問題点を除去した公式を紹介する(2.1). \S 2.2 では拡散行列に対する新しい係数決定公式 (2.2) を示し，その有効性を

簡単な例を通して議論する．最後に \S 3 において Lorenz system と R\"ossler systemの数値解の平面射影データに対して，旧公式と新公式の結果を比較することで新公式の有効性を数値的に検証する．

2

統計的係数決定公式

2.1

ドリフトベクトルの計算式に関する注意点

ドリフトベクトルの計算においては，時間スケ–, すなゎちデータのサンプル間隔$\triangleright$

ム$t$と空間スケールすなわちセルのサイズ$\triangle$

x のバランスが大切になる．すなゎち，42
平均移動距離より比較的大きいビンサイズを取るときに (1.1) は直感と異なる結果を与えることがある．また，ドリフトの計算式から対象の時系列にふさわしいムオと$\triangle x$の取り方がわかる．これを順を追って説明する．

まず用語を定義する．時系列データがあるセルに入ってから最初に出るまでをそのセルにおける “一連のイベンド’ とよぶことにする．セル$X_{t+k\triangle t}$

かつで出た場合，すなわち$X_{t+k\triangle t}\not\in\overline{x}_{m}$$X_{t-\triangle t}\not\in\overline{x}_{m}$$\overline{x}_{m}$に $X_{t}$ で入ってかつ $X_{t+j\Delta t}\in\overline{x}_{m}(j=0,1, \ldots, k-1)$のとき，その一連のイベントは長さが であるという．
$k$

Example 2.1
二次元相空間 $(N=2)$ 上のセル $\overline{x}_{m}=[a, a+\Delta x$ )$\cross[b, b+\triangle x$$X_{t_{1}+j\Delta t}=(a+ \frac{j\triangle x}{4}, b+\frac{j\triangle x}{8})(j=0,1, \ldots, 4)$) の近傍を,$X_{t_{2}+j\Delta t}=(a+ \frac{j\triangle x}{4}, b+\frac{(j+4)\triangle x}{8})(j=0,1, \ldots, 4)$というように長さ 4 の 2 つの一連のイベントが$\overline{x}_{m}$$\overline{x}_{m}$を通過するとする．ここでセルに滞在したデータ点は合計 8 点で $\mathbb{T}_{m}=\{t_{1}+j\triangle l, t_{2}+j\triangle t(i=0,1,2,3)\}$ であることに注意する．このときドリフトベクトルは$\frac{1}{8}\sum_{t\in \mathbb{T}}\frac{1}{8}\frac{1}{8}\frac{1}{8}\frac{1}{8}\frac{1}{8}t\in T(\frac{}{}\frac{x_{t_{2}+4\triangle t}^{1}-x_{t_{2}}^{j--0}x_{t_{1}+j\Delta t_{+\sum_{1}^{3}}}^{1}}{\Delta t})=\frac{\triangle x}{4\triangle t’}\sum^{(}\frac{}{\frac{X}{}+}\frac{X_{t_{2}+(j+1)\triangle t}^{1}-X_{t_{2}+j\triangle t}^{1}}{\triangle t})(\frac{\sum X_{t_{1}^{\frac{X}{}}}^{1}\sum_{1}^{3}j_{m}^{m}X_{t_{1}}^{1}j_{-}^{-0}=03}{}\frac{X_{t_{1}+j\Delta t_{+\sum_{x_{t_{2}+4\triangle t}^{1}-x_{t_{2}^{1}}^{j--0}}^{3}}}^{2}}{\triangle t})=\frac{\triangle x}{8\triangle t}(\frac{X_{t_{1}+(j+1)\triangle t}^{1}-X_{t_{1}+(j+1)\Delta t}^{2}-t+\Delta t_{t}t4\Delta t_{\triangle t}t+\Delta t12\triangle\iota_{-X_{t}^{2^{+}}}^{-X_{t_{1}}^{1^{\triangle t}}}\triangle^{-X^{1}}}{+4\Delta tt_{1},\triangle\iota^{-x_{+}}\triangle t}\frac{X_{t_{2}+(j+1)\triangle t}^{2}-X_{t_{1}+j\triangle t}^{2}}{\triangle t})$$\ovalbox{\tt\small REJECT} A^{1}(\overline{x})A^{2}(\overline{x})$と求められる．ドリフトの計算式では和の線形性のために，単位時間$(\Delta t)$

あたりに進む速度 (すなわち一連のイベントの長さの逆数) は結果に関係せず，一連のイベントにおける最初の点から最後の点までの各成分の (符合込みの) 直線距離だ

けが計算に利いていることがわかる．この例のように一つのセル内をほぼ平行にほ43
ぼ同速度で軌道が通過するような決定論的な時系列データの場合では特に問題は起きない．

Example 2.2

次に，同じセル内で軌道が交差するとき，例えば$X_{t_{1}+j\triangle t}=(a+ \frac{j\triangle x}{4}, b+\frac{\triangle x}{2})(j=0,1, \ldots, 4)$,$X_{t_{2}+j\triangle t}=(a+ \frac{\triangle x}{2}, b+\frac{j\triangle x}{8})(j=0,1, \ldots, 8)$というような 2 つの一連のイベントが$j\triangle t(j=0,1,2,3)$, $t_{2}+k\triangle t(k=0,1,$$\overline{x}_{m}$を通過するときを考える $(\mathbb{T}_{m}=\{t_{1}+$
この場合では，時亥」$7$$\ovalbox{\tt\small REJECT}$$\ldots,$に平行に通過するときの一連のイベントの方が$t_{2}$

から$y$$t_{1}$

から

$x$

軸

軸に平行に通過するときより 2 倍速い．一方で，ドリフトの計算結果は$\ovalbox{\tt\small REJECT} A^{2}(\overline{x})A^{1}(\overline{x})$$\frac{1}{12}\frac{}{}\frac{1}{12}\frac{x_{22}221}{}+\sum_{j=0}^{\frac{1}{12}}\frac{X_{t_{2}+(j+1)\triangle t}^{2}-X_{t_{1}+j\triangle t}^{2}}{\triangle t})\frac{4\triangle tX}{}\frac{1}{12}\frac{X_{t_{1+}}^{1}t\in \mathbb{T}_{m}(\sum_{j=0}^{3}\sum}{X_{t_{2}+8\triangle t}^{1^{=0}}t\in \mathbb{T}_{m}(_{j}\sum_{1}^{3}\sum}\frac{1}{12}\frac{\frac{}{t}X_{t_{1}+(j+1)\triangle t}^{1}-X_{t_{1}+j\Delta t}^{1}t+\triangle t1\triangle^{-X_{t}^{1}}}{\triangle\iota^{-X_{t_{2}}^{-x_{\triangle}}}\triangle^{-X_{t_{1}}}X_{t_{1}+(j+1)\triangle t}-X_{t_{1}+j\triangle t}t+\triangle ttt^{=\frac{\triangle x}{12\triangle t’}}\triangle t=\frac{t_{\triangle x}}{12\triangle t}\triangle t}+\sum^{7}\frac{X_{t_{2}+(j+1)\triangle t}^{1}-X_{t_{2}+j\triangle t}^{1}}{\triangle t})\frac{1}{12}j=07$となり，ドリフトの

$x$

成分も

$y$

成分も同じ値という実際の時系列データから受ける直感とは異なる結果が出される．やはりセル内で進んだ距離のみが計算式に寄与するため，元の時系列データの速さは計算には加味されないことがわかる．このような同一のセル内での軌道交差は，本研究で対象とするような確率論的ノイズにょって摂動を受ける軌道や高次元力学の軌道を低次元相空間に射影した場合で起こり得る．したがって，時系列データの平均移動距離より大きいセルサイズを用いるときは適切にドリフトが計算されているか注意が必要となる．
Example 2.3

44
また，このドリフトの計算ではセルサイズ x より小さい変動は算定に加味してい$\triangle$

ない．例えば，$X_{t+j\triangle t}=(a+ \frac{j\triangle x}{4}, b+\frac{(2-|j-2|)\triangle x}{4})(j=0,1, \ldots, 4)$という同じセル内で軌道が大きく方向を変えるような一連のイベントを考える．このときは，$\{\begin{array}{l}A^{1}(\overline{x})=\frac{\Delta x}{4\triangle t}A^{2}(\overline{x})=0\end{array}$となり，ドリフトの計算値はセルのサイズが最小の解像度となっていることがわかる．すなわち得られる確率微分方程式モデルはそのスケールでの粗視化がされているものであり，したがって$\triangle x$

としては着目する現象の最小スケールを取るべきである．

とはいえ，ムオ当たりの平均移動距離は，実験や測定の都合上 (もしくは別の目
的があって) 着目する現象の最小スケールであるとは限らず，比較的小さい場合がある．この場合は前述の Example 2.2 で挙げた問題点が出てくるため，それに合わせて計算手法を工夫する必要がある．ここでは 3 つの手法を挙げることにする．
1. 時系列データの平均移動距離がセルサイズと同等になるくらいにデータを等間隔に間引く．すなわち，時系列データを $\mathbb{T}’=\{t=\tau j\triangle t|j=1, 2, . . . , [K/\tau]\}$と$\tau$

毎にサブサンプルするのである．同一セルに 2 点連続して滞在しないよ

うに出来れば Example 2.2 のような問題は起こらない．また，Example 2.3の空間解像度の限界からデータの間引きによる時間解像度の劣化は評価には大きく影響しない．ただし，各セルにおけるアンサンブル点はほぼ$1/\tau$となるので，この手法はデータ数が十分ある場合に有効である．
2. ドリフトベクトルの計算式において，次のステップのデータとの差を計算すと同等になるくらいの時間サンプ
るのではなく，平均移動距離がやはり$\triangle x$

ル幅を用いて計算する．すなわち，を取り，$\{t\in \mathbb{T}_{m}|t+\tau\triangle t\leq K\triangle t\}$$\langle X_{t+\tau\triangle t}^{i}-X_{t}^{i}\rangle\approx\triangle x$となるような$\tau\in \mathbb{N}$が空でない限り(2.1)$\overline{A}^{i}(\overline{x}_{m})=\frac{\langle X_{t+\tau\Delta t}^{i}-X_{t}^{i}\rangle}{\tau\triangle t}=\frac{1}{\#\mathbb{T}_{m}-\tau}t+\tau\triangle t\leq K\triangle t\sum_{t\in \mathbb{T}_{m}}\frac{X_{t+\tau\Delta t}^{i}-X_{t}^{i}}{\tau\triangle t}$45とおこなう．データは間引かずに差の計算だけ$\tau$

ステップ先のデータとおこ

なうため，各セルにおけるアンサンブル数は時間サンプル幅の変更前とほぼ変わらない点が上記 1 とは異なる．この手法でも Example 2.2 のような影響を減らすことができる．

3. また，各一連のイベントの平均に対するアンサンブル平均を用いて計算することもできる．すなわち$\overline{x}_{m}$に滞在するデータ点が$X_{t_{j}+k\triangle t}(j=1,2, \ldots, J, k=0,1, \ldots, k_{j}-1)$という長さがそれぞれ $k_{j}(i=1,2, \ldots, J)$ である$J$

本の一連のイベントで成り立っているとき，$\overline{A}^{i}(X_{m})=\frac{1}{J}\sum_{j=1}^{J}\frac{1}{k_{j}}\sum_{k=0}^{k_{j}-1}\frac{X_{t_{j}+(k+1)\Delta t}^{i}-X_{t_{j}+k\triangle t}^{i}}{\Delta t}$

によってドリフトを求めるのである．これも Example 2.2 のような問題点は解消されている．上記 2 つの手法のような時系列の時間間隔を一定に保つ計算法とは異なって各イベントに着目するため，相空間の場所によって平均移動距離が大きく変化する場合にも使える．

いずれにせよ，平均移動速度$\approx\triangle x/\triangle t$となるように時系列のサンプル間隔とセルサイズを取るのが理想である．

しばしば時系列データの時間解像度は対象の実験現象や測定手法に依存するものであり，こちらが自由に設定できるとは限らない．仮に時間解像度に合わせて
空間解像度を取るとすると，着目したい空間スケールより小さいセルサイズになる
場合があり，そうすると解析の目的とは合致しなくなる．また，気象データなど繰
り返し観測を行うことができない場合では計算に必要なだけのデータ長を必ずしも取れないため，各セルにおけるアンサンブルメンバーを多く確保するために粗い
セルサイズで相空間を分割せねばならないこともある．このようにデータセットの性質に由来するセルサイズの問題は様々であるが，それらを打破してより適切なドリフトベクトルを計算するには上記の 2 の手法が良いと思う．セルサイズは着目す
る空間解像度を採用し，平均移動速度が$\triangle x/\triangle t$と釣り合うように$\tau$

を採用すれば，

Example 2.2 などの影響は小さくなり，かつアンサンブルメンバーの数も減らさず
に済む．

46
拡散行列の新計算公式

2.2(1.4) による係数の計算公式ではドリフトの強い領域で誤差が大きくなることがわ

かっている．特に誤差が顕著に現れる拡散行列に対してその評価式からドリフトの影響を除去するには，$\overline{B}_{new}^{ij}(\overline{x}_{m})=\frac{\langle(X_{t+\tau\Delta t}^{i}-X_{t}^{i})(X_{t+\tau\Delta t}^{j}-X_{t}^{j})\rangle-\langle X_{t+\tau\Delta t}^{i}-X_{t}^{i}\rangle\langle X_{t+\tau\Delta t}^{j}-X_{t}^{j}\rangle}{2\tau\triangle t}$(2.2)と，二次変動から一次変動に関する二次式を除去する計算式を用いるのが良い．ただし，$\tau$

は，\S 2.1 で注意した，$|X_{t+\tau\triangle t}-X_{t}|\approx\triangle x/\triangle t$となるような$\tau$

である．これの有効性を具体的な確率微分方程式の解を例にとって説明する．

Example 2.4
幾何ブラウン運動(2.3)$\{\begin{array}{l}dX_{t}=aX_{t}dt+bX_{t}dW_{t},X_{0}=x\end{array}$の場合に，ドリフトと拡散係数の近似式として，それぞれ(2.4)$\{\begin{array}{l}\frac{E_{x}[X_{t}-x]}{t},\frac{E_{x}[(X_{t}-x)^{2}]-E_{x}[X_{t}-x]^{2}}{2t}\end{array}$を用いることが有効であることを示そう．ここで $a,$
し，

$W_{t}$$b\in \mathbb{R}$を非ランダムの定数とを Wiener 過程とする．(2.3) のドリフト係数は $A(x)=ax$ で拡散係数は$B(x)=(bx)^{2}/2$ である．この解は変数分離や伊藤の公式を用いることで計算することができ，$X_{t}=xe^{(a-\frac{b^{2}}{2})t+bW_{t}}$で与えられる．
ここで， $E_{x}[\exp(bW_{t})]=\exp(b^{2}t/2)$ より，解の一次変動は$\frac{E_{x}[X_{t}-x]}{t}=\frac{x(e^{at}-1)}{t}=ax+\mathcal{O}(t)(tarrow 0)$となる．誤差項$\mathcal{O}(t)$(2.5)には $x$ も含まれているため， $Eoe[X_{t}-x]/t$ はドリフト $ax$ との誤差が原点から遠いほど大きくなる．しかし実際は$\mathcal{O}(t)=ax\sum_{n=1}^{\infty}\frac{(at)^{n}}{(n+1)!}$47であるため，ドリフト $ax$ との相対誤差は相空間の位置にょらず一様なので問題ない．
一方，二次変動を計算すると$\frac{E_{x}[(X_{t}-x)^{2}]}{2t}=\frac{x^{2}}{2t}E_{x}[e^{(2a-b^{2})t+2bW_{t}}-2e^{(a-\frac{b^{2}}{2})t+bW_{t}}+1]$$= \frac{x^{2}\{e^{(2a+b^{2})t}-2e^{at}+1\}}{2t}=\frac{x^{2}e^{2at}(e^{b^{2}t}-1)+E_{x}[X_{t}-x]^{2}}{2t}$$= \frac{(bx)^{2}}{2}+\frac{t}{2}(\frac{E_{x}[X_{t}-x]}{t})^{2}+\mathcal{O}(t)(tarrow 0)$$= \frac{(bx)^{2}}{2}+\frac{A(x)^{2}}{2}t+\mathcal{O}(t)(tarrow 0)$を得る．(2.6) の第 2 項は$\mathcal{O}(t)$(2.6)に含まれる項ではあるが， $t>0$ である限りは拡散係数の計算式に一次変動の二乗が陽に現れた形となり，拡散係数の値
より
概ね過剰評価となることがわかる．また，その第 2 項はドリフトの大きい領域で誤$(bx)^{2}/2$差への寄与が強くなるため，誤差が の大きさだけではコントロールされずに相空$t$間の場所によって変わるというのは好ましくない．したがって，二次変動から一次変動の二乗を引いておくと，$\frac{E_{x}[(X_{t}-x)^{2}]-E_{x}[X_{t}-x]^{2}}{2t}=\frac{x^{2}e^{2at}(e^{b^{2}t}-1)}{2l}=\frac{(bx)^{2}}{2}+\mathcal{O}(t)(tarrow 0)$(2.7)となり，その誤差項はドリフトの大きさの多寡にょらず，空間一様に $t>0$ の大きさに支配されるものになっている．
以上により，(2.4) がドリフトと拡散係数の近似として有効であることがゎかった．

注童この例からわかるように，(2.5) や (2.7) という $tarrow 0$ の極限を取らない量は，
相対誤差は空間一様に取れるものの，概ねドリフトの強さのパラメータ にょって
絶対誤差・相対誤差の多寡が左右される傾向があることがゎかる．
$a$

Example 2.5$N$

次元の非斉次線形確率微分方程式$\{\begin{array}{l}dX_{t}=(A_{0}+\mathbb{L}X_{t})dt+(\mathbb{S}_{0}+\sigma X_{t})dW_{t},X_{0}=x\end{array}$(2.8)

について，ドリフトベクトルと拡散行列の各成分の近似式として，それぞれ$\{\begin{array}{l}\frac{E_{x}[X_{t}^{i}-x^{i}]}{t},\frac{E_{x}[(X_{t}^{i}-x^{i})(X_{t}^{j}-x^{j})]-E_{x}[X_{t}^{i}-x^{i}]E_{x}[X_{t}^{j}-x^{j}]}{2t}\end{array}$

(2.9)

48
を用いることが有効であることを示す．ここで，の定ベクトル，$\sigma=(\sigma^{ijk})$$\mathbb{L}=(L^{ij})$,$\mathbb{S}_{0}=(S_{0}^{ij})$$1$はそれぞれサイズ $N\cross N,$

は $N\cross M\cross N$ の三階の定テンソルであり，ムである．また，テンソルの縮約はは $N$ 次元

$A_{0}=(A_{0}^{1},A_{0}^{2}, \ldots,A_{0}^{N})^{T}$$A_{0},\mathbb{L},\mathbb{S}_{0},$$N\cross M$$\sigma$の定行列，は全て非ランダ$[ \sigma X_{t}]^{ij}=\sum_{k=1}^{N}\sigma^{ijk}X_{t}^{k}(i=1,2,$$\ldots,$

$N,$

$j=$

, 2, . .. , $M)$ とする．(2.8) の解の陰的公式は 1 階の非斉次常微分方程式系と同様に計算することで$X_{t}=D(t)+Z_{t},$$D(t)= e^{Lt}x+\int_{0}^{t}e^{L(t-\mathcal{S})}A_{0}ds, Z_{t}=\int_{0}^{t}e^{L(t-s)}\mathbb{S}(X_{s})dW_{s},$$\mathbb{S}(X_{t})=\mathbb{S}_{0}+\sigma X_{t}$を得る．この解は決定論的な項 $D(t)$ と確率論的な項 Z$\mathcal{O}$

和として与えられる．こ

こで$\int_{0}^{t}e^{-\mathbb{L}s}\mathbb{S}(X_{s})dW_{s}$はマルチンゲールなので，$E_{x}[Z_{t}]=e^{Lt}E_{x}[\int_{0}^{t}e^{-L_{8}}\mathbb{S}(X_{s})dW_{s}]=0$となり，$\frac{E_{x}[X_{t}-x]}{t}=\frac{D(t)-x}{t}=\frac{(e^{Lt}-\mathbb{I})x}{t}+\frac{1}{t}\int_{0}^{t}e^{L(t-s)}A_{0}ds$$=A_{0}+\mathbb{L}x+\mathcal{O}(t)(tarrow 0)$が従う．この場合でも絶対誤差は$x$に依存するが，$\frac{|E_{x}[X_{t}-x]/t-(A_{0}+\mathbb{L}x)|}{|A_{0}+\mathbb{L}x|}\leq c(A_{0}, \mathbb{L}, t)$となるため，相対誤差は Example 2.4 と同様に空間一様に抑えられる．
また，$(X_{t}^{i}-x^{i})(X_{t}^{j}-x^{j})=(D^{i}(t)-x^{i})(D^{j}(t)-x^{j})+(D^{i}(t)-x^{i})Z_{t}^{j}$$+Z_{t}^{i}(D^{j}(t)-x^{i})+Z_{t}^{i}Z_{t}^{j}$であるため，$\frac{E_{x}[(X_{t}^{i}-x^{i})(X_{t}^{j}-x^{j})]}{2t}=\frac{(D^{i}(t)-x^{i})(D^{j}(t)-x^{j})+E_{x}[Z_{t}^{i}Z_{t}^{j}]}{2t}$(2.10)

49$= \frac{E_{x}[X_{t}^{i}-x^{i}]E_{x}[X_{t}^{j}-x^{j}]+E_{x}[Z_{t}^{i}Z_{t}^{j}]}{2t}$$= \frac{t}{2}\cdot\frac{E_{x}[X_{t}^{i}-x^{i}]}{t}\cdot\frac{E_{x}[X_{t}^{j}-x^{j}]}{t}+\frac{E_{x}[Z_{t}^{i}Z_{t}^{j}]}{2t}(2.11)$が成り立つ．ここで伊藤の公式を用いると，右辺第 2 項はオ
$E_{x}[Z_{t}^{i}Z_{t}^{j}]=\int_{0}$$E_{x}[2B^{ij}(X_{s})]ds$オ$ds\int_{0}^{s}[2\mathbb{L}e^{L(s-\tau)}E_{x}[\mathbb{B}(X_{\tau})]e^{\mathbb{L}^{T}(s-\tau)}]^{ij}d\tau$$+ \int_{0}$オ$+ \int_{0}$$ds\int_{0}^{s}[2\mathbb{L}e^{L(s-\tau)}E_{x}[\mathbb{B}(X_{\mathcal{T}})]e^{\mathbb{L}^{T}(s-\tau)}]^{ji}d\tau$$=2tB^{ij}(x)+ \int_{0}^{t}2E_{x}[B^{ij}(X_{s})-B^{ij}(x)]ds$オ$d_{\mathcal{S}}\int_{0}^{s}[2\mathbb{L}e^{\mathbb{L}(s-\tau)}E_{x}[\mathbb{B}(X_{\tau})]e^{L^{T}(s-\tau)}]^{ij}d\tau$$+ \int_{0}$オ$+ \int_{0}$$ds\int_{0}^{s}[2\mathbb{L}e^{\mathbb{L}(s-\tau)}E_{x}[\mathbb{B}(X_{\tau})]e^{L^{T}(s-\tau)}]^{ji}d\tau$$=2tB^{ij}(x)+\mathcal{O}(t^{2}) (tarrow 0)$(2.12)となる．ただし，$B^{ij}(x)= \frac{1}{2}\sum_{l=1}^{M}S^{il}(x)S^{jl}(x)=\frac{1}{2}\sum_{l=1}^{M}(\sum_{m=1}^{N}S_{0}^{il}+\sigma^{ilm}x^{m})(\sum_{n=1}^{N}S_{0}^{jl}+\sigma^{jln}x^{n})$である．したがって解の二次変動量 (2.11) は，決定論的成分 (右辺第 1 項) と確率論的成分 (右辺第 2 項) の和で表されることがゎかる．だからこの場合も，二次変動量から一次変動量の二次式を引いておくと良い :$\frac{E_{x}[(X_{t}^{i}-x^{i})(X_{t}^{j}-x^{j})]-E_{x}[X_{t}^{i}-x^{i}]E_{x}[X_{t}^{j}-x^{j}]}{2t}=B^{ij}(x)+\mathcal{O}(t)$(2.13)

さらに $E_{x}[X_{t}^{i}]=D^{i}(t)$ と $E_{x}[X_{t}^{i}X_{t}^{j}]=D^{i}(t)D^{j}(t)+E_{x}[Z_{t}^{i}Z_{t}^{j}]$ より［罪罪］ Fij(t),$E_{x}[B^{ij}(X_{t})]=\sum_{l=1}^{M}\sum_{m,n=1}^{N}\sigma^{ilm}\sigma^{jln}E_{x}$$+$$F^{ij}(t)= \sum_{l=1m}^{M}\sum_{n=1}^{N}\sigma^{ilm}\sigma^{jln}D^{m}(t)D^{n}(t)+\sum_{l=1m}^{M}\sum_{=1}^{N}(\sigma^{ilm}S_{0}^{jl}+\sigma^{jlm}S_{0}^{il})D^{m}(t)+B_{0}^{ij}$が従う．これと (2.12) より，$E_{x}[Z_{t}^{i}Z_{t}^{j}]$は$E_{x}[Z_{t}^{i}Z_{t}^{j}]=\int_{0}^{t}\sum_{l=1}^{M}\sum_{m,n=1}^{N}2\sigma^{ilm}\sigma^{jln}E_{x}[Z_{s}^{m}Z_{s}^{n}]ds$

50$+ \int_{0}^{t}ds\int_{0}^{s}\sum_{\iota=1}^{M}\sum_{m,n,p,q=1}^{N}2[\mathbb{L}e^{L(s-\tau)}]^{ip}\sigma^{plm}\sigma^{qln}E_{x}[Z_{\tau}^{m}Z_{\tau}^{n}][e^{L^{T}(s-\tau)}]^{qj}d\tau$$+ \int_{0}^{t}ds\int_{0}^{s}\sum_{\iota=1}^{M}\sum_{m,n,p,q=1}^{N}2[\mathbb{L}e^{L(s-\tau)}]^{jp}\sigma^{plm}\sigma^{qln}E_{x}[Z_{\tau}^{m}Z_{\tau}^{n}][e^{L^{T}(s-\tau)}]^{qi}d\tau$$+G^{ij}(t)$,$G^{ij}(t)= \int_{0}^{t}F^{ij}(s)ds+\int_{0}^{t}ds\int_{0}^{s}[2\mathbb{L}e^{L(s-\tau)}\mathbb{F}(\tau)e^{\mathbb{L}^{T}(s-\tau)}]^{ij}d\tau$$+ \int_{0}^{t}ds\int_{0}^{s}[2\mathbb{L}e^{L(s-\tau)}\mathbb{F}(\tau)e^{L^{T}(s-\tau)}]^{ji}d\tau$という二階の非斉次線形積分方程式を満たすことがわかり，解の評価から$,| \frac{E_{x}[Z_{t}^{i}Z_{t}^{j}]/2t-B^{ij}(x)}{B^{ij}(x)}|\leq c(A_{0}, \mathbb{L}, \mathbb{S}_{0}, \sigma, t)$Nmax.

$i,j=1,2,..$が従う．すなわち，$E_{x}[Z_{t}^{i}Z_{t}^{j}]/2t$の $B^{ij}(x)$ に対する相対誤差は相空間において一様に抑えられる．

以上により，(2.9) がドリフトベクトルと拡散行列の近似式として有効であることがわかった．

Example 2.6

最後に，一般のマルコフ型の確率微分方程式(2.14)$\{\begin{array}{l}dX_{t}=A(X_{t})dt+\mathbb{S}(X_{t})dW_{t},X_{0}=x\end{array}$の場合を考えよう．このときも (2.9) がドリフトベクトルと拡散行列の近似式として有効であることを示す．
解と係数の関係 (1.3) は一般の確率微分方程式に対して成り立ち，その等式は$tarrow 0$への極限を取ることで得られる．したがって，係数の近似計算では解は初期時刻付近だけを考えれば良く， $t\ll 1$ のときに限定してよい． $A(x)$ と$\mathbb{S}(x)$に$C^{1}$

級程度の十分な滑らかさがあれば， を十分小さく取ることで $|X_{t}-x|\ll 1$ も同$t$様に成立するようにできる．このとき，$\{\begin{array}{l}A^{i}(X_{t})\approx A^{i}(x)+\sum_{k=1}^{N}\frac{\partial A^{i}}{\partial x^{k}}(x)(X_{t}^{k}-x^{k}) ,S^{ij}(X_{t})\approx S^{ij}(x)+\sum_{k=1}^{N}\frac{\partial S^{ij}}{\partial x^{k}}(x)(X_{t}^{k}-x^{k})\end{array}$51と係数を一次近似すると，(2.14) は $\hat{Y}_{t}=X_{t}-x$ に対する以下の方程式に線形化される．

$\{\begin{array}{l}d\hat{Y}_{t}^{i}=(A^{i}(x)+\sum_{k=1}^{N}\frac{\partial A^{i}}{\partial x^{k}}(x)\hat{Y}_{t}^{k})dt+\sum_{j=1}^{M}(S^{ij}(x)+\sum_{k=1}^{N}\frac{\partial S^{ij}}{\partial x^{k}}(x)\hat{Y}_{t}^{k})dW_{t}^{j},\hat{Y}_{0}^{i}=0 (i=1,2, \ldots, N)\end{array}$(2.15)
これは$\hat{Y}_{t}$に関する非斉次線形確率微分方程式であるため，係数の近似計算は$A(x)$ , $\mathbb{L}=\nabla A(x)$ ,$X_{t}$と$\hat{Y}_{t}$$\mathbb{S}_{0}=\mathbb{S}(x)$,$\sigma=\nabla \mathbb{S}(x)$$A_{0}=$として Example 2.5 に帰着される．

をそれぞれ (2.14) と (2.15) の解とすると，(2.10) から$\frac{E_{x}[X_{t}^{i}-x^{i}]}{t}\approx\frac{\hat{E}_{0}[\hat{Y}_{t}^{i}]}{t}=A^{i}(x)+\sum_{k=1}^{N}\frac{\partial A^{i}}{\partial x^{k}}(x)\cdot 0+\mathcal{O}(t)=A^{i}(x)+\mathcal{O}(t)(tarrow 0)$が近似的に成り立っ．ただし，$\’{E}_{0}$

は，初期値

$0$

の (2.15) の解が導く確率測度についての期待値である．(2.15) の初期条件により，上式における消えることに注意する．同様に (2.13) から$\hat{Y}_{t}^{k}$の一次の項は$\frac{E_{x}[(X_{t}^{i}-x^{i})(X_{t}^{j}-x^{j})]-E_{x}[X_{t}^{i}-x^{i}]E_{x}[X_{t}^{j}-x^{j}]}{2t}$$\approx\frac{\hat{E}_{0}[\hat{Y}_{t}^{i}\hat{Y}_{t}^{j}]-\hat{E}_{0}[\hat{Y}_{t}^{i}]\^{E}_{0}[\hat{Y}_{t}^{j}]}{2t}$$= \frac{1}{2}\sum_{l=1}^{M}(S^{il}(x)+\sum_{k=1}^{N}\frac{\partial S^{il}}{\partial x^{k}}(x)\cdot 0)(S^{jl}(x)+\sum_{k=1}^{N}\frac{\partial S^{jl}}{\partial x^{k}}(x)\cdot 0)+\mathcal{O}(t)$$=B^{ij}(x)+\mathcal{O}(t)(tarrow 0)$が近似的に成り立っ．したがって，一般のマルコフ型確率微分方程式 (2.14) に対しても， $tarrow 0$ の極限を取ることで得られる確率微分方程式の係数の公式 (1.3) に対する $t>0$ における近似式としては (2.9) が有効である．

以上の議論により，時系列データからその公式を用いて係数を定量的に評価する
ときは，(1.4)2 ではなく (2.2) を用いるのが適切であることがゎかる．ただ，(2.2) は$\frac{\langle(X_{t+\tau\triangle t}^{i}-X_{t}^{i})(X_{t+\tau\triangle t}^{j}-X_{t}^{j})\rangle-(\langle X_{t+\tau\triangle t}^{i}\rangle-\langle X_{t}^{i}\rangle)(\langle X_{t+\tau\triangle t}^{j}\rangle-\langle X_{t}^{j}\rangle)}{2\tau\triangle t}$と表されるため，$X_{t+\tau\Delta t}$と$X_{t}$

が大概同じセルに滞在するような場合，すなゎち

相空間分割セルのサイズと比べて時系列データの移動距離が小さい場合は，既存の計算方法を用いても結果は変わらないというのは事実である．しかし，\S 2. 1$(1.4)_{2}$52でのドリフトベクトルの計算時の注意から，$X_{t+\tau\Delta t}$が$X_{t}$とは異なるセルに属すを取るようにしたのであった．したがって，そのような時系列の取り方をする場合は，適切にドリフトの影響を除去しないと相対誤差が一様に抑えられなるように$\tau$いなどの弊害が生じてくるため，拡散行列の定量的評価としては総じて (2.2) がより適切な計算手法となる．

2.3

拡散行列の統計的計算公式に対する考察

確率微分方程式の本質的なランダム性をつかさどるのが拡散行列であるため，そのノルムを計算すると相空間の中でランダム性の強い部分を調べることができる．非負定値行列$\mathbb{B}$に対しては$\frac{1}{\sqrt{N}}tr\mathbb{B}\leq\Vert \mathbb{B}\Vert\leq tr\mathbb{B}$となるため，ここでは行列ノルムの替わりにとにする．ただし$\Vert B\Vert$$tr\mathbb{B}$を用いてその大きさを評価するこはフロベニウスノルムである．すると，(2.2) の統計的係数

決定公式では，これに対応する量は$tr\overline{\mathbb{B}}_{new}(\overline{x}_{m})=\frac{\langle|X_{t+\tau\triangle t}-X_{t}|^{2}\rangle-|\langle X_{t+\Delta t}-X_{t}\rangle|^{2}}{2\tau\triangle t}$

であり，これはセル$\overline{X}_{m}$から出発する軌道の分散の時間変化度を表している．すな

わちこの量は，連続形として考えると確率微分方程式の拡散行列のノルムの離散近似に対応し，離散形として考えると分割セルを通る軌道の分散の平均変化率が対応しているのである．また，$\overline{\mathbb{B}}_{new}(\overline{x}_{m})$自体は分割セルを通る軌道の共分散行列の時間変化度であり，連続形における拡散行列が確率微分方程式の各成分のノイズのか
かり方の強弱の相関を表すことに対応している．したがって，$\overline{\mathbb{B}}_{new}(\overline{X}_{m})$は連続形と離散形のどちらの視点でも系のランダム性を表現するのに適切であり，この公式によって得られる計算値に対してきちんとした意味付けを与えることができる．

この確率微分方程式モデリングは両者をつなぐ関係性を見出すことができるという意
味でも重要である．
$B^{ij}(x)= \lim_{tarrow 0}\frac{E_{x}[(X_{t}^{i}-x^{i})(X_{t}^{j}-x^{j})]-E_{x}[X_{t}^{i}-x^{i}]E_{x}[X_{t}^{j}-x^{j}]}{2t}$$\approx\frac{\langle(X_{t+\tau\triangle t}^{i}-X_{t}^{i})(X_{t+\tau\Delta t}^{j}-X_{t}^{j})\rangle-\langle X_{t+\tau\Delta t}^{i}-X_{t}^{i}\rangle\langle X_{t+\tau\triangle t}^{j}-X_{t}^{j}\rangle}{2\tau\triangle t}=\overline{B}_{new}^{ij}(\overline{x}_{m})$$trB$は軌道の分散の平均変化率を表すが，評価時間となる

れだけ軌道は分散してしまう．このため， $trB$ を$|A|^{2}$$\tau\Delta t$.

が大きければそで割っておくと，$\tau\triangle t$

の間の

53
変動に対するノイズの相対的寄与を得ることができる．この量はいゎゆる Wayland
test [12] における並進誤差に対応する．Wayland test はスヵラー時系列を遅延座標系を用いて高次元空間に埋め込み，近接軌道群の方向分散を計算して力学の次元の間接推定をする手法であり，短いデータ長の時系列データでも有効であるのが利点である．Wayland test はカオス検定に用いられる手法で用途が異なるものの，

近接軌道群の方向分散を求めるという点では本研究の趣旨に近いが，並進誤差の計
算では方向分散を計算する点の選択性の影響を排除するために幾つかの点における方向分散の中間値や平均値を取るため，時系列データセット 1 つに対して決まる量である．確率微分方程式モデリングは，相空間の場所場所で軌道の分散の強弱を許す (そして実際にそのように考えると時系列データに対するより深い理解につながるという意味で，Wayland test の拡張になっていると言える．スヵラーデータ$)$に対する埋め込み時系列解析の文脈における確率微分方程式モデルを用いた貢献は今後考えるべき問題であろう．

拡散行列の近似式においては，連続形でも離散形でも二次変動量から一次変動量の二次項を除去するのが本稿におけるポイントであるが，これは確率論的な理由
によって生じることではなく，空間離散化 (粗視化) した枠組みで計算することで発生する問題点であることがわかる．例を挙げよう．
$X(t)=(\begin{array}{l}X(t)Y(t)\end{array})=(\begin{array}{ll}cost -sintsint cost\end{array})(\begin{array}{l}xy\end{array})$(2.16)という $x=(x, y)^{T}$ を出発する角速度一定の円周軌道を考える．ここで $x_{0}=(x_{0}, y_{0})^{T}$を中心とする一辺が $\delta>0$ な正方形 $I=[x_{0}-\delta/2, x_{0}+\delta/2$ ) $\cross[y_{0}-\delta/2, yo+\delta/2$ )上から，一様に (2.16) の軌道が出発するとして各変動量を計算すると，$\frac{1}{|I|}\int_{I}(X(t)-x)dx=(\begin{array}{ll}x_{0}(cost-1)- y_{0}costy_{0}(cost-1)x_{0}sint+ \end{array}),$$\frac{1}{|I|}l|X(t)-x|^{2}dx=(1-\cos t)(2|x_{0}|^{2}+\frac{\delta^{2}}{3})$,(2.17)$\frac{1}{|I|}\int_{I}|X(t)-x|^{2}dx-|\frac{1}{|I|}l(X(t)-x)dx|^{2}=(1-\cos t)\frac{\delta^{2}}{3}$を得る．(2.18) は 上の軌道の分散を計算しているため，$I$となることを加味すると$0$$\deltaarrow 0$の極限で $Iarrow\{x\}$に収束することは自然である．一方，(2.17) ではでも $(t>0$ かつ $|x_{0}|>0$ である限り

$)$

値は

$0$

(2.18)$\deltaarrow 0$に収束しない．これは，(2.17) が$I$

上の軌道の 2 次モーメントの平均であることを考えると，やはり自然なことである．

54
これら二つの量は性質の異なる量であり，目的に応じて計算する量を選ばねばならない．本研究で着目しているような系の相空間におけるランダム性の強さには，当然 2 次モーメントでなく軌道の分散度を採用すべきであるから，(2.18) のように二次変動量から一次変動量の二乗を引く量を計算する方が好ましいということがわ
かる．

3

数値的検証

ここでは，3 次元非線形常微分方程式系の Lorenz system と R\"ossler system の数値解をある 2 次元平面に射影した得た時系列データに対して，(2.1) と (2.2) でドリフ

トベクトルと拡散テンソルを計算し，統計的係数決定公式の有効性を数値的に検証する．ここで，カオス時系列を用いる理由は二度と同じ軌道を通らない長い時系列を作成できるためである．また，3 次元空間の解軌道を 2 次元平面に射影すること
で，高次元力学を低次元相空間へ縮約した力学を模擬している．
ここで用いる Lorenz system$x’=-px+py,$ $y’=-xz+rx-y,$ $z’=xy-bz,$ $p=10,$ $r=28,$ $b=8/3$と R\"ossler system$x’=-y-z, y’=x+ay, z’=b+xz-cz, a=0.2, b=0.2, c=5.7$

Figure 2: Lorenz system (左) と R\"ossler system (右) の数値解の軌道．
の数値解は，どちらも初期値 (1.0,1.0,1.0), 時間幅 $1.0\cross 10^{-5}$ を用いて 4 次 RungeKutta 法で時間発展を計算し，初期から $t<100$ までを遷移過程として打ち切り，

55$t=100$ から 100 ステップおきに$10^{8}$個データを採取する．このため，データのサ

ンプル時間間隔は $\triangle t=1.0\cross 10^{-3}$ となる．得られた軌道を Figure 2 に示す．

さらに，得られた 3 次元データに対してそれぞれ主成分分析を行い，上位 2 モー

ドが張る 2 次元部分空間に射影し，これらを解析対象の時系列データとする．射影された軌道は Figure 3 のとおりである．双方とも平面への射影に因る軌道の擬似交差が確認できる．二次元平面上の力学としては，この軌道交差領域では次ステッ
プでの状態の変化が定まらず，二次元の決定論的力学として定式化することができないことを表す．

Figure 3: 主成分分析による上位 2 モードが張る部分空間へ射影された Lorenz system(左) と R\"ossler system (右) の軌道．

Figure 4:$\overline{A}(\overline{x}_{m})$のベクトルフ o ロット: Lorenz system (左) と R\"ossler system (右).

以上の手続きで得られた相空間 (射影された平面) 上の時系列データに対し，(2.1) と (2.2) の統計的係数決定公式を用いて相空間におけるドリフトと拡散の強さを評価する．Lorenz system ではムオあたりの平均移動距離が 0.09 であるため，$\tau=1$とし，相空間の分割セルのサイズを O. 1 と取る．R\"ossler system は平均移動距

離が 0.008 なので， $\tau=10$ として 0.05 の正方形に取る．ドリフトベクトル (Figure

564) は概ね軌道の平均的な挙動を表すことができている．Lorenz system の射影デー
タの場合 (Figure 4 左),軌道が相空間の原点近傍で交差し且つ交差する軌道の数は左右で拮抗しているため，平均的にベクトルは第 2 主成分が負の方向を向く．したがって，左半平面から右半平面へ，またはその逆に移動するためのドリフトベク
トルは得られていない．すなわち，この枠組みではそのような相空間内の移動は 2次元平面上の決定論的力学としては表されないことを示している．また，R\"ossler
system の場合 (Figure 4 右) は軌道の向き (時計回り) と軌道の交差部はきちんと抽出できている．また，元の 3 次元空間では軌道が 3 次元的に他の軌道を乗り越えるときに速く遷移するが，その軌道の速さは射影データにもドリフトベクトルの大きさとして表されている．

$ao$

$2s$

$os$
$lQ$

20
1

30

-a

$/0$

10

$\infty$$\infty$$\alpha, -\infty t0 \prime 0 \infty \infty$

Figure 5: 主成分分析による上位 2 モードが張る部分空間へ射影されたデータから求(下) の常用対数のプロット:Lorenz system (左), R\"ossler
めた $trB$ (上) と$tr\overline{\mathbb{B}}_{new}$system (右). 双方ともに時系列点が 100 以上滞在したセルのみ表示している．
次に拡散の大きさを表す$tr\overline{\mathbb{B}}_{new}$

5 に示す．まず，従来の公式によると従来の公式 $(1.4)_{2}$ で求めた
$tr\overline{\mathbb{B}}$$tr\overline{\mathbb{B}}$の図を Figureの値では，軌道が速い (遅い) 部分で値が大57
きい (小さい)という情報しか得ることができないことが Lorenz system, R\"ossler

system の双方からわかる (Figure 5 左上，右上). これは (2. 11) にあるように，$\overline{\mathbb{B}}$の評価にドリフトベクトルの寄与が大きく出ているためである．一方で，修正された公式での$tr\overline{\mathbb{B}}_{new}$では Lorenz system,

R\"ossler system のどちらも軌道交差の領域

で値が大きくなっており (Figure5 左下，右下), そこでの軌道の分散を捉えられていることがわかる．ただ，$tr\overline{\mathbb{B}}_{new}$にも軌道の速さが影響を及ぼさないゎけではない．

$30\ovalbox{\tt\small REJECT}_{2\mathfrak{B}}^{2}\sim 235-23\fbox{Error::0x0000}-2t2.22..162.\cdot\cdot 05$$\sim 1020\infty t0\triangleleft 0$$\infty$

$0$

10

20

Figure 6: 主成分分析による上位 2 モードが張る部分空間へ射影されたデータから求めたノイズ指標 (下) と従来の公式を用いた$tr\overline{\mathbb{B}}/|\overline{A}|^{2}$(上) の常用対数のプロット:Lorenz system (左), R\"ossler system (右). 双方ともに時系列点が 100 以上滞在したセルのみ表示している．
特に Lorenz system の方が顕著であるが，軌道の速い第 1 主成分の絶対値が大きい
領域では$tr\overline{\mathbb{B}}_{new}$の値が大きい傾向が見え，かつ軌道の遅い第 1 主成分が 付近か$0$

つ第 2 主成分が負の領域では値はとても小さくなっていることがゎかる．軌道の単位時間当たりの移動距離が長い領域では，仮に決定論性が強くとも軌道の分散は大58きくなり得るからである．これは Example2.4 の注意で指摘したように，ドリフト
ではの強さがその評価に関わってくることとも関係している．したがって，$tr\overline{\mathbb{B}}_{new}$

軌道の分散を適切に評価する指標であるといえる一方で，時系列の決定論性と確率論性の分離という点に着目する場合はドリフトベクトルの大きさで規格化する方が良い．

Figure 6 は拡散行列のノルムとドリフトベクトルの大きさの二乗の比を示したものである．これは軌道の変動に対する決定論性と確率論性の強さの指標であり，値が大きければ軌道の変動には確率論的な寄与が大きいことを意味する．従来の公式を用いた$tr\overline{\mathbb{B}}/|A|^{2}-$の場合は (Figure 6 上), 新しい公式を用いたもの (Figure 6 下)と比較すると Lorenz system, R\"ossler system 共に値域が狭く，強弱の傾向を掴み
では，Lorenz system, R\"ossler
づらいことがわかる．Figure 6 (下) の$tr\overline{\mathbb{B}}_{new}/|A|^{2}-$system 双方とも軌道交差の領域で値が大きくなっているが，軌道の速くかつ決定論性の強い領域 (Lorenz では原点から遠い領域，R\"ossler では下半平面の原点から離れた領域や軌道交差直前の領域) では小さくなる傾向が見られる．

$tr\overline{\mathbb{B}}_{new}/|A|^{2}-$では軌道の速さによる影響も規格化されるため．軌道の変動のうち確率論性の寄与が強い領域を検出することができている．特に不安定固定点まわりでその値が小さくなっていることも特筆すべきことである．固定点まわりは軌道速度が遅いものの
その不安定性によって軌道がばらけるため，相対的に確率論性が高い領域となっているとみなせる．

これらの例で示した通り，新しい公式$tr\overline{\mathbb{B}}_{new}$を用いると適切に相空間内の拡散(軌道の分散の平均変化率) の強い領域を抽出することが可能であり，本研究で提案する統計的係数決定公式 (2.1), (2.2) は数学的にも数値的にも機能することがわを合わせて調べることにより，軌道の分散がとかった．また，$tr\overline{\mathbb{B}}_{new}$$tr\overline{\mathbb{B}}_{new}/|A|^{2}-$大きい部分からさらに確率論的寄与が大きい領域を捉えることができる．これにより時系列データの決定論性と確率論性の強弱が説明付けられるため，本稿で扱った
確率微分方程式モデリングの有効性を示している．

4

謝辞

本研究を進めるにあたり，MaeT 数学-気象学連携研究チーム，京都大学大学院理学研究科の坂上貴之教授，京都大学数理解析研究所の山田道夫教授，東北大学原子分子材料科学高等研究機構の小布施祈織助教，京都大学大学院理学研究科の森田英59
俊博士研究員には有意義な議論や厳しくも温かいご指摘をいただきましたことをこの場を借りて御礼申し上げます．また，本研究は文部科学省 H22 年度 「気候変動適応戦略イニシアチブ」気候変動適応研究推進プログラム 「北海道を対象とする総合的ダウンスケール手法の開発と適用」 と日本学術振興会科学研究費助成事業 (挑戦的萌芽研究 :課題番号 25610028) の支援を受けたものです．

References[1] J. BERNER: Linking nonlinearity and non-Gaussianity of planetary wave be-havior by the Fokker-Planck equation. J. Atmos. Sci. 62, 2098-2117 (2003).
[2] G. BRANSTATOR AND J. BERNER: Linear and nonlinear signatures in theplanetary wave dynamics of an AGCM: phase space tendencies. J. Atmos. Sci.

62, 1792-1811 (2003).
[3] A. BHADRA, E. L. IONIDES, K. LANERI, M. PASCUAL, M. BOUMA AND R.
C. DHIMAN: Malaria in northwest India: Data analysis via partially observed

stochastic diﬀerential equation models driven by L\’evy noise. J. $Am$ . Stat. Assoc.
106, 440-451 (2011).
[4] S. DITLEVSEN AND A. DE GAETANO: Mixed eﬀects in stochastic diﬀerential
equation models. REVSTAT 3, 137-153 (2005).

[5] K. HASSELMANN: Stochastic climate models. Part I. Theory. Tellus 28, 473-

485 (1976).
[6] A. F. IVANOV AND A. V. SWISHCHUK: Optimal control of stochastic diﬀerential delay equations with application in economics. Int. J. Qual. Theory Diﬀ.

Eq. Appl. 2, 201-213 (2008).
[7] M. INATSU, N. NAOTO AND H. MUKOUGAWA: Dynamics and predictabilityof extratropical wintertime low-frequency variability examined by a stochastic

diﬀerential equation in a low-dimensional system. J. Atomos. Sci., 70, 939-952(2013).

60[8] E. N. LORENZ: A study of the predictability of a 28-variable atmosphericmodel. Tellus 17, 321-333 (1965).
[9] V. MISRA, W. -B. GONG AND D. TOWSLEY: Stochastic diﬀerential equation
modeling and analysis of TCP-windowsize behavior. J. Atmos. Sci. 62, 1792-

1811 (2003).
[10] 中野直人，稲津將，向川均，楠岡誠一郎: 確率微分方程式を用いた気候モデルについて．京都大学数理解析研究所講究録 1823, 79-96 (2013).
[11] P. SURA, M. NEWMAN, C. PENLAND AND P. D. SARDESHMUKH: Multiplicative noise and non-Gaussianity: A paradigm for atmospheric regimes?
$J.$

Atmos. Sci. 62, 1391-1409 (2005).
[12] R. WAYLAND, D. BROMLEY AND D. PICKETT: Recognizing determinism ina time series. Phys. Rev. Lett. 70, 580-582 (1993).
[13] D. J. WILKINSON: Stochastic modelling for quantitative description of heterogeneous biological systems. Nature 122, 122-133 (2009).
[14] N. ZABARAS AND SETHURAMAN SANKARAN: An information-theoretic approach to stochastic materials modeling. Computing in Science & Engineering
9, 50-59 (2007).