How to do integral in stan

Modeling
1

Hi, could u please share the solution of doing integral in stan?

Specifically, my problem is as follows:
y = E(exp(kx))+b;
where E(exp(kx)) is the integral of exp(kx) in infinite horizon and x~normal(0,1).
k and b are parameters I wanna estimate.

How to do integral in stan? Thanks!

2

any suggestions?

3

Hi !

First, please refrain from re-uping your topic after a so short time, this is considered rude in many forums.

Second, I guess what you want is actually the expectancy of exp(k*x) when x~Normal(0,1). If so, this is analytically solvable, and equals exp(k²/2).

5

Thanks for your reply

I apologize for re-uping. i will keep it in mind
Actually I wanna ask the universal way of doing integral (and also taking expectancy). Because my real problem is not exp(kx), but one that is way more complicated and have no analytical solution. Exp(kx) is just an example.
So could u share some universal way of taking expectancy? thanks!

6

What you want to do is then numeric integration, the general way to do this is to generate a series of x_i distributed according to the Normal law, and for each of those, compute the function you are interested in (here exp(k*x)) and then compute the average of these samples of the function value.

If this is too slow to converge, I would suggest that you have a look at methods like importance sampling that may help you.

7

thanks. yes, the idea is clear, but I wanna know the command and grammar of numeric integration in Stan. It seems that the only way to do integration in Stan is through ODE (maybe I’m wrong). But this is undoable in my case. Besides, I need to do inference based on this integration.

8

thanks. yes, the idea is clear, but I wanna know the command and grammar of numeric integration in Stan. It seems that the only way to do integration in Stan is through ODE (maybe I’m wrong). But this is undoable in my case. Besides, I need to do inference based on this integration.

9

Take a look at: PDA in Stan?

Maybe in the next release or the other we will also have https://github.com/stan-dev/math/pull/566

10

thanks, Mr.randommm. The topic in that thread is about one dimensional integration. But in my case I really need to do 6 dimensions integration. So the ODE method is undoable. Do u know how to do multiple integration in Stan?

12

thanks, Mr.randommm. The topic in that thread is about one dimensional integration. But in my case I really need to do 6 dimensions integration. So the ODE method is undoable. Do u know how to do multiple integration in Stan?

13

If want you really need is 6D numerical integration inside Stan (which is probably the case), than sadly I cannot help you much, but maybe what you want to do is way more simple and can be written directly on Stan. That is, maybe in your problem you want to use MCMC to do numerical integration (which is far far more common) and not do numerical integration inside an MCMC algorithm [since I think I recently saw some confusion on those two matters, maybe it’s worth adding this distinction to the manual once the numerical integration algorithm is merged].

Something like:

model {
x ~ normal(0, 1);
k ~ ...
}

then outside Stan get the expected value from samples obtained for x and k

   for (i in 1:nsim)
     v[i] = exp(k[i] * x[i])

  expected = mean(v)

Well, even if it that’s not what you want maybe it’s worth pointing this out for future readers of this thread.

14

Can you post your exact problem? I think we might be better able to help you if we know what you are trying to do.

15

thank you. but I think my problem is doing integration inside MCMC algorithm.

16

My exact problem is as follows:

y = E{exp(k1x1)/[exp(k1x1)+exp(k2x2)+exp(k3x3)]}+b

where E{ } is expectancy, x1, x2, x3 ~ normal(0,1), k1, k2, k3, b are parameters I wanna estimate through MCMC. k1, k2, k3 are uninformative priors and b follows normal(0,1)

17

and Y?

When you say “E{} is expectancy” do you mean?
E{exp(k1x1)/[exp(k1x1)+exp(k2x2)+exp(k3x3)]} =
\int_0^\infty \int_0^\infty exp(k1x1)/[exp(k1x1)+exp(k2x2)+exp(k3x3)] dx1dx2dx3?

You said you need to do “6 dimensions integration,” but that is only 3.

You will be running into identifiability issues if you are using weekly informative priors on ks.

19

yes, I omit x4, x5, x6

20

You have data y that are expected values of a complicated function of some random x vectors and some unknown parameters

y = E(f(x1,x2,x3…, k1,k2,k3…))

and you know the distribution of x1, x2 etc and you have priors for k1,k2, etc and you want to infer a posterior for k1,k2,k3 and although you’ve somehow observed this expectation… the expectation is also completely un-calculatable in analytical form.

is that the sum of it?

How did you get this y data?

21

If you really have to do this I’d suggest generating a big sample for x1, x2, x3, etc in transformed data, and then in transformed parameters, calculate the exp(k1*x1)/… stuff into a big vector, and then do

expectation = mean(big_vector)

and in model:

y-expectation ~ distrib_for_b(…)