# Using stan as a random number generator

#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

, 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]`.

#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.