Panel Measles Parameter Estimation

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

Published

October 16, 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  |  Scaling

Introduction

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

The panel has 4 units (London, Halesworth, Hastings and Cardiff), and the parameters are mixed:

  • 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. We collapse them to a single shared entry so that each start counts once. Parameters are compared on the estimation scale (log for rates, logit for probabilities and initial-state fractions), which is the scale block-IF2 perturbs on.

The Python and R smoothing spline implementations interpolate the birth rate and population covariates differently, most noticeably for Cardiff. We therefore export the interpolated covariate grid from pomp into ../R_covariates.csv and use it in the pypomp objects, so that both packages search over the same model.

The pypomp run continues past the search into a panel-specific step: mix-and-match across units, then pruning to the single best combination. panelPomp has no equivalent, so the comparison with R uses results_final.csv, which is saved before that step. The mix-and-match result is reported separately.

Show Code 
import numpy as np
import pandas as pd
import checks
from plotnine import ggplot, aes, geom_density, geom_line, facet_wrap, labs

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

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

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"))

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)

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:
    ll_frames.append(py_ll.assign(source="pypomp"))
if r_ll is not None:
    ll_frames.append(r_ll.assign(source="R panelPomp"))
per_unit = (
    pd.concat(ll_frames, ignore_index=True)
    .groupby(["source", "unit", "theta_idx"])["logLik"]
    .agg(logmeanexp=ru.logmeanexp)
    .reset_index()
    if ll_frames
    else pd.DataFrame(columns=["source", "unit", "theta_idx", "logmeanexp"])
)
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>"))
for label, df, path in (
    ("pypomp endpoint log-likelihoods", py_ll, os.path.join(PY_DIR, "pfilter_logliks.csv")),
    ("R panelPomp endpoint log-likelihoods", r_ll, os.path.join(R_DIR, "pfilter_logliks.csv")),
):
    if df is None:
        display(HTML(ru.missing_note(label, path)))

Automated checks

Show Code 
COMPARED = "pypomp vs R panelPomp"
blocks = ["shared"] + sorted(u for u in transformed["unit"].dropna().unique() if u != "shared")
estimates = {
    block: [
        (param, g[g["source"] == "pypomp"]["value"], g[g["source"] == "R panelPomp"]["value"])
        for param, g in transformed[transformed["unit"] == block].groupby("param")
    ]
    for block in blocks
}
endpoint = [
    (unit, g[g["source"] == "pypomp"]["logmeanexp"], g[g["source"] == "R panelPomp"]["logmeanexp"])
    for unit, g in per_unit.groupby("unit")
]
M_TESTS = 2 * (sum(len(s) for s in estimates.values()) + len(endpoint))

rows, details = [], []
for block, samples in estimates.items():
    r, det = checks.group_checks(samples, "Final estimates", f"{block}: {COMPARED}", M_TESTS)
    rows += r
    details.append(det)
r, det = checks.group_checks(endpoint, "Endpoint log-likelihood", COMPARED, M_TESTS)
rows += r
details.append(det)

level = checks.run_level(runs)
rows = checks.gate_smoke(rows, level)
checks.save(rows, level, pypomp_version=runs["pypomp (GPU)"]["meta"].get("pypomp_version"))
display(HTML(checks.table_html(rows, M_TESTS)))
display(HTML(checks.details_html(pd.concat(details, ignore_index=True), "Results for each parameter and unit")))

These checks are computed from the results each time the report is rendered. They only compare summary statistics, so the plots can show problems that the checks miss. Visual review lists what the plots should show.

Check Compared Result Verdict
Final estimates: means shared: pypomp vs R panelPomp largest |z| = 1.66 (logit(cohort)) of 6 PASS
Final estimates: distributions (KS) shared: pypomp vs R panelPomp smallest p = 0.193 (logit(amplitude)) of 6 PASS
Final estimates: means Cardiff: pypomp vs R panelPomp largest |z| = 0.76 (log(iota)) of 7 PASS
Final estimates: distributions (KS) Cardiff: pypomp vs R panelPomp smallest p = 0.193 (log(iota)) of 7 PASS
Final estimates: means Halesworth: pypomp vs R panelPomp largest |z| = 2.02 (log(psi)) of 7 PASS
Final estimates: distributions (KS) Halesworth: pypomp vs R panelPomp smallest p = 0.0187 (log(psi)) of 7 PASS
Final estimates: means Hastings: pypomp vs R panelPomp largest |z| = 1.61 (logit(I_0)) of 7 PASS
Final estimates: distributions (KS) Hastings: pypomp vs R panelPomp smallest p = 0.0806 (logit(E_0)) of 7 PASS
Final estimates: means London: pypomp vs R panelPomp largest |z| = 2.11 (logit(rho)) of 7 PASS
Final estimates: distributions (KS) London: pypomp vs R panelPomp smallest p = 0.0546 (log(iota)) of 7 PASS
Endpoint log-likelihood: means pypomp vs R panelPomp largest |z| = 3.41 (Cardiff) of 4 PASS
Endpoint log-likelihood: distributions (KS) pypomp vs R panelPomp smallest p = 0.00429 (Cardiff) of 4 PASS

PASS when |z| < 3, CHECK when 3 ≤ |z| < 5, and FAIL when |z| ≥ 5. KS p-values use the equivalent thresholds. Both are Bonferroni-adjusted over the 76 tests that receive a verdict.

Results for each parameter and unit
compared name difference SE z KS D KS p
shared: pypomp vs R panelPomp log(R0) -0.01201 0.03353 -0.3581 0.06944 0.3509
shared: pypomp vs R panelPomp log(gamma) 0.00773 0.06688 0.1156 0.05833 0.5733
shared: pypomp vs R panelPomp log(sigma) 0.05516 0.05357 1.03 0.07778 0.2265
shared: pypomp vs R panelPomp log(sigmaSE) -0.03506 0.05212 -0.6728 0.05833 0.5733
shared: pypomp vs R panelPomp logit(amplitude) -0.02122 0.07144 -0.2971 0.08056 0.1934
shared: pypomp vs R panelPomp logit(cohort) -0.3112 0.1877 -1.658 0.05833 0.5733
Cardiff: pypomp vs R panelPomp log(iota) -0.04198 0.05554 -0.7558 0.08056 0.1934
Cardiff: pypomp vs R panelPomp log(psi) -0.00424 0.005953 -0.7123 0.04444 0.8698
Cardiff: pypomp vs R panelPomp logit(E_0) -0.01503 0.08708 -0.1726 0.05556 0.6357
Cardiff: pypomp vs R panelPomp logit(I_0) 0.03045 0.06525 0.4667 0.05 0.7599
Cardiff: pypomp vs R panelPomp logit(R_0) -0.006833 0.03209 -0.2129 0.08056 0.1934
Cardiff: pypomp vs R panelPomp logit(S_0) 0.006848 0.03212 0.2132 0.08056 0.1934
Cardiff: pypomp vs R panelPomp logit(rho) -0.002448 0.00777 -0.3151 0.06111 0.5127
Halesworth: pypomp vs R panelPomp log(iota) -0.02134 0.04103 -0.5202 0.06667 0.4009
Halesworth: pypomp vs R panelPomp log(psi) -0.0581 0.02873 -2.022 0.1139 0.01869
Halesworth: pypomp vs R panelPomp logit(E_0) -0.04574 0.1342 -0.3409 0.05556 0.6357
Halesworth: pypomp vs R panelPomp logit(I_0) -0.05787 0.07964 -0.7267 0.05 0.7599
Halesworth: pypomp vs R panelPomp logit(R_0) 0.001544 0.02846 0.05426 0.06389 0.455
Halesworth: pypomp vs R panelPomp logit(S_0) -0.001101 0.02894 -0.03804 0.06667 0.4009
Halesworth: pypomp vs R panelPomp logit(rho) 0.02311 0.04111 0.562 0.05 0.7599
Hastings: pypomp vs R panelPomp log(iota) -0.02765 0.05035 -0.5492 0.05556 0.6357
Hastings: pypomp vs R panelPomp log(psi) -0.000732 0.007685 -0.09525 0.03611 0.9733
Hastings: pypomp vs R panelPomp logit(E_0) 0.09743 0.0837 1.164 0.09444 0.08059
Hastings: pypomp vs R panelPomp logit(I_0) -0.1178 0.07321 -1.61 0.075 0.2636
Hastings: pypomp vs R panelPomp logit(R_0) -0.02548 0.04065 -0.6269 0.08889 0.1163
Hastings: pypomp vs R panelPomp logit(S_0) 0.02555 0.04067 0.6281 0.08889 0.1163
Hastings: pypomp vs R panelPomp logit(rho) -0.01675 0.02546 -0.658 0.05833 0.5733
London: pypomp vs R panelPomp log(iota) -0.2783 0.1876 -1.484 0.1 0.05459
London: pypomp vs R panelPomp log(psi) -0.008896 0.005438 -1.636 0.06111 0.5127
London: pypomp vs R panelPomp logit(E_0) -0.1101 0.09885 -1.114 0.08611 0.1386
London: pypomp vs R panelPomp logit(I_0) -0.02744 0.059 -0.4651 0.03889 0.9488
London: pypomp vs R panelPomp logit(R_0) -0.0001483 0.02821 -0.005258 0.07222 0.3051
London: pypomp vs R panelPomp logit(S_0) 0.0001799 0.02829 0.00636 0.07222 0.3051
London: pypomp vs R panelPomp logit(rho) 0.0103 0.004883 2.109 0.09167 0.0971
pypomp vs R panelPomp Cardiff 1.62 0.4755 3.407 0.1306 0.004288
pypomp vs R panelPomp Halesworth -0.008586 0.4111 -0.02088 0.06667 0.4009
pypomp vs R panelPomp Hastings 1.207 0.6491 1.86 0.08056 0.1934
pypomp vs R panelPomp London 2.197 1.597 1.375 0.06111 0.5127

Visual review

Summary statistics

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.4265 0.6777 -6.5192 0.2683
Cardiff R panelPomp log(psi) 360 -1.3274 0.0797 -1.5598 -1.1275
Cardiff pypomp log(psi) 360 -1.3316 0.0800 -1.6228 -1.1433
Cardiff R panelPomp logit(E_0) 360 -12.2716 1.2066 -17.7719 -10.1121
Cardiff pypomp logit(E_0) 360 -12.2867 1.1288 -16.2726 -9.9151
Cardiff R panelPomp logit(I_0) 360 -11.9284 0.9089 -15.9739 -10.1714
Cardiff pypomp logit(I_0) 360 -11.8979 0.8405 -14.4467 -10.0712
Cardiff R panelPomp logit(R_0) 360 3.1835 0.4680 -1.3559 4.0635
Cardiff pypomp logit(R_0) 360 3.1767 0.3897 2.3354 4.1338
Cardiff R panelPomp logit(S_0) 360 -3.1841 0.4683 -4.0656 1.3558
Cardiff pypomp logit(S_0) 360 -3.1772 0.3900 -4.1357 -2.3355
Cardiff R panelPomp logit(rho) 360 0.4290 0.1165 0.1532 1.3888
Cardiff pypomp logit(rho) 360 0.4265 0.0904 0.0717 0.6798
Halesworth R panelPomp log(iota) 360 -3.5921 0.5524 -5.2829 -2.3170
Halesworth pypomp log(iota) 360 -3.6134 0.5486 -5.7030 -2.3758
Halesworth R panelPomp log(psi) 360 -0.2210 0.3414 -3.2855 0.5357
Halesworth pypomp log(psi) 360 -0.2791 0.4250 -3.8791 0.5661
Halesworth R panelPomp logit(E_0) 360 -9.9699 1.8703 -16.2223 -6.2764
Halesworth pypomp logit(E_0) 360 -10.0156 1.7279 -14.3299 -6.5564
Halesworth R panelPomp logit(I_0) 360 -8.0779 1.0023 -15.8802 -6.0016
Halesworth pypomp logit(I_0) 360 -8.1358 1.1308 -13.7680 -6.4179
Halesworth R panelPomp logit(R_0) 360 2.9819 0.4057 -0.6082 3.8129
Halesworth pypomp logit(R_0) 360 2.9835 0.3562 2.1444 4.0375
Halesworth R panelPomp logit(S_0) 360 -2.9964 0.4116 -3.8271 0.5992
Halesworth pypomp logit(S_0) 360 -2.9975 0.3634 -4.1365 -2.1592
Halesworth R panelPomp logit(rho) 360 1.1668 0.5406 0.1084 4.3023
Halesworth pypomp logit(rho) 360 1.1899 0.5624 0.0510 5.1775
Hastings R panelPomp log(iota) 360 -1.8849 0.7394 -3.6612 5.6518
Hastings pypomp log(iota) 360 -1.9126 0.6050 -3.9082 0.0440
Hastings R panelPomp log(psi) 360 -0.7913 0.1037 -1.0533 -0.4748
Hastings pypomp log(psi) 360 -0.7920 0.1025 -1.0776 -0.3846
Hastings R panelPomp logit(E_0) 360 -12.7748 1.1655 -18.6498 -9.8697
Hastings pypomp logit(E_0) 360 -12.6774 1.0789 -15.9870 -9.5296
Hastings R panelPomp logit(I_0) 360 -12.2414 0.9971 -15.7087 -9.9430
Hastings pypomp logit(I_0) 360 -12.3592 0.9670 -16.0607 -10.3074
Hastings R panelPomp logit(R_0) 360 3.5908 0.5538 0.9506 5.9753
Hastings pypomp logit(R_0) 360 3.5653 0.5367 2.4700 5.0146
Hastings R panelPomp logit(S_0) 360 -3.5914 0.5542 -5.9783 -0.9506
Hastings pypomp logit(S_0) 360 -3.5658 0.5370 -5.0156 -2.4701
Hastings R panelPomp logit(rho) 360 0.9042 0.3588 0.3281 4.4942
Hastings pypomp logit(rho) 360 0.8875 0.3235 0.0840 2.6131
London R panelPomp log(iota) 360 -0.3036 2.5104 -12.2456 5.5178
London pypomp log(iota) 360 -0.5819 2.5227 -12.2640 3.7206
London R panelPomp log(psi) 360 -2.1882 0.0690 -2.4044 -2.0400
London pypomp log(psi) 360 -2.1971 0.0767 -2.5022 -1.9646
London R panelPomp logit(E_0) 360 -11.0339 1.3433 -16.7741 -9.0569
London pypomp logit(E_0) 360 -11.1440 1.3088 -15.8906 -9.2541
London R panelPomp logit(I_0) 360 -10.7906 0.7941 -13.7053 -8.6269
London pypomp logit(I_0) 360 -10.8180 0.7890 -13.8664 -8.9608
London R panelPomp logit(R_0) 360 3.1025 0.4073 -0.5147 3.8509
London pypomp logit(R_0) 360 3.1024 0.3473 2.2794 3.9730
London R panelPomp logit(S_0) 360 -3.1041 0.4082 -3.8558 0.5137
London pypomp logit(S_0) 360 -3.1039 0.3485 -3.9810 -2.2796
London R panelPomp logit(rho) 360 -0.0349 0.0821 -1.1646 0.1578
London pypomp logit(rho) 360 -0.0246 0.0429 -0.1404 0.1038
shared R panelPomp log(R0) 360 3.3543 0.4566 0.8974 4.3441
shared pypomp log(R0) 360 3.3423 0.4431 2.4717 4.6947
shared R panelPomp log(gamma) 360 4.4236 0.9585 2.9219 12.3555
shared pypomp log(gamma) 360 4.4313 0.8314 2.8912 9.0030
shared R panelPomp log(sigma) 360 4.1788 0.6636 3.0463 6.8606
shared pypomp log(sigma) 360 4.2340 0.7699 2.9746 8.4944
shared R panelPomp log(sigmaSE) 360 -2.9606 0.6861 -4.6905 -1.0842
shared pypomp log(sigmaSE) 360 -2.9957 0.7122 -4.5003 -1.3784
shared R panelPomp logit(amplitude) 360 -1.0288 0.9994 -9.7369 4.6968
shared pypomp logit(amplitude) 360 -1.0500 0.9156 -2.6153 4.2697
shared R panelPomp logit(cohort) 360 -0.2529 2.5117 -9.0060 8.3077
shared pypomp logit(cohort) 360 -0.5642 2.5258 -8.2851 7.9627

Shared parameter distributions

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
        )

Endpoint log-likelihoods

We re-evaluate each start’s final parameter vector with the particle filter and take the logmeanexp over replicate filters, giving one value per start and unit.

Show Code 
if not per_unit.empty:
    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 parameter vector. pypomp can also recombine the best unit-specific blocks across starts and refit, which has no counterpart in panelPomp. 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 -10535.7796 56.5519 London -5211.6519 51.8824 15.1328 178.2723 111.8577 0.0162 0.513 0.1223 1.8194 0.5049 0.1152 0.0510 0.0000 0.0000 0.9489
0 -10535.7796 56.5519 Halesworth -329.9817 0.1723 15.1328 178.2723 111.8577 0.0162 0.513 0.1223 0.0205 0.7741 0.7916 0.0461 0.0004 0.0001 0.9534
0 -10535.7796 56.5519 Hastings -1797.5880 9.5101 15.1328 178.2723 111.8577 0.0162 0.513 0.1223 0.1462 0.6609 0.4336 0.0230 0.0000 0.0000 0.9770
0 -10535.7796 56.5519 Cardiff -3196.5581 20.3926 15.1328 178.2723 111.8577 0.0162 0.513 0.1223 0.4579 0.5992 0.2480 0.0352 0.0000 0.0000 0.9648

Convergence traces

Show Code 
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
        )

Appendix: settings and environment

Show Code 
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.6
Backend / Engine R 4.4.0 JAX 0.11.2
Quant Git Commit e6ac592 9c649ea
Run Timestamp 2026-08-14 17:06:50 2026-09-24 22:24:34
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 61839612