Monte Carlo simulation

A Monte Carlo simulation studies a stochastic model by running it many times. The repeated outcomes are then used to estimate probabilities, expectations, variability and risk.

Core idea. One stochastic trajectory is one possible history. Monte Carlo simulation asks what we learn when we generate many possible histories from the same model.

From one trajectory to many

A Gillespie simulation may produce an outbreak, while another run with exactly the same parameters may become extinct. The model describes a distribution of possible trajectories.

Six independent SIS Gillespie trajectories generated from the same rates. They are correctly drawn as step functions because the infectious count changes only at event times.

The individual paths should not be smoothed. A CTMC remains constant between events and jumps when an infection or recovery occurs.

What Monte Carlo adds

stochastic model→simulate many times→collect outcomes→estimate distribution and risk

If an event \(A\), such as reaching 30 infectious people, occurs in \(M\) of \(n\) independent runs,

\[\boxed{\widehat P(A)=\widehat p=\frac{M}{n}}.\]

If 372 of 1,000 runs satisfy the criterion, the estimated probability is 37.2%.

Why relative frequency estimates probability

For run \(r\), define \(Y_r=1\) if \(A\) occurs and \(Y_r=0\) otherwise. Then

\[\widehat p=\frac1n\sum_{r=1}^{n}Y_r.\]

Because \(E[Y_r]=P(A)\), the average of many independent indicators approaches the underlying probability. This is the law of large numbers.

Estimating an average and variance

If \(Z_r\) is an outcome such as peak infectious count,

\[\bar Z=\frac1n\sum_{r=1}^{n}Z_r,\qquad s^2=\frac1{n-1}\sum_{r=1}^{n}(Z_r-\bar Z)^2.\]

The mean describes the centre; the variance describes the spread.

The distribution can matter more than the mean

Distribution of peak infectious counts from 2,000 independently generated SIS Gillespie simulations. The vertical line marks the Monte Carlo mean peak. Bars are appropriate because the recorded peak count is discrete.

A single mean peak hides the range of plausible epidemic sizes.

Mean trajectory and uncertainty band

At time \(t\), repeated simulations give \(I_1(t),\ldots,I_n(t)\). Their average estimates \(E[I(t)]\), while quantiles show how widely individual outcomes vary.

Monte Carlo meancentral 90% simulation range

Mean infectious count and central 90% range across 2,000 simulations. The mean looks smoother because it averages many step-like paths; the shaded region shows the uncertainty that the mean alone would hide.

The mean trajectory is not the trajectory the epidemic must follow. It is an average across many stochastic histories.

Monte Carlo error

A finite-simulation estimate is itself uncertain. For a probability,

\[\boxed{\operatorname{SE}(\widehat p)\approx\sqrt{\frac{\widehat p(1-\widehat p)}{n}}}.\]

The error decreases approximately as \(1/\sqrt n\), so roughly four times as many runs are needed to halve it.

Watching an estimate converge

running estimatefinal 10,000-run estimateapproximate 95% Monte Carlo interval

Running outbreak-probability estimate from 10,000 independent simulations. Early estimates fluctuate strongly; later fluctuations become much smaller. The reference line and uncertainty band make the convergence clear without artificially smoothing the data.

Why the convergence curve is not perfectly smooth

After \(n\) runs, \(\widehat p_n=M_n/n\). The next run either adds an outbreak or it does not, so the running average changes slightly. These fluctuations are genuine Monte Carlo sampling variation.

A fitted smooth curve would hide the phenomenon we are trying to explain. The correct improvement is more simulations plus a reference estimate and uncertainty band.

Approximate confidence interval

A simple large-sample 95% interval is

\[\widehat p\pm1.96\sqrt{\frac{\widehat p(1-\widehat p)}{n}}.\]

For \(\widehat p=0.40\) and \(n=1000\), this is approximately \((0.370,0.430)\). For small samples or probabilities near 0 or 1, more suitable binomial intervals may be preferred.

Gillespie and Monte Carlo are different

Gillespie	Monte Carlo
generates one CTMC trajectory	uses many stochastic runs
samples event times and event types	estimates probabilities and distributions
one possible history	uncertainty across histories

Biological quantities we can estimate

Outcome	Question
extinction indicator	What is the probability infection dies out?
outbreak indicator	What is the probability a chosen outbreak level is reached?
peak infectious count	How uncertain is epidemic severity?
time of peak	When might maximum pressure occur?
capacity exceedance	What is the probability demand exceeds capacity?

Rare events

If an event has probability about \(10^{-4}\), only about one occurrence is expected per 10,000 ordinary runs. Ordinary Monte Carlo is then very noisy and specialised rare-event methods may be needed.

Independent runs and reproducibility

Runs should explore independent random histories. A seed can make an experiment reproducible, but repeated simulations must still use distinct random-number sequences rather than reproduce the same path.

How many simulations?

There is no universal number. Report Monte Carlo uncertainty and check whether the quantity of interest remains materially stable as more simulations are added.

Monte Carlo versus Kolmogorov equations

Kolmogorov equations	Monte Carlo
evolve state probabilities directly	estimate them from repeated trajectories
no sampling error when solved exactly	finite-simulation sampling error
can become enormous for large state spaces	often practical for complex models

Key idea. Monte Carlo does not make an individual stochastic trajectory smooth. It turns many random trajectories into increasingly stable statistical summaries.

What comes next?

The next section studies extinction and outbreak probabilities: the chance that a small transmission chain disappears versus becoming established.