revised NUTS description; fixes #436 by bob-carpenter · Pull Request #966 · stan-dev/docs

bob-carpenter · 2026-07-16T18:34:14Z

This revises the description of the HMC and NUTS algorithms with corrected and added citations.

Submission Checklist

Builds locally
[n/a] New functions marked with <<{ since VERSION }>>
Declare copyright holder and open-source license: see below

Summary

Copyright and Licensing

Please list the copyright holder for the work you are submitting (this will be you or your assignee, such as a university or company):

Simons Foundation

By submitting this pull request, the copyright holder is agreeing to license the submitted work under the following licenses:

Code: BSD 3-clause (https://opensource.org/licenses/BSD-3-Clause)
Documentation: CC BY-ND 4.0 (https://creativecommons.org/licenses/by-nd/4.0/)

jgabry

Thanks for updating this. No major suggestions. I found a few typos and had a couple of questions about the math. I hope I didn't miss anything else!

jgabry · 2026-07-16T20:06:57Z

+\begin{align}
+\rho &\leftarrow
+      \rho \, + \, \dfrac{\epsilon}{2} \, \nabla \, p(\theta \mid y)
 \\[6pt]
-\theta & \leftarrow
-       & \theta \, + \, \epsilon \, M^{-1} \, \rho
+\theta &\leftarrow
+       \theta \, + \, \epsilon \, M^{-1} \, \rho
 \\[6pt]
-\rho & \leftarrow
-     & \rho \, - \, \frac{\epsilon}{2} \frac{\partial V}{\partial \theta}.
-\end{array}
-$$
+\rho &\leftarrow
+     \rho \, - \, \dfrac{\epsilon}{2} \, \nabla \, p(\theta \mid y)
+\end{align}


A couple things I noticed here:

p(theta | y) should be log(p(theta|y)) on lines 150 and 156, right?

I think both momentum half-steps should have the same sign, right? You use + for one and - for the other. I think both should be + here if you use log(p(theta|y)) instead of V(theta).

jgabry · 2026-07-16T20:10:52Z

 $$
+\min\left(
+1,
+\dfrac{p(\theta^*, \rho)}{p(\theta, \rho)}


Suggested change

\dfrac{p(\theta^*, \rho)}{p(\theta, \rho)}

\dfrac{p(\theta^*, \rho^*)}{p(\theta, \rho)}

jgabry · 2026-07-16T20:15:17Z

+T(\rho) = - \log p(\rho  \mid  M).
 $$
-V(\theta) = - \log p(\theta)
+Here, $\log p(\rho \mid M) = \textrm{MultiNormal}(\rho \mid 0, M)$,


I think here you either need to drop log from the LHS or add it to the RHS

jgabry · 2026-07-16T20:19:34Z

-If the inverse metric $M^{-1}$ is a poor estimate of the posterior covariance,
-the step size $\epsilon$ must be kept small to maintain arithmetic
-precision. This would lead to a large $L$ to compensate.
+If the inverse metric $M^{-1}$ is a poor estimate of the local


With the notation you've been using in this chapter, should this be M (not M^{-1}) that's an estimate of the curvature?

jgabry · 2026-07-16T20:20:35Z

+the same amount of work to generate each new roughly independent
+draw.
+
+Warmup times being vast only holds when the posterior is log concave,


Suggested change

Warmup times being vast only holds when the posterior is log concave,

Warmup times being fast only holds when the posterior is log concave,

jgabry · 2026-07-16T20:23:58Z

-where the term
-
+The two terms into which it separates are treated in
+the Hamiltonian dynamics as a a potential energy function defined by


Suggested change

the Hamiltonian dynamics as a a potential energy function defined by

the Hamiltonian dynamics as a potential energy function defined by

jgabry · 2026-07-16T20:24:09Z

-value is selected via multinomial sampling with a bias toward the
-second half of the steps in the trajectory @Betancourt:2016.^[Stan previously used slice sampling along the trajectory, following the original NUTS paper of @Hoffman-Gelman:2014.]
+determined by the parameters drawn in the last iteration.  Like HMC,
+it starts by Gibbs sampling a a random momentum vector with covariance


Suggested change

it starts by Gibbs sampling a a random momentum vector with covariance

it starts by Gibbs sampling a random momentum vector with covariance

jgabry · 2026-07-16T20:24:37Z

-@Betancourt-Girolami:2013 into the notation of @GelmanEtAl:2013.
-
+method that scales well in dimension due to its use of gradients of
+the the log density function being sampled.  It uses an approximate


Suggested change

the the log density function being sampled. It uses an approximate

the log density function being sampled. It uses an approximate

jgabry · 2026-07-16T20:25:03Z

+consisting of the sum of the densities in the new doubling divided by
+the sum of the densities in the current trajectory.  Because these are
+the same length and the Hamiltonian should be preserved, the
+probabability of jumping to the new doubling is close to 1 when the


Suggested change

probabability of jumping to the new doubling is close to 1 when the

probability of jumping to the new doubling is close to 1 when the

jgabry · 2026-07-16T20:25:25Z

+the sum of the densities in the current trajectory.  Because these are
+the same length and the Hamiltonian should be preserved, the
+probabability of jumping to the new doubling is close to 1 when the
+leapfrog algoritm is stable.


Suggested change

leapfrog algoritm is stable.

leapfrog algorithm is stable.

jgabry · 2026-07-16T21:13:19Z

+U-turn arises when an infinitesimal extension of the trajectory has
+the ends moving closer to each other.  It also stops doubling if a
+maximum number of doublings is achieved (maximum 10 by default, for a
+total of 1024 leapfrog steps).


Suggested change

total of 1024 leapfrog steps).

total of 1024 leapfrog steps, including the initial state).

It's really 1023 steps, right? 1024 including the initial state.

revised NUTS description; fixes #436

0494ce1

bob-carpenter requested review from jgabry and mitzimorris July 16, 2026 18:34

jgabry requested changes Jul 16, 2026

View reviewed changes

jgabry reviewed Jul 16, 2026

View reviewed changes

	\dfrac{p(\theta^*, \rho)}{p(\theta, \rho)}
	\dfrac{p(\theta^, \rho^)}{p(\theta, \rho)}

	Warmup times being vast only holds when the posterior is log concave,
	Warmup times being fast only holds when the posterior is log concave,

	the Hamiltonian dynamics as a a potential energy function defined by
	the Hamiltonian dynamics as a potential energy function defined by

	it starts by Gibbs sampling a a random momentum vector with covariance
	it starts by Gibbs sampling a random momentum vector with covariance

	the the log density function being sampled. It uses an approximate
	the log density function being sampled. It uses an approximate

	probabability of jumping to the new doubling is close to 1 when the
	probability of jumping to the new doubling is close to 1 when the

	total of 1024 leapfrog steps).
	total of 1024 leapfrog steps, including the initial state).

Uh oh!

Uh oh!

Conversation

bob-carpenter commented Jul 16, 2026

Submission Checklist

Summary

Copyright and Licensing

Uh oh!

jgabry left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants