Panel Measles Fixed-Parameter Likelihood Validation

Particle filter log-likelihood distribution: 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('loglik')}</div>"))
Likelihood Evaluation  |  Parameter Estimation  |  Timing & Throughput

Introduction

This report evaluates the panel particle filter’s likelihood calculation itself, by running it many times at the He et al. (2010) maximum likelihood estimates. Nothing varies across replicates except the random seed, so any systematic difference between pypomp and panelPomp is a difference between the implementations.

The panel has 4 units — London, Halesworth, Hastings and Cardiff — with a mixed parameter structure: R0, sigma, gamma, sigmaSE, cohort and amplitude are shared across units, while iota, rho, psi, S_0, E_0, I_0 and R_0 are estimated separately per unit. The panel log-likelihood is the sum over units, so the per-unit likelihoods below are the components that have to agree.

The python and R smoothing spline implementations interpolate the birthrate and population covariates differently — noticeably so for Cardiff — which on its own shifts the log-likelihood. Since the point here is to compare the two filter implementations rather than two slightly different models, the final interpolated covariate grid is exported from pomp into ../R_covariates.csv and injected into the pypomp objects (model.align_covariates). Both packages therefore filter the identical model.

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
Pfilter Particles / Unit ($N_{p,eval}$) 5,000 5,000
Evaluation Replicates ($N_{reps}$) 3,600 3,600
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 b22249b 630cc23
Run Timestamp 2026-08-12 13:00:06 2026-09-02 19:56:04
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 57205142 59687803
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 scipy.stats import ks_2samp
from plotnine import ggplot, aes, geom_density, facet_wrap, labs

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

frames = []
for label, d in PLATFORMS.items():
    path = os.path.join(d, "pfilter_logliks.csv")
    df = ru.read_if_exists(path)
    if df is not None:
        df = df[["unit", "replicate", "logLik"]].copy()
        df["source"] = label
        frames.append(df)

combined = pd.concat(frames, ignore_index=True) if frames else pd.DataFrame(
    columns=["unit", "replicate", "logLik", "source"]
)

Summary Statistics

logmeanexp is the log of the mean likelihood and sits above the mean log-likelihood by roughly sd^2/2. It is the quantity panelPomp reports as the likelihood estimate, so it is listed alongside the raw mean.

Show Code 
if combined.empty:
    display(HTML("<p><em>No log-likelihood data available.</em></p>"))
else:
    summary = (
        combined.groupby(["unit", "source"])["logLik"]
        .agg(n="count", mean="mean", sd="std", min="min", max="max", logmeanexp=ru.logmeanexp)
        .reset_index()
        .sort_values(["unit", "source"])
    )
    display(HTML(summary.round(3).to_html(classes="table table-striped table-hover", index=False)))
unit source n mean sd min max logmeanexp
Cardiff R panelPomp 3600 -2381.299 2.219 -2391.755 -2375.465 -2379.522
Cardiff pypomp (GPU) 3600 -2381.251 2.225 -2390.548 -2375.325 -2379.449
Halesworth R panelPomp 3600 -323.290 1.342 -328.727 -319.271 -322.526
Halesworth pypomp (GPU) 3600 -323.290 1.327 -328.748 -319.421 -322.517
Hastings R panelPomp 3600 -1592.498 1.855 -1599.050 -1586.051 -1590.693
Hastings pypomp (GPU) 3600 -1592.446 1.830 -1599.408 -1585.153 -1590.572
London R panelPomp 3600 -3953.893 2.950 -3966.563 -3944.690 -3950.535
London pypomp (GPU) 3600 -3954.011 2.938 -3968.128 -3943.809 -3950.503

Panel log-likelihood

The panel likelihood is the sum of the per-unit likelihoods. Summing the logmeanexp of each unit gives the estimate each package would report for the whole panel.

Show Code 
if combined.empty:
    display(HTML("<p><em>No log-likelihood data available.</em></p>"))
else:
    panel = (
        combined.groupby(["unit", "source"])["logLik"]
        .agg(logmeanexp=ru.logmeanexp)
        .reset_index()
        .groupby("source")["logmeanexp"]
        .sum()
        .reset_index()
        .rename(columns={"logmeanexp": "panel logLik"})
    )
    display(HTML(panel.round(3).to_html(classes="table table-striped table-hover", index=False)))
source panel logLik
R panelPomp -8243.275
pypomp (GPU) -8243.041

Kolmogorov-Smirnov tests against the R baseline

Show Code 
rows = []
if not combined.empty and "R panelPomp" in set(combined["source"]):
    for unit in sorted(combined["unit"].unique()):
        r_ll = combined[(combined["unit"] == unit) & (combined["source"] == "R panelPomp")]["logLik"]
        py_ll = combined[(combined["unit"] == unit) & (combined["source"].str.startswith("pypomp"))]["logLik"]
        if len(py_ll) == 0 or len(r_ll) == 0:
            continue
        d, p = ks_2samp(py_ll, r_ll)
        rows.append({
            "Unit": unit,
            "Comparison": "pypomp vs R panelPomp",
            "Mean difference (nats)": py_ll.mean() - r_ll.mean(),
            "KS D": d,
            "KS p": p,
        })

ks_table = pd.DataFrame(rows)
if ks_table.empty:
    display(HTML("<p><em>Both sources are needed for a KS comparison.</em></p>"))
else:
    display(HTML(ks_table.round(4).to_html(classes="table table-striped table-hover", index=False)))
Unit Comparison Mean difference (nats) KS D KS p
Cardiff pypomp vs R panelPomp 0.0483 0.0167 0.6994
Halesworth pypomp vs R panelPomp 0.0000 0.0161 0.7385
Hastings pypomp vs R panelPomp 0.0517 0.0253 0.2003
London pypomp vs R panelPomp -0.1178 0.0247 0.2213

Density Comparison

Show Code 
if combined.empty:
    display(HTML("<p><em>No log-likelihood data available.</em></p>"))
else:
    display(
        ggplot(combined, aes(x="logLik", fill="source", color="source"))
        + geom_density(alpha=0.35)
        + facet_wrap("~unit", scales="free")
        + labs(
            title="Panel particle filter log-likelihood at the He et al. (2010) estimates",
            subtitle="One density per implementation, per unit",
            x="Log-likelihood",
            y="Density",
            fill="Source",
            color="Source",
        )
        + ru.scale_fill_premium()
        + ru.scale_color_premium()
        + ru.theme_premium
    )

Expectation

The two densities should overlap heavily within each unit. A shift confined to one unit points at that unit’s model or covariates; a shift present in all four points at the filter.

Evaluation Cost

Wall clock for the full block of likelihood evaluations. The two arms run on different hardware, so this is context for the distributions above rather than the timing benchmark — see the timing report for a controlled comparison.

Show Code 
rows = []
for label, d in PLATFORMS.items():
    tm = ru.normalize_timings(ru.read_if_exists(os.path.join(d, "timings.csv")))
    if tm is None:
        continue
    cfg = ru.load_json(os.path.join(d, "latest.json")).get("run_config", {}) or {}
    for _, r in tm.iterrows():
        rows.append({
            "Source": label,
            "Phase": r["phase"],
            "Time (s)": r["time_seconds"],
            "Particles": cfg.get("NP_EVAL", "?"),
            "Replicates": cfg.get("NREPS_EVAL", "?"),
        })

timing_df = pd.DataFrame(rows)
if timing_df.empty:
    display(HTML("<p><em>No timings recorded.</em></p>"))
else:
    display(HTML(timing_df.round(2).to_html(classes="table table-striped table-hover", index=False)))
Source Phase Time (s) Particles Replicates
R panelPomp pfilter 7042.04 5000 3600
pypomp (GPU) pfilter 201.06 5000 3600