Use a R function in .stan file

Modeling
1

My aim is to perform a bayesian analysis of my model that is written in R. I would like to use Stan to do so. More specifically, I would like to use the Rstan package in R.

My model is made of 128 differential equations which depends on 8 parameters \theta. It thus represents a vector dx of length 128. It is resolved with a C++ algorithm using the Rcpp package in R. This gives a vector x of length 128 for each time step. As there are 100 timesteps, this will give overall a matrix with 128 rows and 121 columns. We will next use the 1st, 13th, 25th, 37th, …, 121st elements of the first row of this matrix in order to compute the likelihood l. I will call this vector of length 11 \tilde{\lambda}. It will be compared to a vector of observations y of length 11 by assuming a Poisson distribution. The likelihood is thus:

L(\theta)=P(y_1|\lambda=\tilde{\lambda_1}) \cdots P(y_{11}|\lambda=\tilde{\lambda_{11}})
with P() a Poisson distribution.

I want now to do a bayesian analysis. But to do this, I have to specify my model in the .stan file. The problem is that it would be too complex to rewrite the whole model (the differential equations) in C++ in the .stan file. Do you have ideas how to specify my model in Stan without rewriting it in C++?

2

https://cran.r-project.org/web/packages/rstan/vignettes/external.html

3

Thanks. But, it means that I have to transform my R script that describe the model into a C++ script right?

4

Or a Stan function

5

The problem is that I don’t know how to tranform my R script either in C++ or in Stan. I knew how to export my R function into C++, but it was only within R and Rcpp package…

6

You can’t do anything in Stan without knowing a good bit about the language, the sampling algorithm, and the core ideas behind Bayesian inference.

7

Actually, I do understand the ideas behind Bayesian inference as well as the sampling algorithm. My question is not about that, it’s really about it’s application: What should I do to transform my R function that describes the model in order to be usable by Stan?

8

Port it to the Stan language or to C++. It will be more or less a line-for-line translation because all three languages are pretty similar.

9

The specification of a function in R is not sufficient for that function to be added to the Stan language. As @bgoodri alluded to, Stan utilizes gradient-based algorithms which means we need to be able to compute the Jacobian of partial derivatives for every function in Stan.

Either you have to specify that Jacobian yourself (see for example https://github.com/stan-dev/math/wiki/Adding-a-new-function-with-known-gradients) or you have to write your function using C++ templates so that our powerful automatic differentiation code can calculate the partials automatically.

10

Have you checked out Stan’s diff eq solver? It’s pretty easy to code up this kind of model. We often solve thousands of differential equations per likelihood evaluation in pharmacology models (though those can take several days to conerge).

Charles Margossian has some nice case studies he did for both StanCons on the web site. The basics are also in the manual.