Repeated calibration

Fit independent instrument datasets with one compiled regression model. Each instrument can have a different number of readings. Stable IDs keep each fit's random stream tied to its instrument, so reordering the batch does not change any instrument's result.

This does not pool instruments: every fit is independent and learns nothing from its neighbours. Use examples/site_effects.py when group estimates should share information. Chunked batch dispatch currently retains inputs and fits; see the roadmap for bounded streaming.

Run it with python examples/repeated_calibration.py from the repository root. The blocks below are cells of that one file and share its state, so they assume the blocks above them have run.

import numpy as np

from instrument_calibration import calibration_model

Three instruments with different sample counts

rng = np.random.default_rng(42)
datasets = []
ids = []
for i, rows in enumerate([40, 75, 120]):
    x = np.linspace(-2, 2, rows)
    y = 0.1 * i + (1.0 + 0.05 * i) * x + rng.normal(0, 0.2, rows)
    ids.append(f"instrument-{i}")
    datasets.append({"x": x, "y": y})

for name, dataset in zip(ids, datasets):
    print(f"{name}: {len(dataset['x'])} readings")

instrument-0: 40 readings
instrument-1: 75 readings
instrument-2: 120 readings

One compiled model, many datasets

calibration_model() comes from examples/instrument_calibration.py: an offset, a gain and a noise scale, fitted to y ~ Normal(offset + gain * x, noise). Compiling once and calling sample_batch builds the graph structure once and reuses it for every dataset. Each fit still validates its own binding and lays out its own evaluator buffers, because those depend on the data's shapes: what is saved is the model construction, not the per-fit setup.

errors="collect" keeps one bad dataset from losing the whole batch: the failing ID lands in batch.errors and the rest still return fits.

Fit the batch

batch = calibration_model().sample_batch(
    datasets,
    ids=ids,
    chains=4,
    warmup=1000,
    draws=1000,
    threads=2,
    seed=42,
    errors="collect",
    show_progress=False,
)

for name in batch.ids:
    print(f"===== {name}")
    if name in batch.errors:
        print(batch.errors[name])
    else:
        print(batch.get(name).summary())

===== instrument-0
4 chains × 1000 draws per chain

Parameter        mean      std     hdi_3%    hdi_97%   ess_bulk   ess_tail    r_hat  mcse_mean
──────────────────────────────────────────────────────────────────────────────────────────────
offset         0.0076   0.0272    -0.0438     0.0583       3438       2690   1.0009   0.000468
gain           1.0161   0.0232     0.9726     1.0613       3943       2641   1.0024   0.000371
noise          0.1720   0.0213     0.1324     0.2108       3238       2606   1.0012   0.000380
──────────────────────────────────────────────────────────────────────────────────────────────
Mean accept rate: 0.93  │  Divergences: 0
===== instrument-1
4 chains × 1000 draws per chain

Parameter        mean      std     hdi_3%    hdi_97%   ess_bulk   ess_tail    r_hat  mcse_mean
──────────────────────────────────────────────────────────────────────────────────────────────
offset         0.0737   0.0177     0.0405     0.1078       4043       2687   1.0009   0.000280
gain           1.0202   0.0150     0.9902     1.0465       3903       2972   1.0011   0.000241
noise          0.1513   0.0127     0.1290     0.1760       3685       2957   1.0009   0.000211
──────────────────────────────────────────────────────────────────────────────────────────────
Mean accept rate: 0.91  │  Divergences: 0
===== instrument-2
4 chains × 1000 draws per chain

Parameter        mean      std     hdi_3%    hdi_97%   ess_bulk   ess_tail    r_hat  mcse_mean
──────────────────────────────────────────────────────────────────────────────────────────────
offset         0.2031   0.0182     0.1694     0.2368       3871       2985   0.9999   0.000293
gain           1.1013   0.0157     1.0723     1.1305       3923       3110   1.0004   0.000250
noise          0.2012   0.0131     0.1777     0.2259       3533       2869   1.0000   0.000224
──────────────────────────────────────────────────────────────────────────────────────────────
Mean accept rate: 0.92  │  Divergences: 0

Reading the result

The true offsets are 0.0, 0.1 and 0.2 and the true gains 1.00, 1.05 and 1.10. The instrument with 40 readings has the widest intervals, as it should: fewer readings, less information, and no borrowing from the other two.

Interval width does not fall monotonically with sample size here. Instrument 2 has 120 readings against instrument 1's 75 and still reports slightly wider intervals. The counts are fixed at 40, 75 and 120; the noise scale is not, and instrument 2's came out higher on this simulated draw. Width falls with the count and rises with the noise, so a larger instrument can still be the less precisely measured one.

Check r_hat and ess_bulk per instrument before comparing them. A batch reports diagnostics per fit for exactly this reason -- one instrument can fail to converge while the others are fine.