Dan Foreman-Mackey writes:

I [Mackey] demonstrate how to use a custom C++ function in a Stan model using the Python interface PyStan. This was previously only possible using the R interface RStan (see an example here) so I hacked PyStan to make this possible in Python as well. . . .

I have some existing C++ code that implements my model and its derivatives so I don’t want to have to re-write the model in the Stan language. Furthermore, part of the model requires solving a transcendental equation numerically; it’s not obvious that applying autodiff to an iterative solver is a great idea, but the analytic gradients are trivial to evaluate.

The example that we’ll use is fitting radial velocity observations of an exoplanet. In particular, we’ll fit recent observations of 51 Peg b, the first exoplanet discovered around a main sequence star.

An exoplanet! Cool.

Mackey continues with tons of details. Great stuff.

I have some issues with his Stan model; in particular he uses priors with hard constraints which I think is generally a bad idea. For example, he has a parameter with a uniform (-10, 5) prior. I can’t imagine this is the right thing to do. From basic recommendations it would be better to do something like normal (-2.5, 7.5), but really I have a feeling he could do a lot more regularization here. The uniform prior might work in the particular example that he was using, but in general it would be safer to control the inference a bit more. Mackey’s got a bunch of these uniform priors in his code and I think he should look carefully at all of them.

The real point of Mackey’s post, though, is that he’s hacking Stan to solve his problem. And that’s great.