2

I read here that it is possible (and I interpreted straightforward) to call Stan routines from a C++ program.

I have some complex log-likelihood functions which I have coded up in C++ and really have no idea how I could code them using the Stan language. Is it possible to call the Monte Carlo routines in Stan using the log-likelihood function I have already coded in C++? If so are there any examples of this?

It seems like quite a natural thing to do but I cannot find any examples or pointers as to how to do this.

Community
  • 1
  • 1
Jeff
  • 718
  • 8
  • 20

2 Answers2

2

I think your question is a bit different from the one you linked to. He had a complete Stan program and wanted to drive it from C++, whereas you are asking if you could circumvent writing a Stan program by calling an external C++ function to evaluate the log-likelihood. But that would not get you very far because you still have to pass in the data in a form that Stan can handle, declare to Stan what are the unknown parameters (plus their support), etc. So, I don't think you can (or should) evade learning the Stan language.

But it is fairly easy to expose a C++ function to the Stan language, which essentially just involves adding your my_loglikelihood.hpp file in the right place under ${STAN_HOME}/lib/stan_math_${VERSION}/stan/math/, adding an include statement to the math.hpp file in that subdirectory, and editing ${STAN_HOME}/src/stan/lang/function_signatures.h. At that point, your .stan program could look as simple as data { // declare data like y, X, etc. } parameters { // declare parameters like theta } model { // call y ~ my_logliklihood_log(theta, X) } But I think the real answer to your question is that if you have already written a C++ function to evaluate the log-likelihood, then rewriting it in the Stan language shouldn't take more than a few minutes. The Stan language is very C-like so that it is easier to parse the .stan file into a C++ source file. Here is a Stan function I wrote for the log-likelihood of a conditionally Gaussian outcome in a regression context: functions { /** * Increments the log-posterior with the logarithm of a multivariate normal * likelihood with a scalar standard deviation for all errors * Equivalent to y ~ normal(intercept + X * beta, sigma) but faster * @param beta vector of coefficients (excluding intercept) * @param b precomputed vector of OLS coefficients (excluding intercept) * @param middle matrix (excluding ones) typically precomputed as crossprod(X) * @param intercept scalar (assuming X is centered) * @param ybar precomputed sample mean of the outcome * @param SSR positive precomputed value of the sum of squared OLS residuals * @param sigma positive value for the standard deviation of the errors * @param N integer equal to the number of observations */ void ll_mvn_ols_lp(vector beta, vector b, matrix middle, real intercept, real ybar, real SSR, real sigma, int N) { increment_log_prob( -0.5 * (quad_form_sym(middle, beta - b) + N * square(intercept - ybar) + SSR) / square(sigma) - # 0.91... is log(sqrt(2 * pi())) N * (log(sigma) + 0.91893853320467267) ); } } which is basically just me dumping what could otherwise be C-syntax into the body of a function in the Stan language that is then callable in the model block of a .stan program.

So, in short, I think it would probably be easiest for you to rewrite your C++ function as a Stan function. However, it is possible that your log-likelihood involves something exotic for which there is currently no corresponding Stan syntax. In that case, you could fall back to exposing that C++ function to the Stan language and ideally making pull requests to the math and stan repositories on GitHub under stan-dev so that other people could use it (although then you would also have to write unit-tests, documentation, etc.).

Ben Goodrich
  • 4,870
  • 1
  • 19
  • 18
  • Thank you for the detailed answer, Ben. My problem is that the log-likelihood I have written comprises 100s of lines of C++ code and makes use of a number of C++ classes which I have written - so it feels wasteful to re-write the whole thing into the stan language just so the stan compiler can convert it back to C++. I was hoping there would be a way to write C++ functions for the priors and log-likelihood, and then pass them to the stan HMC function - is this at all possible? – Jeff Aug 07 '15 at 11:16
  • Maybe you could write a .stan program that fills in all the blocks except `model`, parse that into C++ with CmdStan, inherit from the class that gets created, overwrite the `log_prob` method in the inherited class to call your C++ function (which needs to be templated so it can deal with Stan's custom scalar types), overcome whatever unanticipated complications arise, compile the modified C++ source, and hope for the best. But I still think it would be better to take one of the routes mentioned in my main answer. Whatever time is wasted on this Stan project will serve you well on the next one. – Ben Goodrich Aug 07 '15 at 15:10
  • I think I will write the model in the .stan file - it seems a lot more straightforward! Many thanks for all your help, Ben. – Jeff Aug 07 '15 at 15:28
2

Upon further review (you may want to unaccept my previous answer), you could try this: Write a .stan program with a user-defined function in the functions block that has the correct signature (and parses) but basically does nothing. Like this functions { real foo_log(real[] y, vector beta, matrix X, real sigma) { return not_a_number(); // replace this after parsing to C++ } } data { int<lower=1> N; int<lower=1> K; matrix[N,K] X; real y[N]; } parameters { vector[K] beta; real<lower=0> sigma; } model { y ~ foo(beta, X, sigma); // priors here } Then, use CmdStan to compile that model, which will generate a .hpp file as an intermediate step. Edit that .hpp file inside the body of foo_log to call your templated C++ function and also #include the header file(s) where your stuff is defined. Then recompile and execute the binary.

That might actually work for you, but if whatever you are doing is somewhat widely useful, we would love for you to contribute the C++ stuff.

Ben Goodrich
  • 4,870
  • 1
  • 19
  • 18
  • I just about had this working over the weekend but I need to work on other things for a while now :( I am trying an inhomogeneous Poisson process model where the instantaneous rate is rather complex and so is the integral in the likelihood function. It is quite specific so I don't think it would be a worthy contribution, unfortunately. Unless perhaps I can think of a way to generalise the model a bit more - I will have a think. Thanks again. – Jeff Aug 10 '15 at 14:34
  • I got this to work for some basic examples, thanks again for your solution! I was wondering if you think anyone would be interested in some examples and also where I could post it? – Jeff Aug 31 '15 at 21:39
  • Perhaps the GitHub wiki? ( https://github.com/stan-dev/stan/wiki ) I don't know if it will let you create a new page but one of us could push it for you or we could probably change the permissions. – Ben Goodrich Sep 01 '15 at 23:42
  • I've put it all on GitHub now (https://github.com/jeffpollock9/StanExample). The example I made which samples the posterior of mu and sigma using 50 observations of y ~ N(mu, sigma) is around twice as slow as the normal Stan equivalent. I'm guessing because I didn't code any gradients? – Jeff Sep 02 '15 at 11:27