function y = dinvgamma(x, shape, scale)
% pdf of inverse gamma distribution 
% it's based on the dinvgamma function of MCMCpack in R
% x:  scalar 
% shape:  scalar shape parameter
% scale:  scalar scale parameter

alpha = shape ;
beta = scale ;
log_density = alpha * log(beta) - gammaln(alpha) - (alpha + 1) * log(x) - (beta/x) ;
y = exp(log_density) ;

end