拡散現象を媒介するネットワークのプロファイリング～反応拡散過程の統計力学と統計的推論～（８）

今回は、ランジュバン方程式をコンピュータで数値的に解いて、ゆらぎを伴う系の時間発展の様子を具体的に調べてみよう。

擬似乱数による数値積分

形式的には、ランジバン方程式は、時間についての１階の微分方程式である。微分方程式の初期値問題を数値的に解くには、連続的な時間を離散的な時間間隔τで近似（discrete time approximation）し、微分方程式を差分方程式に変換する。つまり、以下の関数fで規定される微分方程式から、関数gで規定される差分方程式を導いて解く。

差分方程式が近似する精度によって、解の収束性や計算時間が異なる。方程式の形や近似の特性に応じたさまざまなアルゴリズムが考案されている。大きく、以下の３つの要因が、アルゴリズムの特性を左右する。

• 関数gの値を得るのに、関数fの値だけでなく微分係数の値まで必要かどうか？
• x(t+τ)の値を得るのに、x(t)だけでなく、x(t-τ)などの過去の時刻の値まで必要かどうか？
• g=0から、四則演算や初等関数のみでx(t+τ)の値を陽に得られず、ニュートン法などの球根法（root finding algorithm）が必要かどうか？

4次のルンゲ・クッタ法（Runge-Kutta method）は、高い精度と少ない計算量を両立させるバランスのとれたアルゴリズムである。微分係数の値は不要、x(t-τ)などの値は不要、x(t+τ)の値を陽に得られる、という特性がある。一方、確率微分方程式に特化したアルゴリズムであるミルシュタイン法（Mistein method）では、微分係数の値が必要になる。編微分方程式を解くアルゴリズムであるクランク・ニコルソン法（Crank-Nicolson method）では、x(t+τ)の値を陽に得られず、τごとに代数方程式を解く必要がある。

このようなアルゴリズムを適用する際に、確率微分方程式に特有なゆらぎの具体的な時系列を与える必要がある。事前に、時系列を生成して準備しておいてもよい。たいていは、微分方程式の数値積分を行いながら、同時に次の時刻のゆらぎを生成する。擬似乱数によって、与えられた統計的な性質を満たすゆらぎを発生させる。いわゆる、モンテ・カルロ・シミュレーションである。

ルンゲ・クッタ法に現れる差分方程式の時間間隔τごとに、擬似乱数を生成する。ガウス型白色雑音では、ボックス・ミュラー法（Box-Muller method）を用いて、一様乱数から正規分布に従う擬似乱数を生成する。τより大きな時間間隔では、相関のない時系列になる。τより小さな時間間隔では、そうならない。この意味で、連続的なゆらぎを離散的な擬似乱数で近似している。

さて、τを小さくすると、数値積分で得た解が、個々のゆらぎで決まる真の解に厳密に一致する性質を、強収束性（strong convergence）という。一方、数値積分で得た解がさまざまなゆらぎで決まる真の解の全体としての統計的な特性と一致する性質を、弱収束性（weak convergence）という。通常、個々のゆらぎを知ることはできないので、弱収束性が必要な要件となる。

時間発展の数値解析

いくつかの具体的な例で、時間発展の様子を見てみよう。

まず、1次元のブラウン運動する粒子の運動を表すウィーナー過程（Wiener process）を記述するランジバン方程式を解いてみる。ウィーナー過程では、粒子の位置の平均は常に0で、変わらない。分散は、経過時間に比例する。多数の擬似乱数の時系列について、ランジバン方程式を解いた結果を以下の図に示す。

横軸は時間、縦軸は位置である。同じ統計的な性質に従っていても、擬似乱数の時系列ごとに位置の変動が随分異なることが分かる。2枚のうちの上の図は、ゆらぎがガウス型白色雑音の場合である。

もうひとつの下の図は、白色雑音であるが、振幅が正規分布ではなく対数正規分布（logarithmic normal distribution）に従う場合である。どちらも、平均は0、分散は1の雑音である。正規分布と異なり、対数正規分布の3次以上のモーメントの値は0にならない。

3次のモーメントは、歪度（わいど, skewness）を表す。歪度が大きいと、確率密度関数が左右非対称になる。正の値と負の値の現れ方の偏りが大きくなる。例えば、たいていは小さな負の値が表れ、時々それらを相殺する大きな正の値が現れる。4次のモーメントは、尖度（せんど、kurtosis）を表す。尖度が大きいと、確率密度関数の裾の部分の寄与が相対的に大きくなる。いわゆる、ファット・テール型の分布になる。時たま、桁違いの大きな変動が起こる。

株価の変動の分布では、歪度、尖度が大きく、正規分布のゆらぎでは想定外となる大きな負の変動が起こることが知られている。

次に、多数の擬似乱数の時系列について、2次元のブラウン運動する粒子の運動を表すランジバン方程式を解いた結果を以下の図に示す。2つの軸は、x方向、y方向それぞれの次元である。初期状態で原点から出発した粒子の軌道を描いた。

同じ方程式に従う無数の粒子の動きを重ね合わせると、空気中を煙が広がるような運動になる。これを流体の連続的な流れと見た近似が、拡散という現象である。

最後に、単一のノードにおいて、SIRモデルで記述できる感染症の感染者数が拡大する様子を以下の図に示す。感染者数の拡大に対応した成長のトレンドに、ブラウン運動のゆらぎを重ね合わせたような時間発展となる。

次回は、メタポピュレーション・ネットワークにおけるSIRモデルを表すランジュバン方程式を解いてみる。

Dr. Yoshiharu Maeno, Social Design Group.