Using Stan HMC as Metropolis-Within-Gibbs step in C++

Hi everyone!
I’m building a C++ library for Bayesian Nonparametric models, e.g. Dirichlet Process Mixtures but not limited to that.
In order to gain efficiency in cases when likelihood and base-measure are not conjugate (e.g. Neal’s Algorithm 8), I was thinking of using HMC instead of plain Metropolis-Hastings.

However, it is quite unclear to me how to proceed.
Is there any tutorial / gist on how to write a valid Stan model in C++ and perform sampling?
More broadly, do you have any suggestions?


sorry your question fell through a bit.

I don’t think so, but you can check out how CmdStan or Rstan does things - it starts by translating the model into a C++ class that has a log_prob() member that calculates the density. Computing the density with Stan’s autodiff var gives you the gradient. The HMC code can be found in the Stan project, at

Generally however, the devs always noted that using NUTS as a subroutine in another sampler usually doesn’t work because the adaptation in NUTS is crucical but very easy to break.

Hope this gives you a starting point.

Thank you @martinmodrak, I will look into that

Actually I was planning on using standard HMC, as NUTS is not a valid Metropolis step and that might hurt the ergodicity of the chain, from my understanding :)