📂확률론

길레스피 확률 시뮬레이션 알고리즘

알고리즘 ¹

$\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}이라 부른다.

---

Input

초기값 $t=0$ 과 $\mathbf{X}_{0} = \mathbf{x}$ 을 받는다.

---

Step 1.

$\sum_{j} a_{j} \left( \mathbf{x}_{t} \right)$ 를 계산한다.

---

Step 2.

지수분포 $\exp \left( \sum_{j} a_{j} \left( \mathbf{x}_{t} \right) \right)$ 에서 $\tau$ 를 샘플링한다. $\sum_{j} a_{j} \left( \mathbf{x}_{t} \right)$ 가 크면 클수록 $\tau$ 의 기대값은 작아지고, 사건은 빈번하게 일어나게 된다.

---

Step 3.

$\displaystyle {{a_{j}} \over {\sum_{j} a_{j} \left( \mathbf{x}_{t} \right)}}$ 의 확률로 $j$ 를 선택한다. $a_{j}$ 가 크면 클수록 사건 $j$ 가 일어날 확률이 높아진다.

---

Step 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\}$ 는 연속시간 마코프체인이므로 Step 1.으로 돌아가 반복할 수 있다. 원하는만큼 데이터를 얻으면 알고리즘을 멈춘다.

---

Output

주어진 $\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 $$ 이다. Step 3에 따라 이들을 정규화^normalize해서 구해지는 확률로 어떤 사건이 일어날지 결정하면 된다.

코드

다음은 위 예시를 그대로 코드로 옮긴 줄리아 코드다. $$ b = 2 \\ d = 1 $$ 위의 파라미터가 주어져 있을 때, 초기 인구 $N=2$ 로 시작한 시뮬레이션을 만 번 반복하고 $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정도로 꽤 정확하게 계산되었음을 확인할 수 있다.