MCMC and its diagnostics

Session 6 will put all twelve Halloway outlets in one model: a rate per outlet, and above them two more parameters that say where the outlets cluster and how tightly. Each outlet's rate depends on those two, and those two depend on every outlet's rate. That mutual dependence is what "interact" means, and it is the end of the closed form. There is no grid fine enough either; session 3 already explained why the grid dies with the second parameter.

Markov chain Monte Carlo sidesteps the whole problem. Instead of evaluating the posterior everywhere, it starts somewhere, proposes a step, and accepts or rejects it by comparing the posterior's height at the two points, which only needs the unnormalised product of prior and likelihood, the thing you can always compute. Run that for long enough and the places the walk visits are distributed as the posterior. The walk's record is a chain, and a chain is a set of draws. Read it as session 3 read its draws: mean, quantiles, probabilities above a threshold.

• Chain. One walk through the posterior, recorded step by step. Each step depends only on the one before it, which is the Markov part. The record is a set of draws.
• Warm-up, tuning. The first stretch of a chain, spent finding the posterior and adjusting the step size. Discarded. Under nutpie in PyMC 6 the default is 400 steps; under PyMC's own NUTS it is 1000.
• NUTS, nutpie. NUTS is the sampler family PyMC uses. nutpie is a separate implementation of it, and when it is installed PyMC 6 uses it by default. pip install pymc[nutpie].
• R-hat. Do the chains agree with each other? A ratio that is 1 when they do and grows as they do not.
• Effective sample size, ESS. How many independent draws the chain is really worth. Consecutive steps of a walk are correlated, so a chain of many draws can carry the information of far fewer.
• Divergence. A step the sampler had to abandon because the posterior's geometry was too sharp to follow. Each one is a place the sampler could not go.
• DataTree. The container a fit comes back in: an xarray.DataTree with a group per kind of output. The posterior is dt["posterior"].dataset. InferenceData no longer exists.
• Rank plot. plot_trace_rank. Every draw across all chains is ranked, and each chain's ranks are histogrammed. Flat means mixed; sloped means a chain is living in one part of the posterior.

A chain's first stretch is spent getting from wherever it started to the region where the posterior lives, and adjusting its step size so that it moves through that region efficiently. Those steps are not draws from the posterior and are thrown away; the word for them is warm-up, or tuning. PyMC 6 with nutpie installed makes nutpie the default NUTS sampler and tunes for 400 steps by default, against 1000 under PyMC's own NUTS. Install with pip install pymc[nutpie]. PyMC 6 needs Python 3.12 or later, and it is built on PyTensor 3.0 with numba as the default backend.

The figure's fourth chain is the case warm-up exists for, and the case it does not always fix. It started far from the posterior, and 400 steps was not enough for it to arrive. Kept draws from a chain that has not arrived are not draws from the posterior; they will pull the mean toward wherever the chain has been, and nothing in that chain's own record will say so.

Run four chains from four different starting points. If they have all found the posterior, each chain's draws will spread over the same region, and the spread between the chains' means will be small compared with the spread within each chain. If one is stuck, or still travelling, the between-chain spread will be large relative to the within-chain spread. R-hat is that comparison as a single number: 1 when the chains agree, larger as they disagree.

The value of the idea is that it needs nothing from you. You do not have to know where the posterior is, or what it looks like, to notice that four walks meant to describe the same thing are describing different things. What R-hat cannot do is tell you that four chains which agree with each other have found the whole posterior rather than the same part of it, which is why it is one of three diagnostics and not the only one.

R-hat asks whether the chains agree. az.rhat(dt). It catches the figure in Part 1: a chain still travelling, or stuck in one region while the others are elsewhere. Effective sample size asks how many independent draws the chains are really worth. az.ess(dt, sample_dims=["chain", "draw"]), or the xarray accessor dt.azstats.ess() after import arviz_stats. A chain moves in small correlated steps, so consecutive draws repeat information; ESS is the count after that repetition is discounted, and it is the number to compare with the 20,000 draws session 3 used, not the raw length of the chain. Divergences are steps the sampler abandoned because the posterior's geometry was too sharp to follow. Each one marks a place the sampler could not go, and a posterior it could not fully visit is a posterior it is not fully reporting.

The plots have new names. plot_trace no longer exists; it split into plot_trace_dist, which shows each parameter's distribution per chain, and plot_trace_rank, which shows the rank histograms in the figure. The rank plot is the one to read for mixing. Every ArviZ 1 plot returns a PlotCollection rather than axes. ArviZ itself is now three libraries under one namespace, arviz-base for data, arviz-stats for diagnostics, arviz-plots for plotting, all reached through import arviz as az.

What pm.sample() hands back is an xarray.DataTree: a tree of groups, one per kind of output, each group an xarray Dataset. The posterior draws are at dt["posterior"].dataset, with dimensions chain and draw, which is why az.ess takes sample_dims=["chain", "draw"]. arviz.InferenceData no longer exists; anything you read that returns one is describing a version that is gone.

Old files still open. az.from_netcdf now wraps open_datatree, so a netCDF or zarr file saved from an earlier fit loads into the new container without conversion. The groups are reached the same way.

R-hat near 1, a healthy effective sample size, no divergences: that says the sampler explored the posterior of this model. The draws are a faithful picture of what this prior and this likelihood, taken together, imply about the rate. It is a statement about the computation, and it is the statement you need before you trust a single number that comes out of it.

It says nothing about whether the model is right. A binomial likelihood with a flat prior on a single rate can be sampled perfectly and still be the wrong model for a chain of twelve outlets that differ; the diagnostics will be green and the estimate for Northgate will still be far from 0.192. Checking the model against the data it is meant to describe is a different job, with different tools, and session 8 owns it. Tonight's rule is to keep the two questions apart: did the sampler work, and is the model any good. A green diagnostic answers the first. Nothing in tonight's toolkit answers the second.