📂確率論

ジレスピ確率シミュレーションアルゴリズム

アルゴリズム ¹

$\mathbf{x}_{t} = i$について、 $$P \left( \mathbf{X}_{t + \tau} = j | \mathbf{x}_{t} = i \right) = a_{j} \left( \mathbf{x}_{t} \right) \qquad , j = 1, \cdots , n$$ 連続マルコフ連鎖 $\left\{ \mathbf{X}_{t} \right\}$の実現をシミュレーションする方法は ジレスピー確率シミュレーションアルゴリズム^{gillespie Stochastic Simulation Algorithm}、短く SSA^{stochastic Simulation Algorithm}と呼ばれている。

---

入力

初期値$t=0$と$\mathbf{X}_{0} = \mathbf{x}$を受け取る。

---

ステップ1.

$\sum_{j} a_{j} \left( \mathbf{x}_{t} \right)$を計算する。

---

ステップ2.

指数分布$\exp \left( \sum_{j} a_{j} \left( \mathbf{x}_{t} \right) \right)$から$\tau$をサンプリングする。$\sum_{j} a_{j} \left( \mathbf{x}_{t} \right)$が大きいほど、$\tau$の期待値は小さくなり、事象は頻繁に発生する。

---

ステップ3.

$\displaystyle {{a_{j}} \over {\sum_{j} a_{j} \left( \mathbf{x}_{t} \right)}}$の確率で$j$を選ぶ。$a_{j}$が大きいほど、事象$j$が発生する確率が高くなる。

---

ステップ4.

$\left( t + \tau, \mathbf{x}_{t + \tau} \right)$を記録し、 $$\begin{align*} t \gets & t + \tau \\ \mathbf{x}_{t} \gets & \mathbf{x}_{t + \tau} \end{align*}$$ のように更新する。$\left\{ \mathbf{X}_{t} \right\}$は連続時間マルコフ連鎖なので、ステップ1に戻って繰り返し実行することができる。十分なデータを得たら、アルゴリズムを停止する。

---

出力

与えられた$\left\{ \mathbf{X}_{t} \right\}$の実現を得る。

説明

例

アルゴリズムは式が一般的に書かれていて少し難しそうに見えるけど、実際はとても簡単だ。

簡単な誕生-死滅プロセス^{birth-Death process}を考えてみよう。このプロセスは$N$人の人口から始まり、繁殖係数$b > 0$と死亡係数$d > 0$に従って人口の増減が表される。$b$が$d$より大きいか小さいかで異なる様相が現れるけれども、まず事象が発生する頻度そのものはこれらの絶対値に依存する。生存であれ死であれ、その時点での個体数に比例して頻繁に発生するのが妥当だし、人口数$i$の時の事象が発生する間隔は $$\tau \sim \exp \left( (b+d) i \right)$$ に従う。一度事象が発生する時点を定めたら、次に発生する事象が誕生か死かを選ぶ。$t$時点の人口数が$i$の時、$t + \tau$時点で繁殖するということは$j = i+1$、死ぬということは$j = i-1$を意味する。これらはそれぞれ係数$b, d$と$t$時点の人口数$i$に比例するので、 $$a_{i+1} (x) = b i \\ a_{i-1} (x) = d i$$ ステップ3でこれらを正規化^normalizeして得られる確率でどの事象が発生するかを決定すればいい。

コード

以下は上記の例をそのままコードに転写したJuliaのコードだ。 $$b = 2 \\ d = 1$$ 与えられたパラメータで、初期人口$N=2$から始めてシミュレーションを1万回繰り返し、$t = 2$で平均を計算して理論的な予測と比較する。読みやすさを重視して書いたので、アルゴリズムが例の説明通りに実装されているかを確認しよう。

using Distributions
using Plots

b = 2.0
d = 1.0
θ = b/(b+d)

result = plot(xlims = (0,2), ylims = (0, Inf))

Y = Int64[]
@time for _ in 1:10000
    y = [2]
    t = [0.0]
    while t[end] < 2 && y[end] > 0
        push!(t, t[end] + rand(Exponential(1/((b+d) * y[end]))))
        if rand() < θ
            push!(y, y[end] + 1)
        else
            push!(y, y[end] - 1)
        end
    end
    if t[end] > 2
        push!(Y, y[end-1])
    else
        push!(Y, 0)
    end
    plot!(result, t,y, label = :none, alpha = 0.01, color = :black)
end

plot!(result, 0:0.1:2,2exp.((b-d) .* 0:0.1:2), label = :none, color = :red, linewidth = 4, linestyle = :dot)
result

abs(mean(Y) - 2exp.(2(b-d)))

実行結果を可視化すると上のようになる。理論的な予測値$2\exp(2(b-d)) = 14.78$との差は0.1程度で、かなり正確に計算されたことが確認できる。