How to generate random number from inverse gamma distribution in an extreme conditions

hello everyone

I want to get the sample from the inverse gamma distribution as below

x~IG(54,0.004)

but when I try to use the reciprocal relationship between gamma distribution and inverse gamma distribution, like

first, generate a random number from t~G(54,0.004), then set x=1./t, and the result is:

3.66281673846745 4.15049653026671 5.59965910607058

the matlab code is:

tmp=gamrnd(54,0.004,1,3);
res=1./tmp;

the result seems to have some problems, I also test the sample in R, the result is below

7.432166e-05 6.758658e-05 6.506619e-05

the R code is:

rinvgamma(3,54,0.004)

Is this the internal algorithm problem with Matlab in handling such extreme condition?

Answer 1

[Shep Bryan]

Accoding to wikipedia, if X ~ InvGamma(a,b) then 1/X ~ Gamma(a,b), therefore to sample from the inverse gamma you simply take the inverse of a random variable sampled from the gamma distribution.

Answer 2

[Matteo Iacopini]

I don't you whether you have already solved the issue, however you can sample from an Inverse Gamma by drawing values from a Gamma and then takinng the inverse. Note that if X ~ Ga(a, b), where a= SHAPE and b= RATE (i.e., inverse of SCALE), then Y = 1/X ~ IG(a, b).

The point here is that MATLAB implements the SHAPE-SCALE parametrization of the Gamma distribution, not the SHAPE-RATE (which is the one you are looking for, see the relationship above).

Therefore, all you need to do to sample from IG(a,b) is

x = gamrnd(a, 1./b);
y = 1./x;

the first line generates from a Gamma with SHAPE a and RATE b, the second inverts to obtain Inverse Gamma random variates.