Comparing moment solutions with simulations

Closure quality must be assessed against the stochastic model, not inferred from a smooth ODE solution.

Why comparison is necessary

The closed ODE has numerical error, closure error and possibly model mismatch. Tight solver tolerances address only numerical integration. CTMC simulations provide an independent benchmark for the original stochastic model.

Matched comparison design

Use the same \(N,I_0,\beta,\gamma\).
Compare at the same times.
Use enough independent simulations to reduce Monte Carlo error.
Compare both mean and variance.
Report absolute or relative discrepancies.

Interactive Python laboratory

Interactive PythonClosed moments versus CTMC

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

N, I0 = 100, 5
beta, gamma = 0.25, 0.15
times = np.linspace(0, 30, 61)

def rhs(t, y):
    m1, m2 = y
    m3 = 3 * m1 * m2 - 2 * m1 ** 3
    return [
        (beta - gamma) * m1 - beta * m2 / N,
        (2 * (beta - gamma) - beta / N) * m2
        - 2 * beta * m3 / N
        + (beta + gamma) * m1
    ]

sol = solve_ivp(
    rhs, (0, 30), [I0, I0 ** 2],
    t_eval=times, rtol=1e-9, atol=1e-11
)
moment_mean, moment_second = sol.y
moment_variance = np.maximum(
    moment_second - moment_mean ** 2, 0
)

rng = np.random.default_rng(290)
M = 6000

def ctmc_path():
    I, time = I0, 0.0
    event_times, values = [0.0], [I]
    while I > 0 and time < times[-1]:
        b = beta * (N - I) * I / N
        d = gamma * I
        wait = rng.exponential(1 / (b + d))
        if time + wait > times[-1]:
            break
        time += wait
        I += 1 if rng.random() < b / (b + d) else -1
        event_times.append(time)
        values.append(I)
    pos = np.searchsorted(event_times, times, side="right") - 1
    return np.array(values)[pos]

paths = np.array([ctmc_path() for _ in range(M)])
simulation_mean = paths.mean(axis=0)
simulation_variance = paths.var(axis=0, ddof=1)

mean_error = moment_mean - simulation_mean
variance_error = moment_variance - simulation_variance

selected = [0, 5, 10, 20, 30]
idx = [int(t * 2) for t in selected]
table = pd.DataFrame({
    "time": selected,
    "moment mean": moment_mean[idx],
    "simulation mean": simulation_mean[idx],
    "mean difference": mean_error[idx],
    "moment variance": moment_variance[idx],
    "simulation variance": simulation_variance[idx],
    "variance difference": variance_error[idx]
})
print(table.round(3).to_string(index=False))

plt.figure(figsize=(9.5, 5.1))
plt.plot(times, moment_mean, color="#174c34",
         linewidth=2, label="Closed-moment mean")
plt.plot(times, simulation_mean, "--",
         color="#39815e", label="CTMC sample mean")
plt.xlabel("Time")
plt.ylabel("Expected infectious population")
plt.title("First moment: closure versus simulation")
plt.legend()
plt.grid(alpha=0.25)
plt.tight_layout()

plt.figure(figsize=(9.5, 5.1))
plt.plot(times, moment_variance, color="#b42318",
         linewidth=2, label="Closed-moment variance")
plt.plot(times, simulation_variance, "--",
         color="#d97706", label="CTMC sample variance")
plt.xlabel("Time")
plt.ylabel("Variance of infectious population")
plt.title("Variance: closure discrepancy is visible")
plt.legend()
plt.grid(alpha=0.25)
plt.tight_layout()

Output

Run the code to see the result.

Interpret discrepancies

Small differences can arise from finite Monte Carlo error. Persistent differences larger than that sampling error suggest closure error. Agreement for the mean does not guarantee agreement for variance, tails or extinction probability.

Validation is regime-specific

A closure that performs well for one population size or initial state may fail near extinction, under stronger nonlinearity or for longer times. Validation must cover the biological regimes in which conclusions will be used.

What this lesson adds

You can now design a matched validation experiment, compare closed first and second moments with CTMC sample moments and interpret disagreement as a combination of Monte Carlo and closure error.