Search code examples
rstatisticssimulationstanrstan

Simulations of Exponential Random Variables in Stan (RStan Package/Interface)

I'm trying to do simulations of exponential random variables using RStan code.

The R code for the simulated exponential random variables looks like this:

A <- rexp(1000, 4)
B <- rexp(1000, lambda2)
C <- rexp(1000, lambda3)
D <- rexp(1000, lambda4)
E <- rexp(1000, lambda5)
F <- rexp(1000, lambda6)

A + B = AB
C + D = CD

mean(AB)
mean(CD)

As you can see, it is very straightforward for me to simulate exponential random variables in R. Take A <- rexp(1000, 4) as an example: It generates a random sample of 1000 outcomes, where lambda = 4. I can then do statistical analysis on these simulated values (finding means, etc.).

Now I want to do the same in using RStan. Using RStudio's R notebook feature, one can "insert" different types of code, one of which is Stan:

enter image description here

I have the following Stan and R code:

{stan output.var="exponential"}
generated quantities{
  real total;
  real number;

  number = 0.0;
  total = 0.0;
  while(total < 1000) {
    total += exponential_rng(4);
    number += 1.0;
  }
  number -= 1.0;
}

{r}
simul2 <- sampling(exponential, algorithm="Fixed_param")

{r}
print(simul2, pars=c("number"), digits = 5)

However, when I execute this code (specifically, the print command), I get the following output:

enter image description here

I don't understand why I am getting such absurd statistics (4000 mean!?). My goal is to be able to get the same statistics as had I done the analysis in R, as I demonstrated above; in other words, my goal is to get the same values as when If I had done A <- rexp(1000, 4), in this specific case, and X <- rexp(1000, lambda) more generally.

Obviously, my Stan code is incorrect, so I would greatly appreciate it people could please take the time to explain the correct way to use it.

Solution

  • That is telling you the posterior distribution of number, which is just your counter. Your Stan code is doing something different than your R code, but it is not so difficult to do the analogue of your R code if you pass N and the lambdas as data:

    data {
  int<lower=1> N;
  real<lower=0> lambda2;
  real<lower=0> lambda3;
  real<lower=0> lambda4;
}
generated quantities {
  real mean_AB;
  real mean_CD;
  {
    vector[N] A;
    vector[N] B;
    vector[N] C;
    vector[N] D;
    for (n in 1:N) {
      A[n] = exponential_rng(4);
      B[n] = exponential_rng(lambda2);
      C[n] = exponential_rng(lambda3);
      D[n] = exponential_rng(lambda4);
    }
    mean_AB = mean(A + B);
    mean_CD = mean(C + D);
  }
}