Hi,

suppose I have a convolution of an exponential random variable X with unknown parameter \lambda and a normal random variable Y with unknown mean \mu (the latter for the sake of simplicity with known unit variance). The two unknown parameters should be estimated from i.i.d data.

Now, I know that in this case the convolution can be written in an explicit form. However, I am interested in learning something about the implementation of integrals in STAN. Hence, my goal is to implement this model in STAN without usage of the closed form expression; only using the expression:

f(z|\lambda, \mu)=\int_0^{\infty}\lambda e^{-\lambda v}\frac{1}{\sqrt{2\pi}}e^{-0.5*(z-v-\mu)^2}\mathrm{d}v

I started off by writing a function which approximates the integral via importance sampling, i.e. drawing N exponential random variables v_1, \dots v_N with parameter \lambda and then computing the average of \frac{1}{\sqrt{2\pi}}e^{-0.5*(z-v_i-\mu)^2}. However, I realized that I can only accomplish this in STAN by declaring the function with the rng suffix. But if I follow this declaration, I can only use this function in the generated quantities block.

According to the manual:

â€śFunctions that call random number generating functions in their bodies must have a name that ends in _rngâ€¦user-defined functions with names that end in _rng may be used only in the generated quantities block and transformed data block, or within the bodies of user-defined functions ending in _rngâ€ť

So my question: How is it possible to implement such a density in STAN? (via importance sampling or numerical integration in general) Is there perhaps a way to write this function in another language (R, Python etc.) and use this within STAN?