Measles Parameter Estimation

Distribution of IF2 parameter estimates: R pomp 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 IF2 parameter estimates in pypomp and R pomp. 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.

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

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 pomp": os.path.join("results", "R"),
    "pypomp (GPU)": os.path.join("results", "gpu"),
}
runs = ru.load_timing_data(PLATFORMS)

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

py_results = ru.read_if_exists(os.path.join(PY_DIR, "results.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.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)
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>"))

Automated checks

Show Code 
units = sorted(transformed["unit"].dropna().unique())
by_unit = {
    unit: [
        (param, g[g["source"] == "pypomp"]["value"], g[g["source"] == "R pomp"]["value"])
        for param, g in transformed[transformed["unit"] == unit].groupby("param")
    ]
    for unit in units
}
M_TESTS = 2 * sum(len(s) for s in by_unit.values())

rows, details = [], []
for unit, samples in by_unit.items():
    r, det = checks.group_checks(samples, "Final estimates", f"{unit}: pypomp vs R pomp", M_TESTS)
    rows += r
    details.append(det)
if not rows:
    rows = checks.missing_rows(["Final estimates"], "pypomp vs R pomp")

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)))
if details:
    display(HTML(checks.details_html(pd.concat(details, ignore_index=True))))

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 Halesworth: pypomp vs R pomp largest |z| = 1.19 (logit(amplitude)) of 13 PASS
Final estimates: distributions (KS) Halesworth: pypomp vs R pomp smallest p = 0.0665 (log(sigmaSE)) of 13 PASS
Final estimates: means London: pypomp vs R pomp largest |z| = 1.34 (log(gamma)) of 13 PASS
Final estimates: distributions (KS) London: pypomp vs R pomp smallest p = 0.226 (log(sigma)) of 13 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 52 tests that receive a verdict.

Results for each parameter
compared name difference SE z KS D KS p
Halesworth: pypomp vs R pomp log(R0) -0.05559 0.07226 -0.7693 0.08056 0.1934
Halesworth: pypomp vs R pomp log(gamma) 0.03859 0.07959 0.4849 0.05278 0.6985
Halesworth: pypomp vs R pomp log(iota) -0.04232 0.1003 -0.4221 0.05 0.7599
Halesworth: pypomp vs R pomp log(psi) 0.03209 0.0387 0.8291 0.04444 0.8698
Halesworth: pypomp vs R pomp log(sigma) -0.1351 0.147 -0.9195 0.05278 0.6985
Halesworth: pypomp vs R pomp log(sigmaSE) -0.1064 0.09564 -1.112 0.09722 0.06652
Halesworth: pypomp vs R pomp logit(E_0) -0.09595 0.1493 -0.6428 0.05278 0.6985
Halesworth: pypomp vs R pomp logit(I_0) 0.116 0.1083 1.072 0.05833 0.5733
Halesworth: pypomp vs R pomp logit(R_0) -0.001817 0.03287 -0.05528 0.06389 0.455
Halesworth: pypomp vs R pomp logit(S_0) 0.0003802 0.0335 0.01135 0.05833 0.5733
Halesworth: pypomp vs R pomp logit(amplitude) 0.1154 0.09666 1.194 0.07222 0.3051
Halesworth: pypomp vs R pomp logit(cohort) -0.1487 0.1301 -1.143 0.07222 0.3051
Halesworth: pypomp vs R pomp logit(rho) 0.01325 0.0536 0.2473 0.05278 0.6985
London: pypomp vs R pomp log(R0) 0.02821 0.03817 0.7392 0.075 0.2636
London: pypomp vs R pomp log(gamma) -0.06164 0.04597 -1.341 0.05 0.7599
London: pypomp vs R pomp log(iota) -0.2198 0.1848 -1.189 0.07222 0.3051
London: pypomp vs R pomp log(psi) -0.005474 0.004812 -1.138 0.06944 0.3509
London: pypomp vs R pomp log(sigma) -0.04746 0.03937 -1.205 0.07778 0.2265
London: pypomp vs R pomp log(sigmaSE) 0.03032 0.03817 0.7944 0.06944 0.3509
London: pypomp vs R pomp logit(E_0) 0.07211 0.09343 0.7718 0.06944 0.3509
London: pypomp vs R pomp logit(I_0) -0.008401 0.05268 -0.1595 0.05556 0.6357
London: pypomp vs R pomp logit(R_0) 0.0127 0.05201 0.2442 0.03889 0.9488
London: pypomp vs R pomp logit(S_0) -0.0128 0.05205 -0.246 0.03889 0.9488
London: pypomp vs R pomp logit(amplitude) 0.05732 0.05084 1.127 0.06389 0.455
London: pypomp vs R pomp logit(cohort) 0.1187 0.09203 1.29 0.04722 0.8177
London: pypomp vs R pomp logit(rho) -0.01733 0.01918 -0.9038 0.05833 0.5733

Visual review

Summary statistics

Show Code 
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
Halesworth R pomp log(R0) 360 4.4550 0.9914 0.4908 8.6587
Halesworth pypomp log(R0) 360 4.3994 0.9471 0.0801 10.7455
Halesworth R pomp log(gamma) 360 3.9725 1.0828 3.0847 9.1137
Halesworth pypomp log(gamma) 360 4.0111 1.0527 3.0578 9.5296
Halesworth R pomp log(iota) 360 -3.5470 1.3708 -6.7343 9.4817
Halesworth pypomp log(iota) 360 -3.5894 1.3190 -6.2156 7.6940
Halesworth R pomp log(psi) 360 -0.5213 0.5587 -3.4946 3.4984
Halesworth pypomp log(psi) 360 -0.4892 0.4764 -3.2414 3.0007
Halesworth R pomp log(sigma) 360 5.1906 2.0003 2.7814 13.2651
Halesworth pypomp log(sigma) 360 5.0554 1.9427 2.9482 14.0600
Halesworth R pomp log(sigmaSE) 360 -1.4320 1.2795 -10.0608 0.5000
Halesworth pypomp log(sigmaSE) 360 -1.5383 1.2868 -9.2556 0.5137
Halesworth R pomp logit(E_0) 360 -10.0764 2.0268 -15.2387 -6.3634
Halesworth pypomp logit(E_0) 360 -10.1724 1.9779 -15.1612 -5.9337
Halesworth R pomp logit(I_0) 360 -8.2763 1.5068 -15.1470 -5.2978
Halesworth pypomp logit(I_0) 360 -8.1603 1.3961 -17.0937 -4.2302
Halesworth R pomp logit(R_0) 360 3.3767 0.5259 -0.2417 6.2742
Halesworth pypomp logit(R_0) 360 3.3749 0.3352 2.0574 4.8643
Halesworth R pomp logit(S_0) 360 -3.4004 0.5331 -6.3162 0.2414
Halesworth pypomp logit(S_0) 360 -3.4000 0.3462 -4.9345 -2.0630
Halesworth R pomp logit(amplitude) 360 -0.6054 1.3450 -4.8007 3.4891
Halesworth pypomp logit(amplitude) 360 -0.4900 1.2469 -4.5583 3.5687
Halesworth R pomp logit(cohort) 360 -0.8602 1.7321 -5.3148 3.9532
Halesworth pypomp logit(cohort) 360 -1.0089 1.7583 -6.7777 4.8176
Halesworth R pomp logit(rho) 360 1.7254 0.6871 -1.9003 3.9078
Halesworth pypomp logit(rho) 360 1.7387 0.7497 -1.3168 5.5243
London R pomp log(R0) 360 3.5612 0.5191 0.0877 4.3414
London pypomp log(R0) 360 3.5894 0.5050 0.0830 4.2159
London R pomp log(gamma) 360 3.7164 0.6740 2.9109 9.1198
London pypomp log(gamma) 360 3.6548 0.5535 2.8276 8.2414
London R pomp log(iota) 360 0.3916 2.4757 -10.5882 4.7403
London pypomp log(iota) 360 0.1718 2.4838 -7.8134 4.4458
London R pomp log(psi) 360 -2.1469 0.0640 -2.3656 -1.9651
London pypomp log(psi) 360 -2.1524 0.0651 -2.3705 -1.9187
London R pomp log(sigma) 360 3.7334 0.5817 2.8659 10.0207
London pypomp log(sigma) 360 3.6860 0.4688 2.8585 6.9922
London R pomp log(sigmaSE) 360 -2.7533 0.4634 -4.7248 -1.7573
London pypomp log(sigmaSE) 360 -2.7229 0.5565 -9.1746 -1.8539
London R pomp logit(E_0) 360 -10.8300 1.2215 -16.9441 -9.1518
London pypomp logit(E_0) 360 -10.7579 1.2848 -16.9492 -8.7442
London R pomp logit(I_0) 360 -10.1164 0.6618 -13.5654 -8.9395
London pypomp logit(I_0) 360 -10.1248 0.7492 -14.6729 -9.1065
London R pomp logit(R_0) 360 3.2219 0.7039 -3.5770 3.8575
London pypomp logit(R_0) 360 3.2346 0.6916 -3.4592 3.7594
London R pomp logit(S_0) 360 -3.2243 0.7045 -3.8641 3.5765
London pypomp logit(S_0) 360 -3.2371 0.6920 -3.7640 3.4537
London R pomp logit(amplitude) 360 -0.2957 0.7067 -3.0831 1.7683
London pypomp logit(amplitude) 360 -0.2384 0.6566 -2.6862 1.7266
London R pomp logit(cohort) 360 0.3780 1.2452 -5.5017 3.4633
London pypomp logit(cohort) 360 0.4967 1.2241 -6.0856 4.8889
London R pomp logit(rho) 360 -0.0342 0.2948 -0.7942 4.5635
London pypomp logit(rho) 360 -0.0516 0.2132 -1.0644 3.2431

Parameter distributions

Show Code 
for unit in units:
    unit_data = transformed[transformed["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"IF2 parameter estimates: {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
    )

Convergence traces

Show Code 
py_traces = ru.read_if_exists(os.path.join(PY_DIR, "traces.csv.gz"))
r_traces = ru.read_if_exists(os.path.join(R_DIR, "mif_traces.csv.gz"))
if r_traces is not None:
    r_traces = ru.align_pomp_mif_loglik(r_traces, ["unit", "replicate"])

trace_frames = []

if py_traces is not None:
    py_mif = py_traces[py_traces["method"] == "mif"] if "method" in py_traces.columns else py_traces
    py_mif = py_mif.rename(
        columns={"theta_idx": "replicate", "rep": "replicate", "iter": "iteration", "loglik": "logLik"}
    )
    if "replicate" not in py_mif.columns:
        py_mif["replicate"] = 1
    id_vars = [c for c in ["unit", "replicate", "iteration"] if c in py_mif.columns]
    value_vars = [c for c in py_mif.columns if c in ru.LOG_PARAMS + ru.LOGIT_PARAMS + ["logLik"]]
    py_long = py_mif.melt(id_vars=id_vars, value_vars=value_vars, var_name="param", value_name="value")
    py_params = ru.to_estimation_scale(py_long[py_long["param"] != "logLik"])
    py_loglik = py_long[py_long["param"] == "logLik"].copy()
    py_combined = pd.concat([py_params, py_loglik], ignore_index=True)
    py_combined["source"] = "pypomp"
    trace_frames.append(py_combined)

if r_traces is not None:
    r_mif = r_traces.rename(
        columns={"theta_idx": "replicate", "rep": "replicate", "iter": "iteration", "loglik": "logLik"}
    )
    if "replicate" not in r_mif.columns:
        r_mif["replicate"] = 1
    id_vars = [c for c in ["unit", "replicate", "iteration"] if c in r_mif.columns]
    value_vars = [c for c in r_mif.columns if c in ru.LOG_PARAMS + ru.LOGIT_PARAMS + ["logLik"]]
    r_long = r_mif.melt(id_vars=id_vars, value_vars=value_vars, var_name="param", value_name="value")
    r_params = ru.to_estimation_scale(r_long[r_long["param"] != "logLik"])
    r_loglik = r_long[r_long["param"] == "logLik"].copy()
    r_combined = pd.concat([r_params, r_loglik], ignore_index=True)
    r_combined["source"] = "R pomp"
    trace_frames.append(r_combined)

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 pomp 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:
    combined_traces = pd.concat(trace_frames, ignore_index=True)
    combined_traces["group_id"] = combined_traces["source"] + "_" + combined_traces["replicate"].astype(str)

    for unit in sorted(combined_traces["unit"].dropna().unique()):
        unit_data = combined_traces[combined_traces["unit"] == unit]
        display(
            ggplot(unit_data, aes(x="iteration", y="value", group="group_id", color="source"))
            + geom_line(alpha=0.4, size=0.4)
            + facet_wrap("~param", scales="free_y")
            + labs(
                title=f"IF2 iteration traces: {unit}",
                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=False)))
Setting / Parameter R pomp pypomp (GPU)
Algorithmic & Workload Settings
Run Level 4 4
Starting Searches ($N_{starts}$) 360 360
IF2 Iterations ($N_{iter}$) 100 100
IF2 Particles ($N_p$) 5,000 5,000
Unit(s) London, Halesworth (2 units) London, Halesworth (2 units)
Random Seed — 594709947
Software & Environment
Pomp Framework pomp 6.3 pypomp 1.0.6
Backend / Engine R 4.4.0 JAX 0.11.2
Quant Git Commit 8a1411f 9c649ea
Run Timestamp 2026-08-12 21:18:25 2026-09-24 20:22:07
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 57201587 61826877