📂Probability Theory

Gillespie Stochastic Simulation Algorithm

Algorithm ¹

For $\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$$ the method for simulating the realization of a continuous Markov chain $\left\{ \mathbf{X}_{t} \right\}$ is called the Gillespie Stochastic Simulation Algorithm, or simply SSA.

---

Input

Receives initial values $t=0$ and $\mathbf{X}_{0} = \mathbf{x}$.

---

Step 1.

Calculate $\sum_{j} a_{j} \left( \mathbf{x}_{t} \right)$.

---

Step 2.

Sample $\tau$ from the exponential distribution $\exp \left( \sum_{j} a_{j} \left( \mathbf{x}_{t} \right) \right)$. The larger $\sum_{j} a_{j} \left( \mathbf{x}_{t} \right)$ is, the smaller the expected value of $\tau$ becomes, and events occur more frequently.

---

Step 3.

Select $j$ with the probability of $\displaystyle {{a_{j}} \over {\sum_{j} a_{j} \left( \mathbf{x}_{t} \right)}}$. The larger $a_{j}$ is, the higher the probability that event $j$ will occur.

---

Step 4.

Record $\left( t + \tau, \mathbf{x}_{t + \tau} \right)$ and update as $$\begin{align*} t \gets & t + \tau \\ \mathbf{x}_{t} \gets & \mathbf{x}_{t + \tau} \end{align*}$$. Since $\left\{ \mathbf{X}_{t} \right\}$ is a continuous-time Markov chain, one can return to Step 1 and repeat. Stop the algorithm when enough data has been obtained.

---

Output

Obtains the realization of given $\left\{ \mathbf{X}_{t} \right\}$.

Explanation

Example

The algorithm might look complicated with all the general expressions, but it’s actually pretty simple.

Consider a simple Birth-Death Process. This process starts with $N$ individuals and shows the increase or decrease in population according to the birth rate $b > 0$ and the death rate $d > 0$. Whether $b$ is greater or smaller than $d$ will determine the pattern, but the frequency of events itself depends on their absolute values. It is plausible that events, whether births or deaths, happen more frequently in proportion to the current number of individuals, and the interval at which events occur when the population is $i$ follows $$\tau \sim \exp \left( (b+d) i \right)$$. Once the timing of an event is decided, we choose which event, birth or death, will occur. At $t$, if the population is $i$, to breed at $t + \tau$ means $j = i+1$, and to die means $j = i-1$. These are proportional to the rate $b, d$ and the population at $t$, $i$, respectively, so $$a_{i+1} (x) = b i \\ a_{i-1} (x) = d i$$. According to Step 3, decide which event will occur with the probability calculated by normalizing these.

Code

The following is the Julia code that directly translates the above example. $$b = 2 \\ d = 1$$ With the given parameters, repeat the simulation ten thousand times starting with the initial population $N=2$ and compute the average at $t = 2$ to compare with the theoretical prediction. It is written to be easy to read rather than efficient, so let’s check if the algorithm is implemented as described in the example.

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)))

Visualizing the result looks like the above. The difference from the theoretical prediction $2\exp(2(b-d)) = 14.78$ is about 0.1, confirming that it has been calculated quite accurately.