# Using stan as a random number generator

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> 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?

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

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, …) ?

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, ]
  2.5555562  2.5627603  4.5832294  1.6346871  1.4185565  0.9012479  0.2097981  0.8691214  0.4677514
  0.9363055  1.2994146  6.1237233 11.9161933  6.9381447  0.4072198  4.1231350 19.5466629  8.8823631
  0.2421334 57.0499708  0.6370643  1.2703736  1.2794420  0.8445436  0.4417089 58.6737416  0.5681631
  1.2310256  1.2259458  1.4682253  2.6714699  4.0564212  1.3518744  9.2318232 14.8092812  0.6626966
  3.2753341  3.4416742  1.4426440  1.6569463  0.6947641  0.7250604  1.2911643  2.2276262  0.1892672
 16.7677540  3.4273376  0.7333991  0.4308113  0.6837480 20.3680304  1.2686530  2.2916267  2.0116976
  0.3112171  1.9261007  1.4042962  0.9594311  0.9798000 93.0867078  9.0961214  1.5949804  1.8458280
 40.7629444  0.8932531  1.6885196  0.4620151  1.2741189  1.0047026  0.8385737  0.8121868  1.3590331
  0.3512074 21.2992401  0.6695918  0.3317261 13.6922654  0.9992000  2.4554088  1.9495224  0.3835895
 13.1012419 31.7251914  2.1527747  0.3718639  0.9745706  0.9929655  0.8244007  0.9086866  6.6010477
  4.1552772  0.8379916  1.0556177  0.7874264 13.4156752 11.0079549  0.5005472  5.6516660  4.7824275
  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.