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?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].
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, ]
[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.