Using Stan with stand-alone models



I have a model that takes the form of a moderately large C++ program (about 10K lines of source code). The model has a number of unknown parameters, and it takes those parameters and produces a bunch of outputs that I can compare to an observational dataset in the usual way. Clearly I’d like to use a Markov Chain Monte Carlo, but I’d like to do something better than using the crusty old Metropolis-Hastings code I wrote in graduate school and have been carting around from job to job for the last 20 years. Can Stan provide a better solution?

My first impression, after skimming the documentation, is that it cannot, because there is no good way to get derivatives of the model outputs with respect to the parameters. (Retrofitting the necessary data structures for Stan’s autodifferentiation into 10K lines of code seems like an invitation to madness.)

However, while I was poking around learning about Stan I came across this talk:

At about 3:26 in the talk the speaker has a slide titled “Who is using Stan?”, and on that slide he says

Stan is used for fitting climate models, …

That’s exactly what I’m trying to do (the model I described above is in fact a simple climate model). So, my real question is, does anybody know what work the speaker is referring to? Who is using Stan to fit climate models, and are there any examples I could use as a jumping off point to figure out how my colleagues and I might be able to use Stan to fit our model (and finally retire the 20-year old Metropolis code!)


Can’t say exactly what the company is or does, but hopefully it is within the realm of possibility for you to modify your C++ program to add template parameters (if you don’t have them already) so that your posterior kernel evaluation function can be called with Stan’s types. Then, Stan will be able to use autodiff with your modified C++ code and a little bit of documented hacks.


Does “the company” refer to the people who are using Stan with the climate models? I.e., it’s not NCAR or someplace like that?

Anyhow, what it sounds like you’re saying is that retrofitting Stan’s data structures into the model source is pretty much the only way to make it work. Is that an accurate summary?


There may be plenty of people using Stan for climate models that I don’t know of (I’m not in that field) but there is a company using Stan for such things that I do know of but can’t talk about.

Anyway, it may be less work to just write everything in the Stan language (and that usually takes less time than people fear), but I was saying the opposite: that it may be less work to just add

template<typename T>

to a bunch of your existing C++ functions, change the signatures of those functions to use T rather than double, and utilize our mechanism to find and use external C++ code from Stan programs. I think some people in astronomy go this route.

How complicated are we talking here? You must already have a function that evaluates the logarithm of the posterior kernel in order to do M-H. And that function is a method of some ugly C++ class that holds a bunch of data? Or that function calls a lot of custom numerical libraries to do stuff that cannot be done analytically?


Here’s the repository for the model I’m talking about, if anyone reading the discussion is curious (or if they’re looking for a good climate model).

The process of computing the likelihood a little complicated because the model was designed as a stand-alone, and the MC was grafted on as an afterthought. The simplified version is that the core of the model is made up of a bunch of components that do things like tracking the movement of carbon between various pools, simulating diffusion in the ocean, tracking the concentrations of non-CO2 GHGs, and so on. As you surmise, there are a lot of numerical solvers and whatnot involved in that process. At the end, what you have is a bunch of vectors of output data, and the function that calculates the likelihood harvests these and compares them to the observational data.

One further wrinkle is that in order to keep the model modular, each component runs independently at each time step, storing its results in vectors that it owns. If a component needs a value from another component (e.g., at each time step the ocean component will need the CO2 concentration from the atmosphere), then it sends a message to the coupler, which routes that message to the appropriate component, which in turn responds with the requested value. I worry that such indirection will confuse the autodifferentiation scheme. Can it handle the case where a function picks up a stored value that was computed in a different function (typically not even in the same translation unit)?


That’s often the case in a lot of big scientific models.

Stan can use such a thing standalone with a little work on the C++ side if you can also compute derivatives. Otherwise, it has to be rewritten to use Stan’s math library so it can do autodiff.

We have solvers for ODEs, etc., but they may be computationally challenging to compute and/or differentiate at large scales.

Yes, as long as it’s stored as it’s proper autodiff type.


I see. I suspect that our ODEs don’t qualify as “large scale”. We’re not spatially resolved like an Earth System Model, so we’re really only tracking a dozen or so variables in the largest systems of equations, and only at annual resolution.

Thanks for all the information. I think I have enough to float the idea of converting to the rest of the development team. If I can get them to go along with it, I’ll be sure to let you know how it turns out.


You might want to have a look at the BayesianTools package in R which provide various MCMC-like samplers for complex models where the likelihood cannot be written down in a simple function.


So you mean something like ABC where you assume you can sample, but assume you can’t evaluate the log density?

If you have the likelihood, you can use Metropolis or an ensemble method and there are programs that let you plug those in. emcee in Python implements an ensemble method. The problem with both Metropolis and the ensemble methods is that they don’t scale well with dimension.