Why NUTS for QSP Calibration?
Calibrating a nonlinear mechanistic model can produce a posterior with narrow, correlated regions: one parameter may need to rise when another falls to preserve a similar model output. Hamiltonian Monte Carlo (HMC) uses gradients to move through such regions more efficiently than a random walk, but ordinary HMC requires a trajectory length to be chosen in advance. A trajectory that is too short may make little progress; one that is too long may retrace its path. No-U-Turn Sampling (NUTS) addresses this tuning problem by growing a trajectory and stopping when its direction indicates that it has begun to turn back.
NUTS does not eliminate tuning or diagnose every sampling problem. It adapts the leapfrog step size during warmup and chooses a trajectory length dynamically for each iteration, subject to a maximum tree depth. Its output therefore gives useful information about numerical integration and the sampler's effort, but it is not by itself proof that the posterior has been explored adequately.
Hamiltonian Motion and the Stopping Rule
Let \(\theta\) denote the vector of model parameters on the scale used by the sampler, and let \(\pi(\theta\mid y)\) be their posterior density given observations \(y\). For example, positive biological rate constants are often represented through unconstrained transformed parameters. Define the potential energy as the negative log posterior, up to an additive constant, and introduce an auxiliary momentum \(p\) with mass matrix \(M\). The Hamiltonian is
HMC numerically follows Hamiltonian dynamics using leapfrog steps. In exact arithmetic the dynamics preserve \(H\); a finite step size \(\epsilon\) introduces energy error. A Metropolis correction accounts for that error, with acceptance probability \(\min(1,\exp(H_{\mathrm{current}}-H_{\mathrm{proposal}}))\). NUTS uses this same basic motion, but constructs a trajectory in both directions by repeatedly doubling the number of leapfrog steps. It checks whether the growing path has started to turn back and stops expansion when a U-turn is detected.
A common Euclidean-metric U-turn check uses the displacement between the two ends of the trajectory and the endpoint momenta. Write \(\theta^{-}\) and \(\theta^{+}\) for the left and right endpoint positions, and \(p^{-}\) and \(p^{+}\) for their momenta. One form stops if either endpoint momentum points back across the trajectory:
This test is a practical stopping rule, not a claim that all coordinates have reversed direction. Implementations can use refinements, including checks based on accumulated momenta, so details should be read in the context of the software being used. In all cases, a U-turn means the trajectory should not continue expanding as if its added states were still carrying it into new posterior territory.
Worked Example: Applying a One-Dimensional U-Turn Check
Suppose a one-dimensional trajectory has left endpoint \(\theta^{-}=-1.2\), right endpoint \(\theta^{+}=1.0\), endpoint momenta \(p^{-}=0.10\) and \(p^{+}=-0.20\), and \(M=1\). Apply the endpoint-momentum criterion.
The displacement is \(\theta^{+}-\theta^{-}=1.0-(-1.2)=2.2\). Since \(M^{-1}=1\), the two products are
The first product is positive, but the second is negative. The criterion therefore signals a U-turn at the right endpoint, so tree expansion stops. The calculation illustrates the direction check; it does not determine which trajectory state is ultimately selected as the next sample.
Step-Size Adaptation During Warmup
The leapfrog step size \(\epsilon\) controls the distance advanced in each numerical step. A large value can produce substantial Hamiltonian error, lowering acceptance or causing numerical instability. A very small value usually reduces integration error, but requires many steps to cover a given distance. NUTS therefore adapts \(\epsilon\) during warmup toward a user- or software-specified target acceptance statistic. A higher target generally calls for a smaller step size; the exact relationship depends on the posterior geometry and mass matrix.
Many implementations use dual averaging, a stochastic adaptation scheme that updates a running step-size estimate using acceptance information. The adaptation uses the trajectory's acceptance statistic, often an average over its leapfrog proposals, rather than treating one proposal's accept/reject outcome as a sufficient tuning signal. Warmup may also estimate or update a mass matrix, which rescales momentum and helps accommodate posterior correlations. Step size and mass matrix are distinct: changing the metric can change the useful step size.
After warmup, the tuned quantities are normally held fixed for the sampling phase. Continuing to adapt the transition rule while collecting ordinary posterior draws can invalidate standard sampling guarantees unless a method specifically accounts for that adaptation. The number of leapfrog steps, by contrast, can still vary from one NUTS sampling iteration to another because the U-turn rule is applied dynamically.
Worked Example: Step Size and Leapfrog Energy Error
Consider a hypothetical one-dimensional standard normal target, with \(U(q)=q^2/2\), \(M=1\), and initial state \(q=0,\ p=1\). For a leapfrog step of size \(\epsilon=1\), calculate the state and the Metropolis acceptance probability. This simple calculation isolates numerical integration; it is not a full QSP posterior.
The initial energy is \(H_0=U(0)+p^2/2=0+1/2=0.5000\). For leapfrog, the first half-step updates momentum using the gradient \(dU/dq=q\):
The second half-step gives \(p_1=1-\frac{1}{2}(1)(1)=0.5\). The new energy is
The acceptance probability is \(\min(1,\exp(-\Delta H))=\exp(-0.125)\approx0.8825\). This single step has a modest positive energy error, so it is accepted with probability about 88.3%. During warmup, repeated trajectory-level acceptance information—not this single result alone—would inform step-size adaptation. If similar trajectories repeatedly had acceptance below the target, adaptation would generally reduce \(\epsilon\).
Interpreting NUTS Output in a Mechanistic Model
NUTS software commonly reports a step size, tree depth or number of leapfrog steps, an acceptance statistic, and divergences. These quantities describe different aspects of the computation. The step size is a numerical resolution; tree depth reflects how far the algorithm expanded before stopping; the acceptance statistic reflects energy error; and a divergence flags a trajectory with an unusually large integration error under the implementation's criterion. None is a direct measure of biological plausibility.
Worked Example: Reading a Hypothetical Calibration Summary
Suppose a calibration of a nonlinear QSP model reports, after warmup, step size \(0.035\), maximum tree depth 9, and the following counts over 1,000 sampling iterations: 18 iterations reached the maximum depth and 7 had a divergence.
The maximum-depth fraction is \(18/1000=0.018\), or 1.8%; the divergence fraction is \(7/1000=0.007\), or 0.7%. These are descriptive rates for this run, not universal pass/fail thresholds. Reaching maximum depth means the algorithm stopped because of its configured computational limit rather than its ordinary U-turn rule. It can indicate that longer trajectories might be useful, but frequent maximum-depth hits can also reflect difficult geometry or expensive trajectories.
The 7 divergences are a numerical warning: some transitions encountered integration error large enough to trigger the divergence rule. In a mechanistic model, potential causes include strong posterior curvature, poorly scaled or highly correlated parameters, a difficult region near a parameter boundary, or inaccurate gradients from the ODE solver. The summary does not identify which cause is responsible. Inspecting the affected parameter regions and checking solver and gradient accuracy are sensible next steps. The 1.8% maximum-depth rate and 0.7% divergence rate alone do not establish that the remaining draws represent the full posterior.
Assumptions, Diagnostics, and Alternatives
NUTS relies on the same core ingredients as gradient-based HMC. The target posterior must be a well-defined distribution, and its log density should have usable gradients in the coordinates being sampled. The numerical integrator should be stable and sufficiently accurate over the relevant posterior region. For a posterior involving an ODE likelihood, this requires not only a suitable ODE solution but also gradients that are consistent with the computed log likelihood. Solver tolerances that are too loose can make the apparent energy error reflect solver noise rather than the intended posterior geometry.
When parameters have constraints, such as positive rates, sampling commonly occurs on transformed unconstrained scales, with the appropriate transformation contribution included in the target density. A discontinuous likelihood, unreliable derivative, or unaccounted-for Jacobian can undermine the calculation; adaptive trajectory length cannot repair an incorrectly specified posterior. Diagnostics should therefore include checking parameter transformations and prior support, testing gradients and ODE accuracy, and reviewing divergences and maximum-depth behavior. Convergence and effective sample size require additional diagnostics and are not established by this tutorial's output summaries.
NUTS is most useful when gradients are available and the posterior is sufficiently smooth for Hamiltonian motion to be informative. Fixed-length HMC can be a reasonable alternative when trajectory length has been tuned successfully and predictable computational cost is valuable. It avoids NUTS's tree-building overhead but requires choosing the number of leapfrog steps. Random-walk Metropolis does not require gradients and can handle some nonsmooth targets, but may explore correlated, high-dimensional posteriors inefficiently. These methods do not remove the need to define and check the likelihood, prior, and model implementation.
Confirm that priors, likelihood, transformations, and ODE calculations define the intended posterior and provide reliable gradients.
Interpret the tuned step size as a warmup result, and check that tuning is not unintentionally continuing during ordinary sampling.
Examine divergences and maximum-depth hits in context rather than treating a single summary value as proof of success or failure.
Assess whether chains have explored consistently and whether the quantities of interest are estimated precisely; acceptance and tree depth alone cannot answer those questions.
Check Your Understanding
Answer each question, then compare your reasoning with the solutions.
- What problem does the U-turn criterion address that a fixed number of HMC leapfrog steps leaves to be tuned?
- If repeated warmup trajectories have acceptance statistics below the target, what change to step size is generally indicated, and why?
- A sampling run has 20 maximum-depth hits among 2,000 iterations. What fraction hit the limit, and what does that observation mean?
- Why can a divergence in an ODE-based QSP calibration not be interpreted simply as evidence that the biological model is wrong?
Solutions. (1) A fixed-length trajectory can be too short or can continue after it has started turning back; the U-turn rule dynamically stops expansion when the path begins to reverse. (2) Step size should generally decrease: smaller leapfrog steps tend to reduce Hamiltonian integration error and increase acceptance, although the actual adaptation uses repeated trajectory information. (3) The fraction is \(20/2000=0.010\), or 1.0%. Those iterations reached the configured computational limit before ordinary U-turn stopping; this may warrant investigation, but the fraction alone does not establish a sampling failure. (4) A divergence indicates problematic numerical integration of the posterior transition. It can arise from posterior geometry, parameter scaling, inaccurate ODE or gradient calculations, or other computational issues, so it does not by itself distinguish biological misspecification from numerical difficulty.