Example 20 — SAEM Estimation

Script: examples/20_saem_estimation.py

Demonstrates the Stochastic Approximation EM (SAEM) algorithm on a theophylline-like population PK dataset (12 subjects, 1-cmt oral model).

Key concepts

Concept	Description
Two-phase algorithm	Phase 1 (γ=1, stochastic exploration) + Phase 2 (γ_k→0, convergence)
Metropolis-Hastings E-step	Samples individual ETAs per chain per subject
Closed-form Ω M-step	`Q_Ω ← Q_Ω + γ·(SS_Ω − Q_Ω)`; no gradient needed for Ω
OFV convergence history	Phase-2 parameter window stability criterion

Usage

from openpkpd.estimation.saem import SAEMMethod

result = SAEMMethod(
    n_iter_phase1=300,
    n_iter_phase2=200,
    n_chains=3,
    seed=42,
).estimate(population_model, init_params)

print(f"OFV = {result.ofv:.4f}")
print(f"Converged: {result.converged}")

# Plot convergence
import matplotlib.pyplot as plt
plt.plot(result.ofv_history)
plt.xlabel("Iteration")
plt.ylabel("OFV")
plt.title("SAEM convergence")
plt.show()

Constructor options

Parameter	Default	Description
`n_iter_phase1`	300	Phase-1 stochastic exploration iterations
`n_iter_phase2`	200	Phase-2 convergence iterations
`n_chains`	3	MH chains per subject (Rao-Blackwell averaging)
`n_parallel`	1	Parallel subject processing threads
`mh_step_size`	0.3	Initial MH proposal scale (adaptive)
`seed`	42	RNG seed for reproducibility

Convergence criterion

SAEM declares convergence when the relative change in the full parameter vector φ = [θ, lower_triangle(Ω), diag(Σ)] between successive phase-2 windows (_PH2_WINDOW = 50 iterations) falls below _PH2_TOL = 1e-3.

Comparison with FOCE

The script fits the same dataset with both SAEM and FOCE, printing parameter estimates side by side. SAEM is slower per iteration but does not require the first-order Taylor approximation that FOCE/FOCEI use.