Hi,
I was wondering if there has been any progress on this question here.
The idea is being able to use Stan from C++ as easily as it is from R or Python. My use case is that I am doing data acquisition and pre-processing in C++ and writing files to pass the data to cmdstan is a real bottleneck.
Naively, I was thinking that, once the stan file is compiled to a C++ class, there would be a simple way to fit the model to data provided as through Eigen data structures.
Thanks!
Garjola
If you have the C++ file generated from a Stan program, then it is not so hard. The C++ file will have a class with a method that starts like
template <bool propto, bool jacobian, typename T_>
T_ log_prob(Eigen::Matrix<T_,Eigen::Dynamic,1>& params_r,
std::ostream* pstream = 0) const {
so it can be evaluated with a long vector of proposals for the real parameters. I think you pretty much need to just call the constructor, although that requires a stan::io::var_context, and then it is ready to go.
Thanks for your quick answer. I think I’ll have to get familiar with the API before understanding how to use stan::io::var_context.
I have been looking at the code generated from the Bernouilli example and I don’t think I’ll be able to understand how to pass the data to fit the model.
Maybe one of the Stan developers could take a simple example and show how to do it?
Is there link to the Doxygen documentation for the API? I have found the description of the model concept which can be useful for starters, but links to other classes would be helpful.
Thanks again!
There’s some code in the unit tests that might be helpful, specifically:
stan/src/test/unit/services/sample/standalone_gqs_big_test.cpp
also, you could look at some of the cmdstan code:
C++ isn’t a very convenient language for optional arguments and heterogeneous data structures like dictionaries. The var_context type is a very low level interface meant to read data out of flat memory structures. What you’d want is a var_context you could add variables to one at a time in an easy way, which I don’t think we currently have. We’re going to need something like that soon to do round trips and restarts in a structured way.
One of the reasons we’ve been reluctant to dump a lot of time into documenting the C++ API is that it’s been pretty volatile compared to the Stan language and interface APIs.
I understand the issue about C++ convenience (alhough heterogeneous containers and optional are coming to the language and already available in boost).
Wouldn’t it be useful to reuse what has been done for python and R ?
The bridge between python/R and C++ seems to be doing what is missing here (and it is already in C++). If such a C++ API was provided in Stan, maintaining other interfaces should be easier also, since it could hide the evolutions of the internal “volatile” API.
I know I am not helping much here, but I am reluctant to write code that somebody else can write better than me ;)
What do you mean? I think the direction is going the wrong way… what was done in Stan (now CmdStan) was reused for Python and R.
As @mitzimorris has mentioned, there are actually a bunch of examples in the Stan repo’s unit tests of creating classes that implement the var_context interface that are based directly on C++ objects like std::vector. So… I think what you’re looking for is already there, but isn’t packaged in a way that’s really robust and general.
If you provide a lot more detail and a concrete use case, I can provide some more help.
Not intend to hijack the thread, I think a C API for Stan & Math would be more helpful. The immediate benefit would be in extending with a module system as well as interfacing with other languages. It’d be a chore, though, especially when comes to specializing and wrapping all the templated functions and catch all the exceptions before turning them into error codes. But to me it’s only long-term solution to reach out beyond C++.
That would actually be amazing and simplify the integration with other languages! But you’re right: a ton of work. That’s why most platforms have a specialized package for dealing with C++ that usually relies on name mangling.
I am also interested in calling stan from C++, in particular so I can run over say Gigabytes of data inside the likelihood. This is crucial for particle-physics applications.
My context: I can evaluate my posterior and the gradient in C++, how do I call stan to run HMC on my posterior? A tutorial example (github gist?) would be terrific. I have been looking through the suggested unit test but it is not clear where the model is actually defined.
What would be even better is to understand how to call HMC with state; i.e., continuing the chain for N iterations w/o repeating the burn-in.
Yes, a C interface makes it easy for nearly any other language to call stan but I’m not going to ask the devs to create or maintain such an interface. I am content with C++.
I also just talked to @betanalpha and he said it’s not hard to do. But I guess he had not tried to do it himself. :)
Take a look at this test code. It calls the default algorithm when sampling. This is directly in C++.
I take it you mean you have a way to evaluate the joint up to a proportionality constant, not the actual log posterior with all the normalizing constants. If you could do that, the problem is already solved.
To take advantage of your function and gradients evaluation, you’ll need to make it conform to a model concept:
I think the algorithms require less than what’s described, so ask when you get there.
This is straightforward, but start with the above first. There are also tests in c++ to show how to do this, but it’s a little more involved.
Take a look at this test code. It calls the default algorithm when sampling. This is directly in C++.
Did you forget to paste some code here?
I take it you mean you have a way to evaluate the joint up to a proportionality constant, not the actual log posterior with all the normalizing constants. If you could do that, the problem is already solved.
Of course I just have likelihood and prior, not the evidence.
To take advantage of your function and gradients evaluation, you’ll need to make it conform to a model concept:
Thanks for the link to the model concept. So if I have all these functions implemented, I would have access to all the same functionality as a model defined in the stan language? Cool
Sorry about that – typed that on my phone and the link didn’t make it: https://github.com/stan-dev/stan/blob/develop/src/test/unit/services/sample/hmc_nuts_diag_e_adapt_test.cpp
Cool. Sorry about being pedantic, but terminology / notation matters when trying to help.
Yes! If you run that test (or run any program), take a look at the C++ class that’s generated. It follows the model concept.
At some point in the future, we’ll trim down what’s needed, but we haven’t had many people wanting to use the C++ API. Oh yeah… don’t trust the interface to be completely stable. We won’t change everything around, but the interface may change over different versions.
The key is to implement the log_prob function correctly: Home · stan-dev/stan Wiki · GitHub You’ll need to return a stan::math::var that’s constructed in such a way that it has the value and gradients are the returned stan::math::var. One way to do this is to use precomputed_gradients in the math library. There are other ways to do it. If you get stuck, reach out – probably on a new thread on Discourse.
Thank you very much for the quick reply. Now I need to sit and build the prototype. If I get stuck (prior for that quite high), I will open a new thread.
Does anyone have an example of implementations like these, yet?
I have a similar situation. I think it would be incredibly useful to pass a likelihood (and gradient) specified in C++ directly to the MCMC functions.
If I understand correctly, it is not possible to code arbitrary functions of the parameters in the interface languages like Python or R and pass it to Stan. And while it is often possible to do things inside the model block, sometimes it is cumbersome to write and test all functions in there instead as standalone functions first.
Even though it requires less amicable C++ coding, I think that would allow a lot of flexibility for models that aren’t as “nice”, and the machinery is already there, anyway. Thanks.
This thread is a bit different I think – asking how to call the samplers themselves from C++. You can build C++ code into Rstan pretty easily though (and let Rstan orchestrate all the sampling and interfacing stuff).
Info on how to hook up Stan/math compatible C++ with Rstan is here: https://cran.r-project.org/web/packages/rstan/vignettes/external.html
Info on how to do custom gradients with Stan/math C++ is here: https://github.com/stan-dev/math/wiki/Adding-a-new-function-with-known-gradients
But to understand how the custom gradients link in with everything you should read the Stan Math library paper (not trying to just give you an exhaustive list of resources – this is worth the time it takes to read): https://arxiv.org/abs/1509.07164
This is some of my own code that does this if you want more examples: https://cdn.rawgit.com/bbbales2/stancon_2018/718654d2/rus.html
I understand RStan is maybe the main interface to Stan, but I’m not a big fan of R, and barely have any libraries installed to do anything beyond plotting heatmaps.
I had seen some of those links, although I have yet to read the paper on the math library (not sure I’ll have time soon, but it’s on my list).
I don’t think it’s all that different from what other people asked above. Unless you just really like C++ and want to be able to do the exact same thing RStan or PyStan do, the advantage of C++ is you would be able to start using Stan from somewhere downstream from the model specification, and possibly bypass anything it is not prepared to do.
For instance, if it can’t handle numerical solutions because I could provide a function that returns an array with the prediction, and use provide numerical gradients, or whatever other wacky handling of the model you may want to implement, at your own risk of messing it up.
Packages like PyMC3 apparently moved to be more Stan-like than its previous version, which allowed a more generic implementation of your f(\theta) and specification of the likelihood [say p(y|\theta) = Poisson(f(\theta))] - and that would be all you needed to feed to the sampler. But that just won’t do it anymore.
I would otherwise just write down the likelihood and gradients, code the model and sampler from scratch in any one language (whether it is Python, C++, Julia, or whatever) and run it.
I feel like that is overkill, and extra effort to implement a lesser version of what Stan does very well, so ideally I’d be able to join the community and take advantage of all the work you have already done. I’m not a developer or a statistician, so maybe some of this doesn’t really make sense to you. Thanks for the reading list, again.
How about Python?
C++ is my frenemy, but even I wouldn’t want to do statistical modeling using C++.
You can already start adding that kind of thing in C++ now. The modeling language lets you declare signatures in the functions block then link them in externally. They can compute whatever you want as long as it integrates into Stan’s autodiff framework. We’re trying to make this even easier, but we have a ways to go until Stan’s truly extensible.
Thanks, @Bob_Carpenter. I much prefer Python over R (for no concrete reason, probably), and I agree doing statistical modeling in C++ is less than ideal to say the least, but so far I see Stan’s interfaces as essentially data input. That leaves the modeling to be done in the model block or potentially in C++ directly. As it is, this is probably too much work and defeats the point of using Stan as opposed to just coding everything from scratch, and both are probably a lot more work than is warranted for a single model.
I have a bunch of custom functions, including gp kernels, which not only are hard to code in Stan (also with me being new to it), but may not work with the autodiff. In some cases I can write the gradients analytically, so in principle I could get around that.
Using arbitrary C++ functions to compute the likelihood/gradients should make Stan as flexible as anything else without having to implement an infinite number of extensions in Stan itself. It would also be useful to be able to do inference with numerical solutions/gradients – even if it’s not ideal, it’s that or MH.
Looking forward to greater extensibility in Stan anyway. Thanks.
This is possible in Stan though, unless I’m misunderstanding what you’re saying. This thread shows how to do it in CmdStan if that is preferable to Rstan: What is the best way to add a user function that is used in likelihood calculation? - #6 by philipm