function y = rinvgamma(n, shape, scale)
% random variable from inverse gamma distribution 
% it's based on the rinvgamma function of MCMCpack in R
% n:  number of observations
% shape:  scalar shape parameter
% scale:  scalar scale parameter

y = 1/gamrnd(shape, 1/scale, 1, n) ;
% y = 1/gamrnd(shape, scale, 1, n) ;

end