Programming many independent SDE trajectories

Independent SDE paths turn possible histories into an empirical distribution. A two-dimensional NumPy array makes the ensemble structure explicit.

From one path to an ensemble

Suppose \(M\) independent paths are stored in a rectangular array. Row \(r\) is trajectory \(r\), and column \(n\) is time \(t_n\):

\[X_{r,n}\approx X_r(t_n).\]

Unlike CTMC paths, Euler–Maruyama paths already share the chosen numerical grid, so no separate alignment step is needed.

Independence in the code

At every time step, each path needs its own standard normal draw. The generator is created once before simulation. Reinitialising it with the same seed for each trajectory would create identical paths, not independent ones.

Interactive Python laboratory

Interactive PythonMany independent SIS SDE paths

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

N = 500
I0 = 10
beta = 0.25
gamma = 0.15
T = 80
dt = 0.02
number_of_paths = 300
seed = 83

steps = int(round(T / dt))
times = np.linspace(0, T, steps + 1)
rng = np.random.default_rng(seed)

# Rows are paths; columns are time points
I = np.zeros((number_of_paths, steps + 1))
I[:, 0] = I0
correction_count = 0

for n in range(steps):
    current = I[:, n]

infection_rate = beta * (N - current) * current / N
    recovery_rate = gamma * current
    drift = infection_rate - recovery_rate
    diffusion = np.sqrt(np.maximum(
        infection_rate + recovery_rate, 0
    ))

# One independent normal draw for every path
    Z = rng.normal(0, 1, size=number_of_paths)
    proposals = (
        current
        + drift * dt
        + diffusion * np.sqrt(dt) * Z
    )

corrected = np.clip(proposals, 0, N)
    correction_count += np.count_nonzero(
        corrected != proposals
    )
    I[:, n + 1] = corrected

path_summary = pd.DataFrame({
    "path": np.arange(1, number_of_paths + 1),
    "final I": I[:, -1],
    "minimum I": I.min(axis=1),
    "maximum I": I.max(axis=1)
})

print("Ensemble design")
print(f"Array shape: {I.shape}")
print(f"Independent paths: {number_of_paths}")
print(f"Stored time points: {steps + 1}")
print(f"Total boundary corrections: {correction_count}")
print()
print("First 10 path summaries")
print(path_summary.head(10).round(3).to_string(index=False))

plt.figure(figsize=(9.5, 5.2))
for path_number in range(20):
    plt.plot(
        times, I[path_number],
        alpha=0.42, linewidth=0.8
    )

plt.xlabel("Time (days)")
plt.ylabel("Infectious population")
plt.title("Twenty of 300 independent SIS SDE paths")
plt.grid(alpha=0.25)
plt.tight_layout()

plt.figure(figsize=(8.5, 4.9))
plt.hist(I[:, -1], bins=25,
         color="#79ae90", edgecolor="white")
plt.xlabel("Infectious population at day 80")
plt.ylabel("Number of paths")
plt.title("Distribution across all independent paths")
plt.grid(axis="y", alpha=0.25)
plt.tight_layout()

Output

Run the code to see the result.

Vectorisation

current = I[:, n] selects all paths at one time. NumPy then calculates rates and updates for the whole ensemble together. This is usually clearer and faster than a nested Python loop.

Axis	Meaning
`axis=0`	Summarise down paths at each time.
`axis=1`	Summarise across time within each path.

What the displays establish

The twenty path lines illustrate variability but are not the entire experiment. The histogram uses all 300 paths and describes the empirical distribution at one stated time. Neither alone is an uncertainty band; that is the next lesson.

What this lesson adds

You can now store an ensemble as a path-by-time array, generate independent Brownian draws, vectorise Euler–Maruyama, diagnose boundary corrections and distinguish illustrative paths from a full outcome distribution.