Using stan as a random number generator

Modeling
1

As a toy example that might help reinforce my understanding of the data structure in stan, I tried the following code

"data{}

parameters { vector<lower=0>[100] temp_rinvg;
}

model {
temp_rinvg ~ inv_gamma(1,1);
}"

I save the above stan codes as a string in R and name it as ing_test_code,

then I call the above stan codes in R,

ing_test <- stan(model_code = ing_test_code, chains = 1)

What I got was

stan_puzzle
stan_puzzle.jpg977Ă—822 708 KB
, so is each temp_rinvg[i] a vector (length = 100) or a single number?

2

In one sense, temp_rinvg[i] is a “single number,” since in the parameters block, you declared temp_rinvg as a vector parameter with 100 elements. However, each temp_rinvg[i] is a random number that has a range of possible values, where the distribution of values is governed by the inverse gamma distribution. Stan is (more or less) sampling 1000 (post-warmup) possible random values for each temp_rinvg[i].

1 Like
3

Thank you. So if what I want is a random sample from inv_gamma(1,1), I could simply stack the column “mean” together (i.e. 5.76, 5.73, 5.47, 6.92, 10.70, 51.92, …) ?

4

No, that’s the mean of several draws, which will have the same expectation but lower variance than the actual draws. Instead, you want to use the extract() function and take a single draw.

> extract(ing_test)$temp_rinvg[1, ]
  [1]  2.5555562  2.5627603  4.5832294  1.6346871  1.4185565  0.9012479  0.2097981  0.8691214  0.4677514
 [10]  0.9363055  1.2994146  6.1237233 11.9161933  6.9381447  0.4072198  4.1231350 19.5466629  8.8823631
 [19]  0.2421334 57.0499708  0.6370643  1.2703736  1.2794420  0.8445436  0.4417089 58.6737416  0.5681631
 [28]  1.2310256  1.2259458  1.4682253  2.6714699  4.0564212  1.3518744  9.2318232 14.8092812  0.6626966
 [37]  3.2753341  3.4416742  1.4426440  1.6569463  0.6947641  0.7250604  1.2911643  2.2276262  0.1892672
 [46] 16.7677540  3.4273376  0.7333991  0.4308113  0.6837480 20.3680304  1.2686530  2.2916267  2.0116976
 [55]  0.3112171  1.9261007  1.4042962  0.9594311  0.9798000 93.0867078  9.0961214  1.5949804  1.8458280
 [64] 40.7629444  0.8932531  1.6885196  0.4620151  1.2741189  1.0047026  0.8385737  0.8121868  1.3590331
 [73]  0.3512074 21.2992401  0.6695918  0.3317261 13.6922654  0.9992000  2.4554088  1.9495224  0.3835895
 [82] 13.1012419 31.7251914  2.1527747  0.3718639  0.9745706  0.9929655  0.8244007  0.9086866  6.6010477
 [91]  4.1552772  0.8379916  1.0556177  0.7874264 13.4156752 11.0079549  0.5005472  5.6516660  4.7824275
[100]  0.8673477

If all you want is an RNG, it’s much much more efficient to draw them exactly in the generated quantieis block using the inv_gamma_rng function.

The other problem you’re running into is that inv_gamma(1, 1) doesn’t have a mean; see the Wikipedia page on inverse gamma. The problem is that if u ~ gamma(1, 1), then 1 / u ~ inv_gamma(1, 1). It’s easy to get draws from gamma(1, 1) that are near zero, which yield very large values when inverted.