Measles Parameter Estimation

Distribution of IF2 parameter estimates: R pomp vs. pypomp

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  |  IF2 vs IFAD

Introduction

This report compares IF2 parameter estimates. 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 parameter estimates should agree. Traces should also look similar.

---

Benchmark Settings & Environment

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

PLATFORMS = {
    "R pomp": os.path.join("results", "R"),
    "pypomp (GPU)": os.path.join("results", "gpu"),
}

runs = ru.load_timing_data(PLATFORMS)
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.0rc1
Backend / Engine	R 4.4.0	JAX 0.11.1
Quant Git Commit	8a1411f	630cc23
Run Timestamp	2026-08-12 21:18:25	2026-09-02 19:41:37
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	59687786

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

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 = 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()
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 IF2 perturbs on.

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.5412	1.1253	0.0409	11.7302
Halesworth	R pomp	log(gamma)	360	3.9725	1.0828	3.0847	9.1137
Halesworth	pypomp	log(gamma)	360	3.9731	1.1110	2.9834	8.6136
Halesworth	R pomp	log(iota)	360	-3.5470	1.3708	-6.7343	9.4817
Halesworth	pypomp	log(iota)	360	-3.5834	1.4211	-6.4481	7.9643
Halesworth	R pomp	log(psi)	360	-0.5213	0.5587	-3.4946	3.4984
Halesworth	pypomp	log(psi)	360	-0.4916	0.5195	-3.3470	3.6905
Halesworth	R pomp	log(sigma)	360	5.1906	2.0003	2.7814	13.2651
Halesworth	pypomp	log(sigma)	360	5.0485	1.8736	3.0259	11.6979
Halesworth	R pomp	log(sigmaSE)	360	-1.4320	1.2795	-10.0608	0.5000
Halesworth	pypomp	log(sigmaSE)	360	-1.3979	1.2592	-8.8433	0.4695
Halesworth	R pomp	logit(E_0)	360	-10.0764	2.0268	-15.2387	-6.3634
Halesworth	pypomp	logit(E_0)	360	-10.0758	2.0404	-15.8079	-6.0489
Halesworth	R pomp	logit(I_0)	360	-8.2763	1.5068	-15.1470	-5.2978
Halesworth	pypomp	logit(I_0)	360	-8.1963	1.5039	-14.0992	-5.1571
Halesworth	R pomp	logit(R_0)	360	3.3767	0.5259	-0.2417	6.2742
Halesworth	pypomp	logit(R_0)	360	3.4356	0.4547	1.3961	8.4759
Halesworth	R pomp	logit(S_0)	360	-3.4004	0.5331	-6.3162	0.2414
Halesworth	pypomp	logit(S_0)	360	-3.4628	0.4650	-8.5163	-1.4086
Halesworth	R pomp	logit(amplitude)	360	-0.6054	1.3450	-4.8007	3.4891
Halesworth	pypomp	logit(amplitude)	360	-0.5187	1.4165	-6.2462	4.0010
Halesworth	R pomp	logit(cohort)	360	-0.8602	1.7321	-5.3148	3.9532
Halesworth	pypomp	logit(cohort)	360	-0.9107	1.5936	-5.8092	4.1027
Halesworth	R pomp	logit(rho)	360	1.7254	0.6871	-1.9003	3.9078
Halesworth	pypomp	logit(rho)	360	1.7720	0.7383	-2.1261	4.5802
London	R pomp	log(R0)	360	3.5612	0.5191	0.0877	4.3414
London	pypomp	log(R0)	360	3.5719	0.5043	0.2819	4.2245
London	R pomp	log(gamma)	360	3.7164	0.6740	2.9109	9.1198
London	pypomp	log(gamma)	360	3.7037	0.5301	2.8509	7.4577
London	R pomp	log(iota)	360	0.3916	2.4757	-10.5882	4.7403
London	pypomp	log(iota)	360	0.0083	2.6417	-8.3780	4.1399
London	R pomp	log(psi)	360	-2.1469	0.0640	-2.3656	-1.9651
London	pypomp	log(psi)	360	-2.1409	0.0654	-2.3408	-1.8832
London	R pomp	log(sigma)	360	3.7334	0.5817	2.8659	10.0207
London	pypomp	log(sigma)	360	3.6726	0.4518	2.9198	6.0847
London	R pomp	log(sigmaSE)	360	-2.7533	0.4634	-4.7248	-1.7573
London	pypomp	log(sigmaSE)	360	-2.7754	0.6818	-9.7712	-1.8131
London	R pomp	logit(E_0)	360	-10.8300	1.2215	-16.9441	-9.1518
London	pypomp	logit(E_0)	360	-10.7279	1.2458	-16.5107	-9.0007
London	R pomp	logit(I_0)	360	-10.1164	0.6618	-13.5654	-8.9395
London	pypomp	logit(I_0)	360	-10.1365	0.6210	-13.2758	-9.0149
London	R pomp	logit(R_0)	360	3.2219	0.7039	-3.5770	3.8575
London	pypomp	logit(R_0)	360	3.2261	0.7152	-2.8416	3.8892
London	R pomp	logit(S_0)	360	-3.2243	0.7045	-3.8641	3.5765
London	pypomp	logit(S_0)	360	-3.2286	0.7158	-3.8945	2.8414
London	R pomp	logit(amplitude)	360	-0.2957	0.7067	-3.0831	1.7683
London	pypomp	logit(amplitude)	360	-0.2495	0.7606	-3.2549	1.6436
London	R pomp	logit(cohort)	360	0.3780	1.2452	-5.5017	3.4633
London	pypomp	logit(cohort)	360	0.4012	1.2598	-4.2926	4.1309
London	R pomp	logit(rho)	360	-0.0342	0.2948	-0.7942	4.5635
London	pypomp	logit(rho)	360	-0.0380	0.1894	-0.7498	2.9779

---

Parameter Distributions

for unit in sorted(transformed["unit"].dropna().unique()):
    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
    )

Expectation

The pypomp and pomp densities should overlap heavily for each parameter.

---

Convergence Traces

from plotnine import geom_line

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

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
        )

Expectation

The pypomp and pomp traces should overlap heavily for each parameter.