確率微分方程式

このページでは、簡単な確率微分方程式の数値計算について考えてみる。

ゆらぎを伴う微分方程式

ランダムなゆらぎ$\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} $$ とする。

次に、時系列信号の時間のスケールを1/2、時間あたりのエネルギーを2倍にするような変換を施してみる。 つまり、波形を細かくしても、「時間刻み」あたりのエネルギーは同じになるように調整するわけである。 そうするためには、グラフの横軸を1/2に、縦軸を$\sqrt{2}$倍に変形すれば良い（下図）。

こうした操作を繰り返すことで、その自己相関が $$ \left\langle v(t) v(t') \right\rangle = R \delta(t-t') $$ で、ランダムなスパイクが連続するような波形が最終的に得られ、白色雑音（ホワイトノイズ）と呼ばれている。

ウィーナー・ヒンチンの定理により、 時系列信号の自己相関関数をフーリエ変換すると信号のパワースペクトル（振動数あたりのエネルギー密度）が得られるので、これを白色雑音に適用すると、 $$ 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()