I’m aware that Stan does not use slice sampler in the NUTS paper but sample directly from the probabilities. However, u in the paper is also used to check divergence while current Stan code just use the initial energy, as seen by
I believe this destroys time-reversibility as different states in the trajectory have different energies. Not sure how much this matters though as divergence should be very rare.
A simple correction is to require the difference between the min and max energies in a tree has to be smaller than divergence threshold.
parameters {
vector[3] x;
}
model {
x ~ std_normal();
if (x[1] < 0)
target += -1000;
}
generated quantities {
real q = -1.0 + 2.0*std_normal_cdf(x[1]); // q ~ uniform(0,1)
}
Theoretically the model implies that P(x[1] < 0) is negligible and x[1] should have a half-normal distribution.
But a difference of -1000 is exactly at the threshold of divergence detection; whether the trajectory is considered divergent or not depends on the initial point.
One million draws from Stan with with adapt delta=0.5 gives a visibly biased histogram:
Difference of -995 does not see any divergencies; and a difference of -1005 sees a divergence every time the trajectory crosses x[1]=0 (independent of the initial point); in either case Stan recovers the correct distribution.
One should also suspect bias from the fact that about 50% of all transitions are divergent. Even if we fix the divergence check to maintain reversibility, samples with divergencies are unreliable because divergent transitions are a sign of HMC being unable to move across the parameter space.
