Panel Measles Parameter Estimation

Distribution of block-IF2 parameter estimates: R panelPomp vs. pypomp

Published

September 2, 2026

Show Code 
import os
import sys
sys.path.append("..")
import report_utils as ru
from IPython.display import display, HTML

display(HTML(f"<div style='margin-bottom: 20px;'>{ru.nav_bar('estimation')}</div>"))
Likelihood Evaluation  |  Parameter Estimation  |  Timing & Throughput

Introduction

This report compares block-IF2 parameter estimates for the mixed panel measles model. Both implementations run the same number of iterations from the same committed starting points (../starting_parameters.csv) with the same perturbation sizes, so the distribution of final estimates should agree.

The parameter structure is mixed across the panel’s 4 units — London, Halesworth, Hastings and Cardiff:

  • Shared across units: R0, sigma, gamma, sigmaSE, cohort, amplitude.
  • Unit-specific: iota, rho, psi, S_0, E_0, I_0, R_0.

pypomp reports the shared parameters on every unit’s row; they are collapsed here to a single shared entry so a shared density is not drawn four times over.

The python and R smoothing spline implementations interpolate the birthrate and population covariates differently — noticeably so for Cardiff — which on its own shifts the estimates. The final interpolated covariate grid is therefore exported from pomp into ../R_covariates.csv and injected into the pypomp objects, so both packages search over the identical model.

The pypomp run continues past the search into the panel-specific part of the workflow: mix-and-match across units, then pruning to the single best combination. R panelPomp has no equivalent, so the like-for-like comparison below uses results_final.csv, captured before that step. The mix-and-match outcome is reported separately as a pypomp-only result.

Benchmark Settings & Environment

The table below summarizes algorithmic parameters, software environments, and compute hardware recorded in latest.json for each configuration.

Show Code 
PLATFORMS = {
    "R panelPomp": os.path.join("results", "R"),
    "pypomp (GPU)": ru.platform_dir(),
}

runs = ru.load_timing_data(PLATFORMS)
display(HTML(ru.build_settings_comparison_html(runs, is_panel=True)))
Setting / Parameter R panelPomp pypomp (GPU)
Algorithmic & Workload Settings
Run Level 4 4
Starting Searches ($N_{starts}$) 360 360
MPIF Iterations ($N_{iter}$) 100 100
MPIF Particles / Unit ($N_{p,fit}$) 5,000 5,000
Pfilter Particles / Unit ($N_{p,eval}$) 5,000 5,000
Evaluation Replicates ($N_{reps}$) 36 36
Unit(s) London, Halesworth, Hastings, Cardiff (4 units) London, Halesworth, Hastings, Cardiff (4 units)
Shared Parameters — R0, sigma, gamma, sigmaSE, cohort, amplitude (6 params)
Random Seed — 594709947
Software & Environment
Pomp Framework panelPomp 1.7.0.0 (pomp 6.3) pypomp 1.0.0rc1
Backend / Engine R 4.4.0 JAX 0.11.1
Quant Git Commit e6ac592 630cc23
Run Timestamp 2026-08-14 17:06:50 2026-09-02 20:18:19
Hardware & Compute
Compute Device Intel(R) Xeon(R) Gold 6154 CPU @ 3.00GHz (36 cores) NVIDIA RTX PRO 6000 Blackwell Server Edition (1 GPU)
Slurm Partition standard gpu-rtx6000
Slurm Job ID 57274684 59687801
Show Code 
missing = [(label, r["dir"]) for label, r in runs.items() if not r["available"]]
for label, path in missing:
    display(HTML(f"<div class='alert alert-warning'><strong>Missing results for {label}:</strong> Expected at <code>{path}</code>.</div>"))
Show Code 
import numpy as np
import pandas as pd
from plotnine import ggplot, aes, geom_density, facet_wrap, labs

os.environ["JAX_PLATFORMS"] = "cpu"

PY_DIR = ru.platform_dir()
R_DIR = os.path.join("results", "R")

py_results = ru.read_if_exists(os.path.join(PY_DIR, "results_final.csv"))
r_coefs = ru.read_if_exists(os.path.join(R_DIR, "mif_coefs.csv"))
Show Code 
long_frames = []
if py_results is not None:
    long_frames.append(
        ru.split_shared_specific(ru.wide_results_to_long(py_results, "pypomp"))
    )
if r_coefs is not None:
    long_frames.append(ru.r_coefs_to_long(r_coefs))

combined = pd.concat(long_frames, ignore_index=True) if long_frames else pd.DataFrame(
    columns=["unit", "param", "value", "source"]
)
combined = combined[combined["param"].isin(ru.LOG_PARAMS + ru.LOGIT_PARAMS)]
transformed = ru.to_estimation_scale(combined)

Summary Statistics

Reported on the estimation scale (log for rates, logit for probabilities and initial-state fractions), which is the scale block-IF2 perturbs on.

Show Code 
if transformed.empty:
    display(HTML("<p><em>No parameter estimates available.</em></p>"))
else:
    summary = (
        transformed.groupby(["unit", "source", "param"])["value"]
        .agg(n="count", mean="mean", sd="std", min="min", max="max")
        .reset_index()
        .sort_values(["unit", "param", "source"])
    )
    display(HTML(summary.round(4).to_html(classes="table table-striped table-hover", index=False)))
unit source param n mean sd min max
Cardiff R panelPomp log(iota) 360 -1.3846 0.8071 -8.5191 1.3889
Cardiff pypomp log(iota) 360 -1.3865 0.6724 -3.0666 2.1724
Cardiff R panelPomp log(psi) 360 -1.3274 0.0797 -1.5598 -1.1275
Cardiff pypomp log(psi) 360 -1.3289 0.0874 -1.6340 -1.1045
Cardiff R panelPomp logit(E_0) 360 -12.2716 1.2066 -17.7719 -10.1121
Cardiff pypomp logit(E_0) 360 -12.0875 1.1604 -16.7719 -9.9929
Cardiff R panelPomp logit(I_0) 360 -11.9284 0.9089 -15.9739 -10.1714
Cardiff pypomp logit(I_0) 360 -11.9939 0.8656 -14.9956 -10.3695
Cardiff R panelPomp logit(R_0) 360 3.1835 0.4680 -1.3559 4.0635
Cardiff pypomp logit(R_0) 360 3.1413 0.4445 0.0176 4.1495
Cardiff R panelPomp logit(S_0) 360 -3.1841 0.4683 -4.0656 1.3558
Cardiff pypomp logit(S_0) 360 -3.1419 0.4448 -4.1509 -0.0177
Cardiff R panelPomp logit(rho) 360 0.4290 0.1165 0.1532 1.3888
Cardiff pypomp logit(rho) 360 0.4334 0.1091 0.1737 1.0183
Halesworth R panelPomp log(iota) 360 -3.5921 0.5524 -5.2829 -2.3170
Halesworth pypomp log(iota) 360 -3.6110 0.5878 -5.4775 -2.1657
Halesworth R panelPomp log(psi) 360 -0.2210 0.3414 -3.2855 0.5357
Halesworth pypomp log(psi) 360 -0.2513 0.3436 -2.4590 0.4673
Halesworth R panelPomp logit(E_0) 360 -9.9699 1.8703 -16.2223 -6.2764
Halesworth pypomp logit(E_0) 360 -10.0486 1.7792 -15.4198 -6.6043
Halesworth R panelPomp logit(I_0) 360 -8.0779 1.0023 -15.8802 -6.0016
Halesworth pypomp logit(I_0) 360 -8.0852 1.0159 -14.3419 -6.4938
Halesworth R panelPomp logit(R_0) 360 2.9819 0.4057 -0.6082 3.8129
Halesworth pypomp logit(R_0) 360 2.9529 0.4120 -0.2377 3.9218
Halesworth R panelPomp logit(S_0) 360 -2.9964 0.4116 -3.8271 0.5992
Halesworth pypomp logit(S_0) 360 -2.9662 0.4168 -3.9579 0.2376
Halesworth R panelPomp logit(rho) 360 1.1668 0.5406 0.1084 4.3023
Halesworth pypomp logit(rho) 360 1.2369 0.6357 0.0478 4.6529
Hastings R panelPomp log(iota) 360 -1.8849 0.7394 -3.6612 5.6518
Hastings pypomp log(iota) 360 -1.9034 0.6519 -3.7708 0.0863
Hastings R panelPomp log(psi) 360 -0.7913 0.1037 -1.0533 -0.4748
Hastings pypomp log(psi) 360 -0.7867 0.1079 -1.1025 -0.4372
Hastings R panelPomp logit(E_0) 360 -12.7748 1.1655 -18.6498 -9.8697
Hastings pypomp logit(E_0) 360 -12.8158 1.0917 -16.6527 -10.2459
Hastings R panelPomp logit(I_0) 360 -12.2414 0.9971 -15.7087 -9.9430
Hastings pypomp logit(I_0) 360 -12.3074 0.9222 -15.3185 -9.6045
Hastings R panelPomp logit(R_0) 360 3.5908 0.5538 0.9506 5.9753
Hastings pypomp logit(R_0) 360 3.5425 0.5765 0.1325 4.9867
Hastings R panelPomp logit(S_0) 360 -3.5914 0.5542 -5.9783 -0.9506
Hastings pypomp logit(S_0) 360 -3.5430 0.5767 -4.9900 -0.1326
Hastings R panelPomp logit(rho) 360 0.9042 0.3588 0.3281 4.4942
Hastings pypomp logit(rho) 360 0.9499 0.4138 0.2682 3.8125
London R panelPomp log(iota) 360 -0.3036 2.5104 -12.2456 5.5178
London pypomp log(iota) 360 -0.5266 2.4967 -11.1779 4.5633
London R panelPomp log(psi) 360 -2.1882 0.0690 -2.4044 -2.0400
London pypomp log(psi) 360 -2.1926 0.0747 -2.4554 -2.0119
London R panelPomp logit(E_0) 360 -11.0339 1.3433 -16.7741 -9.0569
London pypomp logit(E_0) 360 -10.9341 1.1250 -15.9268 -9.2685
London R panelPomp logit(I_0) 360 -10.7906 0.7941 -13.7053 -8.6269
London pypomp logit(I_0) 360 -10.8029 0.8057 -14.5075 -9.2203
London R panelPomp logit(R_0) 360 3.1025 0.4073 -0.5147 3.8509
London pypomp logit(R_0) 360 3.0676 0.4275 -0.4781 3.9755
London R panelPomp logit(S_0) 360 -3.1041 0.4082 -3.8558 0.5137
London pypomp logit(S_0) 360 -3.0692 0.4284 -3.9829 0.4776
London R panelPomp logit(rho) 360 -0.0349 0.0821 -1.1646 0.1578
London pypomp logit(rho) 360 -0.0243 0.0709 -0.2981 0.9426
shared R panelPomp log(R0) 360 3.3543 0.4566 0.8974 4.3441
shared pypomp log(R0) 360 3.3177 0.4703 0.9048 4.5728
shared R panelPomp log(gamma) 360 4.4236 0.9585 2.9219 12.3555
shared pypomp log(gamma) 360 4.4408 0.9347 3.0492 9.3928
shared R panelPomp log(sigma) 360 4.1788 0.6636 3.0463 6.8606
shared pypomp log(sigma) 360 4.1773 0.6354 2.9454 6.7107
shared R panelPomp log(sigmaSE) 360 -2.9606 0.6861 -4.6905 -1.0842
shared pypomp log(sigmaSE) 360 -2.9796 0.7168 -4.9314 -1.2002
shared R panelPomp logit(amplitude) 360 -1.0288 0.9994 -9.7369 4.6968
shared pypomp logit(amplitude) 360 -1.0263 0.8966 -2.6614 2.6985
shared R panelPomp logit(cohort) 360 -0.2529 2.5117 -9.0060 8.3077
shared pypomp logit(cohort) 360 -0.2546 2.7761 -9.5456 8.2214

Shared Parameter Distributions

The parameters estimated jointly across all four units.

Show Code 
shared = transformed[transformed["unit"] == "shared"]
if shared.empty:
    display(HTML("<p><em>No shared parameter estimates available.</em></p>"))
else:
    display(
        ggplot(shared, aes(x="value", fill="source", color="source"))
        + geom_density(alpha=0.35)
        + facet_wrap("~param", scales="free")
        + labs(
            title="Block-IF2 estimates of the shared parameters",
            subtitle="Estimation scale; one density per implementation",
            x="Parameter value",
            y="Density",
            fill="Source",
            color="Source",
        )
        + ru.scale_fill_premium()
        + ru.scale_color_premium()
        + ru.theme_premium
    )

Unit-Specific Parameter Distributions

Show Code 
specific = transformed[transformed["unit"] != "shared"]
if specific.empty:
    display(HTML("<p><em>No unit-specific parameter estimates available.</em></p>"))
else:
    for unit in sorted(specific["unit"].dropna().unique()):
        unit_data = specific[specific["unit"] == unit]
        display(
            ggplot(unit_data, aes(x="value", fill="source", color="source"))
            + geom_density(alpha=0.35)
            + facet_wrap("~param", scales="free")
            + labs(
                title=f"Block-IF2 estimates of the unit-specific parameters: {unit}",
                subtitle="Estimation scale; one density per implementation",
                x="Parameter value",
                y="Density",
                fill="Source",
                color="Source",
            )
            + ru.scale_fill_premium()
            + ru.scale_color_premium()
            + ru.theme_premium
        )

Expectation

The pypomp and panelPomp densities should overlap heavily for each parameter, both shared and unit-specific.

Log-Likelihood at the Search Endpoints

Each start’s final parameter vector is re-evaluated by particle filter. The distribution over starts shows how consistently the search reaches the same region.

Show Code 
py_ll = ru.read_if_exists(os.path.join(PY_DIR, "pfilter_logliks.csv"))
r_ll = ru.read_if_exists(os.path.join(R_DIR, "pfilter_logliks.csv"))

ll_frames = []
if py_ll is not None:
    py_ll = py_ll.copy()
    py_ll["source"] = "pypomp"
    ll_frames.append(py_ll)
if r_ll is not None:
    r_ll = r_ll.copy()
    r_ll["source"] = "R panelPomp"
    ll_frames.append(r_ll)

missing_ll = []
if py_ll is None:
    missing_ll.append(("pypomp endpoint log-likelihoods", os.path.join(PY_DIR, "pfilter_logliks.csv")))
if r_ll is None:
    missing_ll.append(("R panelPomp endpoint log-likelihoods", os.path.join(R_DIR, "pfilter_logliks.csv")))

for label, path in missing_ll:
    display(HTML(ru.missing_note(label, path)))

if ll_frames:
    combined_ll = pd.concat(ll_frames, ignore_index=True)
    per_unit = (
        combined_ll.groupby(["source", "unit", "theta_idx"])["logLik"].agg(logmeanexp=ru.logmeanexp).reset_index()
        if "theta_idx" in combined_ll.columns
        else combined_ll.groupby(["source", "unit"])["logLik"].agg(logmeanexp=ru.logmeanexp).reset_index()
    )
    display(
        ggplot(per_unit, aes(x="logmeanexp", fill="source", color="source"))
        + geom_density(alpha=0.35)
        + facet_wrap("~unit", scales="free")
        + labs(
            title="Per-unit log-likelihood across search endpoints",
            subtitle="One logmeanexp value per start per implementation",
            x="Log-likelihood",
            y="Density",
            fill="Source",
            color="Source",
        )
        + ru.scale_fill_premium()
        + ru.scale_color_premium()
        + ru.theme_premium
    )

Mix-and-Match

panelPomp fits each start as one coherent parameter vector. pypomp can additionally recombine the best unit-specific blocks across starts and re-fit, which has no counterpart in the R baseline. The table below is the pruned best combination after that step.

Show Code 
pruned = ru.read_if_exists(os.path.join(PY_DIR, "results.csv"))
if pruned is None:
    display(HTML(ru.missing_note("pypomp mix-and-match result",
                                 os.path.join(PY_DIR, "results.csv"))))
else:
    display(HTML(pruned.round(4).to_html(classes="table table-striped table-hover", index=False)))
theta_idx shared logLik shared logLik se unit unit logLik unit logLik se R0 sigma gamma sigmaSE cohort amplitude iota rho psi S_0 E_0 I_0 R_0
0 -10037.0817 8.8164 London -5701.3041 8.7045 18.6382 66.4044 240.4364 0.0229 0.4636 0.1717 1.7311 0.4913 0.1014 0.0350 0.0001 0.0000 0.9649
0 -10037.0817 8.8164 Halesworth -349.6454 0.2334 18.6382 66.4044 240.4364 0.0229 0.4636 0.1717 0.0525 0.7474 0.6346 0.0475 0.0000 0.0007 0.9519
0 -10037.0817 8.8164 Hastings -1591.5398 0.3840 18.6382 66.4044 240.4364 0.0229 0.4636 0.1717 0.0953 0.7627 0.4007 0.0439 0.0000 0.0000 0.9561
0 -10037.0817 8.8164 Cardiff -2394.5925 1.3257 18.6382 66.4044 240.4364 0.0229 0.4636 0.1717 0.2448 0.6076 0.2455 0.0484 0.0000 0.0000 0.9516

Convergence Traces

Show Code 
from plotnine import geom_line

target_cols = ["unit", "theta_idx", "replicate", "iteration", "method", "logLik", "loglik"] + ru.LOG_PARAMS + ru.LOGIT_PARAMS

py_traces = ru.load_traces(os.path.join(PY_DIR, "traces.csv.gz"), "pypomp", usecols=target_cols, max_iters=20)

r_trace_path = os.path.join(R_DIR, "mif_traces.csv.gz")
if not os.path.exists(r_trace_path):
    r_trace_path = os.path.join(R_DIR, "traces.csv.gz")
r_traces = ru.load_traces(r_trace_path, "R panelPomp", usecols=target_cols, max_iters=20)

trace_frames = []
if py_traces is not None:
    if "method" in py_traces.columns:
        py_traces = py_traces[py_traces["method"] == "mif"]
    py_traces = ru.split_traces_shared_specific(py_traces)
    py_long = ru.to_long(py_traces, max_iters=20)
    trace_frames.append(py_long)

if r_traces is not None:
    if "method" in r_traces.columns:
        r_traces = r_traces[r_traces["method"] == "mif"]
    r_traces = ru.split_traces_shared_specific(r_traces)
    r_long = ru.to_long(r_traces, max_iters=20)
    trace_frames.append(r_long)

missing_traces = []
if py_traces is None:
    missing_traces.append(("pypomp traces", os.path.join(PY_DIR, "traces.csv.gz")))
if r_traces is None:
    missing_traces.append(("R panelPomp traces", os.path.join(R_DIR, "mif_traces.csv.gz")))

for label, path in missing_traces:
    display(HTML(ru.missing_note(label, path)))

if trace_frames:
    long = pd.concat(trace_frames, ignore_index=True)
    long = long.dropna(subset=["param_value"])
    long = ru.to_estimation_scale(long, param_col="quantity", value_col="param_value")
    long["group_id"] = long["source"].astype(str) + "_" + long["rep"].astype(str)

    unit_col = "unit" if "unit" in long.columns else None
    groups = sorted(long[unit_col].dropna().unique()) if unit_col else [None]
    for unit in groups:
        data = long[long[unit_col] == unit] if unit_col else long
        if data.empty:
            continue
        data = data.copy()
        data["quantity"] = data["quantity"].astype(str)
        display(
            ggplot(data, aes(x="iter", y="param_value", group="group_id", color="source"))
            + geom_line(alpha=0.2, size=0.4)
            + facet_wrap("~quantity", scales="free_y")
            + labs(
                title=f"Block-IF2 iteration traces{f': {unit}' if unit else ''}",
                x="IF2 iteration",
                y="Value",
                color="Source",
            )
            + ru.scale_color_premium()
            + ru.theme_premium
        )

Expectation

Traces from independent starts should converge toward a common region rather than each settling at its own starting neighbourhood.