I have a problem which requires the implementation of a user defined probability function. However, the function is not analytically available and requires (numerical) integration. My question is: What is the best way to approach this problem in Stan:

approximate the integral via quadrature rule?

call some numerical integration routine (are these available in Stan?)

generate within the function random numbers and approximate the integral via Monte-Carlo methods?

As I am not very familiar with Stan it would be great if there were some example code which demonstrate the implementation of some of these approaches (in an abitrary modeling framework).

I don’t know if I’ve seen an example where option three works well but you can approximate it as a sum if you’re in one dimension and get on with the modeling. Lots of other choices will have a bigger impact (e.g.-number of leapfrog strips per iteration). There is a 1d integrator in the math lib but you’d have to check if it found it’s way into the interfaces, it came up on the forums a while back. The joint model in rstanarm uses quadrature so that’s an option too.

Thank you very much. As far as I know there is no simple 1d integrator directly available.
It would be most useful if you could point me to some (toy) stan code which implements a 1d integration because I would profit the most from a simple example. (using any of the three approaches)