Hamiltonian Monte Carlo, in anger

The posterior is a funnel. Below a certain scale the step size that works in the neck is catastrophically wrong in the mouth, and the sampler either grinds to a halt or sails straight past the region carrying the mass. Non-centred parameterisation is not a trick; it is a change of coordinates that makes the geometry tractable.

Where the sampler stalls

Neal's funnel is the canonical example. Write a hierarchical model in its natural form and the joint density develops a region so narrow that no single step size can explore it.

The denominator is the whole problem. For anything beyond a toy model it has no closed form, which is why we sample instead of integrate. Inline maths works too: the scale parameter $\tau$ controls how tight the neck becomes, and as $\tau \to 0$ the geometry degenerates.

Fig. 1analysis/funnel_plot.py

Diagnosing it

Divergent transitions are the signal. Do not filter them out.

import arviz as az

idata = az.from_cmdstanpy(fit)
print(az.summary(idata, var_names=["mu", "tau"]))

# divergences are stored alongside the draws
n_div = int(idata.sample_stats["diverging"].sum())
print(f"{n_div} divergent transitions")

The fix

Sample a standardised variable and transform. The sampler explores a well-conditioned space, and the pathological geometry disappears.

parameters {
  vector[J] theta_raw;
  real<lower=0> tau;
}
transformed parameters {
  vector[J] theta = mu + tau * theta_raw;
}
model {
  theta_raw ~ std_normal();
  y ~ normal(theta, sigma);
}

Run it again and the divergences vanish. The lesson generalises: when a sampler struggles, suspect the parameterisation before you suspect the sampler.¹