Samples from scaled inverse chisquare distribution

Question

I want to generate sa scaled-inv-chisquared distribution in R. I know geoR have a R function for generating this. But I want to use gamma-distribution to generate this. I think this two are equivalent:

X ~ rinvchisq(100,  df=d, scale=s)
1/X ~ rgamma(100, shape=d/2, scale=2/(d*s))

isn't it? Can there be any numerical problem due this due to extreme values?

Answer 1

More specifically you would need X <- rinvchisq(...) and X <- 1/rgamma(...) (the ~ notation works this way in programs such as WinBUGS, and in statistics notation, but not in R). If you look at the code of geoR::rinvchisq, the relevant part is just

return((df * scale)/rchisq(n, df = df))

so if you have problems taking the reciprocal of very large or small chi-squared deviates you'll be in trouble anyway (although rchisq is internally using .External(C_rchisq, n, df), which falls through to C code, presumably for efficiency in this special case, rather than calling rgamma). If I were you I would go ahead and superimpose densities of some test samples just to make sure I hadn't screwed up the arithmetic or parameterization somewhere ...

For what it's worth there are also rinvgamma() functions in a variety of packages (library(sos); findFn("rinvgamma"))