I’ve been working on adapting Stan to work with a C++ framework I need to use for generating the likelihoods I need to fit my data. Mostly based on reverse-engineering the output I got from a test model I wrote and ran with CmdStan (and examining CmdStan’s
command.hpp), I’ve derived a new class from
stan::model::prob_grad that implements (amongst other things)
log_prob() as well as a
stan::callbacks::writer derivative to store the samples that come out in a format I can work with. I’m running NUTS with adaptation using a call to
stan::services::sample::hmc_nuts_diag_e_adapt(). (I’d be happy to give more technical details if necessary, but I didn’t want to clutter the original post, and I don’t think they’re needed for this discussion.) This appears to at least run Stan successfully.
Now I’m trying to test whether the fitting is actually working using an extremely simple model before I trust it to do anything sensible elsewhere. I’ve run up against a feature I didn’t expect and don’t understand. My model is a quadratic function with a single parameter, y = a x^2, and my likelihood is a unit normal centered around a true value a_0, with y evaluated only at a single data point of x=1, so that L(x;a) = N(a - a_0). For the sake of simplicity, I’m assuming the prior on a is flat, so it contributes nothing to the total LL. This makes
log_prob look like the plot below as a function of a (in this example a_0 = 1.8, selected at random):
When I run Stan, it’s able to find the correct value a_0 with no trouble, but no matter how many samples I run in warmup (I’ve tried as many as 10K) or how large I set the max tree depth (I’ve gone as high as 17), literally 100% of my post-warmup samples saturate the max tree depth.
I’ve inserted lots of debugging code to ensure that the return value of
log_prob is literally what I plot above, so I know Stan’s getting the right LL values. And the LL surface obviously doesn’t have any strange features for the sampler to get stuck in. I’ve read what I can find in the Stan manual and online (including the page about warnings) and I can’t work out whether I should be concerned.
I’m a bit worried, though, that this might indicate I have a problem in my implementation, and I especially worry that sampling won’t work correctly when I try a more complex model. (I’ve experimented with more complex models and had the same behavior, and where the samples don’t match my expected posterior surface, but it’s less clear from those whether it’s the model or my Stan usage that’s at fault.)
Any insight – suggestions for things to check, further reading in the Stan manual or elsewhere that I might have missed, assurance that everything is ok and I don’t need to worry, etc. – would be extremely welcome. Thanks!