Visualization¶

History matching produces a lot to look at — emulator fits, the shrinking plausible region, how close the simulated outputs sit to their targets. This notebook tours the built-in plotting and display API on the stochastic SIR model used throughout these tutorials.

Every plot_* method returns Matplotlib axes and accepts an ax= argument, so the figures render inline and compose into your own layouts. The same figures are written to the engine's output_dir after each wave.

Each figure is explained inline as we go; the plotting API reference lists every plot_* function.

In [1]:

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

import historymatching as hm
import logging; logging.getLogger("historymatching").propagate = False  # keep INFO logs out of the rendered docs
from model import SIR, generate_observed_data

%matplotlib inline
np.random.seed(0)

import numpy as np import pandas as pd import matplotlib.pyplot as plt import historymatching as hm import logging; logging.getLogger("historymatching").propagate = False # keep INFO logs out of the rendered docs from model import SIR, generate_observed_data %matplotlib inline np.random.seed(0)

WARNING: All log messages before absl::InitializeLog() is called are written to STDERR
I0000 00:00:1782248409.050010    3254 cudart_stub.cc:31] Could not find cuda drivers on your machine, GPU will not be used.
I0000 00:00:1782248409.092585    3254 cpu_feature_guard.cc:227] This TensorFlow binary is optimized to use available CPU instructions in performance-critical operations.
To enable the following instructions: AVX2 FMA, in other operations, rebuild TensorFlow with the appropriate compiler flags.

WARNING: All log messages before absl::InitializeLog() is called are written to STDERR
I0000 00:00:1782248410.249398    3254 cudart_stub.cc:31] Could not find cuda drivers on your machine, GPU will not be used.

/home/runner/work/historymatching/historymatching/.venv/lib/python3.13/site-packages/gpflow/versions.py:1: UserWarning: pkg_resources is deprecated as an API. See https://setuptools.pypa.io/en/latest/pkg_resources.html. The pkg_resources package is slated for removal as early as 2025-11-30. Refrain from using this package or pin to Setuptools<81.
  import pkg_resources

Set up and run a calibration¶

We recover the transmission/recovery rates of a synthetic outbreak, then use the result to demonstrate each plot. (See Basic Workflow for the workflow itself in detail.)

In [2]:

beta_true, gamma_true = 1.3, 0.5
population_size, n_seed = 10_000, 100
incidence_obs, _ = generate_observed_data(beta_true=beta_true, gamma_true=gamma_true,
                                          population_size=population_size,
                                          n_seed_infections=n_seed)

def sir_sim(samples: pd.DataFrame) -> pd.DataFrame:
    df = samples.copy()
    if "rand_seed" not in df.columns:
        df["rand_seed"] = np.random.default_rng(0).integers(0, 2**31, size=len(df))
    rows = []
    for _, row in df.iterrows():
        model = SIR(beta=row["beta"], gamma=row["gamma"],
                    s0=population_size - n_seed, i0=n_seed, seed=int(row["rand_seed"]))
        inc = model.get_incidence()
        rec = {f"incidence_{i}": inc[i] for i in range(len(inc))}
        rec["peak_incidence"] = inc.max()
        rec["total_cases"] = inc.sum()
        rows.append(rec)
    return pd.DataFrame(rows)

observations = {
    "peak_incidence": (incidence_obs.max(), 50),
    "total_cases":    (incidence_obs.sum(), 200),
    "incidence_10":   (incidence_obs.iloc[10], 40),
}

engine = hm.HistoryMatching(
    bounds={"beta": (0.5, 3.0), "gamma": (0.1, 1.0)},
    observations=observations,
    function=sir_sim,
    emulator_type="bayes_linear",
    n_samples=300,
    max_iterations=3,
    random_seed=123,
    feature_selection=["peak_incidence", "incidence_10"],
    output_dir=None,  # focus on interactive plots; no files written here
)
results = engine.run()
print(f"Completed {len(results)} waves.")

beta_true, gamma_true = 1.3, 0.5 population_size, n_seed = 10_000, 100 incidence_obs, _ = generate_observed_data(beta_true=beta_true, gamma_true=gamma_true, population_size=population_size, n_seed_infections=n_seed) def sir_sim(samples: pd.DataFrame) -> pd.DataFrame: df = samples.copy() if "rand_seed" not in df.columns: df["rand_seed"] = np.random.default_rng(0).integers(0, 2**31, size=len(df)) rows = [] for _, row in df.iterrows(): model = SIR(beta=row["beta"], gamma=row["gamma"], s0=population_size - n_seed, i0=n_seed, seed=int(row["rand_seed"])) inc = model.get_incidence() rec = {f"incidence_{i}": inc[i] for i in range(len(inc))} rec["peak_incidence"] = inc.max() rec["total_cases"] = inc.sum() rows.append(rec) return pd.DataFrame(rows) observations = { "peak_incidence": (incidence_obs.max(), 50), "total_cases": (incidence_obs.sum(), 200), "incidence_10": (incidence_obs.iloc[10], 40), } engine = hm.HistoryMatching( bounds={"beta": (0.5, 3.0), "gamma": (0.1, 1.0)}, observations=observations, function=sir_sim, emulator_type="bayes_linear", n_samples=300, max_iterations=3, random_seed=123, feature_selection=["peak_incidence", "incidence_10"], output_dir=None, # focus on interactive plots; no files written here ) results = engine.run() print(f"Completed {len(results)} waves.")

Completed 3 waves.

Text summaries¶

engine.get_status_summary() folds the per-wave convergence and the surviving parameter ranges into one report. engine.nroy_summary() returns the same per-parameter information as a tidy DataFrame.

In [3]:

print(engine.get_status_summary())

print(engine.get_status_summary())

=== History Matching Status ===
State: completed
Progress: 3/3 waves
Acceptance rate: 8.7%
Total samples generated: 17,193
Total samples accepted: 1,800
Emulators trained: 6
[OK] Simulator function configured
[OK] All waves completed successfully

In [4]:

engine.nroy_summary()

engine.nroy_summary()

NROY summary
------------------------------------------------
  waves run:     3
  NROY samples:  300
  NROY fraction: 4.708%
  parameter               median   95% interval
  beta                     1.477   [1.115, 1.825]
  gamma                   0.6412   [0.3074, 0.9646]

Out[4]:

'NROY summary\n------------------------------------------------\n  waves run:     3\n  NROY samples:  300\n  NROY fraction: 4.708%\n  parameter               median   95% interval\n  beta                     1.477   [1.115, 1.825]\n  gamma                   0.6412   [0.3074, 0.9646]'

Convergence¶

The NROY fraction is the share of fresh prior samples passing all emulator constraints so far. A steadily falling fraction means the calibration is successfully ruling out parameter space.

In [5]:

engine.plot_convergence();

engine.plot_convergence();

The NROY parameter cloud¶

The headline result: a corner plot of the non-implausible region — marginals on the diagonal, pairwise scatter below. We overlay the known true values with truth=.

In [6]:

truth = {"beta": beta_true, "gamma": gamma_true}
engine.plot_nroy(truth=truth);

truth = {"beta": beta_true, "gamma": gamma_true} engine.plot_nroy(truth=truth);

plot_marginals shows just the one-dimensional posteriors, with the sample median and the true value marked.

In [7]:

engine.plot_marginals(truth=truth);

engine.plot_marginals(truth=truth);

Z-scores vs targets¶

For every target, the distribution of (simulated - target) / target_std across the NROY samples, by wave. Bands inside the green ±threshold region and centred on zero are consistent with the data.

In [8]:

engine.plot_zscores();

engine.plot_zscores();

Constrained directions¶

A PCA view of which combinations of parameters the data constrained most — useful for spotting non-identifiable directions that marginal plots miss.

In [9]:

engine.plot_constrained_dims(n_top=2);

engine.plot_constrained_dims(n_top=2);

Emulator quality¶

The emulators are only trustworthy if they reproduce the simulator. Inspect the fit per feature for the final wave.

In [10]:

final = results[-1]
final.quality_table()

final = results[-1] final.quality_table()

Out[10]:

	r2	mse	n_train
peak_incidence	0.831108	4168.843742	225.0
incidence_10	0.767963	1355.135836	225.0

In [11]:

final.plot_emulator_quality();

final.plot_emulator_quality();

In [12]:

final.plot_predicted_vs_actual("peak_incidence");

final.plot_predicted_vs_actual("peak_incidence");

Ensemble fan plot vs observed data¶

History matching constrains parameters; to check the outputs, re-run the simulator at NROY parameter sets and compare the trajectories to the observed data. hm.HistoryMatching.plot_ensemble_fan is model-agnostic — give it a 2-D array of trajectories and an observed series.

In [13]:

nroy = engine.get_nroy_samples(40, method="lhs")  # unbiased draw
trajectories = np.array([
    SIR(beta=r["beta"], gamma=r["gamma"], s0=population_size - n_seed, i0=n_seed).get_incidence()
    for _, r in nroy.iterrows()
])
hm.HistoryMatching.plot_ensemble_fan(trajectories, observed=incidence_obs.values,
                     xlabel="Day", ylabel="Daily incidence",
                     title="NROY trajectories vs observed");

nroy = engine.get_nroy_samples(40, method="lhs") # unbiased draw trajectories = np.array([ SIR(beta=r["beta"], gamma=r["gamma"], s0=population_size - n_seed, i0=n_seed).get_incidence() for _, r in nroy.iterrows() ]) hm.HistoryMatching.plot_ensemble_fan(trajectories, observed=incidence_obs.values, xlabel="Day", ylabel="Daily incidence", title="NROY trajectories vs observed");

Domain-object views¶

ParameterSpace, ObservationData, and EmulatorBank each have a summary() and (where it makes sense) a plot.

In [14]:

print(engine.parameter_space.summary())
print()
print(engine.observations.summary())
print()
print(engine.emulator_bank.summary())

print(engine.parameter_space.summary()) print() print(engine.observations.summary()) print() print(engine.emulator_bank.summary())

ParameterSpace: 2 parameter(s)
  beta                 [0.5, 3]
  gamma                [0.1, 1]

ObservationData: 3 target(s)
  peak_incidence       mean=1576, std=50
  total_cases          mean=9091, std=200
  incidence_10         mean=285, std=40

EmulatorBank: 6 emulator(s) across 3 wave(s)
  wave 1: peak_incidence, incidence_10
  wave 2: peak_incidence, incidence_10
  wave 3: peak_incidence, incidence_10

In [15]:

fig, axes = plt.subplots(1, 2, figsize=(12, 4))
engine.parameter_space.plot_bounds(reference=engine.parameter_space, ax=axes[0])
engine.observations.plot_targets(ax=axes[1])
plt.tight_layout();

fig, axes = plt.subplots(1, 2, figsize=(12, 4)) engine.parameter_space.plot_bounds(reference=engine.parameter_space, ax=axes[0]) engine.observations.plot_targets(ax=axes[1]) plt.tight_layout();

All of these methods return Matplotlib axes, so you can drop them into your own multi-panel figures (as above) or save them with fig.savefig(...). For the full catalogue of plotting functions, see the plotting API reference.