function[g]=invgampdf(x,a,b)
%
% g=invgampdf(x,a,b)
% evaluates the inverse gamma pdf at value x, with shape a and scale b,
% using Gelman's notation
% x can be a scalar or column vector

% g=((b^a)/gamma(a))*(x.^-(a+1)).*(exp(-b./x));
% yielded Inf and NAN's

% g=exp((a*log(b))-gammaln(a))*(x.^-(a+1)).*(exp(-b./x));

% even better:



g=exp(((a*log(b))-gammaln(a))-(a+1)*log(x)-(b./x));