Methods · how every number on this site is made

Methods, decisions and limits

Where the data come from, what the model assumes, how each of the four methods and the 2026 checks compute their numbers, how those numbers are verified, and where all of it stops being trustworthy.

Data provenance

The data are the table given in the assignment: seven doses of an experimental influenza medication (0 to 6), ten patients at each dose, and how many of them had improved symptoms a week later. The counts are 3, 5, 4, 4, 7, 6 and 8, so 37 of 70 patients improved. The table holds counts only, with no personal information. The study design behind it (randomisation, blinding, how improvement was judged) is not described, so the site treats the data as a teaching example.

The original analysis is the R Markdown file in coursework/, kept byte for byte as submitted. Two scripts derive reference values from it without editing it. scripts/r-reference.R runs the original functions in both versions of the analysis, and scripts/diagnostics-reference.R runs them again with four seeds and computes R-hat, ESS and MCSE with the posterior package, plus quadrature checks for the prior sensitivity analysis.

Model

For dose $d = 0, \dots, 6$ with $y_d$ improvers out of 10 patients, the assignment specifies

y_d \sim \operatorname{Bin}(10, p_d), \qquad \operatorname{logit} p_d = \beta_0 + \beta_1 d, \qquad p(\beta_0, \beta_1) \propto 1.

\begin{gathered} y_d \sim \operatorname{Bin}(10, p_d) \\ \operatorname{logit} p_d = \beta_0 + \beta_1 d \\ p(\beta_0, \beta_1) \propto 1 \end{gathered}

The posterior has two dimensions, and it is proper despite the flat prior because every dose has both improvers and non-improvers. The site shows two versions. “Corrected” is this model. “As submitted (2023)” reproduces what the 2023 code actually fitted, a different and mis-encoded model for each method (errata, DR-001).

Methods

The four roads (2023, unchanged). The Laplace approximation is R’s glm(), a Gaussian at the posterior mode with the inverse Fisher information as its covariance. Metropolis-Hastings proposes from $N(\beta^{(t)}, c^2 \hat\Sigma_{\text{GLM}})$ with c = 2.4/√2 and runs 100,000 iterations, discarding 50,000. HMC takes L = 20 leapfrog steps of ε = 1/L with an identity mass matrix for 10,000 iterations, discarding 1,000. Expectation propagation uses one Gaussian site per observation and stops when the natural parameters change by less than 0.000001 (relative). All four are line-by-line ports of the submitted R functions, including R’s random number generator, so the browser reproduces R’s draws.

The exact posterior (2026). Because the posterior is two-dimensional, it is also computed on a 401 × 1,601 grid centred on the mode (±9 Laplace standard deviations) with the trapezoidal rule. Its moments agree with nested adaptive quadrature in R to about 10⁻¹¹, and its quantiles and P(β₁ > 0) to about 10⁻⁵. This is the yardstick for the method comparison.

Convergence diagnostics (2026). Each sampler runs as four chains with seeds 90125, 90126, 90127, 90128, each from its own rnorm(p) start. Chain 1 is exactly the 2023 chain. R-hat is the rank-normalised, folded, split R-hat of Vehtari, Gelman, Simpson, Carpenter and Bürkner (2021). Bulk ESS, tail ESS (the smaller of the ESS for the 5% and 95% quantiles) and MCSEs for means and quantiles follow the same paper. The code is a line-by-line port of the R package posterior 1.7.0, including Geyer’s initial monotone sequence and the cap that lets ESS exceed the number of draws for antithetic chains.

Posterior predictive checks (2026). For every post-warm-up HMC draw a replicated trial (10 patients per dose) is simulated with R’s generator. Four count statistics use the mid-p value Pr(T_rep > T_obs) + ½ Pr(T_rep = T_obs), and the chi-square discrepancy uses Pr(T(y_rep, θ) ≥ T(y, θ)).

Prior sensitivity (2026). The exact posterior is recomputed under independent $N(0, s^2)$ priors on both coefficients for s from 0.25 to 10, and compared with the flat prior by the shift in posterior means in units of the flat-prior posterior sd.

Method comparison (2026). Each method’s posterior is compared with the four-chain HMC reference and with the exact posterior. The comparison uses the shift in the mean (in reference sd), the 1-Wasserstein distance between marginals, and the KL divergence. For the Gaussian approximations the KL divergence from the exact posterior is computed by quadrature.

Evaluation design

Uncertainty on the Monte Carlo numbers. Every Monte Carlo mean is shown with its MCSE (posterior standard deviation over the square root of the ESS), and digits beyond the MCSE are not printed. Quantiles carry their own MCSE. Posterior predictive p-values are means of autocorrelated indicator series, so their MCSE uses the ESS of the indicator across chains. The chi-square indicator is binary and gets a Wilson interval with that effective n. The count statistics use a mid-p indicator (0, ½ or 1), which is not Bernoulli, so their MCSE is the indicator’s own sd over the square root of its ESS and their interval is p ± 1.96 MCSE. The HMC tuning comparisons are rerun on ten independent seed sets and reported as a median and range. The exact-posterior quantities carry quadrature error only.

Thresholds. A coefficient is flagged when R-hat exceeds 1.01 or either ESS is below 400, the thresholds recommended by Vehtari et al. (2021). They are conventions, not tests, and the page shows the raw values next to every flag.

Comparisons and effect sizes. Methods are compared on the same data, and the two samplers on the same four seeds. Differences are reported as effect sizes (shift in posterior sd, Wasserstein distance in the units of the coefficient, KL divergence), never as a p-value. A distance is read against the noise floor, which is the distance between two independent halves of the reference. The triangle inequality bounds how far a distance to the HMC reference can differ from the distance to the exact posterior.

Verification. The test suite compares the TypeScript against established implementations, never against itself. R-hat, bulk and tail ESS, ESS for quantiles and both MCSEs are compared with posterior on synthetic AR(1) chains with positive and negative correlation, ties and odd lengths, and on all four four-chain runs. Wilson intervals are compared with prop.test, qbeta with R and scipy, the Wasserstein distance with scipy.stats.wasserstein_distance, the Gaussian KL with numerical integration in scipy, and the grid posterior with nested integrate() in R. The numbers quoted in the decision records and the model card are pinned by tests, so the written record cannot drift from the code.

Seeds. Chains use seeds 90125, 90126, 90127, 90128, and predictive replications use seed 90125. All are R’s Mersenne-Twister, so every chain can be reproduced in R with set.seed() and the original functions.

Assumptions

Patients respond independently, and patients at the same dose share one probability of improving.
The log-odds of improvement is linear in dose over 0 to 6.
The prior is flat on both coefficients, as the assignment specified. The prior-sensitivity analysis checks how much this matters.
The counts are taken as given. Nothing is known about how patients were assigned to doses, so the slope describes an association in this table rather than a causal effect of dose.
The four-chain runs reuse the original settings. Chains start from rnorm(p), which is wide relative to this posterior but is not a designed overdispersed start.

Limitations

Seventy patients in seven groups give a wide interval for the slope (about 0.06 to 0.58 on the log-odds scale), and the posterior predictive checks have little power to detect a wrong curve shape.
Both Gaussian approximations miss the right skew of the slope’s posterior, so their interval endpoints sit up to 0.027 (Laplace) and 0.009 (EP) below the exact ones.
The full HMC tuning sweep uses one set of four seeds per setting, so most of its efficiency figures carry Monte Carlo noise of their own (the per-chain spread is shown in the table). Only the four settings DR-002 compares are rerun on ten seed sets, and their ranges show the spread, not a confidence interval.
The prior sensitivity analysis varies one family (independent normals with a common scale). It does not cover heavier-tailed priors or separate scales for the intercept and slope.
Bit-for-bit parity with R is established for the original settings and seeds. Other settings in the explorer run the same code but have no R reference.
The convergence thresholds are conventions. Passing them shows the samplers agree with each other, not that the model is right. The 2023 HMC run passes them while fitting the wrong data.

What I'd change

Build one data object and design matrix and pass it to every method (DR-001).
Test every method first on simulated data with known coefficients.
Run four or more chains and check R-hat, ESS and a posterior predictive check before reading any estimate (DR-001).
Centre dose, decouple ε from L, and use the Laplace covariance as the HMC mass matrix. On this posterior that is about 7 times as efficient as the original settings (6.4 to 7.7 times across ten seed sets, DR-002).
Use Gauss-Hermite quadrature for EP’s tilted moments, treat the prior as a site, and report a skew-aware interval when the endpoints matter (DR-003).
State a weakly informative prior on a standardised dose instead of a flat prior.

Decision records

Each record states the decision first, then the context, the options considered, why, what happened (weak numbers included) and what I would change. Records are never edited after the fact. A changed decision gets a new record that supersedes the old one. The Markdown sources live in docs/decisions/.

Model card

This card covers the statistical model at the centre of the site and the four ways it is fitted (Laplace approximation, Metropolis-Hastings, Hamiltonian Monte Carlo and expectation propagation). It is a small Bayesian model on a teaching dataset, not a trained system deployed on real decisions. The card follows the usual model-card headings so that its limits are written down in one place.

Model details

Author: Sunchuangyu (Rin) Huang. Written for MAST90125 Bayesian Statistical Learning, University of Melbourne, Semester 2 2023, and revived in 2026. The Metropolis, HMC and EP functions were adapted from the MAST90125 lecture code, as credited in the original submission.
Model: $y_d \sim \text{Binomial}(10, p_d)$ with $\operatorname{logit} p_d = \beta_0 + \beta_1 d$ for doses $d = 0, \dots, 6$ , and a flat prior $p(\beta_0, \beta_1) \propto 1$ , as the assignment specified.
Inference: a Laplace approximation at the posterior mode (a standard logistic GLM), random-walk Metropolis-Hastings, HMC with 20 leapfrog steps of 1/20 and an identity mass matrix, and expectation propagation with Gaussian sites. The two-parameter posterior is also computed exactly on a grid, as a reference.
Implementation: the original R functions, ported line by line to TypeScript in web/src/lib/stats/ and run in the browser. With the original seed, the ports reproduce R's draws to about $10^{-11}$ .
Versions: "Corrected" (the default and the subject of this card) and "As submitted (2023)", which keeps the original data-encoding bugs for transparency (DR-001).

Intended use

Teaching and portfolio use. The site shows how four approximate inference methods compare on a posterior small enough to compute exactly, and what a careful Bayesian workflow adds around them: multi-chain convergence diagnostics, posterior predictive checks, prior sensitivity and method comparison.
Reproducing the 2023 coursework faithfully, and documenting what was wrong with it.

Out of scope

Any clinical, dosing or treatment decision. Nothing here is medical advice.
Inference about the real-world effect of a real medication. The counts come from a coursework exercise, and the study design behind them (randomisation, blinding, how improvement was judged) is not described.
Doses outside 0 to 6. The linear-in-dose logit has no support from the data beyond that range.

Data provenance

The data are the table given in the assignment: seven doses (0 to 6), ten patients per dose, and the number whose symptoms had improved after a week, $y = (3, 5, 4, 4, 7, 6, 8)$ , or 37 of 70 overall.
The table contains counts only. There is no personal information, and nothing was collected by me.
The assignment specification itself is not redistributed. The task is paraphrased in the README and on the site.

Evaluation

All results below are for the corrected model. Monte Carlo results use four chains per sampler with set.seed(90125) to set.seed(90128) and the original settings, and every estimate is reported with its Monte Carlo standard error (MCSE). The convergence diagnostics (rank-normalised split R-hat, bulk and tail ESS, MCSE) are ports of R's posterior package, and the test suite checks them against posterior 1.7.0 on the same draws.

Convergence. Every coefficient passes the thresholds of Vehtari et al. (2021), R-hat at most 1.01 and bulk and tail ESS at least 400. The largest R-hat is 1.0002. The slope's bulk ESS is about 26,500 for Metropolis-Hastings (from 200,000 draws) and 87,200 for HMC (from 36,000 draws, larger than the number of draws because successive HMC draws are negatively correlated).

Estimates of the dose slope $\beta_1$ (log-odds per unit dose):

Method	Mean	95% interval	KL from the exact posterior
Exact posterior (quadrature)	0.3143	[0.063, 0.582]	0
HMC, 4 chains	0.3149 ± 0.0005 (MCSE)	[0.062, 0.582]	Monte Carlo error only
Metropolis-Hastings, 4 chains	0.3154 ± 0.0008 (MCSE)	[0.065, 0.585]	Monte Carlo error only
Expectation propagation	0.3142	[0.056, 0.573]	0.0027
Laplace (GLM)	0.3016	[0.048, 0.555]	0.0085

The posterior probability that the slope is positive is 0.993, and the posterior median odds ratio per unit dose is 1.37. The Laplace mean sits 0.10 posterior sd below the HMC mean, and every other method is within 0.006 sd of it. The best Gaussian possible (the exact mean and covariance) has a KL divergence of 0.0027 from the exact posterior, so EP is as good as a Gaussian can be and Laplace is about three times worse (DR-003).

Posterior predictive checks. For each of the 36,000 HMC draws a replicated trial is simulated. The posterior predictive p-values for the total improved, the number of dose-to-dose drops, the counts at doses 0 and 6, and the chi-square discrepancy are 0.50, 0.39, 0.51, 0.38 and 0.74, with MCSEs between 0.0018 and 0.0025. All 14 observed cells fall inside their 90% predictive intervals. These checks find no misfit, but with 70 patients they have little power to find any.

Prior sensitivity. Under independent $N(0, s^2)$ priors on both coefficients, computed exactly, a weakly informative $s = 2.5$ moves the slope by 0.06 posterior sd and $s = 1$ by 0.29 sd. Only a tight $s = 0.25$ changes the conclusion materially. It shrinks the slope to 0.148 and its 95% interval then includes zero. $P(\beta_1 > 0)$ stays between 0.966 and 0.993 across all the priors.

Known failure modes

Mis-encoded data. The 2023 submission fed each method a different broken version of the data. Four-chain R-hat catches the Metropolis bug (R-hat up to 2.76), but the HMC and EP bug passes every convergence check and is caught only by the posterior predictive check (p = 0 for the total improved). See DR-001.
Skewness. Both Gaussian approximations miss the right skew of the slope's posterior, so their interval endpoints sit 0.007 to 0.009 (EP) and 0.014 to 0.027 (Laplace) below the exact ones.
HMC step size. With the identity mass matrix, step sizes of 0.2 or more almost never accept a proposal, because the intercept and slope are strongly correlated (posterior correlation −0.84). See DR-002.
Proposal scale. Random-walk Metropolis depends on its proposal covariance. With the wrong one, as in 2023, acceptance fell to 5.0%.
Small data. Seven binomial counts cannot distinguish the linear-in-dose logit from many curved alternatives, and the slope remains uncertain (95% interval about 0.06 to 0.58).
The flat prior. A flat prior on the logit scale is not "uninformative" about probabilities. It puts most of its mass near 0 and 1. It is used because the assignment specified it, and the prior-sensitivity analysis shows the conclusions do not hinge on it.

Ethical considerations

The dataset has no personal information. The health setting (an influenza medication) makes over-interpretation the main risk, so the site states plainly that the model is a teaching exercise and that its results say nothing about any real treatment.
The site keeps the original mistakes visible rather than presenting only the corrected analysis, so that the record of the coursework stays honest.
The optional AI feature ("Explain these diagnostics") is bring-your-own-key, sees only a rounded numeric summary, labels its output as AI-generated, runs automatic checks on it and records every call and the human decision in a local audit log. Its design is informed by the transparency principles of the Australian Government's policy for the responsible use of AI in government, the EU AI Act and the NIST AI Risk Management Framework. It does not claim compliance with any of them. See the AI use statement on the methods page.

AI use statement

What AI does here. One optional feature, “Explain these diagnostics” on the diagnostics page. It sends a rounded numeric summary of the convergence diagnostics and posterior predictive p-values to a language model and shows back a plain-English reading, a verdict for each sampler and caveats. The site is complete without it.

What it never does. It never computes, changes or selects a number, a chart or a conclusion on the site. Every statistic comes from the deterministic code described above. It never sees the draws, the data table beyond its totals, other pages or anything about you. It never runs without a click, and it never uses a key belonging to this site, because there is none.

Your key. You bring your own Anthropic or OpenAI key. It is kept in this browser’s sessionStorage, or in localStorage only if you tick “remember on this device”, and a “forget key” button removes it. The browser sends it straight to the provider you chose. It is never sent to this site’s server, never logged and never written to the audit log. The default Anthropic model is Claude Haiku 4.5 (claude-haiku-4-5), with Claude Sonnet 5.5 as an option. For OpenAI the model id is editable (default gpt-5-mini).

Data sent to the provider. Fixed instructions and a JSON summary of about 2,800 characters for the corrected run and 5,700 for the 2023 run (the code refuses anything over 8,000). It holds the model description, the thresholds, and for each sampler the seeds, acceptance rates and each coefficient’s mean, sd, MCSE, R-hat, bulk and tail ESS and failed checks, plus the predictive p-values. You can read the exact summary before sending it. The provider’s own terms govern what happens to it there.

Checks and human review. The answer must match a fixed JSON schema, which is validated in the browser. Two automatic checks run before you see it. One flags any number in the answer that does not appear in the summary. The other compares each sampler’s verdict with the threshold rule the page applies and flags any contradiction. Every answer is labelled “AI-generated” and asks you to accept, edit or reject it.

Audit trail. Each call is recorded with its input, output, model, latency, token usage, check results and your decision in this browser’s IndexedDB. The record is written before the request is sent, and if the browser cannot write it the call is not made. Failed and cancelled calls are recorded too, and when the model replies in a form that fails validation its raw reply is kept. You can read the log and export it as JSON or CSV at /ai-log. It never leaves your device unless you export it.

Evaluation. Because the model’s verdicts compete with a deterministic rule, every answer is also a test case. The log page reports, for each provider, model and summary, the share of answers whose verdicts all agree with the rule and the share that quote a number not in the summary, each as k of n with a Wilson 95% interval, and exports that table as CSV. Nothing is pre-filled: the numbers come only from calls made in your browser.

Frameworks. The design is informed by the Australian Government’s policy for the responsible use of AI in government (Digital Transformation Agency), the transparency principles of the EU AI Act and the NIST AI Risk Management Framework. It does not claim compliance with, or certification under, any of them.