run_level = 3
RL = {'SRJF': {'J': 1000, 'J_eval': 1000, 'Nmif': 250, 'n_starts': 21, 'pf_reps': 20}, 'SIRJPF2': {'J': 1000, 'J_eval': 1000, 'Nmif': 250, 'n_starts': 21, 'pf_reps': 10}}Tutorial for PanelPomp Data Analysis using Pypomp: Supplement for “Mechanistic models for panel data: Analysis of ecological experiments with four interacting species”
This tutorial guides readers through implementation, using the Python package pypomp, of the methodology presented by the article, Mechanistic models for panel data: Analysis of ecological experiments with four interacting species, https://arxiv.org/abs/2506.04508.
Introduction
Ecological experiments often yield panel time series data, wherein multiple trajectories are observed across replicated experimental units. These replicates may exhibit correlation through shared parameters governing the underlying stochastic processes. This tutorial demonstrates the construction, estimation, and validation of mechanistic models for panel data using the partially observed Markov process (POMP) framework. Specifically, we focus on PanelPOMP models, which extend single-unit POMP models to collections of related but independent stochastic processes linked through common parameters.
We analyze panel time series data from a controlled mesocosm experiment examining the population dynamics of two freshwater zooplankton species (Daphnia dentifera and D. lumholtzi), an algal food source (Ankistrodesmus falcatus), and a fungal parasite historically identified as Metschnikowia bicuspidata and subsequently described as Australozyma monospora (Lachance et al. 2025). The experiment, conducted by Searle et al. (2016), was designed to investigate how interspecific competition between the native D. dentifera and invasive D. lumholtzi is modified by their different susceptibility to the parasite.
The panel iterated filter (PIF) algorithm (Bretó, Ionides, and King 2020) enables plug-and-play likelihood-based inference for general PanelPOMP models. We employ both standard PIF and its marginalized variant, MPIF, to maximize likelihood functions. MPIF is used below for the SRJF model containing a unit-specific parameter, whereas the all-shared SIRJPF2 model uses standard PIF. Profile likelihood confidence intervals are constructed using the Monte Carlo Adjusted Profile (MCAP) method (Ionides et al. 2017; Ning, Ionides, and Ritov 2021) to account for Monte Carlo uncertainty in likelihood evaluations. General PanelPOMP notation and algorithmic details are provided in Section S2 of the Supplement.
We also evaluate competing model specifications using Akaike’s Information Criterion, \(\text{AIC} = 2p - 2\ell(\hat{\theta})\), where \(p\) is the number of estimated parameters and \(\ell(\hat{\theta})\) is the maximized log-likelihood (Akaike 1974). AIC balances goodness of fit against model complexity.
Finally, we evaluate models through multiple diagnostic approaches. Unit-level log-likelihood contributions help identify experimental units that are poorly described by a shared parameter vector, while simulation-based and convergence diagnostics assess numerical and model adequacy.
We demonstrate the complete PanelPOMP workflow using two models of increasing complexity. Section 1 presents the SRJF model (Susceptible-Removed-Juvenile-Food), which describes single-species dynamics of D. dentifera in the absence of parasites. Section 2 extends this framework to the SIRJPF2 model (Susceptible-Infected-Removed-Juvenile-Parasite-Food-2species), incorporating parasite transmission, disease dynamics, and interspecific competition. Across the tutorial, we demonstrate parameter estimation via PIF and, for a model containing unit-specific parameters, MPIF. We then demonstrate uncertainty quantification via MCAP, model comparison via AIC, and model validation through simulation-based diagnostics.
GPU execution and numerical precision
This tutorial runs on a CUDA GPU and does not fall back to CPU computation. The calculations shown here use an NVIDIA A100. Float64 is required because the negative-binomial likelihood contains differences of large gammaln terms that may lose substantial accuracy in float32. The double_precision switch is therefore True by default. Setting DAPHNIA_DOUBLE_PRECISION=0 selects float32 for the complete document, including estimation, profiles, and diagnostics. A full float32 and float64 comparison must use separate render processes because JAX compilation and arrays retain the precision state in which they were created. The float32 artifact is diagnostic only. The final and subsequent standard renders use float64.
The active run uses run_level = 2. The object RL contains the corresponding particle counts, MIF iterations, independent starts, and likelihood-evaluation replicates. Level 1 is a rapid prototype, level 2 is the tutorial analysis, and level 3 uses the production settings from the R and regular Python tutorials.
Numerical output throughout this document is reported to four significant digits. Log likelihoods and their differences are the one exception: they are of order \(10^2\) to \(10^3\) and are compared at the scale of their Monte Carlo standard errors, so they are given to two decimal places, which preserves those comparisons.
The following calculation evaluates the same negative-binomial log likelihood in float32 and float64 on the active GPU. The difference is numerical rather than statistical. A material discrepancy indicates that float32 is unsuitable for this likelihood calculation.
Complete tutorial precision: float64
precision log_likelihood difference_from_float64
float32 20.74 36.17
float64 -15.43 0Section 1: SRJF Model
Experimental Design and Data Structure
The experimental treatment consists of \(U = 10\) replicated mesocosms containing D. dentifera populations in the absence of parasites. Each unit \(u \in 1{:}U\) was observed at \(N_u = N = 10\) time points, \(t_{u,1{:}N}\), with \(t_{u,n} = 5n + 2\) days for \(n \in 1{:}N\), yielding a balanced panel design. At each observation time, population densities were quantified by microscopic enumeration, producing observations \(y^*_{u,n}\) comprising adult and juvenile counts. As described by Searle et al. (2016), this single-species, parasite-free treatment establishes a baseline for population dynamics prior to introducing parasite-mediated interactions.
The PanelPOMP data structure comprises:
- Units: \(u \in 1{:}10\) independent mesocosms.
- Observation times: \(t_{u,n} = 5n + 2\) days for \(n \in 1{:}10\).
- Observed data: \(y^*_{u,n} = N^n_{S,u,n}\), the count of adult susceptible D. dentifera in a one-litre sample.
- Latent process: \(X_u(t) = (S^n_u(t),\, J^n_u(t),\, F_u(t))^\top\) for \(t \in [t_{u,0}, t_{u,N}]\).
Here \(S^n_u(t)\) denotes adult susceptible density (individuals/litre), \(J^n_u(t)\) denotes juvenile density (individuals/litre), and \(F_u(t)\) represents algal food resource density (\(10^6\) cells/litre). The superscript \(n\) indicates the native species (D. dentifera). The algal resource is unobserved but hypothesised to mediate population dynamics through resource limitation.
Figure 1 displays the observed adult and juvenile densities \(\{y^*_{u,n}\}\) across all units. Adult densities typically attain maximum values during days 20-30 before declining, with substantial inter-unit variability. The temporal structure of population growth and decline motivates an age-structured modeling framework incorporating juvenile maturation and resource-dependent reproduction.
Show data processing and plotting code
# Load the dent-only sheet. The Excel file is in reverse chronological order
# within each rep (day 10 → day 1), so we sort (rep, day) ascending. See
# quality_reports/audits/DATA_SCHEMA.md for the full schema.
xls = pd.ExcelFile('../data/Mesocosmdata.xls')
srjf_raw = xls.parse('dent-only treatments').iloc[0:100].copy()
srjf_raw = srjf_raw.sort_values(['rep', 'day']).reset_index(drop=True)
srjf_raw['day'] = (srjf_raw['day'] - 1) * 5 + 7
srjf_data = srjf_raw[['rep', 'day', 'dent.adult', 'dent.juv']].copy()
trials = ['A', 'B', 'C', 'D', 'E', 'F', 'G', 'H', 'I', 'J']
trial_to_mesocosm = {t: f"Mesocosm {i + 1}" for i, t in enumerate(trials)}
srjf_data['Mesocosm'] = srjf_data['rep'].map(trial_to_mesocosm)
srjf_data['Mesocosm'] = pd.Categorical(
srjf_data['Mesocosm'],
categories=[f"Mesocosm {i + 1}" for i in range(10)],
ordered=True,
)
srjf_data['dent.adult.plot'] = srjf_data['dent.adult']
srjf_data['dent.juv.plot'] = srjf_data['dent.juv']
fig, axes = plt.subplots(
2, 10,
sharex=True,
sharey='row',
figsize=(10, 4),
gridspec_kw={'hspace': 0.15, 'wspace': 0.25},
)
adult_color = PALETTE["adult"]
juv_color = PALETTE["juvenile"]
line_kwargs = dict(linewidth=0.8, linestyle='-')
for j, label in enumerate([f"Mesocosm {i + 1}" for i in range(10)]):
sub = srjf_data[srjf_data['Mesocosm'] == label]
ax_top = axes[0, j]
ax_bot = axes[1, j]
ax_top.plot(sub['day'], sub['dent.adult.plot'], color=adult_color, **line_kwargs)
ax_bot.plot(sub['day'], sub['dent.juv.plot'], color=juv_color, **line_kwargs)
ax_top.set_title(
label.replace('Mesocosm ', 'Mesocosm-'),
fontsize=FONTS["panel_title"],
)
ax_top.set_yscale('function', functions=(np.sqrt, np.square))
ax_bot.set_yscale('function', functions=(np.sqrt, np.square))
for ax in (ax_top, ax_bot):
ax.set_xlim(0, 52)
ax.set_xticks([0, 25, 50])
ax.tick_params(axis='both', labelsize=FONTS["tick"])
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
if j != 0:
ax.tick_params(axis='y', which='both', labelleft=False)
if j == 0:
ax_top.set_ylabel(
'Adult density\n(ind./L)',
fontsize=FONTS["panel_title"],
)
ax_bot.set_ylabel(
'Juvenile density\n(ind./L)',
fontsize=FONTS["panel_title"],
)
ax_bot.set_xlabel('Day', fontsize=FONTS["panel_title"])
fig.align_ylabels(axes[:, 0])
fig.subplots_adjust(left=0.07, right=0.99, top=0.92, bottom=0.13)
plt.show()
Mechanistic Model
Model Specification
The SRJF model describes the dynamics of susceptible adults (\(S^n_u\)), juveniles (\(J^n_u\)), and the latent algal food resource (\(F_u\)), where the superscript \(n\) denotes the native species D. dentifera. Dead or removed individuals transition to the unobserved compartment \(R\).
The latent process satisfies the following system of stochastic differential equations:
\[\begin{aligned} dS^n_u(t) &= \lambda^n_J J^n_u(t)\,dt - (\theta^n_S + \delta) S^n_u(t)\,dt + S^n_u(t)\,d\zeta^n_{S,u}(t), \\ dJ^n_u(t) &= r^n f^n_S F_u(t) S^n_u(t)\,dt - (\theta^n_J + \delta + \lambda^n_J) J^n_u(t)\,dt + J^n_u(t)\,d\zeta^n_{J,u}(t), \\ dF_u(t) &= - f^n_S F_u(t)\,\big(S^n_u(t) + \xi_J J^n_u(t)\big)\,dt - \delta F_u(t)\,dt + \mu\,dt + F_u(t)\,d\zeta_{F,u}(t), \end{aligned}\]
where the stochastic increments are Gaussian white noise:
\[d\zeta^n_{S,u}(t) \sim \mathcal{N}\!\big(0,(\sigma^n_S)^2\,dt\big), \quad d\zeta^n_{J,u}(t) \sim \mathcal{N}\!\big(0,(\sigma^n_J)^2\,dt\big), \quad d\zeta_{F,u}(t) \sim \mathcal{N}\!\big(0,\sigma^2_F\,dt\big).\]
These equations correspond to the SRJF formulation in Section S9 of the Supplement, specialised to the single-species case.
Biological Mechanisms
The SRJF model encodes three coupled ecological processes. Reproduction and maturation are resource-limited: adults produce juveniles at rate \(r^n f^n_S F_u(t) S^n_u(t)\), where the birth efficiency \(r^n\) converts ingested algae into offspring and \(f^n_S\) is the adult filtration rate. Juveniles transition to the adult class at rate \(\lambda^n_J\). Resource competition is exploitative: both adults and juveniles filter algae from the medium, with juveniles consuming at a rate scaled by \(\xi_J\), and algae are replenished at constant rate \(\mu\) following the Searle et al. (2016) protocol. Mortality arises from natural causes (rates \(\theta^n_S\), \(\theta^n_J\)) and from experimental sampling at rate \(\delta\).
The complete flow diagram illustrating these transitions is provided in Figure S20 of the Supplement to the main text.
Parameter Specification
The full parameter table for the SRJF model is given in Tables S10 and S11 of the Supplement (for D. dentifera and D. lumholtzi, respectively). The parameters used here are:
- \(r^n\): birth rate coefficient (juveniles \(\cdot\) individual\(^{-1}\) \(\cdot\) cell\(^{-1}\) \(\cdot\) day\(^{-1}\)).
- \(f^n_S\): adult filtration rate (L \(\cdot\) individual\(^{-1} \cdot\) day\(^{-1}\)).
- \(\xi_J\): juvenile filtration ratio relative to adults (dimensionless).
- \(\lambda^n_J\): juvenile maturation rate (day\(^{-1}\)).
- \(\theta^n_S\), \(\theta^n_J\): natural mortality rates for adults and juveniles (day\(^{-1}\)).
- \(\delta\): sampling/dilution rate (day\(^{-1}\)).
- \(\mu\): algal replenishment rate (\(10^6\) cells \(\cdot\) L\(^{-1} \cdot\) day\(^{-1}\)).
- \(\sigma^n_S\), \(\sigma^n_J\), \(\sigma_F\): process noise standard deviations.
- \(\tau^n_S\): measurement overdispersion parameter (negative binomial).
Several parameters are fixed from the experimental protocol. These are \(\delta = 0.013\) day\(^{-1}\) for one-litre sampling out of a fixed mesocosm volume, \(\mu = 0.37 \times 10^6\) cells \(\cdot\) L\(^{-1} \cdot\) day\(^{-1}\) for twice-weekly algal supplementation, \(\lambda^n_J = 0.1\) day\(^{-1}\) for an approximate 10-day maturation period, and \(\xi_J = 1\) for similar filtration efficiency in juveniles and adults (Ebert 2005; Searle et al. 2016).
As discussed in Section S9 of the Supplement, profile-likelihood identifiability analysis indicates that the adult susceptible noise (\(\sigma^n_S\)) can be fixed at zero without substantively affecting model fit, whereas \(\sigma^n_J\) and \(\sigma_F\) remain statistically resolved. Complete parameter estimates and MCAP confidence intervals are reported in Tables S10 and S11.
PanelPOMP Implementation Framework
The next several subsections define the model’s plug-and-play components: rproc, rinit, dmeas, rmeas, and partrans.
Process Model Simulator
The rproc component simulates the latent SRJF dynamics. This plug-and-play interface requires only the ability to simulate one Euler step from the process model, never to evaluate transition densities. We discretize the continuous-time SDE with step \(\Delta t = 0.25\) days. This is small enough to control discretization bias while keeping the particle-filter and iterated-filtering inner loops practical.
The callable below implements one Euler-Maruyama step. Independent Gaussian innovations with variances \((\sigma^n_S)^2 \Delta t\), \((\sigma^n_J)^2 \Delta t\), and \(\sigma^2_F \Delta t\) are combined with the deterministic increments of the SRJF equations. Following Ebert (2005), the fixed experimental constants \(\delta = 0.013\), \(\mu = 0.37\), \(\lambda^n_J = 0.1\), and \(\xi_J = 1\) are written directly into the function body.
After updating \(S^n_u\), \(J^n_u\), and \(F_u\), the simulator applies the same numerical guard used for the production fits. Any state below zero or above its declared upper bound is reset to zero and increments an error_count accumulator. The measurement model converts any positive error_count into a \(-150\) log-likelihood penalty for that observation. This reset-and-penalise rule is part of the fitted transition model. Replacing it with upper clipping would define a different model. The non-negative quantity \(T^n_S = |S^n_u|\) is exposed for the measurement model.
Show Euler-Maruyama process simulator
# State variables tracked by the process model. `error_count` is reset at every
# observation time via Pypomp's `accumvars` mechanism.
STATENAMES = ["Sn", "Jn", "F", "T_Sn", "error_count"]
# Canonical parameter ordering used throughout this section.
PARAM_NAMES = ["rn", "f_Sn", "theta_Sn", "theta_Jn",
"sigSn", "sigJn", "sigF", "k_Sn"]
def srjf_rproc(X_, theta_, key, covars, t, dt):
"""One Euler-Maruyama step for the SRJF model."""
Sn, Jn, F = X_["Sn"], X_["Jn"], X_["F"]
error_count = X_["error_count"]
sigSn, sigJn, sigF = theta_["sigSn"], theta_["sigJn"], theta_["sigF"]
theta_Sn, theta_Jn = theta_["theta_Sn"], theta_["theta_Jn"]
rn, f_Sn = theta_["rn"], theta_["f_Sn"]
# Fixed experimental constants
delta = 0.013 # sampling/dilution rate (day^-1)
mu_food = 0.37 # algal replenishment (10^6 cells L^-1 day^-1)
lambda_J = 0.1 # juvenile maturation rate (day^-1)
xi_J = 1.0 # juvenile filtration ratio
# Independent Gaussian innovations with SD * sqrt(dt)
k1, k2, k3 = jax.random.split(key, 3)
sqdt = jnp.sqrt(dt)
noiSn = sigSn * sqdt * jax.random.normal(k1)
noiJn = sigJn * sqdt * jax.random.normal(k2)
noiF = sigF * sqdt * jax.random.normal(k3)
# Deterministic + stochastic increments
Sn_term = (lambda_J * Jn * dt
- theta_Sn * Sn * dt
- delta * Sn * dt
+ Sn * noiSn)
Jn_term = (rn * f_Sn * F * Sn * dt
- lambda_J * Jn * dt
- theta_Jn * Jn * dt
- delta * Jn * dt
+ Jn * noiJn)
F_term = (-f_Sn * F * (Sn + xi_J * Jn) * dt
- delta * F * dt
+ mu_food * dt
+ F * noiF)
Sn_new = Sn + Sn_term
Jn_new = Jn + Jn_term
F_new = F + F_term
# Production reset rule, written with jnp.where for JIT compatibility.
Sn_violated = (Sn_new < 0.0) | (Sn_new > 1e5)
Jn_violated = (Jn_new < 0.0) | (Jn_new > 1e5)
F_violated = (F_new < 0.0) | (F_new > 1e20)
Sn_new = jnp.where(Sn_violated, 0.0, Sn_new)
Jn_new = jnp.where(Jn_violated, 0.0, Jn_new)
F_new = jnp.where(F_violated, 0.0, F_new)
error_count_new = (error_count
+ jnp.where(Sn_violated, 1.0, 0.0)
+ jnp.where(Jn_violated, 0.001, 0.0)
+ jnp.where(F_violated, 1000.0, 0.0))
# Non-negative observable consumed by the measurement model.
T_Sn_new = jnp.abs(Sn_new)
return {
"Sn": Sn_new,
"Jn": Jn_new,
"F": F_new,
"T_Sn": T_Sn_new,
"error_count": error_count_new,
}The error_count accumulator implements the soft-constraint mechanism just described: parameter values that frequently produce biologically implausible trajectories accumulate penalties that manifest as reduced likelihood values during optimisation.
Initial State Specification
Particle filters require initial states that are biologically consistent with the experimental protocol. Following the executable R specification, we set \(t_{u,0} = 0\), giving a seven-day pre-observation integration window before the first sample at day \(7\). At \(t_0\), \(S^n_u(t_0) = 3\) ind./L (45 adults released into a 15 L mesocosm), \(J^n_u(t_0) = 0\), and \(F_u(t_0) = 16.667 \times 10^6\) cells/L (an initial \(2.5 \times 10^8\) cells in a 15 L volume). The accumulator error_count and the measurement variable \(T^n_S\) are initialised to zero.
Show initial-state specification
def srjf_rinit(theta_, key, covars, t0):
"""Deterministic initial conditions, identical across units."""
return {
"Sn": jnp.array(3.0),
"Jn": jnp.array(0.0),
"F": jnp.array(16.667),
"T_Sn": jnp.array(0.0),
"error_count": jnp.array(0.0),
}These initial conditions are shared across units. Between-unit variability in the latent state arises solely from the stochastic dynamics in srjf_rproc.
Measurement Model
The measurement model connects the latent adult density \(S^n_u(t_n)\) to the observed count \(N^n_{S,u,n}\) via a negative binomial distribution: \[N^n_{S,u,n} \mid S^n_u(t_n) \sim \text{NBinomial}\big(S^n_u(t_n), \tau^n_S\big),\] parameterised with mean \(\mu = S^n_u(t_n)\) and variance \(\mu + \mu^2/\tau^n_S\). The dispersion \(\tau^n_S > 0\) quantifies overdispersion, with smaller values indicating greater variance relative to the mean. This specification accommodates the substantial count variability common in ecological sampling, arising from both demographic stochasticity and measurement error.
The measurement density evaluates to \[f_{Y_{u,n}|X_{u,n}}(y^*_{u,n} \mid x_{u,n}; \theta) = f_{\text{NBinomial}}\big(y^*_{u,n}; S^n_u(t_n), \tau^n_S\big),\] corresponding to equation (8) of the main text. Trajectories flagged by error_count > 0 have their log-likelihood replaced by the constant penalty \(-150\), which removes them from consideration in the particle filter while avoiding the numerical instabilities associated with \(-\infty\).
We use \(\tau^n_S\) directly as the size parameter, called k_Sn in code, of the NB mean-dispersion parameterisation. The log-pmf is computed via gammaln because jax.scipy.stats.nbinom.logpmf uses the \((n, p)\) parameterisation and expects integer counts, both awkward inside a JIT-compiled particle filter.
Show measurement log-density (negative binomial)
def srjf_dmeas(Y_, X_, theta_, covars, t):
"""Negative binomial log-pmf of dent.adult given latent S^n."""
y = Y_["dentadult"]
mu = jnp.maximum(X_["T_Sn"], 1e-10)
size = jnp.maximum(theta_["k_Sn"], 1e-10)
error_count = X_["error_count"]
ll = (jax.scipy.special.gammaln(y + size)
- jax.scipy.special.gammaln(size)
- jax.scipy.special.gammaln(y + 1.0)
+ size * jnp.log(size / (size + mu))
+ y * jnp.log(mu / (size + mu)))
# Soft penalty: replace ll on biologically implausible trajectories.
return jnp.where(error_count > 0.0, -150.0, ll)The srjf_rmeas component simulates \(\tilde{N}^n_{S,u,n} \sim \text{NBinomial}(S^n_u(t_n), \tau^n_S)\) as a Gamma-Poisson mixture, which JAX supports natively.
Show measurement sampler
def srjf_rmeas(X_, theta_, key, covars, t):
"""Sample dent.adult ~ NBinomial(mean = T_Sn, size = k_Sn)."""
mu = jnp.maximum(X_["T_Sn"], 1e-10)
size = jnp.maximum(theta_["k_Sn"], 1e-10)
k1, k2 = jax.random.split(key)
scale = mu / size
gamma_sample = jax.random.gamma(k1, size) * scale
obs = jax.random.poisson(k2, gamma_sample)
# Pypomp 1.0.0 requires rmeas to return a mapping keyed by the observation
# names of `ys`, rather than the positionally ordered array accepted by
# 0.4.x. Naming the observation removes the ordering assumption.
return {"dentadult": obs.astype(float)}To ensure biological plausibility and improve numerical stability during optimisation, log transformations map the constrained space \(\Theta^+ = \{\theta : \theta_i > 0\}\) to \(\mathbb{R}^p\), so the iterated filter can apply Gaussian perturbations without violating positivity: \[\tilde{\theta}_i = \log(\theta_i) \quad \text{for} \quad \theta_i \in \{r^n, f^n_S, \theta^n_S, \theta^n_J, \sigma^n_J, \sigma_F, \tau^n_S\}.\] During iterated filtering, the perturbation is applied to \(\tilde{\theta}\) via additive Gaussian noise. Parameters are back-transformed to the natural scale via \(\theta_i = \exp(\tilde{\theta}_i)\) before evaluating srjf_rproc or srjf_dmeas. The scheme is scale invariant because relative steps scale with parameter magnitude, and it prevents numerical underflow as parameters approach zero.
Because \(\sigma^n_S \equiv 0\) (see the parameter discussion above), it is not log-transformed. The remaining seven parameters are mapped through jnp.log / jnp.exp.
Show parameter transformations (log)
# Parameters that are strictly positive and benefit from log-transform.
_LOG_PARAMS = ("rn", "f_Sn", "theta_Sn", "theta_Jn",
"sigJn", "sigF", "k_Sn")
def srjf_to_est(theta):
"""Natural -> estimation scale (log for positive params)."""
out = {**theta}
for name in _LOG_PARAMS:
out[name] = jnp.log(jnp.maximum(theta[name], 1e-30))
# sigSn is fixed at 0; pass through untransformed.
out["sigSn"] = theta["sigSn"]
return out
def srjf_from_est(theta):
"""Estimation -> natural scale."""
out = {**theta}
for name in _LOG_PARAMS:
out[name] = jnp.exp(theta[name])
out["sigSn"] = theta["sigSn"]
return out
srjf_par_trans = pp.ParTrans(to_est=srjf_to_est, from_est=srjf_from_est)We construct one Pomp object per replicate \(u \in \{A, B, \ldots, J\}\). Each unit shares the same model components (rinit, rproc, dmeasure, rmeasure, partrans) and differs only in its observed time series. The sampling grid has ten observations at days \(7, 12, \ldots, 52\) in every replicate. The simulation initial time \(t_0 = 0\) leaves a seven-day pre-observation integration window before the first likelihood evaluation. The integration step \(\Delta t = 0.25\) days is passed to pp.Pomp via dt=0.25, and error_count is reset to zero at every observation time via the accumvars argument.
Following the model selection analysis in Section S9 of the Supplement (Tables S10 to S12), all SRJF parameters are shared across units: \(\theta = \phi\) with \(\psi_u = \emptyset\) for all \(u\).
Code
# Shared SRJF starting vector.
srjf_shared_theta = {
"rn": 1.535539e+03,
"f_Sn": 1.306857e-04,
"theta_Sn": 6.353239e-01,
"theta_Jn": 1.376217e-03,
"sigSn": 0.0,
"sigJn": 3.018772e-01,
"sigF": 8.658726e-07,
"k_Sn": 1.417860e+01,
}
unit_names = ['A', 'B', 'C', 'D', 'E', 'F', 'G', 'H', 'I', 'J']
t0_srjf = 0.0 # First observation is at day 7.
srjf_pomp_dict = {}
for u in unit_names:
sub = (srjf_data[srjf_data['rep'] == u][['day', 'dent.adult']]
.rename(columns={'dent.adult': 'dentadult'})
.sort_values('day'))
ys_u = sub.set_index('day')[['dentadult']].astype(float)
srjf_pomp_dict[u] = pp.Pomp(
ys=ys_u,
theta=pp.PompParameters(srjf_shared_theta),
statenames=STATENAMES,
t0=t0_srjf,
rinit=srjf_rinit,
rproc=srjf_rproc,
dmeas=srjf_dmeas,
rmeas=srjf_rmeas,
par_trans=srjf_par_trans,
dt=0.25,
accumvars=("error_count",),
)In Pypomp, the PanelPomp constructor expects a theta dictionary with 'shared' and 'unit_specific' DataFrames. With an all-shared parameterisation there are no unit-specific entries, but the PanelPomp validator reads unit names from the columns of unit_specific. We therefore pass an empty DataFrame whose columns are the unit names (zero rows) to register the panel structure. This is an edge case of the all-shared parameterisation.
Show PanelPomp object construction
# Shared parameters: phi (all units use identical values).
srjf_shared_df = pd.DataFrame(
{"shared": list(srjf_shared_theta.values())},
index=list(srjf_shared_theta.keys()),
)
# Empty unit_specific DataFrame whose columns register the unit names.
# Pypomp reads unit names from `unit_specific.columns`; passing `None`
# causes the validator to return an empty list and refuse the panel.
srjf_unit_specific_df = pd.DataFrame(index=[], columns=unit_names)
panelfood = pp.PanelPomp(
pomp_dict=srjf_pomp_dict,
theta=pp.PanelParameters(
{"shared": srjf_shared_df, "unit_specific": srjf_unit_specific_df}
),
)
print(f"Number of units: {len(panelfood.unit_objects)}")
print(f"Shared parameters: {panelfood.canonical_shared_param_names}")Number of units: 10
Shared parameters: ['rn', 'f_Sn', 'theta_Sn', 'theta_Jn', 'sigSn', 'sigJn', 'sigF', 'k_Sn']The resulting panelfood object represents the complete PanelPOMP model with \(U = 10\) units, each evolving under SRJF dynamics with common parameter vector \(\phi\). Eight model entries are carried, but sigSn = 0 is fixed, so the fitted all-shared dimension is \(p = 7\).
Parameter Estimation via Panel Iterated Filtering
Panel iterated filtering (PIF) extends IF2 to panel data by coupling a coordinated random walk in parameter space with simultaneous filtering of all \(U\) units (Bretó, Ionides, and King 2020). At iteration \(m \in 1{:}M\), each particle \(j\) carries a parameter value \(\phi^{(m)}_j\) that is perturbed before each observation by a Gaussian random walk with increments whose variance is multiplied by a geometric cooling factor, \(a^{2m/50}\), with \(a < 1\) being the cooling fraction in 50 iterations. Pseudocode for PIF and its marginalised variant MPIF appears as Section S2 of the Supplement.
The compute budget is set by the RUN_LEVELS table in the global setup chunk. Its three levels match the particle counts, likelihood replicates, and MIF iterations in the R and regular Python tutorials. Run level 2 is used here.
Code
# All MIF/pfilter calls in this section read their compute budget from RL,
# the run-level dictionary set in the global-setup chunk. Editing run_level
# at the top of the document propagates here without further changes.
J_mif = RL_SRJF["J"]
M_mif = RL_SRJF["Nmif"]
Nstarts = RL_SRJF["n_starts"]
J_pf_eval = RL_SRJF["J_eval"]
Npf_reps = RL_SRJF["pf_reps"]
print(f"MIF compute budget: J={J_mif}, M={M_mif}, Nstarts={Nstarts}")
print(f"Final pfilter eval: J={J_pf_eval}, reps={Npf_reps}")MIF compute budget: J=1000, M=250, Nstarts=21
Final pfilter eval: J=1000, reps=20The random-walk perturbation \(\sigma_{\mathrm{rw}}\) controls parameter exploration during MIF. Larger values promote exploration but slow convergence and can prevent settling near a local maximum. Smaller values promote local refinement but risk premature convergence. Because srjf_par_trans log-transforms every positive parameter, the perturbation is applied on the log scale, so a single value \(\sigma_{\mathrm{rw}} = 0.02\) produces proportionally comparable steps for parameters spanning many orders of magnitude (e.g., \(r^n \sim 10^3\) vs. \(f^n_S \sim 10^{-4}\)). We pass an RWSigma object listing the seven log-perturbed parameters. The parameter sigSn is held at zero (Section S9 of the Supplement) so its random-walk standard deviation is also zero.
Code
dent_rw_sd = 0.02
# RWSigma must list every canonical parameter (Pypomp checks set equality).
# Setting sigSn=0 keeps it fixed at its initial value during MIF.
srjf_rw_sd = pp.RWSigma(
sigmas={
"rn": dent_rw_sd,
"f_Sn": dent_rw_sd,
"theta_Sn": dent_rw_sd,
"theta_Jn": dent_rw_sd,
"sigSn": 0.0,
"sigJn": dent_rw_sd,
"sigF": dent_rw_sd,
"k_Sn": dent_rw_sd,
},
init_names=[],
).geometric_cooling(0.7)Maximum likelihood estimation for PanelPOMP models is a non-convex problem with multiple local maxima. We compare three labeled starting-point qualities, spaced so that the groups differ on the likelihood surface and not merely in the parameter vector. Each start is the R tutorial vector perturbed on the log scale: medium starts use standard deviation 0.45, poor starts 1.6, and extreme starts 3.0. An extreme start moves a typical parameter by roughly a factor of twenty. The three values were chosen by measurement rather than by eye: twelve starts drawn at each, evaluated by particle filter before any MIF, had median panel log likelihoods of about \(-568\), \(-1{,}865\) and \(-10{,}400\), so each group sits about half a decade of magnitude below the one before it. The extreme group is included to ask whether this compute budget recovers a start that far away, and it is allowed to answer no. A comparison in which every group converges shows only that the search is robust to the displacements tried; it takes a group that struggles to locate where the method’s limit lies. Whether the extreme group is that case is a question for the run, not an assumption: the tables and figures below report what happened, and the text does not presume the answer. Every chain uses the same MIF hyperparameters and a distinct random-number stream. The R tutorial vector itself is not one of these starts; it is evaluated on its own in the chunk below, so the candidate MLE is the product of the displaced searches alone. The run with the highest independently evaluated final likelihood is taken as that candidate.
In Pypomp, this multi-start search is expressed by replicating the panel-parameter object before calling mif: each replicate of panel.theta becomes one independent chain, and a single call to panel.mif(...) advances all chains in parallel under JAX. After MIF, a final particle filter on the terminal estimates yields a low-variance evaluation of \(\log L(\hat\phi^{(M)})\).
The current Pypomp API supports a genuinely all-shared parameterization. The fixed parameter sigSn remains in the shared block with random-walk standard deviation zero. It is neither estimated nor used as an artificial unit-specific placeholder.
Code
# Build a genuinely all-shared parameter object. sigSn is carried as a fixed
# model entry but excluded from the free-parameter count.
shared_keys_srjf = list(srjf_shared_theta)
srjf_panel_shared_df = pd.DataFrame(
{"shared": [srjf_shared_theta[k] for k in shared_keys_srjf]},
index=shared_keys_srjf,
)
srjf_panel_unit_df = pd.DataFrame(index=[], columns=unit_names)
srjf_start_qualities = start_quality_schedule(Nstarts)
srjf_panel_theta_multi, srjf_start_qualities = make_quality_panel_starts(
srjf_panel_shared_df,
srjf_panel_unit_df,
seed=100,
fixed=("sigSn",),
qualities=srjf_start_qualities,
)
def _evaluate_srjf_start(budget=None, key_seed=0):
"""Evaluate the R tutorial starting vector before MIF."""
budget = RL_SRJF if budget is None else budget
panel = pp.PanelPomp(
pomp_dict=srjf_pomp_dict,
theta=pp.PanelParameters(panelfood.theta),
)
panel.pfilter(
J=budget["J_eval"], reps=budget["pf_reps"], key=jax.random.key(key_seed),
)
return np.asarray(panel.results_history[-1].logLiks.values)
def _run_srjf_mif(budget=None, starts=None, key_seed=100):
"""Run MIF + final pfilter and return what the figures and tables need."""
budget = RL_SRJF if budget is None else budget
starts = srjf_panel_theta_multi if starts is None else starts
panel = pp.PanelPomp(
pomp_dict=srjf_pomp_dict, theta=pp.PanelParameters(starts),
)
panel.pfilter(
J=budget["J_eval"], reps=budget["pf_reps"], key=jax.random.key(key_seed),
)
initial_ll = np.asarray(panel.results_history[-1].logLiks.values)
panel.mif(
J=budget["J"], M=budget["Nmif"], rw_sd=srjf_rw_sd,
block=False,
key=jax.random.key(key_seed + 1),
)
panel.pfilter(
J=budget["J_eval"], reps=budget["pf_reps"], key=jax.random.key(key_seed + 102),
)
return {
"initial_ll": initial_ll,
"ll_arr": np.asarray(panel.results_history[-1].logLiks.values),
"traces": panel.traces(),
"final_theta": panel.theta.params(as_list=True),
}
def _finite_argmax(values):
"""Return the best finite index, or stop if no result is usable."""
values = np.asarray(values, dtype=float)
finite = np.flatnonzero(np.isfinite(values))
if not len(finite):
raise ValueError("Cannot select a best result: all values are non-finite.")
return int(finite[np.argmax(values[finite])])Saving computations with bake and stew
bake and stew save computations so they can be reused. Both execute Python code supplied as a string when no matching archive exists. On later calls, they load the saved results if the code, declared dependencies, and seed settings match. bake returns one result; stew restores the variables assigned by its code block and returns their names. Here, bake saves the initial particle-filter likelihood array, while stew saves the MIF results used by the tables and figures. We use compatibility="R" to select the tested R-style archive behavior; the code and random draws remain Python/JAX.
In each call, file names the archive, dir selects its directory, and dependson=CACHE_METADATA tracks the tutorial source, data, software versions, precision, and computation settings. The helpers retain their explicit JAX keys. Because timing is enabled, R-mode bake returns an ArchiveValue: .value contains the likelihood array and .system_time["elapsed"] records computation time in seconds.
srjf_start_archive = pp.bake(
file="srjf_start_pfilter.bake.pkl",
expr="_evaluate_srjf_start()",
dir=CACHE_DIR,
dependson=CACHE_METADATA,
compatibility="R",
)
srjf_start_ll_arr = srjf_start_archive.value
TIMINGS["srjf_start_pfilter"] = srjf_start_archive.system_time["elapsed"]stew restores the following assigned variables into the notebook workspace. The mif_out bundle includes initial and final likelihood evaluations; the additional names give direct access to the final likelihood array, convergence traces, and parameter estimates. Its timing metadata is restored as _system_time. For both functions, a loaded timing describes the original computation, not the time spent loading the file.
SRJF_STEW_EXPR = """
mif_out = _run_srjf_mif()
mif_ll_arr = mif_out["ll_arr"]
mif_traces = mif_out["traces"]
srjf_final_theta = mif_out["final_theta"]
"""
srjf_saved_names = pp.stew(
file="srjf_mif_search.stew.pkl",
expr=SRJF_STEW_EXPR,
dir=CACHE_DIR,
dependson=CACHE_METADATA,
compatibility="R",
)
TIMINGS["srjf_mif_search"] = _system_time["elapsed"]The document also enables Quarto’s execution cache with execute.cache: true. That cache can reuse notebook outputs without executing Python; bake and stew check their archives when Python runs. An initial computation is still required. After changing external data, the installed package, or environment settings, use quarto render daphnia_tut_pypomp_advanced.qmd --cache-refresh so the dependency checks execute. To also recompute every saved analysis, set DAPHNIA_FORCE_RECOMPUTE=1 for that command. The existing provenance checks conservatively invalidate analysis archives whenever the QMD changes. See the installation and cache instructions for the required package revision and jupyter-cache dependency.
Measuring computation, saving, and reuse
An archive saves time by avoiding repeated computation; it does not accelerate the particle filter or MIF itself. The first call computes the result and writes it to disk, which adds overhead. A matching later call reads that result instead. The following experiment measures all three paths: direct computation without an archive, computation plus saving to a new archive, and loading that archive. It also checks that loading returns the complete saved result and does not call the computation helper.
The demonstration uses the same SRJF model and helpers as the analysis, with a separate run-level-1 budget: 50 particles, two MIF iterations, three starts, and two likelihood replicates. This keeps the demonstration small even when the main analysis uses run level 3. We time one warm-up call separately, including any compilation still needed. We then report medians over three direct/save/load trials with identical inputs and explicit JAX keys. Each save uses a new temporary file, and each load uses a fresh workspace. The timer covers the complete call, including dependency hashing, serialization, and file I/O. Results are materialized as NumPy arrays and pandas tables before returning, and the timer also waits for any remaining JAX array results.
Code
import sys
import tempfile
def _archive_max_error(expected, actual):
"""Check every structure/label/dtype; return the largest numeric difference."""
assert type(expected) is type(actual), "Restored object type changed"
if isinstance(expected, pd.DataFrame):
pd.testing.assert_index_equal(expected.index, actual.index, exact=True)
pd.testing.assert_index_equal(expected.columns, actual.columns, exact=True)
assert expected.dtypes.equals(actual.dtypes), "Table dtypes changed"
return max((_archive_max_error(expected[c].to_numpy(), actual[c].to_numpy())
for c in expected.columns), default=0.0)
if isinstance(expected, np.ndarray):
assert expected.shape == actual.shape, "Array shape changed"
assert expected.dtype == actual.dtype, "Array dtype changed"
if expected.dtype.kind not in "biufc":
np.testing.assert_array_equal(expected, actual)
return 0.0
same = (expected == actual) | (np.isnan(expected) & np.isnan(actual))
if np.all(same):
return 0.0
if not np.all(np.isfinite(expected[~same]) & np.isfinite(actual[~same])):
return float("inf")
difference = np.subtract(expected[~same], actual[~same],
dtype=np.result_type(expected.dtype, np.float64))
return float(np.max(np.abs(difference)))
if isinstance(expected, dict):
assert expected.keys() == actual.keys(), "Saved variable/parameter names changed"
return max((_archive_max_error(expected[k], actual[k]) for k in expected),
default=0.0)
if isinstance(expected, (list, tuple)):
assert len(expected) == len(actual), "Number of saved items changed"
return max((_archive_max_error(a, b) for a, b in zip(expected, actual)),
default=0.0)
assert expected == actual, "Saved nonnumeric value changed"
return 0.0
def _archive_timed(call):
started = time.perf_counter()
result = call()
jax.block_until_ready(result)
return result, time.perf_counter() - started
DEMO_STEW_EXPR = SRJF_STEW_EXPR.replace("_run_srjf_mif()", "_counted_compute()")
SRJF_SAVED_NAMES = {"mif_out", "mif_ll_arr", "mif_traces", "srjf_final_theta"}
def _demo_archive(operation, filename, directory, compute, dependencies):
# A new namespace on every call prevents old variables from masking a miss.
scope = {"_counted_compute": compute}
options = dict(file=filename, dir=directory, envir=scope,
dependson=dependencies, compatibility="R")
if operation == "bake":
return pp.bake(expr="_counted_compute()", **options).value
names = pp.stew(expr=DEMO_STEW_EXPR, **options)
return {name: scope[name] for name in names}Code
DEMO_BUDGET = dict(RUN_LEVELS[1])
demo_starts, _ = make_quality_panel_starts(
srjf_panel_shared_df, srjf_panel_unit_df, seed=100, fixed=("sigSn",),
qualities=start_quality_schedule(DEMO_BUDGET["n_starts"]),
)
demo_timing_rows, demo_check_rows, demo_invalidation_rows = [], [], []
# Only these scratch archives are removed afterwards; analysis archives persist.
with tempfile.TemporaryDirectory(prefix="archive-demo-", dir=CACHE_DIR) as demo_dir:
for operation in ("bake", "stew"):
demo_seed = 0 if operation == "bake" else 100
helper_calls = 0
def _counted_compute():
global helper_calls
helper_calls += 1
if operation == "bake":
return _evaluate_srjf_start(budget=DEMO_BUDGET, key_seed=demo_seed)
return _run_srjf_mif(budget=DEMO_BUDGET, starts=demo_starts,
key_seed=demo_seed)
dependencies = dict(tutorial=CACHE_METADATA, budget=DEMO_BUDGET,
start_seed=100, key_seed=demo_seed)
_, warmup_seconds = _archive_timed(_counted_compute)
elapsed = {path: [] for path in ("direct", "save", "load")}
calls = {path: 0 for path in elapsed}
direct_error = 0.0
for trial in range(3):
filename = f"{operation}-{trial}.pkl"
# 1. No bake/stew: run the actual computation each time.
before = helper_calls
direct, seconds = _archive_timed(_counted_compute)
elapsed["direct"].append(seconds)
calls["direct"] += helper_calls - before
# 2. New archive: compute the same result and save it.
before = helper_calls
saved, seconds = _archive_timed(lambda: _demo_archive(
operation, filename, demo_dir, _counted_compute, dependencies))
elapsed["save"].append(seconds)
calls["save"] += helper_calls - before
# 3. Same archive and dependencies: restore into a fresh workspace.
before = helper_calls
loaded, seconds = _archive_timed(lambda: _demo_archive(
operation, filename, demo_dir, _counted_compute, dependencies))
elapsed["load"].append(seconds)
calls["load"] += helper_calls - before
if operation == "stew":
assert saved.keys() == SRJF_SAVED_NAMES
assert _archive_max_error(saved, loaded) == 0.0
direct_error = max(direct_error, _archive_max_error(
direct, saved if operation == "bake" else saved["mif_out"]))
assert calls == {"direct": 3, "save": 3, "load": 0}
med = {path: float(np.median(values)) for path, values in elapsed.items()}
demo_timing_rows.append({
"archive": operation, "warm-up (s)": warmup_seconds,
"direct (s)": med["direct"], "compute + save (s)": med["save"],
"load (s)": med["load"], "direct / load": med["direct"] / med["load"],
})
demo_check_rows.append({
"archive": operation, "saved = loaded": True,
"max |direct - saved|": direct_error,
"helper calls: direct/save/load": f"{calls['direct']}/{calls['save']}/{calls['load']}",
})
# Change a real input to the helper and declare it in dependson.
# Reuse the same filename and expression: the dependency change must miss.
demo_seed += 1
dependencies = dict(dependencies, key_seed=demo_seed)
before = helper_calls
changed = _demo_archive(
operation, filename, demo_dir, _counted_compute, dependencies)
changed_calls = helper_calls - before
before = helper_calls
reloaded = _demo_archive(
operation, filename, demo_dir, _counted_compute, dependencies)
reload_calls = helper_calls - before
assert (changed_calls, reload_calls) == (1, 0)
assert _archive_max_error(changed, reloaded) == 0.0
changed_direct = _counted_compute()
demo_invalidation_rows.append({
"archive": operation, "new JAX key seed": demo_seed,
"helper calls: changed/repeat": f"{changed_calls}/{reload_calls}",
"output changed": _archive_max_error(saved, changed) > 0.0,
"max |new direct - saved|": _archive_max_error(
changed_direct, changed if operation == "bake" else changed["mif_out"]),
"saved = reloaded": True,
})
print(f"Archive demonstration: {jax.default_backend()}, "
f"{'float64' if double_precision else 'float32'}, budget={DEMO_BUDGET}")
print("\nWall times in seconds; medians of three trials after one warm-up:")
print(table_text(pd.DataFrame(demo_timing_rows)))
print("\nComplete output checks (across all three trials):")
print(table_text(pd.DataFrame(demo_check_rows)))
for row in demo_timing_rows:
speedup = row["direct / load"]
if speedup > 1:
print(f"{row['archive']}: loading was {speedup:.4g} times faster than "
"direct computation in this run.")
else:
print(f"{row['archive']}: loading did not save time in this run "
f"(direct/load = {speedup:.4g}).")Archive demonstration: gpu, float64, budget={'J': 50, 'J_eval': 50, 'Nmif': 2, 'n_starts': 3, 'pf_reps': 2}
Wall times in seconds; medians of three trials after one warm-up:
archive warm-up (s) direct (s) compute + save (s) load (s) direct / load
bake 3.906 0.1455 0.1428 0.001838 79.16
stew 34.66 35.18 34.46 0.002701 1.302e+04
Complete output checks (across all three trials):
archive saved = loaded max |direct - saved| helper calls: direct/save/load
bake True 0 3/3/0
stew True 0 3/3/0
bake: loading was 79.16 times faster than direct computation in this run.
stew: loading was 1.302e+04 times faster than direct computation in this run.The equality checks cover every initial and final likelihood entry, every trace row, and every final shared and unit-specific parameter table, including names, shapes, index/column labels, and dtypes. Saved-versus-loaded equality is required exactly. The separate maximum direct-versus-saved difference measures reproducibility of a new computation with the same JAX keys; zero means numerical equality, while a nonzero value is reported rather than concealed. Float32 and float64 renders perform these checks separately; results across precisions need not be equal. These checks establish faithful reuse of results; likelihood accuracy and optimizer convergence are assessed in the scientific diagnostics below.
The ratio direct / load measures the benefit of reuse for this workload on the current machine. The first save still pays the computation cost, and small results or slow storage can make loading less beneficial. These repeated in-process loads can benefit from the operating system’s file cache; they do not measure a cold filesystem or the full run-level-3 speedup. Quarto may reuse the displayed measurements on later renders. Use --cache-refresh to measure them again; the demonstration always creates new scratch archives whenever its Python cells execute.
Changing an input invalidates its archive
The experiment also changes the JAX key seed used by each real computation and updates that value in dependson, while retaining the same filename and expression. The first call must run the helper once, and a repeated call must run it zero times. We compare the new saved result against a direct computation with the new key, then verify exact restoration again. Archive functions do not automatically discover changes inside referenced helpers or external data: declaring those inputs is essential. The main analysis includes the QMD and data hashes in CACHE_METADATA for this reason.
archive new JAX key seed helper calls: changed/repeat output changed max |new direct - saved| saved = reloaded
bake 1 1/0 True 0 True
stew 101 1/0 True 0 TrueRestoring the full analysis in a fresh Python process
The small benchmark above is supplemented by a reload of the actual analysis archives created earlier, at the selected run level. The child process imports Pypomp but has no SRJF model or computation helpers; a cache miss would therefore fail. It uses the CPU backend to load the saved NumPy arrays and pandas tables without reserving additional GPU memory; their stored precision is preserved, and no SRJF computation runs in the child. Its code is embedded here, so no separate Python script is required. It returns both the bake likelihood array and all four stew variables for complete equality checks against the current analysis. The times below measure the two archive calls inside that process; the separately reported round trip also includes Python startup, imports, and result transfer.
Code
fresh_process_code = r'''
import pickle
import sys
import time
import pypomp as pp
request = pickle.load(sys.stdin.buffer)
options = dict(dir=request["directory"], dependson=request["dependencies"],
compatibility="R")
started = time.perf_counter()
initial = pp.bake(file="srjf_start_pfilter.bake.pkl",
expr="_evaluate_srjf_start()", envir={}, **options).value
bake_seconds = time.perf_counter() - started
scope = {}
started = time.perf_counter()
names = pp.stew(file="srjf_mif_search.stew.pkl", expr=request["stew_expr"],
envir=scope, **options)
stew_seconds = time.perf_counter() - started
pickle.dump(dict(initial=initial, names=names,
outputs={name: scope[name] for name in names},
bake_seconds=bake_seconds, stew_seconds=stew_seconds),
sys.stdout.buffer)
'''
request_bytes = pickle.dumps(dict(directory=str(CACHE_DIR.resolve()),
dependencies=CACHE_METADATA,
stew_expr=SRJF_STEW_EXPR))
started = time.perf_counter()
child = subprocess.run([sys.executable, "-c", fresh_process_code],
input=request_bytes, capture_output=True,
env={**os.environ, "JAX_PLATFORMS": "cpu"})
fresh_process_seconds = time.perf_counter() - started
if child.returncode:
raise RuntimeError("Fresh-process archive reload failed:\n" + child.stderr.decode())
restored = pickle.loads(child.stdout)
assert set(restored["names"]) == SRJF_SAVED_NAMES
assert _archive_max_error(srjf_start_ll_arr, restored["initial"]) == 0.0
analysis_outputs = dict(mif_out=mif_out, mif_ll_arr=mif_ll_arr,
mif_traces=mif_traces, srjf_final_theta=srjf_final_theta)
assert _archive_max_error(analysis_outputs, restored["outputs"]) == 0.0
print(f"Full analysis archives: run level {run_level}, "
f"{'float64' if double_precision else 'float32'}")
print(table_text(pd.DataFrame([
{"archive": "bake", "load (s)": restored["bake_seconds"],
"all outputs restored exactly": True, "computation helpers available": False},
{"archive": "stew", "load (s)": restored["stew_seconds"],
"all outputs restored exactly": True, "computation helpers available": False},
])))
print(f"bake: likelihood array {restored['initial'].shape}, "
f"dtype={restored['initial'].dtype}")
print("stew: " + ", ".join(restored["names"]))
print(f"Fresh-process round trip: {fresh_process_seconds:.4g} s")
del restored, analysis_outputs, child, request_bytesFull analysis archives: run level 3, float64
archive load (s) all outputs restored exactly computation helpers available
bake 0.001831 True False
stew 0.01231 True False
bake: likelihood array (1, 10, 20), dtype=float64
stew: mif_ll_arr, mif_out, mif_traces, srjf_final_theta
Fresh-process round trip: 7.813 sCode
(
ll_per_unit_all,
ll_per_unit_se_all,
panel_ll_all,
panel_ll_se_all,
) = summarize_panel_loglik(srjf_start_ll_arr)
_, _, srjf_initial_ll, srjf_initial_se = summarize_panel_loglik(
mif_out["initial_ll"]
)
unit_ll_per_start, unit_se_per_start, panel_ll_per_start, panel_se_per_start = (
summarize_panel_loglik(mif_ll_arr)
)
best_idx = _finite_argmax(panel_ll_per_start)
best_ll = float(panel_ll_per_start[best_idx])
best_ll_se = float(panel_se_per_start[best_idx])
unit_ll_best = unit_ll_per_start[best_idx]
best_theta_dict = srjf_final_theta[best_idx]
mle_shared_df = best_theta_dict["shared"]
srjf_mif_candidate_theta = {**srjf_shared_theta}
for name in mle_shared_df.index:
srjf_mif_candidate_theta[name] = float(mle_shared_df.loc[name, "shared"])
srjf_mif_candidate_ll = best_ll
srjf_mif_candidate_se = best_ll_se
srjf_start_ll = float(panel_ll_all[0])
srjf_start_se = float(panel_ll_se_all[0])
mif_summary = pd.DataFrame({
"start": np.arange(Nstarts),
"quality": srjf_start_qualities,
"initial_logLik": srjf_initial_ll,
"panel_logLik": panel_ll_per_start,
"MCSE": panel_se_per_start,
"improvement": panel_ll_per_start - srjf_initial_ll,
})
srjf_quality_summary = (
mif_summary.groupby("quality", sort=False)
.agg(
starts=("start", "count"),
mean_initial_logLik=("initial_logLik", "mean"),
mean_terminal_logLik=("panel_logLik", "mean"),
best_terminal_logLik=("panel_logLik", "max"),
mean_improvement=("improvement", "mean"),
)
.reset_index()
)
agreement_tol = np.maximum(
1.0,
2.0 * np.sqrt(
panel_se_per_start**2 + srjf_mif_candidate_se**2
),
)
srjf_near_best_count = int(np.sum(
srjf_mif_candidate_ll - panel_ll_per_start <= agreement_tol
))
srjf_mle_theta = dict(srjf_mif_candidate_theta)
srjf_change_from_start = srjf_mif_candidate_ll - srjf_start_ll
srjf_start_combined_se = float(np.sqrt(
srjf_mif_candidate_se**2 + srjf_start_se**2
))
srjf_mle_panel_theta = pp.PanelParameters({
"shared": pd.DataFrame(
{"shared": [srjf_mle_theta[k] for k in shared_keys_srjf]},
index=shared_keys_srjf,
),
"unit_specific": pd.DataFrame(index=[], columns=unit_names),
})
print(f"MIF + final pfilter wall: {sf(TIMINGS['srjf_mif_search'])}s")
print(table_text(mif_summary))
print("\nSummary by starting-point quality:")
print(table_text(srjf_quality_summary))
print(f"Best MIF start: {best_idx}, panel logLik = "
f"{ll(srjf_mif_candidate_ll)} (MCSE {sf(srjf_mif_candidate_se)})")
print(
"Starts within the stated MCSE-based tolerance of the best: "
f"{srjf_near_best_count} of {Nstarts}"
)
print(
"Best MIF minus R starting-vector logLik: "
f"{ll(srjf_change_from_start, signed=True)} "
f"(combined MCSE {sf(srjf_start_combined_se)})"
)
srjf_reference_ll = -498.7808
srjf_reference_se = 0.324
srjf_reference_gap_se = float(np.sqrt(
srjf_mif_candidate_se**2 + srjf_reference_se**2
))
srjf_reference_difference = srjf_mif_candidate_ll - srjf_reference_ll
print(
"Best MIF minus stored reference logLik: "
f"{ll(srjf_reference_difference, signed=True)} "
f"(combined MCSE {sf(srjf_reference_gap_se)})"
)MIF + final pfilter wall: 179.6s
start quality initial_logLik panel_logLik MCSE improvement
0 medium -634.96 -498.28 0.6324 136.67
1 poor -3979.59 -501.04 0.6481 3478.55
2 extreme -688.16 -498.93 0.2792 189.22
3 medium -538.79 -502.43 0.6191 36.37
4 poor -525.97 -501.53 0.2391 24.44
5 extreme -816.07 -518.91 0.0002676 297.16
6 medium -595.52 -500.09 0.1191 95.43
7 poor -1544.65 -577.02 0.2065 967.63
8 extreme -315733.42 -518.93 0.0002043 315214.49
9 medium -520.41 -499.65 0.1777 20.75
10 poor -581.16 -499.77 0.1872 81.39
11 extreme -1091.42 -518.90 0.0001381 572.51
12 medium -542.48 -499.12 0.3744 43.37
13 poor -584.10 -499.98 0.3308 84.12
14 extreme -14634.93 -503.12 0.3332 14131.81
15 medium -517.02 -498.83 0.19 18.19
16 poor -14150.89 -500.56 0.7324 13650.33
17 extreme -2465.75 -499.77 0.4987 1965.98
18 medium -591.33 -499.61 0.3979 91.72
19 poor -1011.82 -499.88 0.3093 511.94
20 extreme -778.77 -500.63 0.2953 278.15
Summary by starting-point quality:
quality starts mean_initial_logLik mean_terminal_logLik best_terminal_logLik mean_improvement
medium 7 -562.93 -499.72 -498.28 63.21
poor 7 -3196.88 -511.40 -499.77 2685.48
extreme 7 -48029.79 -508.46 -498.93 47521.33
Best MIF start: 0, panel logLik = -498.28 (MCSE 0.6324)
Starts within the stated MCSE-based tolerance of the best: 6 of 21
Best MIF minus R starting-vector logLik: +2.40 (combined MCSE 0.6857)
Best MIF minus stored reference logLik: +0.50 (combined MCSE 0.7106)All three start qualities use identical particles, iterations, perturbations, cooling, and final likelihood evaluation. Only the initial parameter displacement changes. The table reports the initial likelihood, terminal likelihood, Monte Carlo standard error, and improvement for every start, followed by a group summary. The best finite terminal parameter vector is used below. Two details of the call are useful:
- Cooling is configured on
RWSigma. Calling.geometric_cooling(0.7)makes 0.7 the multiplicative reduction in perturbation standard deviation accumulated over fifty iterations. block=Falseselects standard PIF. The MPIF implementation withblock=Trueis demonstrated in Diagnostic 2 for the mixed SRJF model containing unit-specific adult mortality.
The unit-specific evidence diagnostic below illustrates a mixed parameterisation in which a subset of parameters varies by unit. The all-shared model just estimated serves as the parsimonious baseline against which the unit-specific variant is compared via AIC.
Diagnostic 1: Parameter Scaling Verification
Particle filters are sensitive to parameter scaling, especially when parameters span many orders of magnitude. To confirm that the constructed object is correctly scaled, we run a quick particle filter at the development run level and inspect the per-unit log-likelihoods.
Code
ll_arr = srjf_start_ll_arr # (theta_idx, unit, replicate)
ll_per_unit_all, ll_per_unit_se_all, panel_ll_all, panel_ll_se_all = (
summarize_panel_loglik(ll_arr)
)
ll_per_unit = ll_per_unit_all[0]
print(f"panel pfilter wall: {sf(TIMINGS['srjf_start_pfilter'])}s "
f"(J={RL_SRJF['J_eval']}, reps={RL_SRJF['pf_reps']})")
print("Per-unit log-likelihood (log-mean-exp over replicates):")
for u, v in zip(unit_names, ll_per_unit):
print(f" {u}: {ll(v):>10s}")
print(f"Total panel logLik: {ll(panel_ll_all[0])} "
f"(MCSE {sf(panel_ll_se_all[0])})")panel pfilter wall: 12.54s (J=1000, reps=20)
Per-unit log-likelihood (log-mean-exp over replicates):
A: -53.03
B: -54.30
C: -50.75
D: -48.74
E: -46.42
F: -48.14
G: -49.92
H: -46.26
I: -51.69
J: -51.41
Total panel logLik: -500.68 (MCSE 0.265)A correctly scaled panel produces stable per-unit log likelihoods on the order of \(-50\) for SRJF data and a substantial effective sample size (ESS) at every observation time. Deliberately mis-scaled parameter vectors with rates perturbed by several orders of magnitude yield degenerate filters whose ESS collapses toward unity. The mis-scaled counter-example below quantifies this contrast.
If the pfilter output looks pathological (very negative per-unit log-likelihoods with high Monte Carlo variance), four checks usually localise the problem: (i) parameter units in the simulator (day\(^{-1}\), ind./L, etc.) are consistent with those of the data, (ii) dimensional and dimensionless parameters have biologically plausible magnitudes, (iii) srjf_par_trans correctly applies log transforms to all positive-valued parameters, and (iv) srjf_rinit returns initial densities consistent with the experimental protocol. This advice is based on similar obstacles that we encountered and resolved while developing this data analysis. Mis-scaled starting points can prevent iterated filtering from locating the MLE even with a correct algorithm, because the optimisation can become trapped in regions where numerical instability precludes accurate likelihood evaluation.
The numerical pathologies of parameter mis-specification are most easily seen empirically. We compare particle-filter performance between the well-scaled panelfood (the initial-guess panel from the smoke check above) and a deliberately corrupted panelfood_wrong whose adult filtration rate \(f^n_S\) has been multiplied by \(100\). The comparison uses initial parameters rather than an MLE because the mis-scaling pathology appears regardless of starting point: a mis-specified filtration rate destabilises the filter even at otherwise-reasonable parameter values.
Code
# Build a panel with a deliberately mis-scaled f_Sn (100x too large).
wrong_theta = dict(srjf_shared_theta)
wrong_theta["f_Sn"] = wrong_theta["f_Sn"] * 100.0
wrong_shared_keys = [k for k in wrong_theta if k != "sigSn"]
wrong_shared_df = pd.DataFrame(
{"shared": [wrong_theta[k] for k in wrong_shared_keys]},
index=wrong_shared_keys,
)
wrong_unit_df = pd.DataFrame(
{u: [0.0] for u in unit_names},
index=["sigSn"],
)
panelfood_wrong = pp.PanelPomp(
pomp_dict=srjf_pomp_dict,
theta=pp.PanelParameters(
{"shared": wrong_shared_df, "unit_specific": wrong_unit_df}
),
)
# Replicated pfilter at the well-scaled initial guess (same `panelfood` the
# smoke check used above) and at the mis-scaled variant.
def _run_srjf_scaling():
panel_correct = pp.PanelPomp(
pomp_dict=srjf_pomp_dict,
theta=pp.PanelParameters(panelfood.theta),
)
panel_wrong = pp.PanelPomp(
pomp_dict=srjf_pomp_dict,
theta=pp.PanelParameters(panelfood_wrong.theta),
)
panel_correct.pfilter(
J=RL_SRJF["J_eval"], reps=RL_SRJF["pf_reps"],
key=jax.random.key(801)
)
correct = np.asarray(panel_correct.results_history[-1].logLiks.values)
panel_wrong.pfilter(
J=RL_SRJF["J_eval"], reps=RL_SRJF["pf_reps"],
key=jax.random.key(801)
)
wrong = np.asarray(panel_wrong.results_history[-1].logLiks.values)
return {"correct": correct, "wrong": wrong}
scaling_out = cached("srjf_scaling_diagnostic", _run_srjf_scaling)
ll_correct = scaling_out["correct"] # (1, U, reps)
ll_wrong = scaling_out["wrong"]
unit_ll_correct_all, unit_se_correct_all, panel_ll_correct_all, panel_se_correct_all = (
summarize_panel_loglik(ll_correct)
)
unit_ll_wrong_all, unit_se_wrong_all, panel_ll_wrong_all, panel_se_wrong_all = (
summarize_panel_loglik(ll_wrong)
)
unit_ll_correct = unit_ll_correct_all[0]
unit_se_correct = unit_se_correct_all[0]
panel_ll_correct = float(panel_ll_correct_all[0])
panel_se_correct = float(panel_se_correct_all[0])
unit_ll_wrong = unit_ll_wrong_all[0]
unit_se_wrong = unit_se_wrong_all[0]
panel_ll_wrong = float(panel_ll_wrong_all[0])
panel_se_wrong = float(panel_se_wrong_all[0])
print("Well-scaled initial parameters:")
print(f" Panel log-likelihood: {ll(panel_ll_correct)} "
f"(MCSE: {sf(panel_se_correct)})")
print(f" Unit-level mean range: "
f"[{ll(np.nanmin(unit_ll_correct))}, "
f"{ll(np.nanmax(unit_ll_correct))}]")
print()
print("Mis-scaled parameters (f_Sn x 100):")
print(f" Panel log-likelihood: {ll(panel_ll_wrong)} "
f"(MCSE: {sf(panel_se_wrong)})")
print(f" Unit-level mean range: "
f"[{ll(np.nanmin(unit_ll_wrong))}, "
f"{ll(np.nanmax(unit_ll_wrong))}]")
print()
print(f"Log-likelihood difference: "
f"{ll(panel_ll_correct - panel_ll_wrong)}")Well-scaled initial parameters:
Panel log-likelihood: -499.97 (MCSE: 0.4259)
Unit-level mean range: [-54.33, -46.23]
Mis-scaled parameters (f_Sn x 100):
Panel log-likelihood: -1417.59 (MCSE: 287.1)
Unit-level mean range: [-266.13, -99.17]
Log-likelihood difference: 917.61The well-scaled initial parameters yield a finite panel log-likelihood with small Monte Carlo standard error, indicating stable filter performance. Mis-scaling generates severely degraded likelihood values: an \(f^n_S\) that is too large drives the algal resource \(F_u\) rapidly toward zero, starving the Daphnia population and producing simulated trajectories far below the observed adult densities. The mismatch surfaces as much more negative per-unit log-likelihoods and as inflated Monte Carlo variance, because the particle weights become dominated by a small fraction of trajectories.
Figure 2 plots the unit-level contributions side by side on a single shared vertical axis, so the two cases are compared on one scale rather than each within its own.
Code
# Connected markers on one shared vertical axis. The earlier version drew
# bars measured from each panel's own mean on independent axes, which made
# two very different quantities look comparable.
fig, axes = plt.subplots(1, 2, figsize=(10, 4), sharey=True)
unit_means_correct = unit_ll_correct
unit_means_wrong = unit_ll_wrong
positions = np.arange(1, len(unit_names) + 1)
finite_correct = np.isfinite(unit_means_correct)
baseline_correct = (
float(np.mean(unit_means_correct[finite_correct]))
if np.any(finite_correct)
else np.nan
)
finite_wrong = np.isfinite(unit_means_wrong)
baseline_wrong = (
np.mean(unit_means_wrong[finite_wrong])
if np.any(finite_wrong)
else -150.0
)
# A -Inf unit has no position on a log-likelihood axis. It is drawn just
# below the lowest finite unit and labelled, so it is visible without
# setting the axis limits from a value that is not on the scale. The offset
# is measured against the range of the shared axis, not against the
# mis-scaled panel alone: on one axis a fraction of that panel's own spread
# can be too small to separate the marker from the line.
wrong_plot = unit_means_wrong.copy()
if np.any(~finite_wrong) and np.any(finite_wrong):
shared_finite = np.concatenate([
unit_means_correct[np.isfinite(unit_means_correct)],
unit_means_wrong[finite_wrong],
])
shared_span = float(np.ptp(shared_finite)) if shared_finite.size else 0.0
wrong_plot[~finite_wrong] = (
float(np.min(unit_means_wrong[finite_wrong]))
- max(1.0, 0.05 * shared_span)
)
elif np.any(~finite_wrong):
wrong_plot[~finite_wrong] = -150.0
for ax, values, baseline, colour, title in (
(axes[0], unit_means_correct, baseline_correct, PALETTE["adult"],
"Well-scaled initial parameters"),
(axes[1], wrong_plot, baseline_wrong, PALETTE["warn"],
"Mis-scaled parameters ($f^n_S \\times 100$)"),
):
ax.plot(positions, values, marker="o", linestyle="-", linewidth=1.2,
markersize=5, color=colour)
ax.axhline(baseline, color=PALETTE["fit"], linestyle="--", linewidth=1.5,
label="mean over finite units")
ax.set_xticks(positions)
ax.set_xticklabels(unit_names)
ax.set_xlabel("Unit", fontsize=FONTS["axis_label"])
ax.set_title(title, fontsize=FONTS["panel_title"])
ax.tick_params(axis="both", labelsize=FONTS["tick"])
# One legend, on the panel whose data sits away from it. The dashed line means
# the same thing in both, and a second copy would collide with the -Inf labels.
axes[0].legend(fontsize=FONTS["legend"], frameon=False, loc="lower right")
axes[0].set_ylabel("Unit log-likelihood", fontsize=FONTS["axis_label"])
for index in np.flatnonzero(~finite_wrong):
axes[1].annotate(
"-Inf", (positions[index], wrong_plot[index]),
xytext=(0, 6), textcoords="offset points", ha="center",
fontsize=FONTS["annotation"], color=PALETTE["warn"],
)
fig.tight_layout()
Under correctly scaled parameters, unit-level log likelihoods exhibit moderate variation reflecting differences among observed trajectories. On the shared axis that variation is barely visible, which is the point: it is small compared with the cost of mis-scaling. Under mis-scaling every unit is far worse, by about \(140\) log units, and the degradation is uneven. The measurement model applies its fixed \(-150\) penalty once per violating observation. A unit whose trajectory violates a state bound at many times is therefore pushed lower than one with fewer violations. The spread across the mis-scaled panel counts constraint violations and is not a measure of how well the units are otherwise described. If every observation in every unit violated a bound, the panel log likelihood would reach \(-150 \times U \times N = -15{,}000\).
The well-scaled panel filters cleanly and mis-scaling is detectable rather than silent, so we proceed to maximum-likelihood estimation. The next subsection sets the algorithmic budget and perturbation scheme that the MIF call inherits.
Diagnostic 2: Evidence for Unit-Specific Parameterization
To assess whether unit-specific parameterisation of adult mortality \(\theta^n_S\) is statistically justified, we follow a three-step workflow: (i) examine the unit-level log-likelihood contributions under the all-shared MLE in Figure 3, (ii) re-estimate the model with \(\theta^n_{S,u}\) allowed to vary by unit, and (iii) compare the two via AIC. Substantial heterogeneity in unit-level log-likelihoods under the all-shared model is a necessary but not sufficient signal of true between-unit heterogeneity. Only the AIC comparison can adjudicate whether the apparent gap is real or absorbed by Monte Carlo noise.
Code
# Replicated pfilter at the all-shared MLE (best replicate from MIF above).
def _run_srjf_mle_pfilter():
panel = pp.PanelPomp(
pomp_dict=srjf_pomp_dict,
theta=srjf_mle_panel_theta,
)
panel.pfilter(J=J_pf_eval, reps=Npf_reps, key=jax.random.key(701))
return np.asarray(panel.results_history[-1].logLiks.values)[0]
unit_ll_at_mle = cached("srjf_mle_pfilter", _run_srjf_mle_pfilter) # (U, reps)
unit_ll_means = logmeanexp(unit_ll_at_mle, axis=1, ignore_nan=True)
unit_ll_ses = logmeanexp_se(unit_ll_at_mle, axis=1, ignore_nan=True)
unit_ll_baseline = float(np.nanmean(unit_ll_means))
fig, ax = plt.subplots(figsize=(8, 4))
ax.bar(
np.arange(len(unit_names)),
unit_ll_means - unit_ll_baseline, # height = deviation from the mean
bottom=unit_ll_baseline, # base = the panel mean
yerr=2 * unit_ll_ses,
capsize=3,
color=PALETTE["adult"],
)
ax.axhline(unit_ll_baseline, linestyle="--", color=PALETTE["fit"],
linewidth=1.5, label="panel mean")
ax.set_xticks(np.arange(len(unit_names)))
ax.set_xticklabels(unit_names)
ax.set_xlabel("Unit", fontsize=FONTS["axis_label"])
ax.set_ylabel("Unit log-likelihood", fontsize=FONTS["axis_label"])
ax.legend(loc="lower right", fontsize=FONTS["legend"], frameon=False)
ax.set_title("Unit-Level Likelihood Contributions (all-shared MLE)",
fontsize=FONTS["panel_title"])
ax.tick_params(axis="both", labelsize=FONTS["tick"])
fig.tight_layout()
These unit contributions are descriptive: different observed trajectories can have different likelihood magnitudes even when a shared parameterization is correct. The formal check below instead fits a nested alternative in which the adult mortality rate \(\theta^n_{S,u}\) differs across units and verifies that the optimization at least recovers the embedded null before calculating AIC.
Code
# Unit-specific variant: theta_Sn now lives in the unit_specific block with
# its initial MLE value replicated across all 10 units. Every other parameter
# (including the fixed sigSn) stays in the shared block; sigSn's rw_sd is
# zero so it never moves. p counts only *estimated* parameters: six shared
# (rn, f_Sn, theta_Jn, sigJn, sigF, k_Sn) plus ten unit-specific theta_Sn.
spec_keys = ["theta_Sn"]
shared_keys_uspec = [k for k in srjf_shared_theta if k not in spec_keys]
estimated_shared_uspec = [k for k in shared_keys_uspec if k != "sigSn"]
shared_uspec_df = pd.DataFrame(
{"shared": [srjf_mle_theta[k] for k in shared_keys_uspec]},
index=shared_keys_uspec,
)
unit_uspec_df = pd.DataFrame(
{u: [srjf_mle_theta[k] for k in spec_keys] for u in unit_names},
index=spec_keys,
)
panel_theta_uspec, _ = make_quality_panel_starts(
shared_uspec_df,
unit_uspec_df,
seed=300,
fixed=("sigSn",),
qualities=srjf_start_qualities,
)
# rw_sd: same magnitudes; theta_Sn now perturbs per-unit because it is
# unit-specific in panel_theta_uspec, not because we changed RWSigma.
srjf_rw_sd_uspec = pp.RWSigma(
sigmas={
"rn": dent_rw_sd,
"f_Sn": dent_rw_sd,
"theta_Sn": dent_rw_sd,
"theta_Jn": dent_rw_sd,
"sigSn": 0.0,
"sigJn": dent_rw_sd,
"sigF": dent_rw_sd,
"k_Sn": dent_rw_sd,
},
init_names=[],
).geometric_cooling(0.7)
def _evaluate_srjf_embedded_null():
"""Evaluate the all-shared MLE represented inside the specific model."""
panel = pp.PanelPomp(
pomp_dict=srjf_pomp_dict,
theta=panel_theta_uspec.subset(0),
)
panel.pfilter(
J=J_pf_eval, reps=Npf_reps, key=jax.random.key(202),
)
return np.asarray(panel.results_history[-1].logLiks.values)
embedded_ll_arr = cached(
"srjf_uspec_embedded_null", _evaluate_srjf_embedded_null,
)
_, _, embedded_panel_ll, embedded_panel_se = summarize_panel_loglik(
embedded_ll_arr
)
embedded_ll = float(embedded_panel_ll[0])
embedded_se = float(embedded_panel_se[0])
embedded_combined_mcse = float(np.sqrt(best_ll_se**2 + embedded_se**2))
embedded_difference = embedded_ll - best_ll
def _run_srjf_uspec():
panel = pp.PanelPomp(
pomp_dict=srjf_pomp_dict,
theta=pp.PanelParameters(panel_theta_uspec),
)
panel.mif(
J=J_mif, M=M_mif, rw_sd=srjf_rw_sd_uspec, block=True,
key=jax.random.key(303),
)
panel.pfilter(J=J_pf_eval, reps=Npf_reps, key=jax.random.key(404))
return {
"ll": np.asarray(panel.results_history[-1].logLiks.values),
"final_theta": panel.theta.params(as_list=True),
"traces": panel.traces(),
}
uspec_out = cached("srjf_uspec_mif", _run_srjf_uspec)
uspec_ll = uspec_out["ll"] # (Nstarts, U, reps)
uspec_unit_ll, uspec_unit_se, uspec_ll_per_start, uspec_se_per_start = (
summarize_panel_loglik(uspec_ll)
)
uspec_best = _finite_argmax(uspec_ll_per_start)
ll_uspec_best = float(uspec_ll_per_start[uspec_best])
uspec_se = float(uspec_se_per_start[uspec_best])
uspec_best_theta = uspec_out["final_theta"][uspec_best]
theta_Sn_by_unit = (
uspec_best_theta["unit_specific"].loc["theta_Sn", unit_names].astype(float)
)
theta_Sn_cv = float(theta_Sn_by_unit.std(ddof=1) / theta_Sn_by_unit.mean())
# AIC bookkeeping. p counts only *estimated* parameters (sigSn is held at 0).
p_shared = 7
p_specific = len(estimated_shared_uspec) + len(spec_keys) * len(unit_names)
combined_mcse = float(np.sqrt(best_ll_se**2 + uspec_se**2))
aic_shared = 2 * p_shared - 2 * best_ll
aic_specific = 2 * p_specific - 2 * ll_uspec_best
delta_aic = aic_specific - aic_shared
aic_table = pd.DataFrame({
"model": ["all-shared", "unit-specific theta_Sn"],
"p": [p_shared, p_specific],
"logLik": [best_ll, ll_uspec_best],
"AIC": [aic_shared, aic_specific],
})
print(f"unit-specific MIF + pfilter wall: {sf(TIMINGS['srjf_uspec_mif'])}s")
print(table_text(aic_table))
print(
f"Embedded-null representation: {ll(embedded_ll)} "
f"(MCSE {sf(embedded_se)}); difference from all-shared fit: "
f"{ll(embedded_difference, signed=True)} "
f"(combined MCSE {sf(embedded_combined_mcse)})"
)
print(f"Delta AIC (unit-specific - all-shared): {ll(delta_aic, signed=True)}")
print(f"Likelihood improvement: {ll(ll_uspec_best - best_ll, signed=True)}")
print(f"Combined likelihood MCSE: {sf(combined_mcse)}")
print(f"theta_Sn unit-specific coefficient of variation: {sf(theta_Sn_cv)}")unit-specific MIF + pfilter wall: 570.4s
model p logLik AIC
all-shared 7 -498.28 1010.56
unit-specific theta_Sn 16 -490.97 1013.94
Embedded-null representation: -546.69 (MCSE 0.3856); difference from all-shared fit: -48.41 (combined MCSE 0.7407)
Delta AIC (unit-specific - all-shared): +3.37
Likelihood improvement: +7.31
Combined likelihood MCSE: 0.6379
theta_Sn unit-specific coefficient of variation: 0.6256Replacing one shared parameter with \(U\) unit-specific values adds \(U-1\) free parameters. Thus, in this 10-unit panel, making \(\theta^n_S\) unit-specific adds \(10-1=9\) parameters. A negative \(\Delta\text{AIC}\) of substantial magnitude indicates that the likelihood improvement outweighs this penalty. A positive or near-zero \(\Delta\text{AIC}\) favours the all-shared model on grounds of parsimony. The embedded-null difference and combined Monte Carlo standard error are printed beside the AIC table because disagreement between equivalent parameterisations or a lower likelihood under the larger model indicates incomplete numerical optimisation.
The shared-versus-unit-specific decision must integrate multiple lines of evidence beyond AIC alone. The coefficient of variation of the estimated \(\hat\theta^n_{S,u}\) across units is a complementary diagnostic. Estimates that cluster tightly around their mean indicate that the unit-specific model is recovering essentially identical values for every unit (supporting the all-shared specification). Estimates that span many orders of magnitude without interpretable pattern indicate that the model is fitting idiosyncratic noise (also supporting the parsimonious choice). Genuine parameter heterogeneity manifests as moderate, biologically plausible between-unit variation that aligns with documented experimental conditions.
The displayed table always reports the two likelihood estimates, parameter counts, and corresponding AIC values. These are numerical outputs of the current run rather than an automated decision about model adequacy. Production SRJF results are reported in Section S9 and Table S12 of the Supplement.
Diagnostic 3: MIF searches visualization
Each MIF run produces a trace of log-likelihood values across iterations. Figure 4 compares the three starting-point qualities under the same algorithmic settings. Individual chains are noisy and overplot each other, so the figure shows one smoothed median curve per group. Its vertical axis is logarithmic in the magnitude of the log likelihood, under the increasing map \(y \mapsto -\log_{10}(-y)\); a better likelihood is therefore still higher on the axis, but a group sitting tens of thousands of log units below no longer compresses the converging curves onto the top border, which is what a linear axis and an interquartile band around such a group both did. The accompanying table quantifies convergence speed from the same group medians. Every group is timed against one shared target, the best final median any group reached less 10 log units, so the column answers the same question for each row: how many iterations until this group came within 10 log units of the best outcome? A group already inside the target at its first evaluated iteration reports that iteration, and a group that never arrives reports NaN. Timing each group against its own gain instead would not compare groups: a group starting near the maximum has almost no gain to make, and a group with a large gain reaches most of it early while still finishing far below the others.
Code
# Restrict to the shared-parameter rows; the unit-specific rows hold sigSn,
# which is fixed at 0 by rw_sd and therefore uninformative. Group medians
# summarise the chains without overplotting their raw traces.
srjf_quality_curves = mif_trace_curves(mif_traces, srjf_start_qualities)
# Groups that were started, and the subset that produced any finite
# likelihood. A group with none cannot be drawn, but it must still be named:
# vanishing from the figure is exactly how a failed start would be hidden.
srjf_quality_groups = [q for q in QUALITY_ORDER if q in set(srjf_start_qualities)]
srjf_groups_present = [
q for q in srjf_quality_groups
if q in set(srjf_quality_curves["group"])
]
srjf_groups_absent = [q for q in srjf_quality_groups
if q not in srjf_groups_present]
fig = plot_quality_curves(
srjf_quality_curves, srjf_groups_present,
"SRJF MIF convergence",
)
Code
srjf_speed_table, srjf_speed_target = convergence_speed_table(
srjf_quality_curves, srjf_quality_groups
)
print(
"Convergence speed by starting-point quality. Target is the best final "
f"median less {sf(SPEED_TOLERANCE)} log units = {ll(srjf_speed_target)}:"
)
print(table_text(srjf_speed_table))
print(describe_arrivals(srjf_speed_table))
if srjf_groups_absent:
print(
"No finite MIF trace at any iteration for: "
+ ", ".join(srjf_groups_absent)
+ ". These starts are shown in the table above as NaN and are absent "
"from the figure because they have no plottable value."
)Convergence speed by starting-point quality. Target is the best final median less 10 log units = -510.38:
group first_evaluated_logLik final_median_logLik median_gain iterations_to_target
medium -521.02 -500.38 20.63 2
poor -783.17 -501.23 281.93 5
extreme -692.38 -506.47 185.92 71
Reached the target, soonest first: medium at iteration 2, poor at iteration 5, extreme at iteration 71. Every group reached the target, so this budget recovered even the most displaced starts.The initial and terminal likelihood table quantifies improvement, while the curves show when that improvement occurs. A poor start that rises toward the terminal region demonstrates recovery from a deliberately displaced parameter vector. The extreme group is displaced by a log-scale standard deviation of 3.0 and so begins far below the others: twelve such starts, evaluated before any MIF, had a median panel log likelihood near \(-10{,}400\) against \(-568\) for medium and \(-1{,}865\) for poor, and the worst of them was below \(-300{,}000\). The fixed \(-150\) per violating observation would floor this ten-unit panel at \(-15{,}000\), so most of those values are not the constraint penalty at all: they are ordinary negative-binomial densities evaluated where the model badly misfits the counts, on a surface that is steep rather than flat. A start can still land on the penalty, and one that does has no gradient for MIF to follow; that is one of the ways an extreme start fails, not the only one. Whether 320 iterations at these perturbation sizes are enough to climb back from it is the question the group is there to answer, and a group that does not is the intended demonstration of the method’s limit rather than a defect in the run. Whether the same budget that moves the poor starts also recovers the extreme ones is settled by the table above, not here; the sentence printed beneath it states which groups came within 10 log units of the best final median and which did not. Persistent separation of this kind is what indicates that additional starts, more iterations, or a broader search are needed; it is not evidence that the algorithm is wrong. Production SRJF results are reported in Section S9 of the Supplement.
Trace patterns also guide tuning. Erratic trajectories can indicate that \(J\) is too small relative to the state dimension and measurement model. Remedies include increasing \(J\), modestly enlarging rw_sd, or slowing the cooling schedule by raising \(a\) from 0.7 to 0.8. Sudden log-likelihood drops warrant checking parameter scales, state bounds, and the log transform in srjf_par_trans. Large differences between the best and median terminal values suggest that more independent searches are needed.
Diagnostic 4: Monte Carlo Adjusted Profile
Following Ionides et al. (2017) and Ning, Ionides, and Ritov (2021), the Monte Carlo Adjusted Profile (MCAP) method provides confidence intervals when the likelihood is evaluated and maximised by Monte Carlo algorithms. A smoothed estimate of the profile likelihood reduces Monte Carlo error, quantifies it, and adjusts the confidence-interval cutoff so that nominal coverage is maintained. MCAP is especially useful in high-dimensional problems such as PanelPOMP inference, where it is impractical to spend enough compute to make Monte Carlo error negligible. Section S1 of the Supplement summarises the algorithm and its theoretical guarantees.
The interpretation of an MCAP result depends on the observed profile. A sharply curved profile with the unrestricted MIF estimate near the smoothed maximum differs from a nearly flat profile or one whose maximum lies on a grid boundary. Following the R tutorial, we examine the adult mortality rate \(\theta^n_S\) and the juvenile mortality rate \(\theta^n_J\) to illustrate these contrasting identifiability patterns.
Pypomp’s pp.mcap(parameter, loglik, level=0.95, span=0.75) accepts a one-dimensional parameter grid and the corresponding profile log likelihoods. It returns the MCAP estimate, confidence interval, standard errors, and fitted profile curves.
For each focal parameter we take nprof values spaced logarithmically over the explored window, and at every focal value we start nprof independent MIF searches whose remaining free parameters are drawn log-uniformly across a box spanning a decade in each direction around the estimate. The default window for the focal parameter is that same decade box; \(\theta^n_J\) overrides it with the explicit range given below, because a decade proved too narrow to contain either endpoint of its interval. The focal parameter is held fixed by setting its random-walk standard deviation to zero, every search is evaluated by a multi-replicate pfilter, and the profile point is the best of the searches at that focal value. This dispersed multi-start design is what keeps the profile away from the fixed \(-150\) constraint penalty: a focal value whose starts are all clustered in one region can be trapped where the measurement model penalises every trajectory, whereas starts spread across the box will normally include a feasible one.
Only 95% confidence intervals are reported. The first example, \(\theta^n_S\), illustrates a well-identified parameter with an interior maximum and visible curvature. The second, \(\theta^n_J\), illustrates a weakly identified parameter whose interval may reach the evaluated range. The figures and concise summaries are shown without raw profile tables or smoothing-sensitivity tables.
Code
MCAP_LEVEL = 0.95
def _mcap_checked(parameter, loglik, span):
"""Run a 95% MCAP fit and return any numerical warnings."""
with warnings.catch_warnings(record=True) as caught:
warnings.simplefilter("always")
result = pp.mcap(
parameter=parameter,
loglik=loglik,
level=MCAP_LEVEL,
span=span,
)
warning_messages = sorted({str(item.message) for item in caught})
return result, warning_messages
def _mcap_ci_array(result):
"""Return MCAP interval endpoints as floats, using NaN for open ends."""
return np.asarray([
np.nan if value is None else float(value) for value in result.ci
])
def _emit_figure(fig, label, caption):
"""Embed a numbered Matplotlib figure directly in the rendered Markdown."""
buffer = io.BytesIO()
fig.savefig(buffer, format="png", dpi=150, bbox_inches="tight")
plt.close(fig)
encoded = base64.b64encode(buffer.getvalue()).decode("ascii")
print(f"{{#{label}}}")
def _theta_payload(theta_dict):
"""Construct one all-shared PanelParameters payload from a dictionary.
The unit-specific block is an empty DataFrame whose *columns* carry the
unit names. Pypomp reads a panel's unit names from those columns, and
PanelPomp refuses a theta whose unit names do not match its pomp_dict, so
`None` cannot be used here even though the model is all-shared. The two
representations must never be mixed within one list of starts.
"""
return {
"shared": pd.DataFrame(
{"shared": [theta_dict[k] for k in theta_dict]},
index=list(theta_dict),
),
"unit_specific": pd.DataFrame(index=[], columns=unit_names),
}
# ---- Profile machinery ---------------------------------------------------
# R draws every profile start log-uniformly across the full parameter box and
# keeps the best start at each focal value. That dispersed design, rather than
# any screening of the output, is what keeps the R profile clear of the fixed
# -150 constraint penalty.
PROFILE_BOX_FACTOR = 10.0 # box is [theta / 10, theta * 10]
PROFILE_RW_SD = 0.05 # random-walk sd for the profile searches
PROFILE_BATCH = 50 # starts per mif call; R used one mif2 per row
PROFILE_INFORMATIVE_DROP = 50.0
PROFILE_IMPLAUSIBLE_GAIN = 10.0
PROFILE_MIN_POINTS = 7
def _profile_implausible(ll_grid, unrestricted_loglik):
"""Points that beat the unrestricted maximum by an impossible margin.
A constrained maximum cannot legitimately exceed the unconstrained one, so
a point far above it is a numerical artefact rather than a better fit: a
widely dispersed start can reach a region where the particle filter
degenerates and reports an inflated likelihood. Because the profile takes
the best of many dispersed starts, it actively selects for that failure,
so the check must be applied before any maximum is taken.
"""
ll = np.asarray(ll_grid, dtype=float)
return np.isfinite(ll) & (ll > unrestricted_loglik + PROFILE_IMPLAUSIBLE_GAIN)
def _profile_informative(ll_grid, unrestricted_loglik):
"""Points usable for the MCAP fit.
Implausibly high points are removed first, and only then is the low-side
cutoff measured from the best surviving point. Anchoring on the raw
maximum instead would let a single artefact define the reference and
exclude the entire genuine profile.
"""
ll = np.asarray(ll_grid, dtype=float)
usable = np.isfinite(ll) & ~_profile_implausible(ll, unrestricted_loglik)
if not np.any(usable):
return usable
return usable & (ll > np.nanmax(ll[usable]) - PROFILE_INFORMATIVE_DROP)
def _profile_mask(ll_grid, unrestricted_loglik, label):
"""Return the MCAP fitting mask, reporting what it removes and why."""
ll = np.asarray(ll_grid, dtype=float)
implausible = _profile_implausible(ll, unrestricted_loglik)
keep = _profile_informative(ll, unrestricted_loglik)
n_high = int(np.sum(implausible))
if n_high:
print(f" {label}: {n_high} point(s) exceed the unrestricted "
f"log-likelihood by more than {PROFILE_IMPLAUSIBLE_GAIN:.0f} "
"units. These are treated as filter artefacts, not as a better "
"fit, and are excluded from the smoothing.")
n_low = int(np.sum(~keep)) - n_high
if n_low > 0:
print(f" {label}: {n_low} further point(s) lie more than "
f"{PROFILE_INFORMATIVE_DROP:.0f} units below the best usable "
"point and are excluded from the smoothing.")
if int(np.sum(keep)) < PROFILE_MIN_POINTS:
print(f" {label}: fewer than {PROFILE_MIN_POINTS} usable points "
"remain. Fitting on all finite points so that the render "
"completes; inspect the resulting fit with particular care.")
keep = np.isfinite(ll)
return keep
def generate_parameter_profile(prof_name, nprof, seed, focal_log_range=None):
"""Build a randomised profile design.
`nprof` focal values spanning [theta/10, theta*10] on the log scale, or
`focal_log_range` when the focal parameter needs a wider window than a
decade either side of the estimate; at each, `nprof` starts with every
other free parameter drawn log-uniformly across the decade box. sigSn is
held at zero. The free-parameter box is unaffected by `focal_log_range`:
widening the focal grid asks where the profile falls away, not where the
nuisance parameters live.
"""
rng = np.random.default_rng(seed)
lb = {k: v / PROFILE_BOX_FACTOR for k, v in srjf_mle_theta.items()}
ub = {k: v * PROFILE_BOX_FACTOR for k, v in srjf_mle_theta.items()}
free = [k for k in srjf_mle_theta if k not in ("sigSn", prof_name)]
focal_lo, focal_hi = (
(np.log(lb[prof_name]), np.log(ub[prof_name]))
if focal_log_range is None else focal_log_range
)
focal_values = np.exp(np.linspace(focal_lo, focal_hi, nprof))
rows = []
for value in focal_values:
for _ in range(nprof):
row = dict(srjf_mle_theta)
for k in free:
row[k] = float(np.exp(
rng.uniform(np.log(lb[k]), np.log(ub[k]))
))
row[prof_name] = float(value)
row["sigSn"] = 0.0
rows.append(row)
return pd.DataFrame(rows)
def generate_sd(x, profile_name):
"""R: generate_sd(). Uniform perturbation x; sigSn and the focal frozen."""
sd_list = {k: x for k in srjf_mle_theta}
sd_list["sigSn"] = 0.0
sd_list[profile_name] = 0.0
return sd_list
def run_profile(prof_name, nprof, seed, focal_log_range=None):
"""Fit every row of the profile design, then keep the best log-likelihood
at each focal value.
Pypomp vectorises MIF over the start dimension, so rows are processed in
batches of PROFILE_BATCH rather than one at a time. The batch size trades
GPU memory against wall time and does not change the method.
"""
parameter_shared = generate_parameter_profile(
prof_name, nprof, seed, focal_log_range=focal_log_range,
)
dent_rw_sd_first = generate_sd(PROFILE_RW_SD, prof_name)
rw_sd_prof = pp.RWSigma(
sigmas=dent_rw_sd_first, init_names=[]
).geometric_cooling(0.7)
loglik, mcse = [], []
for i in range(0, len(parameter_shared), PROFILE_BATCH):
rows = parameter_shared.iloc[i:i + PROFILE_BATCH]
panel = pp.PanelPomp(
pomp_dict=srjf_pomp_dict,
theta=pp.PanelParameters(
[_theta_payload(r) for r in rows.to_dict("records")]
),
)
panel.mif(
J=J_mif, M=M_mif, rw_sd=rw_sd_prof, block=False,
key=jax.random.key(seed + i),
)
panel.pfilter(
J=J_pf_eval, reps=Npf_reps, key=jax.random.key(seed + i + 1),
)
_, _, ll, ll_se = summarize_panel_loglik(
np.asarray(panel.results_history[-1].logLiks.values)
)
loglik.extend(np.asarray(ll, dtype=float))
mcse.extend(np.asarray(ll_se, dtype=float))
final_params = parameter_shared.copy()
final_params["loglik"] = loglik
final_params["mcse"] = mcse
subset = (
final_params
.loc[final_params.groupby(prof_name)["loglik"].idxmax()]
.sort_values(prof_name)
.reset_index(drop=True)
)
subset[f"log_{prof_name}"] = np.log(subset[prof_name])
return final_params, subset
# Profile: theta_Sn --------------------------------------------------------
name_str = "theta_Sn"
focal_Sn = name_str
nprof_Sn = 20 # focal values, and starts per focal value
final_params_Sn, subset_data_theta_Sn = cached(
"srjf_profile_theta_Sn",
lambda: run_profile(name_str, nprof_Sn, seed=1000),
)
log_grid_Sn = subset_data_theta_Sn["log_theta_Sn"].to_numpy()
ll_grid_Sn = subset_data_theta_Sn["loglik"].to_numpy()
se_grid_Sn = subset_data_theta_Sn["mcse"].to_numpy()
grid_Sn = np.exp(log_grid_Sn)
focal_mle_Sn = float(srjf_mle_theta[focal_Sn])
keep_Sn = _profile_mask(ll_grid_Sn, best_ll, "theta_Sn")
mcap_Sn, mcap_Sn_warnings = _mcap_checked(
parameter=log_grid_Sn[keep_Sn], loglik=ll_grid_Sn[keep_Sn],
span=0.6,
)
if mcap_Sn_warnings:
print(" Captured MCAP numerical diagnostics:")
for message in mcap_Sn_warnings:
print(f" - {message}")
print(f"theta_Sn profile wall: {sf(TIMINGS['srjf_profile_theta_Sn'])}s "
f"(nprof={nprof_Sn}, starts={nprof_Sn}, J={J_mif}, M={M_mif})")
mif_log_Sn = float(np.log(srjf_mle_theta[focal_Sn]))
# Summarise the profile that is actually fitted and displayed. Points already
# classified as filter artefacts are not used to define its maximum.
kept_idx_Sn = np.flatnonzero(keep_Sn)
best_kept_Sn = _finite_argmax(ll_grid_Sn[keep_Sn])
raw_max_idx_Sn = int(kept_idx_Sn[best_kept_Sn])
profile_max_gap_Sn = float(ll_grid_Sn[raw_max_idx_Sn] - best_ll)
profile_gap_se_Sn = float(np.sqrt(
se_grid_Sn[raw_max_idx_Sn] ** 2 + best_ll_se ** 2
))
mcap_Sn_ci = _mcap_ci_array(mcap_Sn)
print("95% MCAP summary for theta_Sn:")
print(f" MLE (log / natural): {sf(mcap_Sn.mle)} / {sf(np.exp(mcap_Sn.mle))}")
print(f" 95% CI (log): {interval(mcap_Sn_ci)}")
if np.all(np.isfinite(mcap_Sn_ci)):
print(f" 95% CI (natural): {interval(np.exp(mcap_Sn_ci))}")
print(
" se_stat / se_mc / se_total: "
f"{sf(mcap_Sn.se_stat)} / {sf(mcap_Sn.se_mc)} / {sf(mcap_Sn.se_total)}"
)
print(
" raw profile maximum minus unrestricted logLik: "
f"{ll(profile_max_gap_Sn, signed=True)} "
f"(combined MCSE {sf(profile_gap_se_Sn)})"
)
print(f" unrestricted MIF estimate (log): {sf(mif_log_Sn)}") theta_Sn: 3 further point(s) lie more than 50 units below the best usable point and are excluded from the smoothing.
theta_Sn profile wall: 1541s (nprof=20, starts=20, J=1000, M=250)
95% MCAP summary for theta_Sn:
MLE (log / natural): -0.2127 / 0.8084
95% CI (log): (-0.764, 0.4239)
95% CI (natural): (0.4658, 1.528)
se_stat / se_mc / se_total: 0.3093 / 0.03844 / 0.3117
raw profile maximum minus unrestricted logLik: -0.46 (combined MCSE 0.7687)
unrestricted MIF estimate (log): -0.3053Code
fig, ax = plt.subplots(figsize=(7, 4.5))
# Raw profile points.
ax.errorbar(
log_grid_Sn[keep_Sn], ll_grid_Sn[keep_Sn], yerr=2 * se_grid_Sn[keep_Sn],
fmt="o", color=PALETTE["data"], ecolor=PALETTE["quad"],
capsize=2, zorder=3, label="profile logLik ± 2 MCSE",
)
# Smoothed and quadratic curves on the dense grid returned by mcap().
fit_x = mcap_Sn.fit["parameter"]
ax.plot(fit_x, mcap_Sn.fit["smoothed"], color=PALETTE["fit"], linewidth=1.6,
label="LOWESS smoothed")
ax.plot(fit_x, mcap_Sn.fit["quadratic"], color=PALETTE["quad"], linestyle="--",
linewidth=1.2, label="local quadratic")
# CI and MLE markers.
if mcap_Sn.ci[0] is not None:
ax.axvline(mcap_Sn.ci[0], color=PALETTE["ci"], linestyle="--",
linewidth=1.0)
if mcap_Sn.ci[1] is not None:
ax.axvline(mcap_Sn.ci[1], color=PALETTE["ci"], linestyle="--",
linewidth=1.0)
ax.axvline(mcap_Sn.mle, color=PALETTE["mle"], linewidth=1.4, label="MCAP MLE")
ax.axvline(np.log(srjf_mle_theta["theta_Sn"]), color=PALETTE["mif"],
linestyle=":", linewidth=1.4, label="MIF point estimate")
# Cutoff line: max(smoothed) - delta/2.
cutoff = float(np.nanmax(mcap_Sn.fit["smoothed"])) - 0.5 * mcap_Sn.delta
ax.axhline(cutoff, color=PALETTE["ci"], linestyle=":", linewidth=0.8,
label="max - $\\delta/2$")
# Keep the vertical range on the informative points, so that a single
# constraint-penalty point cannot flatten the profile into a line.
_ll_keep = ll_grid_Sn[keep_Sn]
if _ll_keep.size and np.any(np.isfinite(_ll_keep)):
_top = float(np.nanmax(_ll_keep))
# delta is NaN when the local quadratic has non-positive curvature, which
# is the normal outcome for a flat profile, so it cannot be trusted here.
_cands = [float(_top - np.nanmin(_ll_keep)), 1.0]
if np.isfinite(mcap_Sn.delta):
_cands.append(float(mcap_Sn.delta))
_span = max(_cands)
ax.set_ylim(_top - 1.5 * _span, _top + 0.5 * _span)
for _x in log_grid_Sn[~keep_Sn]:
ax.annotate("excluded", (_x, _top - 1.45 * _span), rotation=90,
ha="center", va="bottom", fontsize=FONTS["annotation"],
color=PALETTE["warn"])
ax.set_xlabel(r"$\log(\theta^n_S)$", fontsize=FONTS["axis_label"])
ax.set_ylabel("panel log-likelihood", fontsize=FONTS["axis_label"])
ax.set_title(r"MCAP profile for $\theta^n_S$", fontsize=FONTS["panel_title"])
ax.tick_params(axis="both", labelsize=FONTS["tick"])
ax.legend(
loc="upper center",
bbox_to_anchor=(0.5, -0.18),
fontsize=FONTS["legend"],
ncol=3,
frameon=False,
)
fig.tight_layout()
fig.subplots_adjust(bottom=0.30)
_emit_figure(
fig,
label="fig-srjf-mcap-theta-Sn",
caption=(
"MCAP profile log-likelihood for $\\theta^n_S$ on the log "
"scale. Points and error bars are independently evaluated profile "
"log-likelihoods with $\\pm 2$ MCSE. Curves and vertical markers show "
"the MCAP smoothing, interval, and unrestricted estimate. Grid points "
"that exceed the unrestricted log-likelihood, or that fall far below "
"the best usable point, are filter artefacts rather than likelihood "
"evaluations. They are excluded from the fit and marked at the axis."
),
)
The fitted profile has an interior maximum and clear curvature. This is the well-identified example. The 95% interval and its Monte Carlo and statistical standard errors provide the uncertainty summary.
MCAP does not always yield well-defined confidence intervals: the result depends on both information in the data and Monte Carlo optimisation. Juveniles are not directly observed in the fitted SRJF likelihood, so \(\theta^n_J\) is informed indirectly through adult recruitment and deserves particular scrutiny. A nearly flat profile, a maximum at a grid endpoint, an interval touching the evaluated range, non-positive fitted curvature, or a large difference from the unrestricted fit can all be read directly from the output below.
Code
# Profile: theta_Jn --------------------------------------------------------
name_str = "theta_Jn"
focal_Jn = name_str
nprof_Jn = 50 # focal values, and starts per focal value
# theta_Jn is explored over an explicit window rather than a decade either
# side of the estimate. In the previous run the grid ran from log -10.06 to
# -5.46 and the smoothed profile peaked at its lower edge, so neither endpoint
# of the interval could be located inside the evaluated range. The window
# below covers ten log units and reaches theta_Jn = 1, which is four orders of
# magnitude above the unrestricted estimate: wide enough for the profile to
# fall away on the right if it is going to.
THETA_JN_LOG_RANGE = (-10.0, 0.0)
final_params_Jn, subset_data_theta_Jn = cached(
"srjf_profile_theta_Jn",
lambda: run_profile(
name_str, nprof_Jn, seed=2000, focal_log_range=THETA_JN_LOG_RANGE,
),
)
log_grid_Jn = subset_data_theta_Jn["log_theta_Jn"].to_numpy()
ll_grid_Jn = subset_data_theta_Jn["loglik"].to_numpy()
se_grid_Jn = subset_data_theta_Jn["mcse"].to_numpy()
grid_Jn = np.exp(log_grid_Jn)
focal_mle_Jn = float(srjf_mle_theta[focal_Jn])
keep_Jn = _profile_mask(ll_grid_Jn, best_ll, "theta_Jn")
mcap_Jn, mcap_Jn_warnings = _mcap_checked(
parameter=log_grid_Jn[keep_Jn], loglik=ll_grid_Jn[keep_Jn],
span=0.95,
)
if mcap_Jn_warnings:
print(" Captured MCAP numerical diagnostics:")
for message in mcap_Jn_warnings:
print(f" - {message}")
print(f"theta_Jn profile wall: {sf(TIMINGS['srjf_profile_theta_Jn'])}s "
f"(nprof={nprof_Jn}, starts={nprof_Jn}, J={J_mif}, M={M_mif})")
# As for theta_Sn: assess the fitted profile, not points already excluded.
kept_idx_Jn = np.flatnonzero(keep_Jn)
best_kept_Jn = _finite_argmax(ll_grid_Jn[keep_Jn])
raw_max_idx_Jn = int(kept_idx_Jn[best_kept_Jn])
mcap_Jn_combined_se = float(np.sqrt(
se_grid_Jn[raw_max_idx_Jn] ** 2 + best_ll_se ** 2
))
profile_max_gap_Jn = float(ll_grid_Jn[raw_max_idx_Jn] - best_ll)
mcap_Jn_ci = _mcap_ci_array(mcap_Jn)
print("95% MCAP summary for theta_Jn:")
print(f" MLE (log / natural): {sf(mcap_Jn.mle)} / {sf(np.exp(mcap_Jn.mle))}")
print(f" 95% CI (log): {interval(mcap_Jn_ci)}")
if np.all(np.isfinite(mcap_Jn_ci)):
print(f" 95% CI (natural): {interval(np.exp(mcap_Jn_ci))}")
print(
" se_stat / se_mc / se_total: "
f"{sf(mcap_Jn.se_stat)} / {sf(mcap_Jn.se_mc)} / {sf(mcap_Jn.se_total)}"
)
print(
" raw profile maximum minus unrestricted logLik: "
f"{ll(profile_max_gap_Jn, signed=True)} "
f"(combined MCSE {sf(mcap_Jn_combined_se)})"
)
print(f" unrestricted MIF estimate (log): {sf(np.log(focal_mle_Jn))}")
print(f" explored range (log): {interval(THETA_JN_LOG_RANGE)}")theta_Jn profile wall: 7506s (nprof=50, starts=50, J=1000, M=250)
95% MCAP summary for theta_Jn:
MLE (log / natural): -7.017 / 0.0008965
95% CI (log): (-10, -3.453)
95% CI (natural): (4.54e-05, 0.03164)
se_stat / se_mc / se_total: 2.011 / 0.1596 / 2.018
raw profile maximum minus unrestricted logLik: +0.27 (combined MCSE 0.7638)
unrestricted MIF estimate (log): -5.735
explored range (log): (-10, 0)Code
fig, ax = plt.subplots(figsize=(7, 4.5))
ax.errorbar(
log_grid_Jn[keep_Jn], ll_grid_Jn[keep_Jn], yerr=2 * se_grid_Jn[keep_Jn],
fmt="o", color=PALETTE["data"], ecolor=PALETTE["quad"],
capsize=2, zorder=3, label="profile logLik ± 2 MCSE",
)
fit_x_J = mcap_Jn.fit["parameter"]
ax.plot(fit_x_J, mcap_Jn.fit["smoothed"], color=PALETTE["fit"], linewidth=1.6,
label="LOWESS smoothed")
ax.plot(fit_x_J, mcap_Jn.fit["quadratic"], color=PALETTE["quad"],
linestyle="--", linewidth=1.2, label="local quadratic")
if mcap_Jn.ci[0] is not None:
ax.axvline(mcap_Jn.ci[0], color=PALETTE["ci"], linestyle="--",
linewidth=1.0)
if mcap_Jn.ci[1] is not None:
ax.axvline(mcap_Jn.ci[1], color=PALETTE["ci"], linestyle="--",
linewidth=1.0)
ax.axvline(mcap_Jn.mle, color=PALETTE["mle"], linewidth=1.4, label="MCAP MLE")
ax.axvline(np.log(srjf_mle_theta["theta_Jn"]), color=PALETTE["mif"],
linestyle=":", linewidth=1.4, label="MIF point estimate")
cutoff_J = float(np.nanmax(mcap_Jn.fit["smoothed"])) - 0.5 * mcap_Jn.delta
ax.axhline(cutoff_J, color=PALETTE["ci"], linestyle=":", linewidth=0.8,
label="max - $\\delta/2$")
_ll_keep_J = ll_grid_Jn[keep_Jn]
if _ll_keep_J.size and np.any(np.isfinite(_ll_keep_J)):
_top_J = float(np.nanmax(_ll_keep_J))
# delta is NaN for a flat profile; see the theta_Sn figure above.
_cands_J = [float(_top_J - np.nanmin(_ll_keep_J)), 1.0]
if np.isfinite(mcap_Jn.delta):
_cands_J.append(float(mcap_Jn.delta))
_span_J = max(_cands_J)
ax.set_ylim(_top_J - 1.5 * _span_J, _top_J + 0.5 * _span_J)
for _x in log_grid_Jn[~keep_Jn]:
ax.annotate("excluded", (_x, _top_J - 1.45 * _span_J), rotation=90,
ha="center", va="bottom", fontsize=FONTS["annotation"],
color=PALETTE["warn"])
ax.set_xlabel(r"$\log(\theta^n_J)$", fontsize=FONTS["axis_label"])
ax.set_ylabel("panel log-likelihood", fontsize=FONTS["axis_label"])
ax.set_title(r"MCAP profile for $\theta^n_J$", fontsize=FONTS["panel_title"])
ax.tick_params(axis="both", labelsize=FONTS["tick"])
ax.legend(
loc="upper center",
bbox_to_anchor=(0.5, -0.18),
fontsize=FONTS["legend"],
ncol=3,
frameon=False,
)
fig.tight_layout()
fig.subplots_adjust(bottom=0.30)
_emit_figure(
fig,
label="fig-srjf-mcap-theta-Jn",
caption=(
"MCAP profile for $\\theta^n_J$ on the log scale. Points "
"and error bars are independently evaluated profile log-likelihoods "
"with $\\pm 2$ MCSE. Curves and vertical markers show the MCAP "
"smoothing, interval, and unrestricted estimate."
),
)
This is the weakly identified example. If the 95% interval reaches the explored range, it should not be interpreted as a resolved finite endpoint. This behavior is scientifically plausible because juveniles are not directly observed in the fitted SRJF likelihood. The grid here spans ten log units, from \(\log\theta^n_J = -10\) to \(0\), rather than the decade either side of the estimate used for \(\theta^n_S\). The narrower grid used previously ran from \(-10.06\) to \(-5.46\) and its smoothed profile peaked at the lower edge, which is a boundary effect and not an estimate; widening the window is what allows a boundary effect and a genuinely flat profile to be told apart. An interval that still reaches an endpoint of this wider range is evidence for the flat profile rather than for the grid.
Section 1 wall times at run_level = 3 (s)
srjf_start_pfilter 12.54
srjf_mif_search 179.6
srjf_scaling_diagnostic 17.61
srjf_mle_pfilter 9.248
srjf_uspec_embedded_null 7.34
srjf_uspec_mif 570.4
srjf_profile_theta_Sn 1541
srjf_profile_theta_Jn 7506Section 2: SIRJPF2 Model
Experimental Design and Data Structure
The SIRJPF2 model extends the SRJF baseline of Section 1 to the full ecological complexity of the Searle et al. (2016) mesocosm experiment: native D. dentifera and invasive D. lumholtzi compete for a shared algal food resource (Ankistrodesmus falcatus) while both are exposed to A. monospora. The treatment consists of \(U = 8\) replicated mesocosms, each initialised with both host species in the presence of parasite inoculum, sampled at \(N = 10\) time points \(t_{u,n} = 5n + 2\) days for \(n \in 1{:}10\). Susceptible and infected adult densities are recorded for both species. Juvenile counts and parasite-spore / algal-resource densities remain latent.
As in the R tutorial, this section presents the all-shared SIRJPF2 model and its standard panel iterated-filtering analysis. The unit-specific MPIF example remains in the SRJF section above.
The PanelPOMP data structure comprises:
- Units: \(u \in 1{:}8\) independent mesocosms with two-species competition (replicate labels \(K, L, M, N, O, P, Q, S\)).
- Observation times: \(t_{u,n} = 5n + 2\) days for \(n \in 1{:}10\).
- Observed data: \(y^*_{u,n} = (N^n_{S,u,n}, N^n_{I,u,n}, N^l_{S,u,n}, N^l_{I,u,n})^\top\), the four-dimensional count vector of susceptible and infected adults for the native (\(n\)) and invasive (\(l\)) species.
- Latent process: \(X_u(t) = (S^n_u(t), I^n_u(t), J^n_u(t), S^l_u(t), I^l_u(t), J^l_u(t), F_u(t), P_u(t))^\top\) for \(t \in [t_{u,0}, t_{u,N}]\).
For each species \(k \in \{n, l\}\) we track susceptible adult density \(S^k_u(t)\), infected adult density \(I^k_u(t)\), and juvenile density \(J^k_u(t)\). These compartments are coupled through the shared parasite spore pool \(P_u(t)\) (in units of \(10^3\) spores/L) and the algal food resource \(F_u(t)\) (\(10^6\) cells per litre). Dead or removed individuals leave the dynamics into an absorbing \(R^k_u(t)\) compartment that is not modelled explicitly.
The acronym SIRJPF2 encodes this structure: Susceptible, Infected, and Removed disease compartments, Juveniles for age structure, Parasite spores for the environmental transmission stage, Food for bottom-up regulation, and the trailing 2 signals the two-species extension. We read sheet 3 ('both species combined') of the Excel file and slice iloc[90:170]. As documented in quality_reports/audits/DATA_SCHEMA.md, the column 'dent.ephip ' carries a trailing space that we strip after read. The calendar-day conversion (day - 1) * 5 + 7 matches Section 1.
Show data load and plotting code
# Load sheet 3 ("both species combined"). The 'dent.ephip ' column has a
# trailing space; strip after read. Column-name and slicing details are in
# quality_reports/audits/DATA_SCHEMA.md.
xls = pd.ExcelFile('../data/Mesocosmdata.xls')
sirjpf_raw = xls.parse('both species combined').iloc[90:170].copy()
sirjpf_raw.columns = sirjpf_raw.columns.str.strip()
sirjpf_raw['day'] = (sirjpf_raw['day'] - 1) * 5 + 7
sirjpf_data = (sirjpf_raw[['rep', 'day', 'dent.adult', 'dent.inf',
'lum.adult', 'lum.adult.inf']]
.sort_values(['rep', 'day'])
.reset_index(drop=True))
sirjpf_unit_names = ['K', 'L', 'M', 'N', 'O', 'P', 'Q', 'S']
sirjpf_trial_to_mesocosm = {t: f"Mesocosm {i+1}"
for i, t in enumerate(sirjpf_unit_names)}
sirjpf_data['Mesocosm'] = sirjpf_data['rep'].map(sirjpf_trial_to_mesocosm)
sirjpf_data['Mesocosm'] = pd.Categorical(
sirjpf_data['Mesocosm'],
categories=[f"Mesocosm {i+1}" for i in range(8)],
ordered=True,
)
# Four-row faceted plot. `sharey='row'` keeps each observation channel on a
# single comparable scale across the eight mesocosms, with y-axis tick labels
# only in the left column.
fig, axes = plt.subplots(
4, 8,
sharex=True,
sharey='row',
figsize=(10, 6),
gridspec_kw={'hspace': 0.18, 'wspace': 0.20},
)
obs_cols = ['dent.adult', 'dent.inf', 'lum.adult', 'lum.adult.inf']
row_titles = [r"Native $S^n$", r"Native $I^n$",
r"Invasive $S^l$", r"Invasive $I^l$"]
row_colors = [PALETTE["adult"], PALETTE["infected"],
PALETTE["lum_adult"], PALETTE["lum_inf"]]
row_styles = ['-', '-', '--', '--']
sqrt_scale = dict(value='function', functions=(np.sqrt, np.square))
for j, label in enumerate([f"Mesocosm {i+1}" for i in range(8)]):
sub = sirjpf_data[sirjpf_data['Mesocosm'] == label]
for i, col in enumerate(obs_cols):
ax = axes[i, j]
ax.plot(sub['day'], sub[col],
color=row_colors[i], linestyle=row_styles[i], linewidth=0.7)
ax.set_yscale(**sqrt_scale)
ax.yaxis.set_major_locator(plt.MaxNLocator(nbins=3))
ax.set_xlim(0, 52)
ax.set_xticks([0, 25, 50])
ax.tick_params(axis='both', labelsize=FONTS["tick"])
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
if i == 0:
ax.set_title(label.replace('Mesocosm ', 'Mesocosm-'),
fontsize=FONTS["panel_title"])
if j == 0:
ax.set_ylabel(row_titles[i] + '\n(ind./L)',
fontsize=FONTS["panel_title"])
if i == 3:
ax.set_xlabel('Day', fontsize=FONTS["panel_title"])
fig.align_ylabels(axes[:, 0])
fig.subplots_adjust(left=0.08, right=0.99, top=0.94, bottom=0.08)
Figure 7 displays the four observed density streams. Consistent with the competitive dominance reported by Searle et al. (2016) under shared parasitism, the native species typically attains higher peak densities than the invasive species, with infection prevalence varying considerably between mesocosms. These species-specific trajectories, coupled with cross-unit variability in parasite dynamics, motivate a mechanistic model that explicitly represents competitive and parasitic interactions.
Mechanistic Model
Model Specification
The SIRJPF2 latent process is an eight-dimensional system of stochastic differential equations describing two interacting host species, a shared parasite spore pool, and a shared algal resource. For each species \(k \in \{n, l\}\) the state vector contains susceptible adults \(S^k_u\), infected adults \(I^k_u\), and juveniles \(J^k_u\). These populations interact through \(P_u\) and \(F_u\), which are common to both species within each unit:
\[\begin{aligned} dS^{k}_{u}(t) &= \lambda^{k}_{J}\, J^{k}_{u}(t)\,dt - \big\{\theta^{k}_{S} + p^{k} f^{k}_{S} P_{u}(t) + \delta\big\}\, S^{k}_{u}(t)\,dt + S^{k}_{u}(t)\, d\zeta^{k}_{S,u}(t),\\ dI^{k}_{u}(t) &= p^{k} f^{k}_{S}\, S^{k}_{u}(t)\, P_{u}(t)\,dt - \big\{\theta^{k}_{I} + \delta\big\}\, I^{k}_{u}(t)\,dt + I^{k}_{u}(t)\, d\zeta^{k}_{I,u}(t),\\ dJ^{k}_{u}(t) &= r^{k} f^{k}_{S}\, F_{u}(t)\, S^{k}_{u}(t)\,dt - \big\{\theta^{k}_{J} + \delta + \lambda^{k}_{J}\big\}\, J^{k}_{u}(t)\,dt + J^{k}_{u}(t)\, d\zeta^{k}_{J,u}(t),\\ dP_{u}(t) &= \sum_{k\in\{n,l\}}\!\!\Big(\beta^{k} \theta^{k}_{I}\, I^{k}_{u}(t) - f^{k}_{S}\,\big\{S^{k}_{u}(t) + \xi I^{k}_{u}(t)\big\}\, P_{u}(t)\Big)\,dt - \theta_{P} P_{u}(t)\,dt - \delta P_{u}(t)\,dt + P_{u}(t)\, d\zeta_{P,u}(t),\\ dF_{u}(t) &= -\sum_{k\in\{n,l\}}\!\!f^{k}_{S}\, F_{u}(t)\,\big\{S^{k}_{u}(t) + \xi I^{k}_{u}(t) + J^{k}_{u}(t)\big\}\, dt - \delta F_{u}(t)\,dt + \mu\, dt + F_{u}(t)\, d\zeta_{F,u}(t), \end{aligned}\]
with stochastic increments
\[\begin{aligned} d\zeta^{k}_{S,u}(t) &\sim \mathcal{N}\!\big(0,(\sigma^{k}_{S})^{2} dt\big),\quad d\zeta^{k}_{I,u}(t) \sim \mathcal{N}\!\big(0,(\sigma^{k}_{I})^{2} dt\big),\quad d\zeta^{k}_{J,u}(t) \sim \mathcal{N}\!\big(0,(\sigma^{k}_{J})^{2} dt\big),\\ d\zeta_{F,u}(t) &\sim \mathcal{N}\!\big(0,\sigma_{F}^{2} dt\big),\quad d\zeta_{P,u}(t) \sim \mathcal{N}\!\big(0,\sigma_{P}^{2} dt\big). \end{aligned}\]
These equations correspond to the SIRJPF2 specification in Section S6 of the Supplement. The spore-production coefficient \(\beta^{k}\) is fixed at \(30\) for both species, \(\delta = 0.013\) day\(^{-1}\), \(\mu = 0.37 \times 10^{6}\) cells \(\cdot\) L\(^{-1}\) day\(^{-1}\), and \(\lambda^{k}_{J} = 0.1\) day\(^{-1}\) for both species. The day-4 inoculum adds 25 internal parasite units: 25,000 spores/L, equivalently 25 spores/mL.
Biological Mechanisms
The SIRJPF2 model couples four ecological and epidemiological processes. Reproduction is resource-limited: susceptible adults produce juveniles at rate \(r^{k} f^{k}_{S} F_{u}(t) S^{k}_{u}(t)\), where the birth efficiency \(r^{k}\) converts ingested algae into offspring and the filtration rate \(f^{k}_{S}\) scales encounter with food. Juveniles mature into the susceptible adult class at rate \(\lambda^{k}_{J}\). Mortality arises from natural causes (rates \(\theta^{k}_{S}\), \(\theta^{k}_{I}\), \(\theta^{k}_{J}\)) and from experimental sampling at rate \(\delta\).
Parasite transmission is environmentally mediated mass action: incidence is proportional to susceptible density and spore density, with coefficient \(p^{k} f^{k}_{S}\). Infected adults experience elevated mortality \(\theta^{k}_{I} > \theta^{k}_{S}\) and, following Searle et al. (2016), do not reproduce. Upon death, an infected individual releases \(\beta^{k} = 30\) internal spore units into the environmental pool. The shared spore pool couples the epidemiological dynamics of the two species.
Resource competition is exploitative: \(F_{u}(t)\) is depleted by all host classes of both species, with the infected-class scaling \(\xi\) representing reduced filtration in diseased individuals. Replenishment at constant rate \(\mu\) reflects the twice-weekly algal supplementation in the Searle et al. (2016) protocol. The complete model and parameter table appear in Section S6 of the Supplement.
Parameter Estimation via Panel Iterated Filtering
Marginalized panel iterated filtering (MPIF) is designed for panel models containing unit-specific as well as shared parameters. In Pypomp, block=True selects MPIF and changes the resampling treatment of unit-specific parameter particles across panel units. This use of “block” should not be confused with temporal block sampling or spatial block particle filters. Those terms refer to different algorithms. The marginalized-Bayes-map interpretation of MPIF is developed by Wheeler, Abkemeier, and Ionides (2025).
The SIRJPF2 specification used here contains only shared parameters. Consequently, it provides no unit-specific parameter particles for MPIF to marginalize, and a block=True versus block=False comparison is not informative for this model. We therefore use standard panel iterated filtering (block=False) for SIRJPF2 and present the MPIF implementation in the unit-specific SRJF example above.
Below we present the complete Pypomp implementation of standard panel iterated filtering for the SIRJPF2 model. It reuses the SRJF construction pattern and adds an eight-state Euler-Maruyama process simulator, deterministic initial states, a four-dimensional negative-binomial measurement model, and log transformations for positive parameters.
Show 8-state SDE simulator with parasite dynamics
# Eight latent states + four non-negative observables (T_*) + error_count.
SIRJPF_STATENAMES = [
"Sn", "In", "Jn", "Si", "Ii", "Ji", "F", "P",
"T_Sn", "T_In", "T_Si", "T_Ii", "error_count",
]
def sirjpf_rproc(X_, theta_, key, covars, t, dt):
"""One Euler-Maruyama step of the SIRJPF2 SDE system."""
Sn, Jn, In = X_["Sn"], X_["Jn"], X_["In"]
Si, Ji, Ii = X_["Si"], X_["Ji"], X_["Ii"]
F, P = X_["F"], X_["P"]
error_count = X_["error_count"]
sigSn, sigIn = theta_["sigSn"], theta_["sigIn"]
sigSi, sigIi = theta_["sigSi"], theta_["sigIi"]
sigJn, sigJi = theta_["sigJn"], theta_["sigJi"]
sigF, sigP = theta_["sigF"], theta_["sigP"]
theta_Sn, theta_In = theta_["theta_Sn"], theta_["theta_In"]
theta_Si, theta_Ii = theta_["theta_Si"], theta_["theta_Ii"]
theta_Jn, theta_Ji = theta_["theta_Jn"], theta_["theta_Ji"]
theta_P = theta_["theta_P"]
f_Sn, f_Si = theta_["f_Sn"], theta_["f_Si"]
rn, ri = theta_["rn"], theta_["ri"]
probn, probi = theta_["probn"], theta_["probi"]
xi = theta_["xi"]
# Fixed experimental constants
delta = 0.013 # sampling/dilution rate
mu_food = 0.37 # algal replenishment
lambda_J = 0.1 # juvenile maturation rate (both species)
# Eight independent Gaussian innovations
keys = jax.random.split(key, 8)
sqdt = jnp.sqrt(dt)
noiSn = sigSn * sqdt * jax.random.normal(keys[0])
noiIn = sigIn * sqdt * jax.random.normal(keys[1])
noiSi = sigSi * sqdt * jax.random.normal(keys[2])
noiIi = sigIi * sqdt * jax.random.normal(keys[3])
noiJn = sigJn * sqdt * jax.random.normal(keys[4])
noiJi = sigJi * sqdt * jax.random.normal(keys[5])
noiF = sigF * sqdt * jax.random.normal(keys[6])
noiP = sigP * sqdt * jax.random.normal(keys[7])
# Native species (n)
Sn_term = (lambda_J * Jn * dt
- theta_Sn * Sn * dt
- probn * f_Sn * Sn * P * dt
- delta * Sn * dt + Sn * noiSn)
Jn_term = (rn * f_Sn * F * Sn * dt
- lambda_J * Jn * dt
- theta_Jn * Jn * dt
- delta * Jn * dt + Jn * noiJn)
In_term = (probn * f_Sn * Sn * P * dt
- theta_In * In * dt
- delta * In * dt + In * noiIn)
# Invasive species (l, written `i` in code)
Si_term = (lambda_J * Ji * dt
- theta_Si * Si * dt
- probi * f_Si * Si * P * dt
- delta * Si * dt + Si * noiSi)
Ji_term = (ri * f_Si * F * Si * dt
- lambda_J * Ji * dt
- theta_Ji * Ji * dt
- delta * Ji * dt + Ji * noiJi)
Ii_term = (probi * f_Si * Si * P * dt
- theta_Ii * Ii * dt
- delta * Ii * dt + Ii * noiIi)
# Shared resources
F_term = (- f_Sn * F * (Sn + xi * In + Jn) * dt
- f_Si * F * (Si + xi * Ii + Ji) * dt
- delta * F * dt
+ mu_food * dt + F * noiF)
P_term = (30.0 * theta_In * In * dt
+ 30.0 * theta_Ii * Ii * dt
- f_Sn * (Sn + xi * In) * P * dt
- f_Si * (Si + xi * Ii) * P * dt
- theta_P * P * dt
- delta * P * dt + P * noiP)
Sn_new = Sn + Sn_term
In_new = In + In_term
Jn_new = Jn + Jn_term
Si_new = Si + Si_term
Ii_new = Ii + Ii_term
Ji_new = Ji + Ji_term
F_new = F + F_term
P_new = P + P_term
# Add the inoculum exactly once on the transition beginning at day 4.
inoculate = (t <= 4.0) & ((t + dt) > 4.0)
P_new = P_new + jnp.where(inoculate, 25.0, 0.0)
# Production reset rule with weighted error_count contributions.
def viol(x, hi):
return (x < 0.0) | (x > hi)
eps = 0.0
eps += jnp.where(viol(Sn_new, 1e5), 1.0, 0.0)
eps += jnp.where(viol(Si_new, 1e5), 1.0e6, 0.0)
eps += jnp.where(viol(F_new, 1e20), 1.0e3, 0.0)
eps += jnp.where(viol(In_new, 1e5), 1.0e-3, 0.0)
eps += jnp.where(viol(Ii_new, 1e5), 1.0e-9, 0.0)
eps += jnp.where(viol(Jn_new, 1e5), 1.0e-3, 0.0)
eps += jnp.where(viol(Ji_new, 1e5), 1.0e-9, 0.0)
eps += jnp.where(viol(P_new, 1e20) & (t > 3.9), 1.0e-6, 0.0)
Sn_new = jnp.where(viol(Sn_new, 1e5), 0.0, Sn_new)
In_new = jnp.where(viol(In_new, 1e5), 0.0, In_new)
Jn_new = jnp.where(viol(Jn_new, 1e5), 0.0, Jn_new)
Si_new = jnp.where(viol(Si_new, 1e5), 0.0, Si_new)
Ii_new = jnp.where(viol(Ii_new, 1e5), 0.0, Ii_new)
Ji_new = jnp.where(viol(Ji_new, 1e5), 0.0, Ji_new)
F_new = jnp.where(viol(F_new, 1e20), 0.0, F_new)
P_new = jnp.where(viol(P_new, 1e20) & (t > 3.9), 0.0, P_new)
return {
"Sn": Sn_new, "In": In_new, "Jn": Jn_new,
"Si": Si_new, "Ii": Ii_new, "Ji": Ji_new,
"F": F_new, "P": P_new,
"T_Sn": jnp.abs(Sn_new), "T_In": jnp.abs(In_new),
"T_Si": jnp.abs(Si_new), "T_Ii": jnp.abs(Ii_new),
"error_count": error_count + eps,
}The simulator returns an updated state dictionary rather than mutating in place. Invalid states are reset to zero with jnp.where, matching the production transition rule while remaining JIT-compatible. The error_count accumulator is reset at every observation time via accumvars, so its measurement penalty is local to an observation interval.
Show initial-state specification
def sirjpf_rinit(theta_, key, covars, t0):
"""Deterministic initial conditions; t0 = 1."""
return {
"Sn": jnp.array(2.333), # 35 / 15 L
"In": jnp.array(0.0),
"Jn": jnp.array(0.0),
"Si": jnp.array(0.667), # 10 / 15 L
"Ii": jnp.array(0.0),
"Ji": jnp.array(0.0),
"F": jnp.array(16.667), # 250e6 cells / 15 L
"P": jnp.array(0.0),
"T_Sn": jnp.array(0.0),
"T_In": jnp.array(0.0),
"T_Si": jnp.array(0.0),
"T_Ii": jnp.array(0.0),
"error_count": jnp.array(0.0),
}The four-dimensional measurement density is the product of four independent negative-binomial likelihoods, one per observed count, with a per-stream overdispersion parameter (\(k_{S}^{n}\), \(k_{S}^{l}\), \(k_{I}^{n}\), \(k_{I}^{l}\)). As in Section 1, the log-pmf is computed via gammaln because jax.scipy.stats.nbinom.logpmf expects integer counts and the \((n, p)\) parameterisation, both awkward inside a JIT-compiled particle filter.
Show 4-D negative-binomial log-density
def _sirjpf_nb_logpmf(y, mu, size):
mu = jnp.maximum(mu, 1e-10)
size = jnp.maximum(size, 1e-10)
return (jax.scipy.special.gammaln(y + size)
- jax.scipy.special.gammaln(size)
- jax.scipy.special.gammaln(y + 1.0)
+ size * jnp.log(size / (size + mu))
+ y * jnp.log(mu / (size + mu)))
def sirjpf_dmeas(Y_, X_, theta_, covars, t):
"""4-D NB log-pmf: dent.adult, dent.inf, lum.adult, lum.adult.inf."""
error_count = X_["error_count"]
ll_dent_adult = _sirjpf_nb_logpmf(Y_["dentadult"], X_["T_Sn"], theta_["k_Sn"])
ll_dent_inf = _sirjpf_nb_logpmf(Y_["dentinf"], X_["T_In"], theta_["k_In"])
ll_lum_adult = _sirjpf_nb_logpmf(Y_["lumadult"], X_["T_Si"], theta_["k_Si"])
ll_lum_inf = _sirjpf_nb_logpmf(Y_["luminf"], X_["T_Ii"], theta_["k_Ii"])
ll = ll_dent_adult + ll_dent_inf + ll_lum_adult + ll_lum_inf
# Soft penalty when the simulator violated a state bound during the step.
return jnp.where(error_count > 0.0, -150.0, ll)Show 4-D measurement sampler
def _sirjpf_nb_sample(key, mu, size):
mu = jnp.maximum(mu, 1e-10)
size = jnp.maximum(size, 1e-10)
k1, k2 = jax.random.split(key)
scale = mu / size
g = jax.random.gamma(k1, size) * scale
return jax.random.poisson(k2, g)
def sirjpf_rmeas(X_, theta_, key, covars, t):
"""Simulate the four-dimensional NB observation."""
keys = jax.random.split(key, 4)
y_dent_adult = _sirjpf_nb_sample(keys[0], X_["T_Sn"], theta_["k_Sn"])
y_dent_inf = _sirjpf_nb_sample(keys[1], X_["T_In"], theta_["k_In"])
y_lum_adult = _sirjpf_nb_sample(keys[2], X_["T_Si"], theta_["k_Si"])
y_lum_inf = _sirjpf_nb_sample(keys[3], X_["T_Ii"], theta_["k_Ii"])
# Keyed by observation name, as Pypomp 1.0.0 requires; see `srjf_rmeas`.
return {
"dentadult": y_dent_adult.astype(float),
"dentinf": y_dent_inf.astype(float),
"lumadult": y_lum_adult.astype(float),
"luminf": y_lum_inf.astype(float),
}The parameter transformation log-transforms every strictly positive parameter, leaving the two fixed-at-zero process-noise terms \(\sigma^{n}_{S} = \sigma^{l}_{S} = 0\) untouched. As in Section 1, MIF perturbation is applied on the log scale for the perturbed parameters and the two fixed nuisance terms have zero random-walk standard deviation.
Show parameter transformations
# Every positive parameter is log-transformed. sigSn and sigSi are fixed at 0
# in MIF (rw_sd = 0) and therefore pass through to_est / from_est unchanged.
_SIRJPF_LOG_PARAMS = (
"rn", "ri", "f_Sn", "f_Si", "probn", "probi", "xi",
"theta_Sn", "theta_Si", "theta_In", "theta_Ii",
"theta_Jn", "theta_Ji", "theta_P",
"sigIn", "sigIi", "sigJn", "sigJi", "sigF", "sigP",
"k_Sn", "k_Si", "k_In", "k_Ii",
)
def sirjpf_to_est(theta):
"""Natural -> estimation scale."""
out = {**theta}
for n in _SIRJPF_LOG_PARAMS:
out[n] = jnp.log(jnp.maximum(theta[n], 1e-30))
out["sigSn"] = theta["sigSn"]
out["sigSi"] = theta["sigSi"]
return out
def sirjpf_from_est(theta):
"""Estimation -> natural scale."""
out = {**theta}
for n in _SIRJPF_LOG_PARAMS:
out[n] = jnp.exp(theta[n])
out["sigSn"] = theta["sigSn"]
out["sigSi"] = theta["sigSi"]
return out
sirjpf_par_trans = pp.ParTrans(to_est=sirjpf_to_est, from_est=sirjpf_from_est)The shared starting parameter vector matches the R tutorial. We set \(t_{0} = 1\), leaving a six-day pre-observation integration window during which the parasite inoculum is added at \(t = 4\).
Code
# Shared SIRJPF2 starting vector.
sirjpf_shared_theta = {
"ri": 1.307600e+04, "rn": 5.904676e+01,
"f_Si": 1.838259e-05, "f_Sn": 1.105668e-03,
"probi": 3.110083e+01, "probn": 2.565626e-01,
"xi": 2.865620e+01,
"theta_Sn": 1.479834e-01, "theta_Si": 3.186040e-02,
"theta_Ii": 3.531879e-01, "theta_In": 5.489315e-01,
"theta_P": 2.024991e-02,
"theta_Ji": 1.299562e-04, "theta_Jn": 1.532613e-04,
"sigSn": 0.0, "sigSi": 0.0,
"sigIn": 3.063207e-04, "sigIi": 2.208698e-02,
"sigJi": 2.727418e-01, "sigJn": 2.836891e-01,
"sigF": 1.551729e-01, "sigP": 2.385890e-01,
"k_Ii": 1.241092e+00, "k_In": 1.005756e+00,
"k_Si": 4.715556e+00, "k_Sn": 4.282648e+00,
}
sirjpf_unit_names_full = ['K', 'L', 'M', 'N', 'O', 'P', 'Q', 'S']
t0_sirjpf = 1.0 # 6-day pre-observation window; parasite inoculum at t = 4.
# Build one Pomp object per replicate.
sirjpf_pomp_dict = {}
for u in sirjpf_unit_names_full:
sub = (sirjpf_data[sirjpf_data['rep'] == u]
[['day', 'dent.adult', 'dent.inf', 'lum.adult', 'lum.adult.inf']]
.rename(columns={
'dent.adult': 'dentadult',
'dent.inf': 'dentinf',
'lum.adult': 'lumadult',
'lum.adult.inf': 'luminf',
})
.sort_values('day'))
ys_u = (sub.set_index('day')
[['dentadult', 'dentinf', 'lumadult', 'luminf']]
.astype(float))
sirjpf_pomp_dict[u] = pp.Pomp(
ys=ys_u,
theta=pp.PompParameters(sirjpf_shared_theta),
statenames=SIRJPF_STATENAMES,
t0=t0_sirjpf,
rinit=sirjpf_rinit,
rproc=sirjpf_rproc,
dmeas=sirjpf_dmeas,
rmeas=sirjpf_rmeas,
par_trans=sirjpf_par_trans,
dt=0.25,
accumvars=("error_count",),
)
# All-shared parameterisation.
sirjpf_shared_df = pd.DataFrame(
{"shared": list(sirjpf_shared_theta.values())},
index=list(sirjpf_shared_theta.keys()),
)
sirjpf_unit_specific_df = pd.DataFrame(
index=[], columns=sirjpf_unit_names_full,
)
panelfood_sirjpf = pp.PanelPomp(
pomp_dict=sirjpf_pomp_dict,
theta=pp.PanelParameters(
{"shared": sirjpf_shared_df,
"unit_specific": sirjpf_unit_specific_df}
),
)
print(f"Number of units: {len(panelfood_sirjpf.unit_objects)}")
print(f"Shared parameters: {len(panelfood_sirjpf.canonical_shared_param_names)}")Number of units: 8
Shared parameters: 26A smoke check confirms the constructed object is correctly scaled.
Code
def _run_sirjpf_smoke():
panel = pp.PanelPomp(
pomp_dict=sirjpf_pomp_dict,
theta=pp.PanelParameters(panelfood_sirjpf.theta),
)
panel.pfilter(
J=RL_SIRJPF["J_eval"], reps=RL_SIRJPF["pf_reps"],
key=jax.random.key(0)
)
return np.asarray(panel.results_history[-1].logLiks.values)
sirjpf_ll0 = cached("sirjpf_smoke_pfilter", _run_sirjpf_smoke) # (theta, U, reps)
sirjpf_unit_ll0_all, sirjpf_unit_se0_all, sirjpf_panel_ll0_all, sirjpf_panel_se0_all = (
summarize_panel_loglik(sirjpf_ll0)
)
sirjpf_panel_ll0 = float(sirjpf_panel_ll0_all[0])
sirjpf_panel_se0 = float(sirjpf_panel_se0_all[0])
print(f"panel pfilter wall: {sf(TIMINGS['sirjpf_smoke_pfilter'])}s "
f"(J={RL_SIRJPF['J_eval']}, reps={RL_SIRJPF['pf_reps']})")
print(f"Initial panel logLik: {ll(sirjpf_panel_ll0)} "
f"(MCSE {sf(sirjpf_panel_se0)})")
print("Per-unit logLik (log-mean-exp over replicates):")
for u, v in zip(sirjpf_unit_names_full,
sirjpf_unit_ll0_all[0]):
print(f" {u}: {ll(v):>10s}")panel pfilter wall: 8.138s (J=1000, reps=10)
Initial panel logLik: -1081.31 (MCSE 8.358)
Per-unit logLik (log-mean-exp over replicates):
K: -128.34
L: -121.27
M: -133.08
N: -141.45
O: -111.95
P: -136.54
Q: -150.24
S: -158.44We refine the all-shared R tutorial starting vector with standard PIF. The two fixed process-noise terms remain shared and have random-walk standard deviation zero. The other 24 parameters are perturbed on the log scale with \(\sigma_{\mathrm{rw}} = 0.05\). Medium, poor, and extreme starting vectors are evaluated with the same first-stage settings. The strongest distinct initial and terminal candidates then receive the same refinement and high-particle likelihood evaluation.
Code
sirjpf_shared_keys_est = list(sirjpf_shared_theta)
sirjpf_panel_shared_df = pd.DataFrame(
{"shared": [sirjpf_shared_theta[k] for k in sirjpf_shared_keys_est]},
index=sirjpf_shared_keys_est,
)
sirjpf_panel_unit_df = pd.DataFrame(index=[], columns=sirjpf_unit_names_full)
sirjpf_start_qualities = start_quality_schedule(RL_SIRJPF["n_starts"])
sirjpf_panel_theta_multi, sirjpf_start_qualities = make_quality_panel_starts(
sirjpf_panel_shared_df,
sirjpf_panel_unit_df,
seed=500,
fixed=("sigSn", "sigSi"),
qualities=sirjpf_start_qualities,
)
# Random walk: 0.05 on every estimated parameter (log scale), 0 on fixed ones.
sirjpf_rw_value = 0.05
sirjpf_rw_sigmas = {k: sirjpf_rw_value for k in sirjpf_shared_theta}
sirjpf_rw_sigmas["sigSn"] = 0.0
sirjpf_rw_sigmas["sigSi"] = 0.0
sirjpf_rw_sd = pp.RWSigma(
sigmas=sirjpf_rw_sigmas, init_names=[]
).geometric_cooling(0.7)
sirjpf_refine_rw_sd = pp.RWSigma(
sigmas=sirjpf_rw_sigmas, init_names=[]
).geometric_cooling(0.8)
def _copy_panel_payloads(theta):
"""Deep-copy PanelParameters payloads for reproducible candidate pools.
`params(as_list=True)` normalises an all-shared panel's unit block to
None, but PanelPomp reads unit names from the columns of that block and
refuses a theta whose unit names do not match its pomp_dict. The empty
frame is therefore rebuilt from the source object's own unit names rather
than propagated as None, so a copied pool can be used to construct a new
panel.
"""
theta_unit_names = list(theta.get_unit_names())
copied = []
for payload in theta.params(as_list=True):
unit_payload = payload["unit_specific"]
copied.append({
"shared": payload["shared"].copy(deep=True),
"unit_specific": (
pd.DataFrame(index=[], columns=theta_unit_names)
if unit_payload is None
else unit_payload.copy(deep=True)
),
})
return copied
def _panel_payload_signature(payload):
"""Numerical signature used to avoid refining duplicate candidates."""
values = []
for block in ("shared", "unit_specific"):
frame = payload[block]
if frame is not None and frame.size:
values.extend(np.asarray(frame, dtype=float).ravel().tolist())
return tuple(np.round(np.asarray(values, dtype=float), 10))
def _select_distinct_payloads(payloads, scores, count):
"""Select the strongest numerically distinct parameter payloads."""
order = np.argsort(np.where(np.isfinite(scores), scores, -np.inf))[::-1]
selected = []
signatures = set()
for idx in order:
signature = _panel_payload_signature(payloads[int(idx)])
if signature in signatures:
continue
selected.append(payloads[int(idx)])
signatures.add(signature)
if len(selected) == count:
break
if not selected:
selected.append(payloads[0])
return selected
def _run_sirjpf_mif():
panel = pp.PanelPomp(
pomp_dict=sirjpf_pomp_dict, theta=sirjpf_panel_theta_multi,
)
initial_payloads = _copy_panel_payloads(panel.theta)
panel.pfilter(
J=RL_SIRJPF["J_eval"], reps=RL_SIRJPF["pf_reps"],
key=jax.random.key(601)
)
initial_ll = np.asarray(panel.results_history[-1].logLiks.values)
panel.mif(
J=RL_SIRJPF["J"], M=RL_SIRJPF["Nmif"], rw_sd=sirjpf_rw_sd,
block=False,
key=jax.random.key(501),
)
terminal_payloads = _copy_panel_payloads(panel.theta)
panel.pfilter(
J=RL_SIRJPF["J_eval"], reps=RL_SIRJPF["pf_reps"],
key=jax.random.key(601)
)
terminal_ll = np.asarray(panel.results_history[-1].logLiks.values)
_, _, initial_panel_ll, _ = summarize_panel_loglik(initial_ll)
_, _, terminal_panel_ll, _ = summarize_panel_loglik(terminal_ll)
stage1_payloads = initial_payloads + terminal_payloads
stage1_scores = np.concatenate([initial_panel_ll, terminal_panel_ll])
refine_payloads = _select_distinct_payloads(
stage1_payloads,
stage1_scores,
count=min(3, len(stage1_payloads)),
)
refine_panel = pp.PanelPomp(
pomp_dict=sirjpf_pomp_dict,
theta=pp.PanelParameters(refine_payloads),
)
refine_initial_payloads = _copy_panel_payloads(refine_panel.theta)
refine_panel.pfilter(
J=RL_SIRJPF["J_eval"], reps=RL_SIRJPF["pf_reps"],
key=jax.random.key(601)
)
refine_initial_ll = np.asarray(
refine_panel.results_history[-1].logLiks.values
)
refine_panel.mif(
J=RL_SIRJPF["J"],
M=RL_SIRJPF["Nmif"],
rw_sd=sirjpf_refine_rw_sd,
block=False,
key=jax.random.key(701),
)
refine_terminal_payloads = _copy_panel_payloads(refine_panel.theta)
refine_panel.pfilter(
J=RL_SIRJPF["J_eval"], reps=RL_SIRJPF["pf_reps"],
key=jax.random.key(601)
)
refine_terminal_ll = np.asarray(
refine_panel.results_history[-1].logLiks.values
)
candidate_ll = np.concatenate(
[initial_ll, terminal_ll, refine_initial_ll, refine_terminal_ll],
axis=0,
)
candidate_theta = (
initial_payloads
+ terminal_payloads
+ refine_initial_payloads
+ refine_terminal_payloads
)
candidate_source = (
[
f"initial {quality} {i}"
for i, quality in enumerate(sirjpf_start_qualities)
]
+ [
f"stage-1 terminal {quality} {i}"
for i, quality in enumerate(sirjpf_start_qualities)
]
+ [f"refinement start {i}" for i in range(len(refine_initial_payloads))]
+ [
f"refinement terminal {i}"
for i in range(len(refine_terminal_payloads))
]
)
return {
"candidate_ll": candidate_ll,
"initial_ll": initial_ll,
"stage1_terminal_ll": terminal_ll,
"candidate_theta": candidate_theta,
"candidate_source": candidate_source,
"refinement_ll": refine_terminal_ll,
"refinement_theta": refine_terminal_payloads,
"stage1_traces": panel.traces(),
"refinement_traces": refine_panel.traces(),
}
sirjpf_mif_out = cached("sirjpf_mif", _run_sirjpf_mif)
sirjpf_mif_ll = sirjpf_mif_out["candidate_ll"]
(
sirjpf_unit_ll_per_start,
sirjpf_unit_se_per_start,
sirjpf_panel_ll_per_start,
sirjpf_panel_se_per_start,
) = summarize_panel_loglik(sirjpf_mif_ll)
sirjpf_best_idx = _finite_argmax(sirjpf_panel_ll_per_start)
sirjpf_best_ll = float(sirjpf_panel_ll_per_start[sirjpf_best_idx])
sirjpf_best_se = float(sirjpf_panel_se_per_start[sirjpf_best_idx])
sirjpf_unit_ll_best = sirjpf_unit_ll_per_start[sirjpf_best_idx]
sirjpf_unit_se_best = sirjpf_unit_se_per_start[sirjpf_best_idx]
sirjpf_best_theta_dict = sirjpf_mif_out["candidate_theta"][sirjpf_best_idx]
sirjpf_selected_shared_df = sirjpf_best_theta_dict["shared"]
sirjpf_selected_theta = {**sirjpf_shared_theta}
for name in sirjpf_selected_shared_df.index:
sirjpf_selected_theta[name] = float(
sirjpf_selected_shared_df.loc[name, "shared"]
)
# The R tutorial starting vector is no longer one of the MIF starts: every
# start is displaced, so none of the candidates is that vector. Its likelihood
# is the one evaluated on its own above, under the same J_eval and replicate
# count, which is what makes the difference below meaningful.
sirjpf_production_ll = sirjpf_panel_ll0
sirjpf_production_se = sirjpf_panel_se0
sirjpf_production_gap_se = float(np.sqrt(
sirjpf_best_se ** 2 + sirjpf_production_se ** 2
))
_, _, sirjpf_stage1_initial_ll, _ = summarize_panel_loglik(
sirjpf_mif_out["initial_ll"]
)
_, _, sirjpf_stage1_terminal_ll, sirjpf_stage1_terminal_se = (
summarize_panel_loglik(sirjpf_mif_out["stage1_terminal_ll"])
)
sirjpf_start_table = pd.DataFrame({
"start": np.arange(len(sirjpf_start_qualities)),
"quality": sirjpf_start_qualities,
"initial_logLik": sirjpf_stage1_initial_ll,
"terminal_logLik": sirjpf_stage1_terminal_ll,
"MCSE": sirjpf_stage1_terminal_se,
"improvement": sirjpf_stage1_terminal_ll - sirjpf_stage1_initial_ll,
})
sirjpf_quality_summary = (
sirjpf_start_table.groupby("quality", sort=False)
.agg(
starts=("start", "count"),
mean_initial_logLik=("initial_logLik", "mean"),
mean_terminal_logLik=("terminal_logLik", "mean"),
best_terminal_logLik=("terminal_logLik", "max"),
mean_improvement=("improvement", "mean"),
)
.reset_index()
)
print("\nMatched all-shared starts:")
print(table_text(sirjpf_start_table))
print("\nSummary by starting-point quality:")
print(table_text(sirjpf_quality_summary))
sirjpf_reference_ll = -880.5596
sirjpf_reference_se = 0.232
sirjpf_reference_gap_se = np.sqrt(
sirjpf_best_se**2 + sirjpf_reference_se**2
)
sirjpf_reference_difference = sirjpf_best_ll - sirjpf_reference_ll
(
_,
_,
sirjpf_refinement_panel_ll,
sirjpf_refinement_panel_se,
) = summarize_panel_loglik(sirjpf_mif_out["refinement_ll"])
sirjpf_refinement_agreement_count = int(np.sum(
sirjpf_best_ll - sirjpf_refinement_panel_ll
<= np.maximum(
1.0,
2.0 * np.sqrt(
sirjpf_refinement_panel_se**2 + sirjpf_best_se**2
),
)
))
sirjpf_refinement_tolerance = np.maximum(
1.0,
2.0 * np.sqrt(sirjpf_refinement_panel_se**2 + sirjpf_best_se**2),
)
sirjpf_refinement_table = pd.DataFrame({
"refinement": np.arange(len(sirjpf_refinement_panel_ll)),
"terminal_logLik": sirjpf_refinement_panel_ll,
"MCSE": sirjpf_refinement_panel_se,
"terminal_minus_selected": sirjpf_refinement_panel_ll - sirjpf_best_ll,
"comparison_tolerance": sirjpf_refinement_tolerance,
})
print(
"Selected minus stored manuscript-reference logLik: "
f"{ll(sirjpf_reference_difference, signed=True)} "
f"(combined MCSE {sf(sirjpf_reference_gap_se)})"
)
print(
"Refinement terminals within the stated MCSE-based tolerance of the "
f"selected candidate: {sirjpf_refinement_agreement_count} of "
f"{len(sirjpf_refinement_panel_ll)}"
)
print(table_text(sirjpf_refinement_table))
Matched all-shared starts:
start quality initial_logLik terminal_logLik MCSE improvement
0 medium -1293.21 -886.51 0.2795 406.70
1 poor -1798.97 -884.43 0.2845 914.53
2 extreme -12000.00 -5410.00 0.2126 6590.00
3 medium -1022.01 -890.28 0.2859 131.72
4 poor -3505.55 -884.66 0.2771 2620.89
5 extreme -11574.90 -900.54 0.7012 10674.36
6 medium -985.17 -888.68 0.1617 96.49
7 poor -4379.71 -886.60 0.3019 3493.12
8 extreme -3890.40 -899.91 0.3462 2990.49
9 medium -1004.29 -883.90 0.2551 120.39
10 poor -2529.51 -887.96 0.1723 1641.55
11 extreme -5289.22 -882.81 0.3019 4406.41
12 medium -1382.83 -886.67 0.254 496.17
13 poor -1899.86 -883.55 0.364 1016.31
14 extreme -12000.00 -887.27 0.2501 11112.73
15 medium -1218.42 -885.75 0.3347 332.66
16 poor -1655.89 -885.67 0.2342 770.22
17 extreme -11044.99 -906.23 0.2653 10138.76
18 medium -1283.29 -885.44 0.2296 397.85
19 poor -3918.87 -891.65 0.2332 3027.22
20 extreme -7171.71 -902.20 0.3439 6269.51
Summary by starting-point quality:
quality starts mean_initial_logLik mean_terminal_logLik best_terminal_logLik mean_improvement
medium 7 -1169.89 -886.75 -883.90 283.14
poor 7 -2812.62 -886.36 -883.55 1926.26
extreme 7 -8995.89 -1541.28 -882.81 7454.61
Selected minus stored manuscript-reference logLik: -2.25 (combined MCSE 0.3807)
Refinement terminals within the stated MCSE-based tolerance of the selected candidate: 0 of 3
refinement terminal_logLik MCSE terminal_minus_selected comparison_tolerance
0 -893.56 0.3194 -10.75 1.00
1 -891.85 0.256 -9.04 1.00
2 -887.02 0.378 -4.21 1.00The SIRJPF2 searches use standard PIF because the fitted specification contains only shared parameters. The three start qualities use the same algorithmic settings, including the extreme group described in Section 1. Terminal parameter estimates are evaluated by replicated particle filters, and the search with the largest estimated panel log likelihood is retained for subsequent inference.
Code
sirjpf_quality_curves = mif_trace_curves(
sirjpf_mif_out["stage1_traces"], sirjpf_start_qualities
)
sirjpf_quality_groups = [
q for q in QUALITY_ORDER if q in set(sirjpf_start_qualities)
]
sirjpf_groups_present = [
q for q in sirjpf_quality_groups
if q in set(sirjpf_quality_curves["group"])
]
sirjpf_groups_absent = [q for q in sirjpf_quality_groups
if q not in sirjpf_groups_present]
fig = plot_quality_curves(
sirjpf_quality_curves, sirjpf_groups_present,
"SIRJPF2 MIF convergence",
)
Code
sirjpf_quality_speed, sirjpf_speed_target = convergence_speed_table(
sirjpf_quality_curves, sirjpf_quality_groups
)
print(
"Convergence speed by starting-point quality. Target is the best final "
f"median less {sf(SPEED_TOLERANCE)} log units = {ll(sirjpf_speed_target)}:"
)
print(table_text(sirjpf_quality_speed))
print(describe_arrivals(sirjpf_quality_speed))
if sirjpf_groups_absent:
print(
"No finite MIF trace at any iteration for: "
+ ", ".join(sirjpf_groups_absent)
+ ". These starts are shown in the table above as NaN and are absent "
"from the figure because they have no plottable value."
)Convergence speed by starting-point quality. Target is the best final median less 10 log units = -897.22:
group first_evaluated_logLik final_median_logLik median_gain iterations_to_target
medium -971.94 -889.19 82.75 62
poor -1246.57 -887.22 359.35 105
extreme -8214.46 -902.71 7311.75 NaN
Reached the target, soonest first: medium at iteration 62, poor at iteration 105. Did not come within 10 log units of the best final median: extreme.The start table measures the improvement from each initial vector, the curves show when the likelihood increases, and the speed table states how many iterations each group needed to come within 10 log units of the best final median. Persistent separation between start groups would indicate that additional search effort is needed before interpreting the selected estimate.
Fit diagnostic: per-unit log-likelihood breakdown
The following figure shows how the eight observed trajectories contribute to the panel total at the selected candidate. It is descriptive rather than a heterogeneity test: Monte Carlo error quantifies computation, not the biological sampling variation between mesocosms.
Code
fig, ax = plt.subplots(figsize=(8, 4))
sirjpf_unit_ll_baseline = np.nanmean(sirjpf_unit_ll_best)
ax.bar(
sirjpf_unit_names_full,
sirjpf_unit_ll_best - sirjpf_unit_ll_baseline,
bottom=sirjpf_unit_ll_baseline,
yerr=2 * sirjpf_unit_se_best,
capsize=3,
color=PALETTE["adult"],
)
ax.axhline(sirjpf_unit_ll_baseline, color=PALETTE["fit"],
linestyle="--", linewidth=1.5, label="panel mean")
ax.set_xlabel("Unit", fontsize=FONTS["axis_label"])
ax.set_ylabel("Unit log-likelihood", fontsize=FONTS["axis_label"])
ax.set_title("SIRJPF2 unit-level log-likelihood at selected estimate",
fontsize=FONTS["panel_title"])
ax.tick_params(axis="both", labelsize=FONTS["tick"])
ax.legend(loc="lower right", fontsize=FONTS["legend"])
fig.tight_layout()
_emit_figure(
fig,
label="fig-sirjpf-mif-unit-decomposition",
caption=(
"Per-unit log-likelihood contributions for the SIRJPF2 panel under "
"the selected estimate. Error bars show $\\pm 2$ Monte Carlo "
"standard errors. The dashed line is the mean unit contribution."
),
)
Matched PIF/MPIF comparison with unit-specific infected mortality
PIF and MPIF differ only when estimated unit-specific parameters are present. We therefore compare them under the SIRJPF2 alternative in which native and invasive infected-adult mortality (theta_In, theta_Ii) vary by mesocosm. Both methods begin from the same medium-quality parameter vectors. Particle counts, iterations, perturbations, cooling, MIF keys, and final-evaluation keys are matched. The trace comparison measures convergence speed during the first stage, and the final table compares independently evaluated terminal likelihoods and Monte Carlo standard errors.
Code
sirjpf_comparison_base = sirjpf_selected_theta
sirjpf_specific_names = ["theta_In", "theta_Ii"]
sirjpf_specific_shared_names = [
k for k in sirjpf_comparison_base if k not in sirjpf_specific_names
]
sirjpf_specific_shared_df = pd.DataFrame(
{"shared": [
sirjpf_comparison_base[k] for k in sirjpf_specific_shared_names
]},
index=sirjpf_specific_shared_names,
)
sirjpf_specific_unit_df = pd.DataFrame(
{
unit: [
sirjpf_comparison_base[name] for name in sirjpf_specific_names
]
for unit in sirjpf_unit_names_full
},
index=sirjpf_specific_names,
)
sirjpf_specific_quality = ["medium"] * RL_SIRJPF["n_starts"]
sirjpf_specific_starts, sirjpf_specific_quality = make_quality_panel_starts(
sirjpf_specific_shared_df,
sirjpf_specific_unit_df,
seed=1500,
fixed=("sigSn", "sigSi"),
qualities=sirjpf_specific_quality,
)
sirjpf_specific_rw = pp.RWSigma(
sigmas=sirjpf_rw_sigmas,
init_names=[],
).geometric_cooling(0.7)
sirjpf_specific_refine_rw = pp.RWSigma(
sigmas=sirjpf_rw_sigmas,
init_names=[],
).geometric_cooling(0.8)
def _run_sirjpf_specific_stage1(block):
panel = pp.PanelPomp(
pomp_dict=sirjpf_pomp_dict,
theta=pp.PanelParameters(sirjpf_specific_starts),
)
initial_theta = _copy_panel_payloads(panel.theta)
panel.pfilter(
J=RL_SIRJPF["J_eval"],
reps=RL_SIRJPF["pf_reps"],
key=jax.random.key(1602),
)
initial_ll = np.asarray(panel.results_history[-1].logLiks.values)
panel.mif(
J=RL_SIRJPF["J"],
M=RL_SIRJPF["Nmif"],
rw_sd=sirjpf_specific_rw,
block=block,
key=jax.random.key(1601),
)
panel.pfilter(
J=RL_SIRJPF["J_eval"],
reps=RL_SIRJPF["pf_reps"],
key=jax.random.key(1602),
)
return {
"initial_ll": initial_ll,
"initial_theta": initial_theta,
"terminal_ll": np.asarray(panel.results_history[-1].logLiks.values),
"terminal_theta": _copy_panel_payloads(panel.theta),
"traces": panel.traces(),
}
sirjpf_pif_stage1 = cached(
"sirjpf_specific_pif_stage1",
lambda: _run_sirjpf_specific_stage1(block=False),
)
sirjpf_mpif_stage1 = cached(
"sirjpf_specific_mpif_stage1",
lambda: _run_sirjpf_specific_stage1(block=True),
)
_, _, specific_initial_panel_ll, _ = summarize_panel_loglik(
sirjpf_pif_stage1["initial_ll"]
)
_, _, specific_pif_stage1_ll, _ = summarize_panel_loglik(
sirjpf_pif_stage1["terminal_ll"]
)
_, _, specific_mpif_stage1_ll, _ = summarize_panel_loglik(
sirjpf_mpif_stage1["terminal_ll"]
)
# mif_trace_curves() drops iterations with no finite likelihood before it
# takes quantiles. Iteration 0 has none, and summarising it made the first
# median NaN, which propagated into the speed statistic below.
pif_stage1_curve = mif_trace_curves(sirjpf_pif_stage1["traces"]).assign(
group="PIF"
)
mpif_stage1_curve = mif_trace_curves(sirjpf_mpif_stage1["traces"]).assign(
group="MPIF"
)
# One table, so both methods are timed against the same target. Timing each
# method against its own gain would let the slower method look faster.
sirjpf_method_speed, sirjpf_method_target = convergence_speed_table(
pd.concat([pif_stage1_curve, mpif_stage1_curve], ignore_index=True),
["PIF", "MPIF"],
)
sirjpf_method_speed = sirjpf_method_speed.rename(columns={"group": "method"})
sirjpf_method_speed.insert(1, "start_quality", "medium")
print(
"Matched medium-start convergence summary. Target is the better final "
f"median less {sf(SPEED_TOLERANCE)} log units = "
f"{ll(sirjpf_method_target)}:"
)
print(table_text(sirjpf_method_speed))
sirjpf_specific_pool_payloads = (
sirjpf_pif_stage1["initial_theta"]
+ sirjpf_pif_stage1["terminal_theta"]
+ sirjpf_mpif_stage1["terminal_theta"]
)
sirjpf_specific_pool_scores = np.concatenate([
specific_initial_panel_ll,
specific_pif_stage1_ll,
specific_mpif_stage1_ll,
])
sirjpf_specific_refine_payloads = _select_distinct_payloads(
sirjpf_specific_pool_payloads,
sirjpf_specific_pool_scores,
count=min(3, len(sirjpf_specific_pool_payloads)),
)
sirjpf_specific_refine_starts = pp.PanelParameters(
sirjpf_specific_refine_payloads
)
def _run_sirjpf_specific_refine(block):
"""Refine one common pooled start set using PIF or MPIF."""
panel = pp.PanelPomp(
pomp_dict=sirjpf_pomp_dict,
theta=pp.PanelParameters(sirjpf_specific_refine_starts),
)
panel.mif(
J=RL_SIRJPF["J"],
M=RL_SIRJPF["Nmif"],
rw_sd=sirjpf_specific_refine_rw,
block=block,
key=jax.random.key(1701),
)
panel.pfilter(
J=RL_SIRJPF["J_eval"],
reps=RL_SIRJPF["pf_reps"],
key=jax.random.key(1702),
)
return {
"ll": np.asarray(panel.results_history[-1].logLiks.values),
"theta": _copy_panel_payloads(panel.theta),
"traces": panel.traces(),
}
sirjpf_pif_out = cached(
"sirjpf_specific_pif_refine",
lambda: _run_sirjpf_specific_refine(block=False),
)
sirjpf_mpif_out = cached(
"sirjpf_specific_mpif_refine",
lambda: _run_sirjpf_specific_refine(block=True),
)
def _method_summary(output):
unit_ll, unit_se, panel_ll, panel_se = summarize_panel_loglik(output["ll"])
best = _finite_argmax(panel_ll)
best_ll = float(panel_ll[best])
best_se = float(panel_se[best])
near_best_count = int(np.sum(
best_ll - panel_ll
<= np.maximum(
1.0,
2.0 * np.sqrt(panel_se**2 + best_se**2),
)
))
return {
"unit_ll": unit_ll,
"unit_se": unit_se,
"panel_ll": panel_ll,
"panel_se": panel_se,
"best": best,
"best_ll": best_ll,
"best_se": best_se,
"near_best_count": near_best_count,
}
sirjpf_pif = _method_summary(sirjpf_pif_out)
sirjpf_mpif = _method_summary(sirjpf_mpif_out)
sirjpf_method_combined_se = float(np.sqrt(
sirjpf_pif["best_se"] ** 2 + sirjpf_mpif["best_se"] ** 2
))
sirjpf_method_difference = sirjpf_pif["best_ll"] - sirjpf_mpif["best_ll"]
sirjpf_method_table = pd.DataFrame({
"method": ["PIF", "MPIF"],
"block": [False, True],
"free_parameters": [38, 38],
"best_logLik": [sirjpf_pif["best_ll"], sirjpf_mpif["best_ll"]],
"MCSE": [sirjpf_pif["best_se"], sirjpf_mpif["best_se"]],
"starts_within_tolerance": [
sirjpf_pif["near_best_count"],
sirjpf_mpif["near_best_count"],
],
"wall_seconds": [
TIMINGS["sirjpf_specific_pif_stage1"]
+ TIMINGS["sirjpf_specific_pif_refine"],
TIMINGS["sirjpf_specific_mpif_stage1"]
+ TIMINGS["sirjpf_specific_mpif_refine"],
],
})
print(table_text(sirjpf_method_table))
print(
"PIF best minus MPIF best logLik: "
f"{ll(sirjpf_method_difference, signed=True)} "
f"(combined MCSE {sf(sirjpf_method_combined_se)})"
)Matched medium-start convergence summary. Target is the better final median less 10 log units = -878.81:
method start_quality first_evaluated_logLik final_median_logLik median_gain iterations_to_target
PIF medium -922.56 -874.05 48.51 151
MPIF medium -922.56 -868.81 53.75 72
method block free_parameters best_logLik MCSE starts_within_tolerance wall_seconds
PIF False 38 -869.53 0.2543 2 88.11
MPIF True 38 -871.59 0.2429 2 86.02
PIF best minus MPIF best logLik: +2.06 (combined MCSE 0.3517)Code
fig, ax = plt.subplots(figsize=(7, 4))
for curve, method, color in (
(pif_stage1_curve, "PIF", PALETTE["adult"]),
(mpif_stage1_curve, "MPIF", PALETTE["infected"]),
):
ax.plot(
curve["iteration"], smooth_trace(curve["median"]), label=method,
color=color, linewidth=1.8,
)
ax.fill_between(
curve["iteration"],
smooth_trace(curve["lower"]),
smooth_trace(curve["upper"]),
color=color, alpha=0.16, linewidth=0,
)
ax.set_xlabel("Iteration", fontsize=FONTS["axis_label"])
ax.set_ylabel("Trace log likelihood", fontsize=FONTS["axis_label"])
ax.set_title("Matched medium-start convergence",
fontsize=FONTS["panel_title"])
ax.legend(frameon=False, fontsize=FONTS["legend"])
ax.tick_params(axis="both", labelsize=FONTS["tick"])
fig.tight_layout()
_emit_figure(
fig,
label="fig-sirjpf-pif-mpif",
caption=(
"Matched PIF and MPIF comparison for unit-specific infected mortality. "
"Median convergence traces and interquartile bands from common "
"medium-quality starts."
),
)
The iterations-to-target column is a descriptive convergence-speed measure against a target shared by both methods. It does not establish that either method reached the global maximum. The final likelihood and MCSE comparison remains necessary because a fast trace can settle at a lower local maximum.
Closing remarks for Section 2
The all-shared and unit-specific calculations report candidate likelihoods, Monte Carlo standard errors, reference differences, nearby terminal starts, and timings. The active analysis uses run level 2. Production SIRJPF2 results are reported in Section S6 of the Supplement.
Section 2 wall times at run_level = 3 (s)
sirjpf_smoke_pfilter 8.138
sirjpf_mif 109.1
sirjpf_specific_pif_stage1 47.64
sirjpf_specific_mpif_stage1 46.19
sirjpf_specific_pif_refine 40.47
sirjpf_specific_mpif_refine 39.83
Akaike, H. 1974. “A New Look at the Statistical Model Identification.” IEEE Transactions on Automatic Control 19 (6): 716–23.
Bretó, Carles, Edward L. Ionides, and Aaron A. King. 2020. “Panel Data Analysis via Mechanistic Models.” Journal of the American Statistical Association 115 (531): 1178–88. https://doi.org/10.1080/01621459.2019.1604367.
Ebert, Dieter. 2005. Ecology, Epidemiology, and Evolution of Parasitism in Daphnia. National Library of Medicine.
Ionides, Edward L, Carles Breto, Joonha Park, Richard A Smith, and Aaron A King. 2017. “Monte Carlo Profile Confidence Intervals for Dynamic Systems.” Journal of the Royal Society Interface 14: 1–10.
Lachance, Marc-André, Carla E. Cáceres, Molly J. Fredericks, Meghan A. Duffy, and Tara E. Stewart Merrill. 2025. “Reviving Élie Metschnikoff’s Monospora: The Obligately Parasitic Yeast Australozyma Monospora Sp. Nov.” FEMS Yeast Research 25: foaf041. https://doi.org/10.1093/femsyr/foaf041.
Ning, N., E. L. Ionides, and Y. Ritov. 2021. “Scalable Monte Carlo Inference and Rescaled Local Asymptotic Normality.” Bernoulli 27: 2532–55.
Searle, Catherine L, Michael H Cortez, Katherine K Hunsberger, Dylan C Grippi, Isabella A Oleksy, Clara L Shaw, Solanus B de la Serna, Chloe L Lash, Kailash L Dhir, and Meghan A Duffy. 2016. “Population Density, Not Host Competence, Drives Patterns of Disease in an Invaded Community.” The American Naturalist 188 (5): 554–66. https://doi.org/10.1086/688402.
Wheeler, Jesse, Aaron J. Abkemeier, and Edward L. Ionides. 2025. “Iterating Marginalized Bayes Maps for Likelihood Maximization.” In Preparation.