Moment-closure approximations and assumptions

Moment closure replaces an unknown higher moment with an explicit approximation based on retained lower moments.

What closure does

A closure replaces an unknown higher moment by a function of retained lower moments:

\[m_3\approx C(m_1,m_2).\]

This converts an infinite hierarchy into a finite ODE system. The resulting equations are approximate even when the original moment identities were exact.

Common second-order closures

Closure	Approximation	Underlying idea
Zero-variance or mean-field	\(m_2\approx m_1^2\)	Distribution concentrated at its mean.
Normal third-moment closure	\(m_3\approx3m_1m_2-2m_1^3\)	Third central moment is zero.
Lognormal closure	\(m_3\approx m_2^3/m_1^3\)	Positive state approximated by a lognormal distribution.

The lognormal expression requires \(m_1>0\). None is universally accurate.

Assumptions must match biology

Near extinction, epidemic distributions can have a point mass at zero and a separate group of large outbreaks. Such a distribution is neither normal nor well represented by a single narrow peak. Closure accuracy can therefore change through time.

Interactive Python laboratory

Use CTMC samples to compare normal and lognormal approximations with the empirical third moment.

Interactive PythonTesting moment closures

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

N, I0 = 100, 5
beta, gamma = 0.25, 0.15
grid = np.linspace(0, 30, 31)
M = 6000
rng = np.random.default_rng(280)

def path_on_grid():
    I, time = I0, 0.0
    event_times, values = [0.0], [I]
    while I > 0 and time < grid[-1]:
        b = beta * (N - I) * I / N
        d = gamma * I
        wait = rng.exponential(1 / (b + d))
        if time + wait > grid[-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, grid, side="right") - 1
    return np.array(values)[pos]

paths = np.array([path_on_grid() for _ in range(M)])
m1 = np.mean(paths, axis=0)
m2 = np.mean(paths ** 2, axis=0)
m3 = np.mean(paths ** 3, axis=0)

normal_m3 = 3 * m1 * m2 - 2 * m1 ** 3
lognormal_m3 = np.where(
    m1 > 0, m2 ** 3 / m1 ** 3, np.nan
)

normal_relative_error = np.abs(normal_m3 - m3) / m3
lognormal_relative_error = np.abs(lognormal_m3 - m3) / m3

selected = [0, 5, 10, 20, 30]
table = pd.DataFrame({
    "time": selected,
    "empirical m3": m3[selected],
    "normal closure": normal_m3[selected],
    "lognormal closure": lognormal_m3[selected],
    "normal relative error": normal_relative_error[selected],
    "lognormal relative error": lognormal_relative_error[selected]
})
print(table.round(4).to_string(index=False))

plt.figure(figsize=(9.5, 5.1))
plt.plot(grid, m3, color="black",
         linewidth=2, label="Empirical m3")
plt.plot(grid, normal_m3, "--",
         label="Normal closure")
plt.plot(grid, lognormal_m3, ":",
         linewidth=2, label="Lognormal closure")
plt.xlabel("Time")
plt.ylabel("Third raw moment")
plt.title("Closure accuracy changes through time")
plt.legend()
plt.grid(alpha=0.25)
plt.tight_layout()

Output

Run the code to see the result.

How to justify a closure

State the formula and distributional assumption.
Check whether it preserves biologically necessary quantities.
Compare closed-moment solutions with simulation moments.
Test several parameter regimes and initial states.
Report where the approximation fails.

What this lesson adds

You can now explain why closure is needed, distinguish common closures, identify their assumptions and evaluate them against simulated higher moments rather than treating closure as an exact identity.