確率微分方程式
このページでは、簡単な確率微分方程式の数値計算について考えてみる。
ゆらぎを伴う微分方程式
ランダムなゆらぎ$\xi(t)$を考え、それによって動かされる粒子の位置$x(t)$が微分方程式 $$ \frac{d x(t)}{dt} = \xi(t) \tag{1} $$ で変化するようなモデルを考えてみよう。 ここで、 ゆらぎの平均は0、すなわち$i$番目のサンプルの揺らぎの時系列を$\xi_i(t)$とすれば $$ \left\langle \xi(t) \right\rangle = \lim_{N\to\infty} \frac{1}{N} \sum_{i=1}^N \xi_i(t) = 0 $$ で、自己相関が $$ \left\langle \xi(t) \xi(t+\tau) \right\rangle = \lim_{N\to\infty} \frac{1}{N} \sum_{i=1}^N \xi_i(t) \xi_i(t+\tau) = R \delta(\tau) $$ であるような場合を考える($\delta(\tau)$はデルタ関数)。 そして、これらの関係は時間の原点には依らずに成り立つ($\xi(t)$は定常的)と仮定する。
この種の確率的(ランダム)な変数を含む微分方程式は確率微分方程式(stochastic differential equation)と呼ばれている。 確率微分方程式を数値的に解くには、ゆらぎの部分の扱い方に、通常の常微分方程式には無かった注意が必要である。
確率微分方程式の差分化
ここで、時間$t$を$\Delta$毎の区間に $$ [0, \Delta),\ [\Delta, 2\Delta),\ [2\Delta, 3\Delta),\cdots,[k \Delta, (k+1)\Delta),\cdots $$ のように等分割して、 $k$番目の区間の揺らぎの累積(積分)を $$ q_k = \int_{k\Delta}^{(k+1)\Delta} \xi(s) ds $$ と定義する。 すると、そのサンプル平均は $$ \left\langle q_k \right\rangle = \lim_{N \to \infty} \frac{1}{N} \sum_{i=1}^{N} \int_{k\Delta}^{(k+1)\Delta} \xi_i(s) ds = \int_{k\Delta}^{(k+1)\Delta} \left\langle \xi(s) \right\rangle ds = 0 $$ 分散は $$ \begin{eqnarray} \left\langle q_k q_{k+\ell} \right\rangle & = & \lim_{N\to\infty} \frac{1}{N} \sum_{i=1}^{N} \left( \int_{k\Delta}^{(k+1)\Delta} \xi_i(s) ds \right) \left( \int_{(k+\ell)\Delta}^{(k+\ell+1)\Delta} \xi_i(t) dt \right) \\ & = & \lim_{N\to\infty} \frac{1}{N} \sum_{i=1}^{N} \left( \int_{k\Delta}^{(k+1)\Delta} \xi_i(s) ds \right) \left( \int_{k \Delta}^{(k+1)\Delta} \xi_i(t + \ell\Delta) dt \right) \\ & = & \int_{k\Delta}^{(k+1)\Delta} \int_{k \Delta}^{(k+1)\Delta} \left[ \frac{1}{N} \lim_{N\to\infty} \sum_{i=1}^{N} \xi_i(s) \xi_i(t + \ell\Delta) \right] ds dt \\ & = & \int_{k\Delta}^{(k+1)\Delta} \int_{k \Delta}^{(k+1)\Delta} \left\langle \xi(s) \xi(t + \ell\Delta) \right\rangle ds dt \\ & = & \int_{k\Delta}^{(k+1)\Delta} \int_{k \Delta}^{(k+1)\Delta} R \delta(t + \ell\Delta -s) dt ds \\ & = & R \int_{k\Delta}^{(k+1)\Delta} \left\{ \int_{k \Delta + \ell \Delta - s}^{(k+1)\Delta + \ell \Delta - s} \delta(u) du \right\} ds \\ \end{eqnarray} $$ となる。 ここで、デルタ関数の積分 $$ \int_{k \Delta + \ell \Delta - s}^{(k+1)\Delta + \ell \Delta - s} \delta(u) du $$ は、$\ell=0$の場合に限り、$k \Delta \le s \le (k+1)\Delta$の範囲で積分区間が$u=0$を「通過」すること、 および $\epsilon \gt 0$ に対して $\int_{0-\epsilon}^{0+\epsilon} \delta(u) du=1$ を考慮すると、 $$ \left\langle q_k q_{k+\ell} \right\rangle = R \Delta \delta_{0,\ell} $$ が得られる($\delta_{0,\ell}$はクロネッカーのデルタ)。
オイラー・丸山の近似
式(1)の両辺を$k\Delta$から$(k+1)\Delta$まで積分すると $$ x(k\Delta + \Delta) - x(k\Delta) = q_k $$ となるが、$x(k\Delta)$をあらためて$x_k$と書けば $$ x_{k+1} = x_k + q_k $$ である。 ここで、分散$R$を持つ乱数を$\xi$とすると、$q_k$の分散はその$\Delta$倍(標準偏差は$\sqrt{\Delta}$倍)であるから、 $$ x_{k+1} = x_k + \sqrt{\Delta} \xi $$ という差分近似によって$\Delta$だけ先の$x(t)$を数値的に見積もることが可能となる(オイラー・丸山法)。
同様に $$ \frac{dx(t)}{dt} = f(x(t)) + \xi(t) $$ に対しては、 $$ x_{k+1} = x_k + f(x_k) \Delta + \sqrt{\Delta} \xi $$ となり、確率的な揺らぎの項は他と時間差分$\Delta$の次数が異なってくる。
解説: 白色ノイズ
$\Delta$毎に変動し、その値が平均0分散$R$で分布するような無限に続くランダムな時系列 $v(t_i)$を考える。 ここで、$i$は何番目の区分か、を表すことにしよう。 そして、$i$番目と$j$番目($i \neq j$)は無相関、すなわち $$ \left\langle v(t_i) v(t_j) \right\rangle = R \delta_{i,j} $$ とする。
一般に波動の時間あたりのエネルギーは振幅の2乗に比例する。
次に、時系列信号の時間のスケールを1/2、時間あたりのエネルギーを2倍にするような変換を施してみる。 つまり、波形を細かくしても、「時間刻み」あたりのエネルギーは同じになるように調整するわけである。 そうするためには、グラフの横軸を1/2に、縦軸を$\sqrt{2}$倍に変形すれば良い(下図)。
こうした操作を繰り返すことで、その自己相関が $$ \left\langle v(t) v(t') \right\rangle = R \delta(t-t') $$ で、ランダムなスパイクが連続するような波形が最終的に得られ、白色雑音(ホワイトノイズ)と呼ばれている。
図はWikipediaへのリンク
ウィーナー・ヒンチンの定理により、 時系列信号の自己相関関数をフーリエ変換すると信号のパワースペクトル(振動数あたりのエネルギー密度)が得られるので、これを白色雑音に適用すると、 $$ S(\omega) = \int_{-\infty}^{\infty} R \delta(\tau) e^{-i \omega \tau} d\tau = R $$ となり、全ての角振動数について、エネルギー密度が一定($R$)となる。 そうすると、全ての振動数についてエネルギーを積分すると無限大となってしまうことになるが、 現実の問題では、髙振動数側のどこかに必ずカットオフが存在するため、信号のエネルギーが発散するようなことは無い。 全ての振動数で強度が一定となることから、光のスペクトルになぞらえられたのが「白色」の由来である。
練習: 確率微分方程式
オイラー・丸山の方法で微分方程式を数値的に解くコードの例を以下に示した。
これを雛形にして、$x(0)=0$を初期条件とした複数の経路を重ね書きするコードを作成してみなさい。
# coding: utf-8 import math import random import matplotlib.pyplot as plt R=1 t=0 dt=0.01 x=0 X=[] T=[] while t<10: X.append(x) T.append(t) x = x + math.sqrt(dt) * random.gauss(0,R) t += dt plt.plot(T,X, color='red', linewidth=1.0) plt.grid(True) plt.xlabel('t') plt.ylabel('x') plt.show()