Sampling from inverse gamma distribution by specifing mean and standard deviation

I am writing a function to draw samples from an inverse gamma distribution by specifing mean and standard deviation.

Here is the code.

function x = inv_gamma_rnd_mean_std(mu, sig)
    % mu is the mean, sig is the standard deviation.
    % shape (a) and scale (b) parameters
    a = mu^2/sig^2+2;
    b = mu*(a-1);
    x = 1/gamrnd(a,1/b);
end

The code looks problematic because my Monte Carlo test returns a standard deviation far from the specified one. For instance, I expected std(tmp) becomes arount 2.0 in the code below, but it is far from 2.0

tmp=zeros(10000,1);
for idx = 1:10000
    tmp(idx) = inv_gamma_rnd_mean_std(0.1,2.0);
end
mean(tmp)

ans = 0.0973

std(tmp)

ans = 0.1528

I could not find what was going wrong in my function. Thank you for your support in advance.

Answer 1

[Torsten]

It works for mu = 0.5, sig = 0.5, e.g. . It seems that your distribution parameter "a" is too close to the critical value of a = 2 to give good results for the random number generation (for a <= 2, the variance is undefined).

rng('default')
% mu is the mean, sig is the standard deviation.
% shape (a) and scale (b) parameters
mu = 0.1;
sig = 2.0;
a = mu^2/sig^2+2;
b = mu*(a-1);
tmp = 1./gamrnd(a,1/b,[10000,1]);
hold on
histogram(tmp,'Normalization','pdf')
xlim([0 1])
x = 0:0.001:5;
f = b^a/gamma(a)*(1./x).^(a+1).*exp(-b./x);
plot(x,f)
hold off

●

[Roy0014]

Thank you, Torsten, this makes sense! I appreciate it.