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?

2 Likes

Hi,
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 https://github.com/stan-dev/stan/tree/develop/src/stan/services/sample

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 :)