Final Project: Monte Carlo Uncertainty, Scaling Laws, and Geometry Stress Tests for Bayesian Samplers

A reproducible PyMC benchmark spanning applied models, repeated-seed uncertainty quantification, high-dimensional scaling, and Neal’s Funnel

1 Introduction

Our original proposal asked a practical question from computational Bayesian statistics: when do more sophisticated MCMC methods such as Hamiltonian Monte Carlo (HMC) and the No-U-Turn Sampler (NUTS) justify their extra computational cost relative to simpler baselines such as Metropolis-Hastings?

Besides comparison of algorithms from a simple Monte Carlo samling task, this project seeks to include more variables in the experiments to study performance of these algorithms from more perspectives. Specifically, this project will conduct following experiments:

1. Repeated-seed Monte Carlo benchmarking, so the benchmark itself has uncertainty quantification rather than only point estimates.
2. Dimension-scaling experiments, so we can estimate how statistical efficiency decays as parameter dimension increases.
3. Exact-posterior accuracy checks on a known Gaussian target, so efficiency can be compared to actual estimation error.
4. Approximate-inference comparisons, adding mean-field and full-rank ADVI to the benchmark.
5. Budget-response curves, so we can study how extra computation changes performance rather than relying on one fixed budget.
6. Centered-versus-non-centered parameterization checks in the hierarchical model.
7. Posterior-geometry stress tests, centered on Neal’s Funnel, to study sampler behavior on pathologically hard targets.
8. Posterior predictive checks to see whether the sampled posterior can recover similar data to what is observed
9. Simulation-based calibration to verify whether the posterior sample produced by algorithms are similar to the actual posterior distribution

The purpose of this project is to explore how sampler performance varies (1) across applied workloads, (2) across random seeds, (3) across dimensions. And how they are sensitive to parametrizations, difficult parameter geometry. Besides performance, this project checks whehter samplers recover real posterior.

This report studies eight related questions:

1. Which samplers perform best on a realistic portfolio of applied Bayesian models?
2. How much uncertainty is there in the benchmark itself when we repeat runs under different random seeds?
3. How quickly does efficiency degrade with increasing parameter dimension?
4. Which exact samplers are most accurate when the true posterior is known?
5. How do approximate methods such as ADVI compare to exact samplers on speed and posterior error?
6. Which methods improve meaningfully when we increase the computational budget?
7. How much do centered versus non-centered parameterizations change the apparent sampler ranking?
8. Which samplers remain reliable when posterior geometry becomes genuinely pathological?

2 Study Design

The notebook is organized as ten connected benchmark blocks:

1. A broad applied battery of five model classes.
2. A repeated-seed uncertainty study over the same applied scenarios.
3. A high-dimensional scaling-law study on correlated Gaussian targets.
4. An exact-posterior accuracy check on a medium-high-dimensional Gaussian target.
5. An approximate-inference comparison using mean-field and full-rank ADVI.
6. A budget-response study over short, medium, and long compute budgets.
7. A hierarchical parameterization study comparing centered and non-centered forms.
8. A dedicated Neal’s Funnel stress test.
9. Posterior predictive checks for the applied models.
10. Simulation-based calibration checks on representative easy and hierarchical targets.

This structure deliberately balances breadth and depth. The applied battery makes the report feel useful to practitioners, while the exact-accuracy, approximation, budget, parameterization, and funnel sections move it much closer to a computational-statistics study.

2.1 Benchmark Assets

All reported HTML results below use the same PyMC sampling configuration: 4 chains, 500 tuning iterations (burn-in), and 1000 retained draws per chain. That gives 4000 post-warmup draws per model-sampler fit.

Block	Purpose	Runs	Seeds	Sampling configuration
Applied battery	Single-seed cross-scenario benchmark for breadth	20	1	4 chains x (500 tune + 1000 draws per chain)
Repeated-seed uncertainty	Estimate benchmark variability and win probabilities	400	20	4 chains x (500 tune + 1000 draws per chain)
Scaling law study	Estimate ESS/sec ~ p^(-alpha) under correlation stress	160	8	4 chains x (500 tune + 1000 draws per chain)
Exact-posterior accuracy	Measure posterior mean, variance, and coverage error against a known Gaussian target	32	8	4 chains x (500 tune + 1000 draws per chain)
Approximate inference	Compare ADVI methods to exact samplers on the same known target	16	8	10000 optimization steps + 4000 posterior draws
Budget-response study	Trace runtime-accuracy tradeoffs on the exact Gaussian target and full-sampler budget diagnostics on harder targets	108	4 for Gaussian; 2 for hard targets	Short: NUTS = 4 chains x (100 tune + 200 draws per chain); ADVI = 500 optimization steps + 4000 posterior draws
Parameterization study	Compare centered and non-centered hierarchical formulations	40	10	4 chains x (500 tune + 1000 draws per chain)
Neal’s Funnel	Stress-test geometry adaptation on a pathological posterior	80	20	4 chains x (500 tune + 1000 draws per chain)

Sampler	Family	Uses gradients
Metropolis-Hastings	Random-walk MH	No
DE-Metropolis	Differential evolution MH	No
Hamiltonian MC	Static HMC	Yes
NUTS	Dynamic HMC	Yes

Scenario	Data source	Likelihood	Observations	Approx. parameters	Why included
Gaussian	Marginal response from scikit-learn diabetes	Normal likelihood	260	2	An easy baseline where simple samplers can still look competitive.
Linear regression	scikit-learn diabetes dataset	Gaussian regression	240	12	A medium-dimensional continuous model where gradients should matter.
Logistic regression	scikit-learn breast cancer dataset	Bernoulli-logit regression	569	9	A non-Gaussian classification task with correlated predictors.
Hierarchical	Synthetic grouped outcomes	Varying-intercept Gaussian model	303	22	A structured multilevel model with shrinkage and uneven group sizes.
Zero-inflated count	Synthetic insurance-style claims	Zero-inflated Poisson	300	22	A hard applied target combining sparse counts, many zeros, and linked latent structure.

3 Applied Portfolio Benchmark

We begin with a broad single-seed battery. This section provides the same kind of applied overview as the previous version of the project, but now it serves as only the first layer of the analysis rather than the whole story.

portfolio_models = {scenario_name: scenarios[scenario_name].builder() for scenario_name in scenario_order}
portfolio_rows = []
portfolio_idatas = {scenario_name: {} for scenario_name in scenario_order}

for scenario_name in scenario_order:
    scenario = scenarios[scenario_name]
    model = portfolio_models[scenario_name]
    for sampler_name in sampler_order:
        idata, metrics = run_model_with_sampler(
            model=model,
            sampler_name=sampler_name,
            summary_vars=scenario.summary_vars,
            sample_config=PORTFOLIO_CONFIG,
            seed=PORTFOLIO_SEED,
        )
        portfolio_idatas[scenario_name][sampler_name] = idata
        portfolio_rows.append(
            {
                "scenario": scenario_name,
                "scenario_short": scenario.short_name,
                "scenario_abbrev": scenario_abbrev[scenario_name],
                "sampler": sampler_name,
                **metrics,
            }
        )

portfolio_results_df = finalize_results_df(pd.DataFrame(portfolio_rows), scenario_order)
portfolio_results_df["runtime_rank"] = portfolio_results_df.groupby("scenario", observed=True)["runtime_s"].rank(method="dense", ascending=True)
portfolio_results_df["ess_rank"] = portfolio_results_df.groupby("scenario", observed=True)["ess_per_second"].rank(method="dense", ascending=False)
portfolio_results_df["stability_rank"] = portfolio_results_df.groupby("scenario", observed=True)["max_rhat"].rank(method="dense", ascending=True)
portfolio_results_df["composite_rank"] = portfolio_results_df[["runtime_rank", "ess_rank", "stability_rank"]].mean(axis=1)

leaderboard_df = (
    portfolio_results_df.groupby("sampler", observed=True)
    .agg(
        mean_runtime_s=("runtime_s", "mean"),
        mean_ess_per_second=("ess_per_second", "mean"),
        mean_max_rhat=("max_rhat", "mean"),
        reliable_scenarios=("reliable_run", "sum"),
        mean_composite_rank=("composite_rank", "mean"),
    )
    .reset_index()
    .sort_values(["mean_composite_rank", "mean_max_rhat", "mean_runtime_s"])
)

scenario_winners_df = pd.DataFrame(
    [
        {
            "Scenario": scenario_name,
            "Fastest": sub.loc[sub["runtime_s"].idxmin(), "sampler"],
            "Best ESS/s": sub.loc[sub["ess_per_second"].idxmax(), "sampler"],
            "Best stability": sub.loc[sub["max_rhat"].idxmin(), "sampler"],
            "Best composite": sub.loc[sub["composite_rank"].idxmin(), "sampler"],
        }
        for scenario_name, sub in portfolio_results_df.groupby("scenario", observed=True)
    ]
)

sampler	mean_runtime_s	mean_ess_per_second	mean_max_rhat	reliable_scenarios	mean_composite_rank
Hamiltonian MC	238.884	36.761	1	5	1.867
NUTS	284.427	15.808	1	5	2.2
Metropolis-Hastings	233.272	2.975	1.238	3	2.4
DE-Metropolis	211.855	3.971	2.182	1	2.4

A compact visual summary of the applied benchmark portfolio.

Single-run portfolio view: ESS per second and maximum r-hat.

Single-run speed-efficiency frontier. Upper-left points are most attractive.

Scenario	Fastest	Best ESS/s	Best stability	Best composite
Gaussian location model	DE-Metropolis	Hamiltonian MC	DE-Metropolis	DE-Metropolis
Bayesian linear regression	DE-Metropolis	Hamiltonian MC	Hamiltonian MC	Hamiltonian MC
Bayesian logistic regression	DE-Metropolis	NUTS	Hamiltonian MC	Hamiltonian MC
Hierarchical regression	DE-Metropolis	NUTS	Hamiltonian MC	Hamiltonian MC
Zero-inflated count regression	DE-Metropolis	NUTS	Hamiltonian MC	Hamiltonian MC

The broad single-seed benchmark already shows a pattern that will keep reappearing: simple random-walk methods often look attractive on raw runtime, but geometry-aware samplers dominate once we care about reliable effective sampling. The problem is that a single run cannot tell us how stable those rankings really are. That is the purpose of the next section.

4 Repeated-Seed Monte Carlo Benchmarking

Single-run benchmarks produce point estimates. This section upgrades the study by running the applied battery repeatedly under different random seeds so the benchmark itself has uncertainty quantification.

For each scenario-sampler pair, we now report:

1. Mean performance over repeated runs.
2. Empirical 95% intervals across repeated runs.
3. Reliability frequency.
4. The probability that a sampler is the best within a scenario.

repeated_models = {scenario_name: scenarios[scenario_name].builder() for scenario_name in scenario_order}
repeated_rows = []

for scenario_name in scenario_order:
    scenario = scenarios[scenario_name]
    model = repeated_models[scenario_name]
    for replicate, seed in enumerate(REPEATED_SEEDS, start=1):
        for sampler_name in sampler_order:
            _, metrics = run_model_with_sampler(
                model=model,
                sampler_name=sampler_name,
                summary_vars=scenario.summary_vars,
                sample_config=REPEATED_CONFIG,
                seed=seed,
            )
            repeated_rows.append(
                {
                    "scenario": scenario_name,
                    "scenario_short": scenario.short_name,
                    "replicate": replicate,
                    "seed": seed,
                    "sampler": sampler_name,
                    **metrics,
                }
            )

repeated_results_df = finalize_results_df(pd.DataFrame(repeated_rows), scenario_order)
repeated_summary_df = summarize_repeated_block(repeated_results_df, scenario_order)

repeated_display_df = repeated_summary_df.assign(
    **{
        "ESS/sec mean [95% interval]": [
            interval_string(row.ess_mean, row.ess_low, row.ess_high, digits=1)
            for row in repeated_summary_df.itertuples()
        ],
        "Runtime mean [95% interval]": [
            interval_string(row.runtime_mean, row.runtime_low, row.runtime_high, digits=2)
            for row in repeated_summary_df.itertuples()
        ],
        "Mean max r-hat": repeated_summary_df["rhat_mean"].round(3),
        "Reliability rate": repeated_summary_df["reliability_rate"].round(2),
        "Pr(best ESS/s)": repeated_summary_df["pr_best_ess"].round(2),
    }
)[
    [
        "scenario_short",
        "sampler",
        "ESS/sec mean [95% interval]",
        "Runtime mean [95% interval]",
        "Mean max r-hat",
        "Reliability rate",
        "Pr(best ESS/s)",
    ]
].rename(columns={"scenario_short": "Scenario", "sampler": "Sampler"})

prob_best_heatmap = repeated_summary_df.pivot(index="scenario_short", columns="sampler", values="pr_best_ess")
reliability_heatmap = repeated_summary_df.pivot(index="scenario_short", columns="sampler", values="reliability_rate")

hardest_applied_scenario = "Zero-inflated count regression"
zip_pairwise_ess = pairwise_probability_matrix(
    repeated_results_df.loc[repeated_results_df["scenario"] == hardest_applied_scenario],
    metric="ess_per_second",
    higher_is_better=True,
)

Scenario	Sampler	ESS/sec mean [95% interval]	Runtime mean [95% interval]	Mean max r-hat	Reliability rate	Pr(best ESS/s)
Gaussian	Metropolis-Hastings	9.5 [1.6, 22.2]	82.29 [18.71, 216.40]	1	1	0
Gaussian	DE-Metropolis	11.5 [2.2, 23.5]	88.06 [21.68, 215.55]	1	1	0
Gaussian	Hamiltonian MC	71.6 [17.9, 154.7]	102.18 [26.15, 218.03]	1	1	1
Gaussian	NUTS	74.1 [14.5, 157.1]	99.11 [21.89, 218.07]	1	1	0
Linear regression	Metropolis-Hastings	2.5 [0.5, 6.8]	119.42 [26.35, 295.44]	2	0	0
Linear regression	DE-Metropolis	0.1 [0.0, 0.3]	99.86 [20.51, 247.72]	3	0	0
Linear regression	Hamiltonian MC	79.1 [15.3, 167.8]	117.44 [30.12, 311.89]	1	1	1
Linear regression	NUTS	52.0 [12.0, 109.6]	114.01 [37.00, 313.44]	1	1	0
Logistic regression	Metropolis-Hastings	5.2 [1.5, 8.0]	62.30 [29.66, 154.05]	1	1	0
Logistic regression	DE-Metropolis	0.1 [0.0, 0.2]	70.21 [27.89, 200.34]	3	0	0
Logistic regression	Hamiltonian MC	70.5 [17.4, 145.6]	77.36 [26.28, 213.45]	1	1	0
Logistic regression	NUTS	87.0 [19.4, 163.5]	75.69 [26.20, 216.86]	1	1	1
Hierarchical	Metropolis-Hastings	8.9 [2.3, 20.9]	105.15 [26.27, 227.70]	1	1	0
Hierarchical	DE-Metropolis	0.2 [0.0, 0.3]	63.52 [19.49, 192.37]	3	0	0
Hierarchical	Hamiltonian MC	94.4 [24.0, 177.8]	90.74 [27.28, 211.48]	1	1	0
Hierarchical	NUTS	101.0 [28.6, 206.5]	103.76 [30.95, 222.63]	1	1	1
Zero-inflated count	Metropolis-Hastings	2.5 [0.6, 4.0]	94.31 [39.74, 247.26]	1	0	0
Zero-inflated count	DE-Metropolis	0.1 [0.0, 0.2]	80.35 [31.71, 215.32]	2	0	0
Zero-inflated count	Hamiltonian MC	35.3 [7.8, 53.7]	86.92 [41.79, 256.31]	1	1	0
Zero-inflated count	NUTS	52.1 [22.5, 82.4]	80.58 [45.39, 189.56]	1	1	1

Repeated-seed uncertainty intervals for ESS per second across the applied scenarios.

Monte Carlo benchmark uncertainty summarized as probabilities and reliability frequencies.

Pairwise win probabilities on ESS/sec for the hardest applied scenario: zero-inflated count regression.

Scenario	Top sampler by mean ESS/s	Pr(best ESS/s)	Reliability rate
Gaussian location model	NUTS	0.2	1
Bayesian linear regression	Hamiltonian MC	0.85	1
Bayesian logistic regression	NUTS	0.65	1
Hierarchical regression	NUTS	0.75	1
Zero-inflated count regression	NUTS	0.75	1

The repeated-seed study changes the interpretation of the benchmark in an important way. Instead of saying only that a sampler “won” on ESS per second, we can now ask how often that result recurs. In a research or industry setting, that is a much more defensible statement than a single-run ranking.

5 Dimension Scaling Laws

The next question is how fast does efficiency decay as the dimension of the target increases?

To study that, we benchmark the samplers on a family of correlated Gaussian targets with fixed condition number \(\kappa(\Sigma)=10^3\) and dimensions

1. \(p = 2\)
2. \(p = 10\)
3. \(p = 50\)
4. \(p = 100\)
5. \(p = 250\)

We then estimate an empirical scaling law of the form

\[ \text{ESS/sec} \propto p^{-\alpha} \]

where smaller \(\alpha\) means slower degradation with dimension.

scaling_rows = []
scaling_models = {}
scaling_kappas = {}

for p in DIMENSION_GRID:
    model, condition_number = build_correlated_gaussian_model(p)
    scaling_models[p] = model
    scaling_kappas[p] = condition_number

for p in DIMENSION_GRID:
    model = scaling_models[p]
    for replicate, seed in enumerate(SCALING_SEEDS, start=1):
        for sampler_name in sampler_order:
            _, metrics = run_model_with_sampler(
                model=model,
                sampler_name=sampler_name,
                summary_vars=["theta"],
                sample_config=SCALING_CONFIG,
                seed=seed,
            )
            scaling_rows.append(
                {
                    "scenario": "Correlated Gaussian scaling",
                    "scenario_short": "Scaling",
                    "p": p,
                    "replicate": replicate,
                    "sampler": sampler_name,
                    "condition_number": scaling_kappas[p],
                    **metrics,
                }
            )

scaling_results_df = finalize_results_df(pd.DataFrame(scaling_rows))

scaling_summary_df = (
    scaling_results_df.groupby(["p", "sampler"], observed=True)
    .agg(
        ess_mean=("ess_per_second", "mean"),
        ess_low=("ess_per_second", lambda x: empirical_interval(x)[0]),
        ess_high=("ess_per_second", lambda x: empirical_interval(x)[1]),
        runtime_mean=("runtime_s", "mean"),
    )
    .reset_index()
)

scaling_alpha_rows = []
for sampler_name, sub in scaling_results_df.groupby("sampler", observed=True):
    mean_by_p = sub.groupby("p", observed=True)["ess_per_second"].mean().sort_index()
    ess_vals = mean_by_p.to_numpy(dtype=float)
    p_vals = mean_by_p.index.to_numpy(dtype=float)
    valid = np.isfinite(ess_vals) & (ess_vals > 0)
    if valid.sum() < 2:
        scaling_alpha_rows.append(
            {"Sampler": sampler_name, "Estimated alpha": np.nan,
             "95% interval low": np.nan, "95% interval high": np.nan, "R^2": np.nan}
        )
        continue
    log_p = np.log(p_vals[valid])
    log_ess = np.log(ess_vals[valid])
    slope, intercept = np.polyfit(log_p, log_ess, 1)
    fitted = slope * log_p + intercept
    r_squared = 1 - np.sum((log_ess - fitted) ** 2) / np.sum((log_ess - np.mean(log_ess)) ** 2)
    low_alpha, high_alpha = bootstrap_scaling_alpha(sub)
    scaling_alpha_rows.append(
        {
            "Sampler": sampler_name,
            "Estimated alpha": -float(slope),
            "95% interval low": low_alpha,
            "95% interval high": high_alpha,
            "R^2": float(r_squared),
        }
    )

scaling_alpha_df = pd.DataFrame(scaling_alpha_rows).sort_values("Estimated alpha")

Dimension p	Condition number
2	1000
10	1000
50	1000
100	1000
250	1000

Empirical scaling of ESS per second with parameter dimension on correlated Gaussian targets.

Sampler	alpha [95% interval]	R^2
NUTS	-0.24 [-0.30, -0.16]	0.759
Metropolis-Hastings	0.08 [-0.03, 0.15]	0.298
Hamiltonian MC	0.23 [0.18, 0.29]	0.793
DE-Metropolis	0.96 [0.87, 1.04]	0.935

NUTS achieves the smallest empirical scaling exponent in this fixed-budget experiment with four chains, 500 tuning iterations, and 1000 retained draws per chain, but its fitted relationship is also noisier than the cleaner monotone decay seen for the other samplers. That is why the report treats the scaling section as empirical evidence rather than a final asymptotic theorem.

6 Exact-Posterior Accuracy Check

Efficiency metrics are informative, but they are not the same as posterior accuracy. On the correlated Gaussian target family we know the exact distribution, so we can go beyond ESS/sec and ask whether the samplers with the best efficiency also achieve the smallest estimation error in posterior moments.

We focus on the \(p=50\) target as a representative medium-high-dimensional case. The exact posterior mean is zero and the exact marginal standard deviations are known from the covariance construction, which lets us measure Monte Carlo error directly.

reference_model, reference_kappa = build_correlated_gaussian_model(REFERENCE_DIM)
reference_chol, _ = make_conditioned_covariance(REFERENCE_DIM)
reference_exact_mean = np.zeros(REFERENCE_DIM)
reference_exact_sd = np.sqrt(np.sum(reference_chol**2, axis=1))
reference_rows = []

for replicate, seed in enumerate(REFERENCE_SEEDS, start=1):
    for sampler_name in sampler_order:
        idata, metrics = run_model_with_sampler(
            model=reference_model,
            sampler_name=sampler_name,
            summary_vars=["theta"],
            sample_config=SCALING_CONFIG,
            seed=seed,
        )
        accuracy_metrics = summarize_exact_gaussian_accuracy(
            idata,
            exact_mean=reference_exact_mean,
            exact_sd=reference_exact_sd,
        )
        reference_rows.append(
            {
                "scenario": f"Exact Gaussian p={REFERENCE_DIM}",
                "scenario_short": "Exact Gaussian",
                "replicate": replicate,
                "seed": seed,
                "sampler": sampler_name,
                "condition_number": reference_kappa,
                **metrics,
                **accuracy_metrics,
            }
        )

reference_results_df = finalize_results_df(pd.DataFrame(reference_rows))
reference_summary_df = (
    reference_results_df.groupby("sampler", observed=True)
    .agg(
        ess_mean=("ess_per_second", "mean"),
        ess_low=("ess_per_second", lambda x: empirical_interval(x)[0]),
        ess_high=("ess_per_second", lambda x: empirical_interval(x)[1]),
        mean_rmse_mean=("mean_rmse", "mean"),
        mean_rmse_low=("mean_rmse", lambda x: empirical_interval(x)[0]),
        mean_rmse_high=("mean_rmse", lambda x: empirical_interval(x)[1]),
        sd_rmse_mean=("sd_rmse", "mean"),
        sd_rmse_low=("sd_rmse", lambda x: empirical_interval(x)[0]),
        sd_rmse_high=("sd_rmse", lambda x: empirical_interval(x)[1]),
        coverage_mean=("coverage_95", "mean"),
        coverage_low=("coverage_95", lambda x: empirical_interval(x)[0]),
        coverage_high=("coverage_95", lambda x: empirical_interval(x)[1]),
        runtime_mean=("runtime_s", "mean"),
        runtime_low=("runtime_s", lambda x: empirical_interval(x)[0]),
        runtime_high=("runtime_s", lambda x: empirical_interval(x)[1]),
    )
    .reset_index()
)

reference_display_df = reference_summary_df.assign(
    **{
        "ESS/sec mean [95% interval]": [
            interval_string(row.ess_mean, row.ess_low, row.ess_high, digits=1)
            for row in reference_summary_df.itertuples()
        ],
        "Posterior mean RMSE [95% interval]": [
            interval_string(row.mean_rmse_mean, row.mean_rmse_low, row.mean_rmse_high, digits=4)
            for row in reference_summary_df.itertuples()
        ],
        "Posterior sd RMSE [95% interval]": [
            interval_string(row.sd_rmse_mean, row.sd_rmse_low, row.sd_rmse_high, digits=4)
            for row in reference_summary_df.itertuples()
        ],
        "95% mean coverage [95% interval]": [
            interval_string(row.coverage_mean, row.coverage_low, row.coverage_high, digits=3)
            for row in reference_summary_df.itertuples()
        ],
    }
)[
    [
        "sampler",
        "ESS/sec mean [95% interval]",
        "Posterior mean RMSE [95% interval]",
        "Posterior sd RMSE [95% interval]",
        "95% mean coverage [95% interval]",
    ]
].rename(columns={"sampler": "Sampler"})

exact_reference_summary_df = reference_summary_df.rename(columns={"sampler": "method"}).assign(
    method_family="Exact MCMC"
)

Sampler	ESS/sec mean [95% interval]	Posterior mean RMSE [95% interval]	Posterior sd RMSE [95% interval]	95% mean coverage [95% interval]
Metropolis-Hastings	0.2 [0.0, 0.7]	2.6596 [2.0108, 3.6122]	1.6991 [1.5786, 1.8214]	1.000 [1.000, 1.000]
DE-Metropolis	0.2 [0.0, 0.7]	0.7951 [0.5898, 0.9715]	9.8659 [9.5175, 10.3712]	1.000 [1.000, 1.000]
Hamiltonian MC	12.1 [4.3, 41.1]	0.3931 [0.3088, 0.4856]	0.1833 [0.1553, 0.2097]	1.000 [1.000, 1.000]
NUTS	33.6 [11.9, 78.6]	0.2025 [0.1675, 0.2851]	0.1642 [0.1411, 0.1791]	1.000 [1.000, 1.000]

Exact-posterior accuracy on the p=50 correlated Gaussian target: higher ESS/sec and lower RMSE are preferable.

On the exact Gaussian target, NUTS achieves the smallest posterior-mean RMSE (0.2025), while NUTS also gives the best recovery of marginal posterior standard deviations (0.1642). By contrast, DE-Metropolis has the weakest variance recovery (9.8659), which shows why ESS/sec alone is not a sufficient notion of sampler quality.

7 Approximate Inference Comparison

The exact-Gaussian target also lets us compare approximate inference methods to exact samplers. Here we add mean-field ADVI and full-rank ADVI. Because these methods optimize a variational objective rather than generate Markov chains, runtime and posterior accuracy are more appropriate comparison metrics than r-hat.

approximation_rows = []

for replicate, seed in enumerate(REFERENCE_SEEDS, start=1):
    for approximation_name in APPROXIMATION_ORDER:
        idata, metrics = run_model_with_approximation(
            model=reference_model,
            approximation_name=approximation_name,
            approx_config=APPROXIMATION_CONFIG,
            seed=seed,
        )
        accuracy_metrics = summarize_exact_gaussian_accuracy(
            idata,
            exact_mean=reference_exact_mean,
            exact_sd=reference_exact_sd,
        )
        approximation_rows.append(
            {
                "scenario": f"Exact Gaussian p={REFERENCE_DIM}",
                "scenario_short": "Exact Gaussian",
                "replicate": replicate,
                "seed": seed,
                "method": approximation_name,
                "method_family": "Approximate inference",
                **metrics,
                **accuracy_metrics,
            }
        )

approximation_results_df = pd.DataFrame(approximation_rows)
approximation_summary_df = (
    approximation_results_df.groupby(["method", "method_family"], observed=True)
    .agg(
        runtime_mean=("runtime_s", "mean"),
        runtime_low=("runtime_s", lambda x: empirical_interval(x)[0]),
        runtime_high=("runtime_s", lambda x: empirical_interval(x)[1]),
        mean_rmse_mean=("mean_rmse", "mean"),
        mean_rmse_low=("mean_rmse", lambda x: empirical_interval(x)[0]),
        mean_rmse_high=("mean_rmse", lambda x: empirical_interval(x)[1]),
        sd_rmse_mean=("sd_rmse", "mean"),
        sd_rmse_low=("sd_rmse", lambda x: empirical_interval(x)[0]),
        sd_rmse_high=("sd_rmse", lambda x: empirical_interval(x)[1]),
        coverage_mean=("coverage_95", "mean"),
        coverage_low=("coverage_95", lambda x: empirical_interval(x)[0]),
        coverage_high=("coverage_95", lambda x: empirical_interval(x)[1]),
        final_vi_loss_mean=("final_vi_loss", "mean"),
        final_vi_loss_low=("final_vi_loss", lambda x: empirical_interval(x)[0]),
        final_vi_loss_high=("final_vi_loss", lambda x: empirical_interval(x)[1]),
    )
    .reset_index()
)

approximation_display_df = approximation_summary_df.assign(
    **{
        "Final VI loss [95% interval]": [
            interval_string(row.final_vi_loss_mean, row.final_vi_loss_low, row.final_vi_loss_high, digits=1)
            for row in approximation_summary_df.itertuples()
        ]
    }
)[["method", "Final VI loss [95% interval]"]].rename(columns={"method": "Method"})

approximate_comparison_summary_df = pd.concat(
    [
        exact_reference_summary_df[
            [
                "method",
                "method_family",
                "runtime_mean",
                "runtime_low",
                "runtime_high",
                "mean_rmse_mean",
                "mean_rmse_low",
                "mean_rmse_high",
                "sd_rmse_mean",
                "sd_rmse_low",
                "sd_rmse_high",
                "coverage_mean",
                "coverage_low",
                "coverage_high",
            ]
        ],
        approximation_summary_df[
            [
                "method",
                "method_family",
                "runtime_mean",
                "runtime_low",
                "runtime_high",
                "mean_rmse_mean",
                "mean_rmse_low",
                "mean_rmse_high",
                "sd_rmse_mean",
                "sd_rmse_low",
                "sd_rmse_high",
                "coverage_mean",
                "coverage_low",
                "coverage_high",
            ]
        ],
    ],
    ignore_index=True,
)

approximate_comparison_display_df = approximate_comparison_summary_df.assign(
    **{
        "Runtime mean [95% interval]": [
            interval_string(row.runtime_mean, row.runtime_low, row.runtime_high, digits=2)
            for row in approximate_comparison_summary_df.itertuples()
        ],
        "Posterior mean RMSE [95% interval]": [
            interval_string(row.mean_rmse_mean, row.mean_rmse_low, row.mean_rmse_high, digits=4)
            for row in approximate_comparison_summary_df.itertuples()
        ],
        "Posterior sd RMSE [95% interval]": [
            interval_string(row.sd_rmse_mean, row.sd_rmse_low, row.sd_rmse_high, digits=4)
            for row in approximate_comparison_summary_df.itertuples()
        ],
        "95% mean coverage [95% interval]": [
            interval_string(row.coverage_mean, row.coverage_low, row.coverage_high, digits=3)
            for row in approximate_comparison_summary_df.itertuples()
        ],
    }
)[
    [
        "method",
        "method_family",
        "Runtime mean [95% interval]",
        "Posterior mean RMSE [95% interval]",
        "Posterior sd RMSE [95% interval]",
        "95% mean coverage [95% interval]",
    ]
].rename(columns={"method": "Method", "method_family": "Family"})

approximate_comparison_results_df = pd.concat(
    [
        reference_results_df.rename(columns={"sampler": "method"}).assign(method_family="Exact MCMC")[
            ["method", "method_family", "coverage_95"]
        ],
        approximation_results_df[["method", "method_family", "coverage_95"]],
    ],
    ignore_index=True,
)

Method	Family	Runtime mean [95% interval]	Posterior mean RMSE [95% interval]	Posterior sd RMSE [95% interval]	95% mean coverage [95% interval]
Metropolis-Hastings	Exact MCMC	148.52 [18.89, 260.30]	2.6596 [2.0108, 3.6122]	1.6991 [1.5786, 1.8214]	1.000 [1.000, 1.000]
DE-Metropolis	Exact MCMC	173.16 [26.61, 261.96]	0.7951 [0.5898, 0.9715]	9.8659 [9.5175, 10.3712]	1.000 [1.000, 1.000]
Hamiltonian MC	Exact MCMC	180.12 [33.83, 268.83]	0.3931 [0.3088, 0.4856]	0.1833 [0.1553, 0.2097]	1.000 [1.000, 1.000]
NUTS	Exact MCMC	165.02 [40.08, 274.38]	0.2025 [0.1675, 0.2851]	0.1642 [0.1411, 0.1791]	1.000 [1.000, 1.000]
Full-rank ADVI	Approximate inference	85.18 [38.52, 257.50]	0.0470 [0.0388, 0.0551]	9.3716 [9.3607, 9.3780]	1.000 [1.000, 1.000]
Mean-field ADVI	Approximate inference	92.90 [32.20, 293.15]	0.0349 [0.0297, 0.0412]	10.1655 [10.1583, 10.1728]	1.000 [1.000, 1.000]

Method	Final VI loss [95% interval]
Full-rank ADVI	25.0 [21.7, 27.7]
Mean-field ADVI	40.5 [30.7, 54.5]

Approximate inference versus exact samplers on the known Gaussian target.

Within the variational family, Mean-field ADVI gives the smallest posterior-mean RMSE (0.0349), and Full-rank ADVI is the fastest approximation (85.18s on average). Within that same family, Full-rank ADVI gives the best variance recovery (9.3716), but both ADVI variants remain much weaker on posterior standard deviations than the best exact methods: the strongest exact mean estimator is NUTS (mean RMSE 0.2025), while NUTS gives the smallest sd RMSE (0.1642).

The results show an asymmetric tradeoff. The ADVI methods produce fast, low-error posterior mean estimates on this Gaussian target, but they perform poorly in recovering posterior uncertainty. This makes them useful for quick posterior screening, initialization, or rough effect ranking, while exact samplers remain the safer choice for final uncertainty reporting.

8 Budget-Response Curves

8.1 Budget-response on the exact Gaussian target

A single computational budget can obscure important behavior. Some methods improve steadily as more computation is allocated, while others plateau quickly or remain biased even when optimization is extended. To study this directly, we compare short, medium, and long budgets on the same exact Gaussian target for NUTS, mean-field ADVI, and full-rank ADVI.

budget_rows = []

for replicate, seed in enumerate(BUDGET_SEEDS, start=1):
    for budget_label in BUDGET_LEVEL_ORDER:
        idata, metrics = run_model_with_sampler(
            model=reference_model,
            sampler_name="NUTS",
            summary_vars=["theta"],
            sample_config=BUDGET_MCMC_CONFIGS[budget_label],
            seed=seed,
        )
        accuracy_metrics = summarize_exact_gaussian_accuracy(
            idata,
            exact_mean=reference_exact_mean,
            exact_sd=reference_exact_sd,
        )
        budget_rows.append(
            {
                "replicate": replicate,
                "seed": seed,
                "budget": budget_label,
                "budget_index": BUDGET_LEVEL_ORDER.index(budget_label) + 1,
                "method": "NUTS",
                "method_family": "Exact MCMC",
                **metrics,
                **accuracy_metrics,
            }
        )

        for approximation_name in APPROXIMATION_ORDER:
            idata, metrics = run_model_with_approximation(
                model=reference_model,
                approximation_name=approximation_name,
                approx_config=BUDGET_APPROXIMATION_CONFIGS[budget_label],
                seed=seed,
            )
            accuracy_metrics = summarize_exact_gaussian_accuracy(
                idata,
                exact_mean=reference_exact_mean,
                exact_sd=reference_exact_sd,
            )
            budget_rows.append(
                {
                    "replicate": replicate,
                    "seed": seed,
                    "budget": budget_label,
                    "budget_index": BUDGET_LEVEL_ORDER.index(budget_label) + 1,
                    "method": approximation_name,
                    "method_family": "Approximate inference",
                    **metrics,
                    **accuracy_metrics,
                }
            )

budget_results_df = pd.DataFrame(budget_rows)
budget_results_df["method"] = pd.Categorical(
    budget_results_df["method"],
    categories=BUDGET_METHODS,
    ordered=True,
)
budget_results_df["budget"] = pd.Categorical(
    budget_results_df["budget"],
    categories=BUDGET_LEVEL_ORDER,
    ordered=True,
)

budget_summary_df = (
    budget_results_df.groupby(["method", "method_family", "budget", "budget_index"], observed=True)
    .agg(
        runtime_mean=("runtime_s", "mean"),
        runtime_low=("runtime_s", lambda x: empirical_interval(x)[0]),
        runtime_high=("runtime_s", lambda x: empirical_interval(x)[1]),
        mean_rmse_mean=("mean_rmse", "mean"),
        mean_rmse_low=("mean_rmse", lambda x: empirical_interval(x)[0]),
        mean_rmse_high=("mean_rmse", lambda x: empirical_interval(x)[1]),
        sd_rmse_mean=("sd_rmse", "mean"),
        sd_rmse_low=("sd_rmse", lambda x: empirical_interval(x)[0]),
        sd_rmse_high=("sd_rmse", lambda x: empirical_interval(x)[1]),
        coverage_mean=("coverage_95", "mean"),
        coverage_low=("coverage_95", lambda x: empirical_interval(x)[0]),
        coverage_high=("coverage_95", lambda x: empirical_interval(x)[1]),
    )
    .reset_index()
    .sort_values(["method", "budget_index"])
)

budget_display_df = budget_summary_df.assign(
    **{
        "Runtime mean [95% interval]": [
            interval_string(row.runtime_mean, row.runtime_low, row.runtime_high, digits=2)
            for row in budget_summary_df.itertuples()
        ],
        "Posterior mean RMSE [95% interval]": [
            interval_string(row.mean_rmse_mean, row.mean_rmse_low, row.mean_rmse_high, digits=4)
            for row in budget_summary_df.itertuples()
        ],
        "Posterior sd RMSE [95% interval]": [
            interval_string(row.sd_rmse_mean, row.sd_rmse_low, row.sd_rmse_high, digits=4)
            for row in budget_summary_df.itertuples()
        ],
        "95% mean coverage [95% interval]": [
            interval_string(row.coverage_mean, row.coverage_low, row.coverage_high, digits=3)
            for row in budget_summary_df.itertuples()
        ],
    }
)[
    [
        "method",
        "budget",
        "Runtime mean [95% interval]",
        "Posterior mean RMSE [95% interval]",
        "Posterior sd RMSE [95% interval]",
        "95% mean coverage [95% interval]",
    ]
].rename(columns={"method": "Method", "budget": "Budget"})

Method	Budget	Runtime mean [95% interval]	Posterior mean RMSE [95% interval]	Posterior sd RMSE [95% interval]	95% mean coverage [95% interval]
NUTS	Short	71.94 [30.45, 113.73]	0.3986 [0.3247, 0.4684]	0.4150 [0.3602, 0.4649]	1.000 [1.000, 1.000]
NUTS	Medium	53.70 [21.85, 110.29]	0.3203 [0.2232, 0.4270]	0.2680 [0.2350, 0.2848]	1.000 [1.000, 1.000]
NUTS	Long	39.21 [37.03, 42.84]	0.1877 [0.1261, 0.2276]	0.1595 [0.1383, 0.1704]	1.000 [1.000, 1.000]
Mean-field ADVI	Short	96.73 [34.35, 159.16]	0.0132 [0.0129, 0.0138]	11.5650 [11.5635, 11.5666]	1.000 [1.000, 1.000]
Mean-field ADVI	Medium	86.37 [45.33, 156.19]	0.0196 [0.0182, 0.0211]	11.3037 [11.3016, 11.3056]	1.000 [1.000, 1.000]
Mean-field ADVI	Long	63.60 [46.20, 102.64]	0.0389 [0.0360, 0.0414]	10.1629 [10.1575, 10.1705]	1.000 [1.000, 1.000]
Full-rank ADVI	Short	95.25 [56.08, 143.72]	0.0240 [0.0230, 0.0249]	10.9390 [10.9364, 10.9424]	1.000 [1.000, 1.000]
Full-rank ADVI	Medium	46.13 [30.21, 68.99]	0.0304 [0.0281, 0.0321]	10.6429 [10.6404, 10.6474]	1.000 [1.000, 1.000]
Full-rank ADVI	Long	44.33 [40.77, 48.96]	0.0517 [0.0471, 0.0539]	9.3636 [9.3577, 9.3708]	1.000 [1.000, 1.000]

Budget-response curves on the exact Gaussian target.

Across increasing budgets, Mean-field ADVI achieves the smallest long-budget posterior-mean RMSE (0.0389), while NUTS gives the best long-budget posterior-sd RMSE (0.1595). For NUTS, posterior-mean RMSE drops from 0.3986 at the short budget to 0.1877 at the long budget, and posterior-sd RMSE drops from 0.4150 to 0.1595. By contrast, Mean-field ADVI stays extremely poor on uncertainty recovery (11.5650 to 10.1629), while Full-rank ADVI improves only modestly on that axis (10.9390 to 9.3636).

NUTS benefits from added compute on both mean and uncertainty error, whereas the variational methods reach very strong posterior means quickly but plateau at much poorer variance recovery. This provides a more substantive conclusion than a single-budget win table.

8.2 Budget-response on harder applied targets and funnel-like geometries

The exact-Gaussian budget-response study above evaluates runtime against posterior accuracy because the Gaussian target has known posterior moments. For the harder applied models and Neal’s Funnel, no analytic posterior is available, so this subsection focuses instead on efficiency and convergence diagnostics. We run the full MCMC sampler portfolio over the same short, medium, and long budget grid used above, keeping the number of chains fixed and varying only the amount of tuning and posterior sampling. This design tests whether additional computation improves effective sampling and reliability, and whether difficult posterior geometry can be resolved by longer chains alone. We use hierarchical regression, zero-inflated count regression, and Neal’s Funnel to represent increasingly challenging posterior structures.

# Run the raw full-sampler budget-response experiments on harder targets.
# Set to False after the experiments have already been run and the raw CSV has been saved.
RUN_HARD_BUDGET = True
HARD_BUDGET_RAW_PATH = "hard_budget_results_raw.csv"

HARD_BUDGET_SEEDS = BUDGET_SEEDS[:2]

HARD_BUDGET_TARGETS = {
    "Hierarchical regression": {
        "target_short": "Hierarchical",
        "target_family": "Hard applied",
        "builder": build_hierarchical_model,
        "summary_vars": scenarios["Hierarchical regression"].summary_vars,
    },
    "Zero-inflated count regression": {
        "target_short": "Zero-inflated count",
        "target_family": "Hard applied",
        "builder": build_zip_model,
        "summary_vars": scenarios["Zero-inflated count regression"].summary_vars,
    },
    "Neal's Funnel": {
        "target_short": "Funnel",
        "target_family": "Geometry stress test",
        "builder": build_funnel_model,
        "summary_vars": ["z", "x"],
    },
}

HARD_BUDGET_TARGET_ORDER = list(HARD_BUDGET_TARGETS.keys())
HARD_BUDGET_SAMPLERS = sampler_order


def failed_hard_budget_row(
    target_name,
    target_info,
    sampler_name,
    budget_label,
    replicate,
    seed,
    error_message,
):
    return {
        "target": target_name,
        "target_short": target_info["target_short"],
        "target_family": target_info["target_family"],
        "replicate": replicate,
        "seed": seed,
        "budget": budget_label,
        "budget_index": BUDGET_LEVEL_ORDER.index(budget_label) + 1,
        "sampler": sampler_name,
        "sampler_family": sampler_specs[sampler_name]["family"],
        "uses_gradients": sampler_specs[sampler_name]["gradients"],
        "runtime_s": np.nan,
        "peak_mem_mib": np.nan,
        "median_ess_bulk": np.nan,
        "min_ess_bulk": np.nan,
        "median_ess_tail": np.nan,
        "ess_per_second": np.nan,
        "tail_ess_per_second": np.nan,
        "max_rhat": np.nan,
        "divergences": np.nan,
        "mean_accept_stat": np.nan,
        "mean_n_steps": np.nan,
        "mean_step_size": np.nan,
        "mean_tree_depth": np.nan,
        "max_treedepth_rate": np.nan,
        "mean_bfmi": np.nan,
        "reliable_run": False,
        "failed_run": True,
        "error_message": str(error_message),
    }


def normalize_failed_flag(series):
    return series.astype(str).str.lower().isin(["true", "1", "yes"])


def completed_hard_budget_keys(df):
    if df.empty:
        return set()

    if "failed_run" not in df.columns:
        return set()

    working = df.copy()
    working["failed_run_bool"] = normalize_failed_flag(working["failed_run"])

    completed = working.loc[~working["failed_run_bool"]].copy()

    if completed.empty:
        return set()

    return set(
        zip(
            completed["target"].astype(str),
            completed["sampler"].astype(str),
            completed["budget"].astype(str),
            completed["seed"].astype(int),
        )
    )


def all_hard_budget_jobs():
    jobs = []

    for target_name in HARD_BUDGET_TARGET_ORDER:
        for replicate, seed in enumerate(HARD_BUDGET_SEEDS, start=1):
            for budget_label in BUDGET_LEVEL_ORDER:
                for sampler_name in HARD_BUDGET_SAMPLERS:
                    key = (
                        target_name,
                        sampler_name,
                        budget_label,
                        int(seed),
                    )

                    jobs.append(
                        {
                            "key": key,
                            "target_name": target_name,
                            "replicate": replicate,
                            "seed": seed,
                            "budget_label": budget_label,
                            "sampler_name": sampler_name,
                        }
                    )

    return jobs


if os.path.exists(HARD_BUDGET_RAW_PATH):
    hard_budget_results_df = pd.read_csv(HARD_BUDGET_RAW_PATH)
else:
    hard_budget_results_df = pd.DataFrame()

if RUN_HARD_BUDGET:
    hard_budget_rows = (
        hard_budget_results_df.to_dict("records")
        if not hard_budget_results_df.empty
        else []
    )

    done_keys = completed_hard_budget_keys(hard_budget_results_df)
    planned_jobs = all_hard_budget_jobs()
    remaining_jobs = [job for job in planned_jobs if job["key"] not in done_keys]

    print(f"Planned jobs: {len(planned_jobs)}")
    print(f"Completed jobs found in CSV: {len(done_keys)}")
    print(f"Remaining jobs to run: {len(remaining_jobs)}")

    for run_id, job in enumerate(remaining_jobs, start=1):
        target_name = job["target_name"]
        replicate = job["replicate"]
        seed = job["seed"]
        budget_label = job["budget_label"]
        sampler_name = job["sampler_name"]
        target_info = HARD_BUDGET_TARGETS[target_name]

        print(
            f"[{run_id}/{len(remaining_jobs)}] "
            f"{target_name} | {sampler_name} | {budget_label} | seed={seed}"
        )

        try:
            model = target_info["builder"]()

            _, metrics = run_model_with_sampler(
                model=model,
                sampler_name=sampler_name,
                summary_vars=target_info["summary_vars"],
                sample_config=BUDGET_MCMC_CONFIGS[budget_label],
                seed=seed,
            )

            row = {
                "target": target_name,
                "target_short": target_info["target_short"],
                "target_family": target_info["target_family"],
                "replicate": replicate,
                "seed": seed,
                "budget": budget_label,
                "budget_index": BUDGET_LEVEL_ORDER.index(budget_label) + 1,
                "sampler": sampler_name,
                "sampler_family": sampler_specs[sampler_name]["family"],
                "uses_gradients": sampler_specs[sampler_name]["gradients"],
                "failed_run": False,
                "error_message": "",
                **metrics,
            }

        except Exception as err:
            row = failed_hard_budget_row(
                target_name=target_name,
                target_info=target_info,
                sampler_name=sampler_name,
                budget_label=budget_label,
                replicate=replicate,
                seed=seed,
                error_message=err,
            )

        hard_budget_rows.append(row)
        hard_budget_results_df = pd.DataFrame(hard_budget_rows)
        hard_budget_results_df.to_csv(HARD_BUDGET_RAW_PATH, index=False)

else:
    hard_budget_results_df = pd.read_csv(HARD_BUDGET_RAW_PATH)

hard_budget_results_df.head()

Scenario	Sampler	Budget	Mean runtime	Mean ESS/sec	Mean tail ESS/sec	Mean max r-hat	Reliability	Mean divergences	Mean BFMI
Hierarchical	Metropolis-Hastings	Short	233.69	1.2	1.2	1.26	0	0	NA
Hierarchical	DE-Metropolis	Short	198.59	0.1	0.3	3.495	0	0	NA
Hierarchical	Hamiltonian MC	Short	236.3	6.7	5.7	1.245	0	0	1.018
Hierarchical	NUTS	Short	47.09	24.9	13.6	1.015	1	0	1.045
Hierarchical	Metropolis-Hastings	Medium	48.87	5.5	5.3	1.125	0	0	NA
Hierarchical	DE-Metropolis	Medium	40.04	0.1	0.4	3.36	0	0	NA
Hierarchical	Hamiltonian MC	Medium	61.57	46	29.7	1	1	0	1.032
Hierarchical	NUTS	Medium	43.56	60.3	34.4	1.005	1	0	1.030
Hierarchical	Metropolis-Hastings	Long	46.46	12	11.7	1.035	1	0	NA
Hierarchical	DE-Metropolis	Long	28.56	0.2	0.5	3.15	0	0	NA
Hierarchical	Hamiltonian MC	Long	237.36	76	47.5	1	1	0	1.016
Hierarchical	NUTS	Long	235.93	102.4	44.4	1	1	0	1.033
Zero-inflated count	Metropolis-Hastings	Short	295.73	0.3	0.5	1.435	0	0	NA
Zero-inflated count	DE-Metropolis	Short	266.2	0.1	0.1	2.85	0	0	NA
Zero-inflated count	Hamiltonian MC	Short	317.61	2	3.7	1.05	0.5	0	1.014
Zero-inflated count	NUTS	Short	313.7	5.3	4.5	1.015	1	0	1.021
Zero-inflated count	Metropolis-Hastings	Medium	62.12	1.2	2.3	1.22	0	0	NA
Zero-inflated count	DE-Metropolis	Medium	50.65	0.1	0.4	2.66	0	0	NA
Zero-inflated count	Hamiltonian MC	Medium	71	17.1	19.8	1.01	1	0	1.035
Zero-inflated count	NUTS	Medium	308.28	12.2	11.2	1	1	0	0.986
Zero-inflated count	Metropolis-Hastings	Long	304.57	1.1	1.7	1.065	0.5	0	NA
Zero-inflated count	DE-Metropolis	Long	270.09	0	0.2	2.19	0	0	NA
Zero-inflated count	Hamiltonian MC	Long	317.91	12.5	15.9	1	1	0	1.017
Zero-inflated count	NUTS	Long	548.43	6.4	5.6	1	1	0	1.017
Funnel	Metropolis-Hastings	Short	36.34	1.3	0.6	1.975	0	0	NA
Funnel	DE-Metropolis	Short	32.53	0.3	1	1.7	0	0	NA
Funnel	Hamiltonian MC	Short	37.31	1.5	1	1.395	0	1	1.058
Funnel	NUTS	Short	38.85	27.6	0.9	1.595	0	1.5	0.121
Funnel	Metropolis-Hastings	Medium	181.18	1.8	0.7	1.605	0	0	NA
Funnel	DE-Metropolis	Medium	179.51	0.3	1.2	1.455	0	0	NA
Funnel	Hamiltonian MC	Medium	192.15	1.4	0.5	1.32	0	1.5	1.579
Funnel	NUTS	Medium	38.33	48.2	5.6	1.295	0	16	0.147
Funnel	Metropolis-Hastings	Long	35.49	2.2	0.9	1.365	0	0	NA
Funnel	DE-Metropolis	Long	29.72	3.1	4.4	1.255	0	0	NA
Funnel	Hamiltonian MC	Long	29.69	4.9	1.2	1.125	0	8.5	1.779
Funnel	NUTS	Long	36.06	84.4	9.5	1.51	0	584.5	0.347

The summary table shows that runtime differences are not always proportional to the nominal budget. This is likely because these models are modest in size and the fixed overhead of PyMC sampling, four-chain multiprocessing, and result collection is substantial. Therefore, runtime alone does not distinguish the samplers very well in this subsection. The more informative columns are ESS/sec, tail ESS/sec, max r-hat, and reliability.

The budget-response plots show a consistent efficiency pattern across the harder targets. As the budget increases from short to medium to long, the gradient-based samplers, especially NUTS and HMC, generally achieve higher ESS/sec and tail ESS/sec. This trend is clearest for the hierarchical and zero-inflated count models. In contrast, the random-walk samplers improve much more slowly: Metropolis-Hastings gains some efficiency with longer budgets, while DE-Metropolis remains close to zero in most settings.

The max r-hat plots show that efficiency alone is not sufficient. The dashed horizontal line marks the reliability threshold, \(\hat{R}=1.05\); values below this line indicate better between-chain agreement and are treated as passing the convergence diagnostic. On the hierarchical target, NUTS is reliable across all budgets, HMC becomes reliable from the medium budget onward, and Metropolis-Hastings only becomes reliable under the long budget. DE-Metropolis remains unreliable throughout. On the zero-inflated count target, HMC and NUTS are reliable across all budgets, while Metropolis-Hastings is only partially reliable at the long budget and DE-Metropolis again fails throughout.

The reliability heatmaps summarize the same diagnostic information as pass rates across the repeated seeds. A value of 1.00 means that all repeated runs for that sampler-budget combination passed the reliability criterion, 0.50 means that only half passed, and 0.00 means that none passed. This makes the sampler-budget pattern easier to read: HMC and NUTS are consistently reliable on the zero-inflated count target, NUTS is consistently reliable on the hierarchical target, and all samplers remain unreliable on Neal’s Funnel.

Neal’s Funnel behaves differently from the two applied targets. Although NUTS has the highest ESS/sec and tail ESS/sec on the funnel target, all samplers fail the reliability criterion across the full budget grid. The max r-hat values remain above the convergence threshold, so the high NUTS efficiency should not be interpreted as successful posterior exploration. Instead, the result suggests that funnel-like geometry is not solved by longer chains alone.

9 Hierarchical Parameterization Study

One important limitation of any sampler benchmark is that performance can be driven as much by parameterization as by the algorithm itself. To study that effect directly, we compare centered and non-centered versions of the hierarchical regression model under the two gradient-based samplers.

parameterization_builders = {
    "Centered": build_hierarchical_model,
    "Non-centered": build_noncentered_hierarchical_model,
}
parameterization_scenarios = [
    "Centered hierarchical regression",
    "Non-centered hierarchical regression",
]
parameterization_rows = []
parameterization_example_traces = {}

for parameterization_name, builder in parameterization_builders.items():
    model = builder()
    for replicate, seed in enumerate(PARAMETERIZATION_SEEDS, start=1):
        for sampler_name in PARAMETERIZATION_SAMPLERS:
            idata, metrics = run_model_with_sampler(
                model=model,
                sampler_name=sampler_name,
                summary_vars=["mu_a", "tau_a", "a", "beta", "sigma"],
                sample_config=REPEATED_CONFIG,
                seed=seed,
            )
            parameterization_rows.append(
                {
                    "scenario": f"{parameterization_name} hierarchical regression",
                    "scenario_short": parameterization_name,
                    "parameterization": parameterization_name,
                    "replicate": replicate,
                    "seed": seed,
                    "sampler": sampler_name,
                    **metrics,
                }
            )
            if replicate == 1:
                parameterization_example_traces[(parameterization_name, sampler_name)] = idata

parameterization_results_df = finalize_results_df(pd.DataFrame(parameterization_rows), parameterization_scenarios)
parameterization_summary_df = summarize_repeated_block(parameterization_results_df, parameterization_scenarios)

parameterization_diag_df = (
    parameterization_results_df.groupby(["scenario", "scenario_short", "sampler"], observed=True)
    .agg(
        bfmi_mean=("mean_bfmi", "mean"),
        bfmi_low=("mean_bfmi", lambda x: empirical_interval(x)[0]),
        bfmi_high=("mean_bfmi", lambda x: empirical_interval(x)[1]),
        tail_ess_mean=("tail_ess_per_second", "mean"),
        tail_ess_low=("tail_ess_per_second", lambda x: empirical_interval(x)[0]),
        tail_ess_high=("tail_ess_per_second", lambda x: empirical_interval(x)[1]),
        accept_mean=("mean_accept_stat", "mean"),
    )
    .reset_index()
)

parameterization_merged_df = parameterization_summary_df.merge(
    parameterization_diag_df,
    on=["scenario", "scenario_short", "sampler"],
    how="left",
)

parameterization_display_df = parameterization_merged_df.assign(
    **{
        "ESS/sec mean [95% interval]": [
            interval_string(row.ess_mean, row.ess_low, row.ess_high, digits=1)
            for row in parameterization_merged_df.itertuples()
        ],
        "Tail ESS/sec mean [95% interval]": [
            interval_string(row.tail_ess_mean, row.tail_ess_low, row.tail_ess_high, digits=1)
            for row in parameterization_merged_df.itertuples()
        ],
        "Mean BFMI [95% interval]": [
            interval_string(row.bfmi_mean, row.bfmi_low, row.bfmi_high, digits=3)
            for row in parameterization_merged_df.itertuples()
        ],
        "Mean accept stat": parameterization_merged_df["accept_mean"].round(3),
        "Mean max r-hat": parameterization_merged_df["rhat_mean"].round(3),
        "Reliability rate": parameterization_merged_df["reliability_rate"].round(2),
    }
)[
    [
        "scenario_short",
        "sampler",
        "ESS/sec mean [95% interval]",
        "Tail ESS/sec mean [95% interval]",
        "Mean BFMI [95% interval]",
        "Mean accept stat",
        "Mean max r-hat",
        "Reliability rate",
    ]
].rename(columns={"scenario_short": "Parameterization", "sampler": "Sampler"})

Parameterization	Sampler	ESS/sec mean [95% interval]	Tail ESS/sec mean [95% interval]	Mean BFMI [95% interval]	Mean accept stat	Mean max r-hat	Reliability rate
Centered	Hamiltonian MC	100.0 [72.0, 118.1]	71.0 [52.1, 85.1]	1.020 [1.004, 1.054]	1	1	1
Centered	NUTS	106.8 [17.8, 160.6]	50.2 [8.3, 75.1]	1.013 [0.991, 1.046]	1	1	1
Non-centered	Hamiltonian MC	147.4 [100.9, 203.3]	27.5 [20.3, 36.7]	0.737 [0.695, 0.769]	1	1	1
Non-centered	NUTS	121.7 [89.3, 169.4]	55.6 [42.3, 74.0]	0.746 [0.719, 0.778]	1	1	1

Repeated-run efficiency and BFMI under centered versus non-centered hierarchical parameterizations.

CSV-based geometry diagnostics for centered versus non-centered hierarchical parameterizations.

Parameterization effects are sampler-specific rather than one-directional in this model. For NUTS, the non-centered form is stronger on mean ESS/sec (121.7 versus 106.8), whereas Hamiltonian MC gains raw ESS/sec under the non-centered form (147.4 versus 100.0). Both parameterizations remain fully reliable here, which means geometry can matter even when standard convergence flags stay quiet.

This section makes the benchmark more informative for real modeling practice. In this example, the non-centered form is not universally better: it improves raw ESS/sec for both static HMC and NUTS, but the centered form keeps BFMI healthier overall. That is exactly why parameterization belongs in a serious sampler comparison instead of being treated as a nuisance detail.

10 Neal’s Funnel Geometry Stress Test

Applied benchmarks are useful, but they do not always expose the geometric pathologies that make one sampler fundamentally different from another. Neal’s Funnel is a classic test case for that reason.

The model is

\[ z \sim \mathcal{N}(0, 3^2), \qquad x_i \mid z \sim \mathcal{N}(0, \exp(z/2)) \]

When \(z\) is very negative, the funnel neck becomes extremely tight. Random-walk proposals can struggle badly, while dynamic HMC methods are designed to adapt much better.

funnel_model = build_funnel_model()
funnel_rows = []
funnel_example_traces = {}

for replicate, seed in enumerate(FUNNEL_SEEDS, start=1):
    for sampler_name in sampler_order:
        idata, metrics = run_model_with_sampler(
            model=funnel_model,
            sampler_name=sampler_name,
            summary_vars=["z", "x"],
            sample_config=FUNNEL_CONFIG,
            seed=seed,
        )
        funnel_rows.append(
            {
                "scenario": "Neal's Funnel",
                "scenario_short": "Funnel",
                "replicate": replicate,
                "seed": seed,
                "sampler": sampler_name,
                **metrics,
            }
        )
        if replicate == 1:
            funnel_example_traces[sampler_name] = idata

funnel_results_df = finalize_results_df(pd.DataFrame(funnel_rows))
funnel_summary_df = summarize_repeated_block(funnel_results_df, ["Neal's Funnel"])
funnel_display_df = funnel_summary_df.assign(
    **{
        "ESS/sec mean [95% interval]": [
            interval_string(row.ess_mean, row.ess_low, row.ess_high, digits=1)
            for row in funnel_summary_df.itertuples()
        ],
        "Runtime mean [95% interval]": [
            interval_string(row.runtime_mean, row.runtime_low, row.runtime_high, digits=2)
            for row in funnel_summary_df.itertuples()
        ],
        "Mean max r-hat": funnel_summary_df["rhat_mean"].round(3),
        "Reliability rate": funnel_summary_df["reliability_rate"].round(2),
        "Pr(best ESS/s)": funnel_summary_df["pr_best_ess"].round(2),
    }
)[
    [
        "sampler",
        "ESS/sec mean [95% interval]",
        "Runtime mean [95% interval]",
        "Mean max r-hat",
        "Reliability rate",
        "Pr(best ESS/s)",
    ]
].rename(columns={"sampler": "Sampler"})

funnel_pairwise_ess = pairwise_probability_matrix(
    funnel_results_df,
    metric="ess_per_second",
    higher_is_better=True,
)

funnel_hmc_diag_df = (
    funnel_results_df.loc[funnel_results_df["sampler"].isin(PARAMETERIZATION_SAMPLERS)]
    .groupby("sampler", observed=True)
    .agg(
        tail_ess_mean=("tail_ess_per_second", "mean"),
        tail_ess_low=("tail_ess_per_second", lambda x: empirical_interval(x)[0]),
        tail_ess_high=("tail_ess_per_second", lambda x: empirical_interval(x)[1]),
        bfmi_mean=("mean_bfmi", "mean"),
        bfmi_low=("mean_bfmi", lambda x: empirical_interval(x)[0]),
        bfmi_high=("mean_bfmi", lambda x: empirical_interval(x)[1]),
        accept_mean=("mean_accept_stat", "mean"),
        n_steps_mean=("mean_n_steps", "mean"),
        tree_depth_mean=("mean_tree_depth", "mean"),
        max_treedepth_rate_mean=("max_treedepth_rate", "mean"),
    )
    .reset_index()
)

funnel_hmc_display_df = funnel_hmc_diag_df.assign(
    **{
        "Tail ESS/sec mean [95% interval]": [
            interval_string(row.tail_ess_mean, row.tail_ess_low, row.tail_ess_high, digits=1)
            for row in funnel_hmc_diag_df.itertuples()
        ],
        "Mean BFMI [95% interval]": [
            interval_string(row.bfmi_mean, row.bfmi_low, row.bfmi_high, digits=3)
            for row in funnel_hmc_diag_df.itertuples()
        ],
        "Mean accept stat": funnel_hmc_diag_df["accept_mean"].round(3),
        "Mean n_steps": funnel_hmc_diag_df["n_steps_mean"].round(1),
        "Mean tree depth": [
            f"{value:.2f}" if np.isfinite(value) else "NA"
            for value in funnel_hmc_diag_df["tree_depth_mean"].to_numpy(dtype=float)
        ],
        "Pr(max treedepth)": [
            f"{value:.2f}" if np.isfinite(value) else "NA"
            for value in funnel_hmc_diag_df["max_treedepth_rate_mean"].to_numpy(dtype=float)
        ],
    }
)[
    [
        "sampler",
        "Tail ESS/sec mean [95% interval]",
        "Mean BFMI [95% interval]",
        "Mean accept stat",
        "Mean n_steps",
        "Mean tree depth",
        "Pr(max treedepth)",
    ]
].rename(columns={"sampler": "Sampler"})

Sampler	ESS/sec mean [95% interval]	Runtime mean [95% interval]	Mean max r-hat	Reliability rate	Pr(best ESS/s)
Metropolis-Hastings	3.8 [0.6, 11.2]	46.91 [23.49, 90.99]	1	0	0
DE-Metropolis	2.2 [0.5, 5.0]	45.74 [23.09, 75.94]	1	0	0
Hamiltonian MC	7.2 [0.9, 15.9]	54.96 [26.96, 108.15]	1	0	0
NUTS	68.8 [9.7, 145.3]	61.72 [27.54, 118.74]	1	0	1

Sampler	Tail ESS/sec mean [95% interval]	Mean BFMI [95% interval]	Mean accept stat	Mean n_steps	Mean tree depth	Pr(max treedepth)
Hamiltonian MC	2.8 [0.7, 5.7]	1.748 [1.062, 1.955]	1	14	NA	NA
NUTS	8.6 [3.3, 19.3]	0.103 [0.068, 0.157]	1	50	4.18	0.00

Neal’s Funnel stress diagnostics by sampler.

Repeated-seed diagnostics for Neal’s Funnel.

Richer HMC-family diagnostics on Neal’s Funnel.

On Neal’s Funnel, NUTS is the clear efficiency leader (68.8 ESS/sec versus 7.2 for Hamiltonian MC), but the richer diagnostics show that the geometry is still hostile: its mean BFMI is only 0.103 and its mean tree depth is 4.18. Thus, the funnel result is not merely a ranking table; it also shows why even the best-performing method remains under geometric stress.

Neal’s Funnel is where the geometric differences between samplers become impossible to ignore. In the applied battery, one might still be tempted to focus on runtime.

At the same time, the funnel results are a reminder that this is still a finite-budget benchmark. Even with four chains, 500 tuning iterations, and 1000 retained draws per chain, Neal’s Funnel remains demanding under the strict reliability threshold used here, while NUTS is clearly the strongest performer by repeated-run ESS/sec. The added BFMI, tail ESS, and tree-depth diagnostics also make the geometry story more mechanistic than a simple speed comparison.

11 Posterior Predictive Checks

Sections 4 through 12 report on sampler efficiency and accuracy. This section and the next add a second dimension to the benchmark—inferential reliability and calibration—by asking not “which sampler is fastest” but “are the posteriors that the winning samplers produce actually trustworthy?” Posterior predictive checks (PPCs) address this question at the model–data level. Conditional on apparently clean sampler diagnostics, such as R-hat and ESS, we ask whether the fitted posterior predictive distribution can reproduce the observed data. If not, sampler rankings may be less meaningful because the benchmark would be comparing algorithms on poorly fitting targets. We draw replicate observations from each fitted posterior and compare them to the real data; if the replicates envelop the observations, the model is at least reproducing the features that the likelihood tries to capture.

Concretely, the posterior predictive density is

\[ p(y_{\text{rep}} \mid y) \;=\; \int p(y_{\text{rep}} \mid \theta)\, p(\theta \mid y)\, d\theta, \]

so a single posterior predictive draw is obtained by first sampling \(\theta \sim p(\theta \mid y)\) from the fitted posterior and then sampling \(y_{\text{rep}} \sim p(y_{\text{rep}} \mid \theta)\) from the likelihood. Repeating this many times yields a sample from \(p(y_{\text{rep}} \mid y)\), which is what pm.sample_posterior_predictive produces and what az.plot_ppc overlays against the observed \(y\). The graphical PPC used here has a quantitative counterpart, the Bayesian posterior predictive p-value of Gelman, Meng & Stern (1996), \(p_B = \Pr\bigl(T(y_{\text{rep}}, \theta) \geq T(y_{\text{obs}}, \theta) \,\big|\, y\bigr)\) for a chosen test statistic \(T\). We retain the graphical version because the question of interest at this stage is whether the whole shape of the response distribution is covered by the replicates—shape mismatches that no single scalar \(T\) would necessarily flag.

We run PPC on all five applied scenarios using the single-seed NUTS posteriors obtained in Section 4.

ppc_idata_by_scenario = {}

for scenario_name in scenario_order:
    scenario = scenarios[scenario_name]
    model = portfolio_models[scenario_name]
    idata = portfolio_idatas[scenario_name]["NUTS"]
    with model:
        ppc = pm.sample_posterior_predictive(
            idata,
            random_seed=PORTFOLIO_SEED,
            progressbar=False,
        )
    idata.extend(ppc)
    ppc_idata_by_scenario[scenario_name] = idata

Posterior predictive checks for the five applied scenarios using NUTS. Points are observed summaries; intervals are posterior predictive 95% ranges.

The PPC panels overlay observed summaries on draws simulated from each fitted model. The hierarchical, logistic, and zero-inflated scenarios show observed summaries that fall well within the posterior predictive envelopes, suggesting that these fits reproduce the main location and dispersion features of the data. The Gaussian and linear regression panels show mild shape mismatch, consistent with observed features that a unimodal Gaussian likelihood may not fully reproduce. However, the observed summaries are still largely contained within the predictive envelopes. This supports the use of these fitted targets for sampler comparison, while also showing why PPCs are useful: they can reveal model-form limitations that R-hat and ESS alone cannot detect.

12 Simulation-Based Calibration

PPC addressed the model–data half of inferential reliability. SBC addresses the sampler half: it asks whether the posterior samples produced by a given sampler are statistically consistent with the true posterior that the sampler is supposed to target, independent of any particular observed dataset. A biased sampler can pass both R-hat diagnostics and PPCs while still producing posterior samples that systematically misplace the true parameter—an error PPC cannot see because PPC never compares the posterior against a known ground truth, and an error r-hat cannot see because r-hat only measures agreement between chains of the same sampler rather than agreement with the true posterior. SBC addresses this blind spot by simulating ground truth from the prior itself and measuring where the true parameter lands within the sampler’s posterior samples. The algorithm is:

1. Draw a true parameter \(\theta^{*}\) from the prior.
2. Simulate data \(y\) from the likelihood given \(\theta^{*}\).
3. Fit the model on \(y\) and record posterior samples \(\theta_{1}, \ldots, \theta_{M}\).
4. Compute the rank of \(\theta^{*}\) within those samples.

Under a correctly specified pipeline those ranks are uniform on \([0, 1]\). Deviations diagnose specific failure modes:

• U-shaped histogram: the posterior is too narrow (the truth lands in the tails too often).
• ∩-shaped histogram: the posterior is too wide.
• Skewed histogram: the posterior mean is biased.

We run SBC on two of the applied scenarios and on all four samplers in the portfolio (Metropolis-Hastings, DE-Metropolis, Hamiltonian MC, NUTS), producing a \(4 \times 2\) design in which each cell is checked by thirty independent rounds. The Gaussian location model is the simplest continuous target and serves as a baseline that every reasonable sampler should pass; the hierarchical regression is the most structurally sensitive applied target and is where we expect sampler-level differences to surface first. Two scenarios (rather than all five applied scenarios) keeps the runtime bounded at \(4 \times 2 \times 30 = 240\) posterior fits while still covering the easy–hard extremes of the portfolio, which is the same design choice used in the original SBC literature (Talts et al., 2018) where calibration is demonstrated on a small number of deliberately chosen representative targets rather than the full model catalog. For each round we use a reduced sampling budget of two chains, 200 tuning iterations, and 500 retained draws per chain.

SBC_CONFIG = {
    "draws": 500,
    "tune": 200,
    "chains": SBC_CHAINS,
    "cores": SBC_CORES,
    "progressbar": False,
    "return_inferencedata": True,
    "compute_convergence_checks": False,
    "idata_kwargs": {"log_likelihood": False},
}
SBC_N_ROUNDS = 30
SBC_SCENARIOS = ["Gaussian location model", "Hierarchical regression"]
SBC_SAMPLERS = sampler_order
SBC_PARAMS = {
    "Gaussian location model": ["mu", "sigma"],
    "Hierarchical regression": ["mu_a", "tau_a", "beta", "sigma"],
}
SBC_SEEDS = {
    "Gaussian location model": [RANDOM_SEED + 3000 + i for i in range(SBC_N_ROUNDS)],
    "Hierarchical regression": [RANDOM_SEED + 3500 + i for i in range(SBC_N_ROUNDS)],
}


def build_gaussian_for_sbc(n_obs):
    coords = {"obs_id": np.arange(n_obs)}
    with pm.Model(coords=coords) as model:
        y_obs = pm.Data("y_obs", np.zeros(n_obs), dims="obs_id")
        mu = pm.Normal("mu", mu=0.0, sigma=1.5)
        sigma = pm.HalfNormal("sigma", sigma=1.0)
        pm.Normal("obs", mu=mu, sigma=sigma, observed=y_obs, dims="obs_id")
    return model


def build_hier_for_sbc(x_vals, group_idx_vals, group_labels_local):
    n_obs = len(x_vals)
    coords = {"obs_id": np.arange(n_obs), "group": group_labels_local}
    with pm.Model(coords=coords) as model:
        x = pm.Data("x", np.asarray(x_vals, dtype=float), dims="obs_id")
        group_idx = pm.Data("group_idx", np.asarray(group_idx_vals, dtype="int64"), dims="obs_id")
        y_obs = pm.Data("y_obs", np.zeros(n_obs), dims="obs_id")
        mu_a = pm.Normal("mu_a", mu=0.0, sigma=1.5)
        tau_a = pm.HalfNormal("tau_a", sigma=1.0)
        a = pm.Normal("a", mu=mu_a, sigma=tau_a, dims="group")
        beta = pm.Normal("beta", mu=0.0, sigma=1.0)
        sigma = pm.HalfNormal("sigma", sigma=1.0)
        mu = a[group_idx] + beta * x
        pm.Normal("obs", mu=mu, sigma=sigma, observed=y_obs, dims="obs_id")
    return model


sbc_models = {
    "Gaussian location model": build_gaussian_for_sbc(gaussian_y.shape[0]),
    "Hierarchical regression": build_hier_for_sbc(x_hier, hier_group_idx, hier_group_labels),
}


def sbc_one_round(scenario_name, seed, sampler_name="NUTS"):
    model = sbc_models[scenario_name]
    with model:
        prior = pm.sample_prior_predictive(draws=1, random_seed=seed)
    y_sim = np.asarray(prior.prior_predictive["obs"]).reshape(-1)
    theta_star = {
        name: float(np.asarray(prior.prior[name]).ravel()[0])
        for name in SBC_PARAMS[scenario_name]
    }
    with model:
        pm.set_data({"y_obs": y_sim})
        idata = pm.sample(
            step=sampler_specs[sampler_name]["factory"](),
            random_seed=seed,
            **SBC_CONFIG,
        )
    ranks = {}
    for param_name, theta in theta_star.items():
        posterior_vals = np.asarray(idata.posterior[param_name]).reshape(-1)
        ranks[param_name] = float((posterior_vals < theta).mean())
    return ranks

sbc_rows = []
for scenario_name in SBC_SCENARIOS:
    for sampler_name in SBC_SAMPLERS:
        for round_idx, seed in enumerate(SBC_SEEDS[scenario_name]):
            ranks = sbc_one_round(scenario_name, seed, sampler_name=sampler_name)
            for param_name, rank_val in ranks.items():
                sbc_rows.append(
                    {
                        "scenario": scenario_name,
                        "scenario_short": scenarios[scenario_name].short_name,
                        "sampler": sampler_name,
                        "param": param_name,
                        "round": round_idx,
                        "seed": seed,
                        "rank": rank_val,
                    }
                )

sbc_results_df = pd.DataFrame(sbc_rows)

Each rank is the empirical quantile of the prior-drawn true value \(\theta^{*}\) within the \(M = 2 \times 500 = 1{,}000\) posterior samples produced by a given sampler from the simulated data. Under a calibrated sampler the collection of ranks across independent rounds is Uniform\((0, 1)\), so each bin of the histogram should contain approximately \(N_{\text{rounds}} / N_{\text{bins}} = 30 / 10 = 3\) observations; this expected count is drawn as a red dashed line for reference. The uniform mean of \(0.5\) and the uniform standard deviation of \(\sqrt{1/12} \approx 0.289\) are the two numerical anchors used in the summary table at the end of the section. Each figure below shows one scenario; rows are samplers and columns are the parameters for which SBC ranks are tracked, so reading a single column top-to-bottom compares the four samplers on the same parameter, while reading a single row left-to-right checks whether one sampler is calibrated across all parameters of that scenario.

SBC rank histograms faceted by sampler (rows) and parameter (columns) for each scenario. Uniform distribution around the red dashed line indicates a calibrated sampler; systematic deviation indicates a miscalibrated posterior that r-hat and ESS cannot detect.

The two figures above expose a clear separation between scenarios and between samplers. In the Gaussian scenario all four samplers produce histograms that hover around the red uniform line with only the bin-to-bin noise expected at thirty rounds; no row shows a systematic U-shape, \(\cap\)-shape, or one-sided skew. In the Hierarchical scenario three of the four samplers (Metropolis-Hastings, Hamiltonian MC, and NUTS) remain close to uniform, but the DE-Metropolis row is visibly different: mu_a and tau_a show pronounced piles at both \(0\) and \(1\) with a hollow middle, beta shows a milder version of the same pattern, and sigma collapses into a single tall spike in the leftmost bin with essentially nothing to the right of \(0.1\). These shapes are the diagnostic signatures listed earlier (U-shape and extreme left skew) and they are the first concrete, sampler-discriminating finding the SBC section is designed to produce. The numerical summary below converts these visual impressions into mean and standard-deviation deviations from the Uniform\((0,1)\) anchors of \(0.5\) and \(0.289\).

In addition to the visual histograms, we report a Kolmogorov–Smirnov test of each rank set against the Uniform\((0,1)\) reference. KS measures the largest vertical gap between the empirical CDF of the 30 ranks and the theoretical Uniform CDF \(F(x) = x\), and returns a p-value for the null hypothesis “ranks are Uniform\((0,1)\).” A p-value below conventional thresholds (e.g. \(0.05\)) constitutes a quantitative rejection of calibration that complements the visual reading of the histograms.

Scenario	Sampler	Parameter	Mean rank	SD of rank	Min rank	Max rank	KS p-value
Gaussian	Metropolis-Hastings	mu	0.532	0.264	0.014	0.983	0.291
Gaussian	Metropolis-Hastings	sigma	0.393	0.3	0.001	0.951	0.0741
Gaussian	DE-Metropolis	mu	0.511	0.243	0.026	0.912	0.3852
Gaussian	DE-Metropolis	sigma	0.385	0.292	0.029	0.998	0.0617
Gaussian	Hamiltonian MC	mu	0.515	0.251	0.019	0.943	0.4483
Gaussian	Hamiltonian MC	sigma	0.408	0.308	0.015	0.996	0.1477
Gaussian	NUTS	mu	0.511	0.255	0.017	0.946	0.6254
Gaussian	NUTS	sigma	0.409	0.304	0.005	0.994	0.1715
Hierarchical	Metropolis-Hastings	mu_a	0.509	0.295	0.041	1	0.7719
Hierarchical	Metropolis-Hastings	tau_a	0.415	0.278	0	1	0.2747
Hierarchical	Metropolis-Hastings	beta	0.552	0.33	0	1	0.5431
Hierarchical	Metropolis-Hastings	sigma	0.454	0.336	0.029	1	0.2078
Hierarchical	DE-Metropolis	mu_a	0.419	0.46	0	1	0
Hierarchical	DE-Metropolis	tau_a	0.434	0.435	0	1	0.0001
Hierarchical	DE-Metropolis	beta	0.405	0.339	0	1	0.1381
Hierarchical	DE-Metropolis	sigma	0.018	0.081	0	0.44	0
Hierarchical	Hamiltonian MC	mu_a	0.5	0.299	0.041	0.986	0.953
Hierarchical	Hamiltonian MC	tau_a	0.423	0.272	0.036	1	0.2404
Hierarchical	Hamiltonian MC	beta	0.56	0.328	0.001	0.991	0.3212
Hierarchical	Hamiltonian MC	sigma	0.456	0.325	0.01	0.999	0.3583
Hierarchical	NUTS	mu_a	0.504	0.292	0.039	0.968	0.9003
Hierarchical	NUTS	tau_a	0.424	0.278	0.032	1	0.2629
Hierarchical	NUTS	beta	0.556	0.33	0.002	0.993	0.3607
Hierarchical	NUTS	sigma	0.456	0.328	0.018	1	0.2459

12.1 Gaussian scenario: every sampler is calibrated

On the Gaussian location model, all four samplers produce mean ranks inside the narrow band \([0.385, 0.532]\) for both parameters and standard deviations between \(0.243\) and \(0.306\). Each of these values is a plausible finite-sample draw from Uniform\((0,1)\) at \(N_{\text{rounds}} = 30\), and the histograms show no systematic shape. The KS test confirms this reading quantitatively: all eight Gaussian rows return p-values above conventional rejection thresholds (lowest is \(0.062\) for DE-Metropolis on sigma), so the Uniform\((0,1)\) null is not rejected for any (sampler, parameter) cell. The small common under-shoot of the sigma mean rank (all four samplers sit at approximately \(0.39\)–\(0.41\) rather than exactly \(0.5\)) is a shared feature across samplers rather than a sampler-specific defect, so the substantive conclusion for the Gaussian scenario is that all four samplers in the portfolio deliver calibrated posteriors on this target. SBC therefore does not discriminate between samplers on the easy baseline, which is the expected outcome: on a simple unimodal Gaussian target, calibration problems are less likely, and the main distinction between samplers is efficiency rather than reliability.

12.2 Hierarchical scenario: DE-Metropolis fails SBC

The hierarchical regression tells a very different story and shows why SBC is especially informative. The Metropolis-Hastings, Hamiltonian MC, and NUTS rows of the summary table are indistinguishable from calibrated pipelines: every mean rank lies in \([0.415, 0.556]\) and every standard deviation lies in \([0.269, 0.336]\), all within sampling noise of the Uniform\((0,1)\) anchors, and every KS p-value sits above \(0.20\) (twelve cells in total, lowest \(0.208\)). These three samplers pass SBC on the hardest target in the SBC design—visually and at the formal \(\alpha = 0.05\) level—which confirms that for the hierarchical regression in this benchmark, r-hat, ESS, and SBC all agree that the posteriors being reported are trustworthy in distribution.

DE-Metropolis is the exception, with clear evidence of calibration failure:

• sigma collapses to the leftmost bin. The mean rank is \(0.018\) and the standard deviation is \(0.081\), against uniform anchors of \(0.5\) and \(0.289\), and the KS p-value rounds to \(0.0000\) (well below \(10^{-10}\)). The histogram has roughly twenty-nine of the thirty rounds concentrated in the first decile and essentially zero counts anywhere above rank \(0.1\). The interpretation is unambiguous: in nearly every round, the true \(\sigma\) drawn from the prior ends up below practically every posterior draw DE-Metropolis produced. DE-Metropolis is systematically returning posteriors for sigma that are shifted far above the truth.
• mu_a and tau_a produce classic U-shaped histograms. The mean ranks are close to \(0.5\) (\(0.419\) and \(0.434\)), but the standard deviations—\(0.460\) and \(0.435\)—are roughly \(50\%\) to \(60\%\) larger than the Uniform\((0,1)\) benchmark of \(0.289\), and the KS test rejects Uniform\((0,1)\) at \(p \approx 0.0000\) for mu_a and \(p = 0.0001\) for tau_a. The U-shape plus inflated SD is the textbook signature that the posterior is too narrow: the true parameter lands in the extreme tails of the sampler’s posterior far more often than a uniform rank distribution allows. DE-Metropolis is therefore under-covering the uncertainty on these two parameters.
• beta shows a milder version of the same pathology. The mean rank of \(0.405\) and SD of \(0.339\) sit closer to uniformity, and the KS p-value of \(0.138\) does not formally reject the Uniform\((0,1)\) null at \(\alpha = 0.05\). The histogram still has visibly heavier tails than middle bins, so the visual signature of the U-shape remains present, but the failure on beta is mild enough that the formal test cannot distinguish it from finite-sample noise at thirty rounds—precisely the kind of borderline call where reporting both the visual diagnostic and the KS p-value is more informative than either alone.

This matters because DE-Metropolis is often used as a relatively inexpensive gradient-free sampler for correlated posteriors, and on speed-only metrics it can look competitive. SBC reveals that this appearance is misleading: on the Hierarchical scenario, DE-Metropolis is producing posteriors whose spread and location do not match the true posterior, and the failure mode (mean rank of \(0.018\) on sigma, standard deviations of \(0.460\) and \(0.435\) on mu_a and tau_a) is the signature of structural bias rather than Monte-Carlo noise that would wash out with more draws. Any downstream inferential claim—credible intervals, posterior predictive summaries, decision-theoretic quantities—that is computed from a DE-Metropolis hierarchical posterior under the budget used here is therefore suspect, and this suspicion is raised specifically by SBC: it would not be raised by ESS/sec, by PPC, or by inter-chain r-hat measured within a DE-Metropolis run.

12.3 Consequences for sampler selection

Combining SBC with the efficiency rankings of the earlier sections produces a sharper recommendation than either measurement could make alone:

• Gaussian-style targets. All four samplers pass SBC, so the choice collapses to efficiency. The ESS/sec and runtime comparisons in Sections 4–6 are directly actionable and the fastest sampler is also a trusted sampler.
• Hierarchical-style targets. DE-Metropolis is ruled out on calibration grounds regardless of its ESS/sec number, because its posteriors are systematically biased in both location (sigma) and spread (mu_a, tau_a). The choice reduces to Metropolis-Hastings, Hamiltonian MC, and NUTS; among these, the efficiency rankings in Sections 4–7 can now be read as trustworthy efficiency rankings rather than merely as “who passes r-hat fastest.” In the benchmark runs of this report, that leaves NUTS and Hamiltonian MC as the dominant recommendation for hierarchical targets.

13 Recommendations and Conclusions

The benchmark supports stronger conclusions because it combines broad applied comparisons with uncertainty quantification, exact-posterior error analysis, approximate-inference comparisons, budget-response curves, parameterization checks, and geometry stress testing.

Use case	Recommended sampler	Reason
Fast exploratory baseline	Metropolis-Hastings or DE-Metropolis	Useful for smoke tests and quick iteration, but not the most trustworthy final answer.
Best average single-run portfolio performance	Hamiltonian MC	Lowest average composite rank in the broad applied battery.
Most reliable repeated-run default	Hamiltonian MC and NUTS	Tied best average reliability frequency across the repeated-seed applied study.
Fast approximate posterior screening	Mean-field ADVI	Best speed-versus-accuracy tradeoff within the variational methods on the exact Gaussian reference target.
Best scaling with dimension	NUTS	Smallest empirical scaling exponent alpha in the correlated Gaussian study, though the fitted relationship is noisier than for the other samplers.
Strongest finite-budget performer on pathological geometry	NUTS	Best repeated-run efficiency on Neal’s Funnel under the 4-chain, 500-tune, 1000-draw-per-chain benchmark setting, though not fully reliable by the convergence criterion.

The report supports eight main conclusions:

1. Single-run benchmarks are informative, but they can overstate certainty. Repeated-seed benchmarking reveals whether a nominal winner actually wins consistently.
2. Random-walk methods remain useful as inexpensive exploratory baselines, but they do not scale well with dimension or pathological geometry.
3. Static HMC performs extremely well on many medium-complexity problems and can win the average-rank competition in broad applied benchmarks.
4. Variational methods can add a useful fast-screening baseline, but the exact-Gaussian comparison shows that excellent posterior-mean accuracy does not automatically translate into good uncertainty recovery.
5. More compute does not help every method in the same way. In the budget-response study, NUTS improves steadily as the budget grows, whereas the ADVI methods remain much weaker on posterior variance even after longer optimization.
6. Parameterization is part of the story, not a side note. In the hierarchical study, centered and non-centered forms change the within-family ranking, so some apparent sampler differences are really geometry differences.
7. NUTS is usually the safest all-around default once we care about reliability under repeated runs and difficult posterior geometry.
8. The PPC and SBC sections further check whether good sampler diagnostics correspond to trustworthy inference.

That last point is important. There is no universally dominant sampler in every metric. What the report does show is something more useful: which samplers remain credible when the problem becomes harder, the benchmark itself becomes uncertain, and posterior geometry stops being friendly.

14 Limitations and Future Work

This report should be interpreted as a reproducible computational study conditional on the benchmark design implemented here, rather than as a universal ranking of samplers. Several limitations matter for interpretation.

14.1 Methodological limitations

1. Every conclusion in the report is conditional on the shared sampling budget of 4 chains, 500 tuning iterations, and 1000 retained draws per chain. A sampler that looks weak at this budget may improve materially at longer horizons, while a sampler with favorable ESS/sec may lose some advantage once startup and adaptation costs are amortized over much longer runs.
2. The benchmark includes a targeted centered-versus-non-centered hierarchical study, but it still does not exhaust the role of parameterization. The hierarchical and funnel sections remain sensitive to prior scaling, transforms, and adaptation behavior, so the comparison is best understood as sampler-plus-implementation rather than abstract algorithm alone.
3. ESS/sec and maximum r-hat are supplemented with exact-posterior RMSE, interval coverage, tail ESS, BFMI, tree-depth diagnostics, and a small variational comparison, but even that collection is incomplete. The report still does not directly measure predictive performance, decision-relevant utility, or error against gold-standard posteriors on every applied model.
4. The target portfolio is broad but not exhaustive. Five applied models, one correlated-Gaussian scaling family, and Neal’s Funnel cover several important regimes, but they do not yet represent multimodality, severe misspecification, discrete latent-state structure, heavy-tailed shrinkage priors, or models whose gradients are especially expensive.
5. Runtime results are hardware-, backend-, and implementation-dependent. The reported wall-clock comparisons reflect this PyMC stack, this machine, and these parallel settings; rankings can shift under JAX-backed inference, different linear algebra libraries, GPU acceleration, or alternative compilation and threading environments.
6. The exact-Gaussian and variational blocks partially address ground-truth accuracy, but only on a target with known analytic structure. Without long reference chains treated as near-gold-standard posteriors on the harder applied models, the report still cannot quantify full posterior error everywhere it benchmarks efficiency or approximation quality.

Taken together, these caveats suggest a more precise interpretation of the benchmark: it identifies which samplers remain credible under this reproducible workflow, not which sampler is universally optimal in all Bayesian computation settings.

14.2 High-value extensions

1. Extend the budget-response analysis to additional hard targets and longer budgets, since this report already includes hierarchical regression, zero-inflated count regression, and Neal’s Funnel.
2. Extend the exact-posterior error analysis beyond the Gaussian target. For a subset of harder applied models, very long reference chains could be used to measure estimation error in posterior means, variances, tail quantiles, and interval coverage.
3. Expand the parameterization study beyond the hierarchical regression example. A stronger version would include centered/non-centered variants for additional multilevel and latent-variable models, especially ones with stronger funnel geometry.
4. Expand the target family to include multimodal posteriors, sparse high-dimensional shrinkage models, latent-variable time-series models, and deliberately misspecified likelihoods. That would improve external validity and expose failure modes not visible in the current portfolio.
5. Broaden the approximate-inference comparison beyond ADVI variants. Methods such as Laplace approximations, Pathfinder-style initialization, stochastic-gradient MCMC, and JAX-backed NUTS would allow the report to study richer speed-quality frontiers rather than only exact-sampler rankings plus a small variational baseline.
6. Link the empirical scaling results more explicitly to theory. A future version of the project could connect the observed scaling exponents to optimal-scaling results for random-walk MH and HMC, spectral-gap intuition, and curvature-based explanations for why geometry-aware methods dominate on difficult targets.

Taken together, the project answers the original proposal at a substantially higher level. It moves from “which sampler looked better on one run?” to “which samplers remain strong under uncertainty, under increasing dimension, and under difficult posterior geometry?”