Botter自主ゼミノート 2.9 ウィーナ過程(2)

やること

https://www.amazon.co.jp/dp/4254209444
を読んで、確率微分方程式による最適化問題を解けるようになることです。

これまでのBotter自主ゼミノート

Botter自主ゼミノート 1.2 確定システムの制御の回顧
Botter自主ゼミノート 2.1 確率過程とは？, 2.2 確率過程の数学的表現
Botter自主ゼミノート 1.2 数式導出
Botter自主ゼミノート 2.3 確率モーメント
Botter自主ゼミノート 2.4 確率過程の分類
Botter自主ゼミノート 2.5 エルゴード性
Botter自主ゼミノート 2.6 確率過程の周波数表現
Botter自主ゼミノート 2.7 マルコフ過程
Botter自主ゼミノート 2.8 正規型確率過程
Botter自主ゼミノート 2.9 ウィーナ過程(1)

2.9 ウィーナ過程

ブラウン運動とウィーナ過程の関係性について

ウィーナ過程は物理現象としてよく知られたブラウン運動を数学的に記述したものです。教科書では、ブラウン運動過程が式(2.58)で与えられるような正規型過程であることを、液体中の粒子の運動方程式を考えることで示しています。

$$p(t,w) = \frac{1}{\sqrt{2\pi t}\sigma}\exp\left(-\frac{w^2}{2\sigma^2 t}\right) \tag{2.58}$$

まず、液体中の粒子は質量$m$の質点とし、その位置を$w(t)$(ただし一次元)とします。この粒子の速度を$u(t)=\frac{dw(t)}{dt}$とします。

すると、ニュートンの第二法則(Wikipedia)より、式(2.70)が成り立ちます。

$$m\frac{du(t)}{dt} = -\xi u(t) + \gamma_0(t)\quad (t \ge 0) \tag{2.70}$$

式(2.70)の右辺第一項は液体による粘性抵抗力で、$\xi$はその係数です。また、$\gamma_0(t)$は粒子に働く外力です。粒子が不規則に動く原因として、外力$\gamma_0(t)$が確率的に不規則に動くものと考えます。外力$\gamma_0(t)$は粒子に衝突する液体分子による力であると考え、白色雑音と仮定します。($\gamma_0(t)$と$\gamma_0(\tau)$ $(t \ne \tau)$は互いに相関を持たないものとします)

式(2.70)のような運動方程式をランジュヴァン方程式といい、ブラウン運動の解析モデルとして用いられています。この式を書き直すと式(2.71)のようになります。

$$\frac{du(t)}{dt} + \beta u(t) = \gamma(t) \tag{2.71}$$

ただし、$\beta = \frac{\xi}{m}, \gamma(t) = \frac{\gamma_0(t)}{m}$とします。また、$\gamma(t)$は平均値と相関関数が式(2.72)のように与えられる確率過程であると仮定します。ここで、$D$は正の定数、$\delta(t)$はディラックのデルタ関数です。これは、液体の分子が粒子に不規則に衝突し、その相関は$\delta(t-s)$と仮定するのが適切という考え方から来ています。

$$\mathcal{E}\{ \gamma(t) \} = 0, \quad \mathcal{E}\{ \gamma(t)\gamma(s)\} = 2D\delta(t-s) \tag{2.72}$$

ディラックのデルタ関数は

$$\delta(x) = \begin{cases} 0 & (x \ne 0)\\ \infty & (x = 0) \end{cases}, \quad \int_{-\infty}^{\infty} \delta(x) dx = 1$$

なので、不規則に衝突する粒子からの力には時間的な相関はまったくないということでしょうか。

式(2.71)は一階線形常微分方程式なので、これを解くと式(2.73)となります。

$$u(t) = u_0 e^{-\beta t} + \int_{0}^{t} e^{-\beta(t-s)} \gamma(s) ds \tag{2.73}$$

粒子の速度$u(t)$は、初期速度$u_0$が粘性抵抗によって弱まりながらも積分される値と、外力$\gamma_0(t)$による加速分がやはり徐々に粘性抵抗によって弱まりながら積分される値で決まっている、ということが分かります。

ラグランジュの定数変化法を使って式(2.71) (以降式Aと呼ぶ)を解いてみます。

$$\frac{du(t)}{dt} + \beta u(t) = \gamma(t) \tag{A}$$

まず式(2.71)の右辺に$0$を代入した斉次方程式の一般解を求めます。

$$\begin{aligned} \frac{du(t)}{dt} + \beta u(t) &= 0 \\ \frac{du(t)}{dt} &= - \beta u(t) \end{aligned}$$

$u(t)$を微分すると$-\beta u(t)$になるということは、

$$u(t) = C e^{-\beta t} \tag{B}$$

になりそうです。

式(B)の$C$を$t$の関数として、式(A)に代入します。

$$\begin{aligned} (C(t) e^{-\beta t}) \frac{d}{dt} + \beta C(t) e^{-\beta t} &= \gamma(t)\\ C'(t) e^{-\beta t} -\beta C(t) e^{-\beta t} + \beta C(t) e^{-\beta t} &= \gamma(t) \\ C'(t)& = e^{\beta t} \gamma(t) \\ C(t) &= \int e^{\beta t} \gamma(t) dt + C\\ \end{aligned} \tag{C}$$

式(C)を式(B)に代入します。このとき式(C)では$t$について積分しているので、変数を$s$に置換しておきます。

$$\begin{aligned} u(t) &= e^{-\beta t} (\int e^{\beta s} \gamma(s) ds + C) \\ &= C e^{-\beta t} + \int e^{-\beta t} e^{\beta s} \gamma(s) ds &= C e^{-\beta t} + \int e^{-\beta(t-s)} \gamma(s) ds \end{aligned}$$

$t = 0$のときの$u(t)$の値を$u_0$とするために、$C$を$u_0$で置き換えると、以下のようになり、式(2.73)と等しくなります。

$$u(t) = u_0 e^{-\beta t} + \int e^{-\beta(t-s)} \gamma(s) ds$$

ここで、$u(t) = \frac{dw(t)}{dt}, w(0) = 0$に留意して式(2.73)を積分します。

$$\begin{aligned} w(t) &= \int_{0}^{t} u(s) ds\\ &= \int_{0}^{t} \left[ u_0 e^{-\beta s} + \int_{0}^{s} e^{-\beta (s-\tau)} \gamma(\tau) d\tau \right] ds \tag{2.74} \end{aligned}$$

式(2.74)を変形して以下の式(2.75)を得ます。

$$\begin{aligned} w(t) - u_0 \int_{0}^{t} e^{-\beta s} ds &= \int_{0}^{t} e^{-\beta s} \left[\int_{0}^{s} e^{\beta \tau} \gamma(\tau) d\tau \right] ds \tag{2.75} \end{aligned}$$

左辺の積分を実行して…

$$\begin{aligned} w(t) - \frac{1}{\beta} u_0 (1 - e^{-\beta t}) &= \int_{0}^{t} e^{-\beta s} \left[\int_{0}^{s} e^{\beta \tau} \gamma(\tau) d\tau \right] ds \end{aligned}$$

ここで教科書では右辺の部分積分を実行して以下の式(2.76)を得るとされています。

$$\begin{aligned} w(t) - \frac{1}{\beta} u_0 (1 - e^{-\beta t}) &= \frac{1}{\beta} \int_{0}^{t} (1 - e^{-\beta (t-s)}) \gamma(s) ds \tag{2.76} \end{aligned}$$

左辺第二項は、粘性抵抗によって減衰し続けている初期速度の$t$までの積分結果で、$t$時点での$w(t)$の平均値に相当します。左辺は白色雑音$\gamma(t)$とは無関係なので、確率変数ではありません。

白色雑音$\gamma(\tau)$が絡んだ部分積分なので、どういう結果になるのかちょっとわからない。と言うことを書いていたら、Twitterに書いたら @d_osh8 さんが教えてくれました。神…ありがとうございます。

まず、$f(s)$を以下のように定義して式(2.75)の右辺を置き換えます。

$$f(s) = \int_{0}^{s} e^{\beta \tau} \gamma(\tau) d\tau$$

すると、式(2.75)は以下のような形になります。

$$\begin{aligned} \int_{0}^{t} e^{-\beta s} \left[ \int_{0}^{s} e^{\beta \tau} \gamma(\tau) d\tau \right] ds &= \int_{0}^{t} e^{-\beta s} f(s) ds\\ \end{aligned}$$

部分積分の公式を使って…

$$\begin{aligned} &= \left[ \left( \int e^{- \beta s} ds \right) \cdot f(s) \right]_{0}^{t} - \int_{0}^{t} \left( \int e^{- \beta s} ds \right) \cdot f'(s) ds \\ \end{aligned}$$

$f(0) = 0$なので、第一項をシンプルにできます。よいしょよいしょ。

$$\begin{aligned} &= \left[ \left( \int e^{- \beta s} ds \right) \cdot f(s) \right]_{0}^{t} - \int_{0}^{t} \left( \int e^{- \beta s} ds\right) \cdot f'(s) ds \\ &= \left( \int e^{- \beta t} ds \right) \cdot f(t) - \int_{0}^{t} \left( \int e^{- \beta s} ds \right) \cdot f'(s) ds \\ &= -\frac{1}{\beta} e^{- \beta t} \cdot f(t) - \int_{0}^{t} \left( \int e^{- \beta s} ds \right) \cdot f'(s) ds \\ \end{aligned}$$

第一項の$f(t)$を、先に出てきた$f(s)$の定義を使って、置換前のもとの形に戻してあげます。これで第一項は完成です。

$$\begin{aligned} &= -\frac{1}{\beta} e^{- \beta t} \cdot \int_{0}^{t} e^{\beta s} \gamma(s) ds - \int_{0}^{t} \left( \int e^{- \beta s} ds \right) \cdot f'(s) ds \\ &= -\frac{1}{\beta} \cdot \int_{0}^{t} e^{-\beta (t - s)} \gamma(s) ds - \int_{0}^{t} \left( \int e^{- \beta s} ds \right) \cdot f'(s) ds \tag{A}\\ \end{aligned}$$

続けて、以下の$f'(s)$の定義を利用して第二項をシンプルにしていきます。ここで再度$f(0) = 0$の特徴を使っています。

$$\begin{aligned} f(s) &= \int_{0}^{s} e^{\beta \tau} \gamma(\tau) d\tau \\ f'(s) &= e^{\beta \s} \gamma(\s) \\ \end{aligned}$$

上記の定義を使って第二項を変形すると…

$$\begin{aligned} &= -\frac{1}{\beta} \cdot \int_{0}^{t} e^{-\beta (t - s)} \gamma(s) ds - \int_{0}^{t} -\frac{1}{\beta} e^{- \beta s} \cdot e^{\beta s} \gamma(s) ds \\ &= -\frac{1}{\beta} \cdot \int_{0}^{t} e^{-\beta (t - s)} \gamma(s) ds + \frac{1}{\beta} \int_{0}^{t} \gamma(s) ds \\ \end{aligned}$$

積分範囲と積分対象の変数が一致しているので、第一項と第二項の積分をくっつけることができます。

$$\begin{aligned} &= \frac{1}{\beta} \left\{\int_{0}^{t} \gamma(s) ds - \int_{0}^{t} e^{-\beta (t - s)} \gamma(s) ds \right\}\\ &= \frac{1}{\beta} \int_{0}^{t} \gamma(s) - e^{-\beta (t - s)} \gamma(s) ds \\ &= \frac{1}{\beta} \int_{0}^{t} \left(1 - e^{-\beta (t - s)} \right) \gamma(s) ds \tag{B} \\ \end{aligned}$$

式(B)は式(2.76)の右辺と一致しています。よかったよかった。

部分積分と言われたとき、てっきり右辺の一番内側を部分積分するのかと思ったら、このように置換するパターンなんですね。一人では気づけなかったと思います。@d_osh8 さん重ねてありがとうございました。

式(2.76)が得られたので、$W(t) := w(t) - \frac{1}{\beta}u_0(1-e^{-\beta t})$の平均値と分散を計算すると…

$$\begin{aligned} \mathcal{E}\{W(t)\} &= \frac{1}{\beta} \int_{0}^{t} \left( 1-e^{-\beta(t-s)}\right) \mathcal{E}\{\gamma(s)\} ds \\ &= \frac{1}{\beta} \int_{0}^{t} \left( 1-e^{-\beta(t-s)}\right) \cdot 0 \: ds \\ &= 0 \end{aligned}$$

$$\begin{aligned} \mathcal{E}\{W^2(t)\} &= \frac{1}{\beta^2} \int_{0}^{t} \int_{0}^{t} \left( 1-e^{-\beta(t-s_1)}\right) \left( 1-e^{-\beta(t-s_2)}\right) \mathcal{E}\{\gamma(s_1)\gamma(s_2)\} ds_1 ds_2 \\ &= \frac{1}{\beta^2} \int_{0}^{t} \int_{0}^{t} \left( 1-e^{-\beta(t-s_1)}\right) \left( 1-e^{-\beta(t-s_2)}\right) 2D \delta(s_1 - s_2) ds_1 ds_2 \\ \end{aligned}$$

$\delta(s_1 - s_2)$は$s_1 = s_2$の時のみ値を持ち、その積分した値は$1$であるので…

$$\begin{aligned} &= \frac{1}{\beta^2} \int_{0}^{t} \left( 1-e^{-\beta(t-s)}\right)^2 2D \: ds \\ &= \frac{2D}{\beta^2} \int_{0}^{t} 1 - 2e^{-\beta t} e^{\beta s} + e^{-2 \beta t} e^{2 \beta s} ds \\ &= \frac{2D}{\beta^2} \left[s - \frac{2}{\beta} e^{-\beta t} e^{\beta s} + \frac{1}{2 \beta} e^{-2 \beta t} e^{2 \beta s} \right]_{0}^{t}\\ &= \frac{2D}{\beta^2} \left(t - \frac{2}{\beta} + \frac{1}{2 \beta} \right) - \left(0 - \frac{2}{\beta} e^{-\beta t} + \frac{1}{2 \beta} e^{-2 \beta t} \right) \\ &= \frac{2D}{\beta^2} \left(t - \frac{3}{2 \beta} + \frac{2}{\beta} e^{-\beta t} - \frac{1}{2 \beta} e^{-2 \beta t} \right) \\ \mathcal{E}\{W^2(t)\} &= \frac{D}{\beta^3} \left(2 \beta t - 3 + 4 e^{-\beta t} - e^{-2 \beta t} \right) \\ \end{aligned}$$

となります。

$w(t) = \int_{0}^{t} u(s) ds$であり、$u(t)$は式(2.73)のように白色雑音の積分和で与えられるので、$w(t)$は中心極限定理によって正規分布します。そのことから、$W(t)$過程の分布は以下のようになります。

$$\begin{aligned} p(t, w) &= \frac{1}{\sqrt{2 \pi \sigma^2}} \exp{-\frac{(w-\mu)^2}{2 \sigma^2}} \\ &= \sqrt{ \frac{\beta^3}{2 \pi D (2 \beta t - 3 + 4 e^{-\beta t} - e^{-2 \beta t} ) } } \exp \left( -\frac{\beta^3 |w - \frac{1}{\beta} u_0(1-e^{-\beta t})|^2}{2D(2\beta t - 3 + 4e^{-\beta t}-e^{-2\beta t}} \right) \tag{2.77}\\ \end{aligned}$$

教科書では上式は$w(t)$過程の分布とされているのですが、多分$W(t)$の間違いではないかな…と思いました。

ここで、$t \gg \frac{1}{\beta}$となる十分長い時刻$t$においては、指数関数$e^{-\beta t}, e^{-2 \beta t}$は微小で、定数項は$2\beta t$に対して無視できます。これに基づいて上式を変形すると…。

$$\begin{aligned} p(t, w) &= \sqrt{ \frac{\beta^2}{4 \pi D t} } \exp \left( -\frac{\beta^2 |w - \frac{1}{\beta} u_0|^2}{4 D t} \right) \tag{2.77}\\ \end{aligned}$$

上式で得られた$W(t)$の分布から、$W(t) \sim \frac{1}{\sqrt{t}}$であることが分かります。であれば$w(t) \sim \frac{1}{\sqrt{t}}$であるので、$\frac{1}{\beta}u_0(1-e^{\beta t})$も$w(t)$と比較すると無視できることが分かります。

これを利用して式(2.77)は式(2.78)に近似できます。

$$\begin{aligned} p(t, w) &\simeq \sqrt{ \frac{\beta^2}{4 \pi D t} } \exp \left( -\frac{\beta^2 w^2}{4 D t} \right) \tag{2.78}\\ \end{aligned}$$

ここで、$D_0 = \frac{D}{\beta^2}$とすると、以下の式(2.79)となります。

$$\begin{aligned} p(t, w) &\simeq \frac{1}{\sqrt{4 \pi D_0 t}} \exp \left( -\frac{w^2}{4 D_0 t} \right) \tag{2.79}\\ \end{aligned}$$

ここで、$\sigma^2 = 2D_0 \: (= \frac{2D}{\beta^2} = \frac{2Dm^2}{\xi^2})$とすると、式(2.79)は前回のノートで出てきた、ウィーナ過程の増分は正規分布するものとする、という定義から求めた確率密度分数式(2.58)と同じになります。

$$p(t,w) = \frac{1}{\sqrt{2\pi t}\sigma}\exp\left(-\frac{w^2}{2\sigma^2 t}\right) \tag{2.58}$$

ブラウン運動の外力$D$の大きさについて

つぎに、外力の大きさ$D$ (あるいは$D_0$)を求めてみます。式(2.73)から、$U(t) := u(t) - u_0 e^{-\beta t}$の平均と分散を求めます。

教科書によると計算は$W(t)$の平均と分散の計算と同様だそうです。

$$\begin{aligned} \mathcal{E}\{U(t)\} &= \int_{0}^{t} e^{-\beta(t-s)} \mathcal{E}\{\gamma(s)\} ds\\ \mathcal{E}\{U(t)\} &= \int_{0}^{t} e^{-\beta(t-s)} \cdot 0 ds\\ \mathcal{E}\{U(t)\} &= 0 \end{aligned}$$

$$\begin{aligned} \mathcal{E}\{U^2(t)\} &= \int_{0}^{t} \int_{0}^{t} e^{-\beta(t-s_1)} e^{-\beta(t-s_2)} \mathcal{E}\{\gamma(s_1) \gamma(s_2)\} ds_1 ds_2\\ &= \int_{0}^{t} \int_{0}^{t} e^{-\beta(t-s_1)} e^{-\beta(t-s_2)} 2D \delta(s_1 - s_2) ds_1 ds_2\\ &= 2D \int_{0}^{t} e^{-2\beta(t-s)} ds\\ &= 2D e^{-2 \beta t} \int_{0}^{t} e^{2\beta s} ds\\ &= 2D e^{-2 \beta t} \frac{1}{2 \beta} \left[ e^{2\beta s} \right]_{0}^{t}\\ &= \frac{D}{\beta} e^{-2 \beta t} (e^{2\beta t} - 1)\\ &= \frac{D}{\beta} (1 - e^{-2 \beta t}) \end{aligned}$$

中心極限定理により、$U(t)$の分布も正規型となるので

$$\begin{aligned} p(t, u) &= \frac{1}{\sqrt{2 \pi \sigma^2}} \exp{-\frac{(u-\mu)^2}{2 \sigma^2}} \\ &= \sqrt{ \frac{\beta}{2 \pi D(1-e^{-2\beta t})} } \exp\left( -\frac{\beta |u-u_0 e^{-\beta t}|^2}{2D(1-e^{-2 \beta t})}\right) \tag{2.80} \end{aligned}$$

$W(t)$の時と同様に$t \to \infty$の安定状態での分布を考えると、

$$lim_{t \to \infty} p(t,u) = \sqrt{\frac{\beta}{2 \pi D}} \exp\left(-\frac{\beta u^2}{2D}\right) \tag{2.81}$$

となる。

教科書によると、統計熱力学で知られていることとして、熱平衡状態における粒子の速度分布は式(2.82)のマックスウェル分布で与えられるそうです。ここで、$k$はボルツマン定数, $T$は絶対温度だそうです。

$$p(t) \propto e^{-\frac{mu^2}{2kT}} \tag{2.82}$$

式(2.81)と式(2.82)は一致しなければいけないので、両者を見比べて、式(2.83)を得ます。

$$\begin{aligned} D &= \frac{\beta k T}{m} = \frac{\xi k T}{m^2}\\ D_0 &= \frac{kT}{\beta m} = \frac{kT}{\xi} \tag{2.83} \end{aligned}$$

もともと$\xi$は粘性抵抗力として導入されていましたが、ストークスの法則により粘性抵抗力は$-6\pi a \eta u$ ($\eta$: 液体の粘性係数)となることから、$\xi = 6\pi a \eta u$で与えられます。

ここまでの計算により、白色雑音$\sigma(t)$の強さを表すものとして導入された係数$D$は、式(2.83)のような物理的意味をもつことが分かりました。また、$D$は実験によって求めることができ、式(2.83)はアインシュタインの関係式と呼ばれています。

ここまでの計算では、外力$\gamma(t)$は時間的に相関のない白色雑音という仮定をしただけで、その分布が正規型である過程はしていないことに注意するようにと教科書にはありました。

ようやっとウィーナ過程の節が終わりました。つかれた…。

時々外部の助けを必要としますが、なんとなく自分で読める距離が伸びてきた気がします！

中心極限定理を使うと、白色雑音から突然正規分布が出てきて驚きます。後で中心極限定理の証明を見ておこう。