inverse-Gamma distribution
Density, distribution function, quantile function, and highest density region calculation for the inverse-Gamma distribution with parameters alpha and beta.
densigamma(x,alpha,beta) pigamma(q,alpha,beta) qigamma(p,alpha,beta) rigamma(n,alpha,beta) igammaHDR(alpha,beta,content=.95,debug=FALSE)
x,q |
vector of quantiles |
p |
vector of probabilities |
n |
number of random samples in |
alpha,beta |
rate and shape parameters of the inverse-Gamma density, both positive |
content |
scalar, 0 < |
debug |
logical; if TRUE, debugging information from the search for the HDR is printed |
The inverse-Gamma density arises frequently in Bayesian analysis of normal data, as the (marginal) conjugate prior for the unknown variance parameter. The inverse-Gamma density for x>0 with parameters α>0 and β>0 is
(beta^alpha)/Gamma(alpha) x^(-alpha-1) exp(-beta/x)
where Γ(x) is the gamma function
Gamma(a) = int_0^infty t^(a-1) exp(-t) dt
and so ensures f(x) integrates to one. The inverse-Gamma density has a mean at beta/(alpha-1) for alpha>1 and has variance beta^2/((alpha-1)^2 (alpha-2)) for alpha>2. The inverse-Gamma density has a unique mode at beta/(alpha+1).
The evaluation of the density,
cumulative distribution function and quantiles is done by
calls to the dgamma, pgamma and igamma
functions, with the arguments appropriately transformed. That
is, note
that if x ~ IG(alpha,beta then 1/x ~ G(alpha,beta).
Highest Density Regions. In general, suppose x has a density f(x), where x \in Θ. Then a highest density region (HDR) for x with content p \in (0,1] is a region (or set of regions) \mathcal{Q} \subseteq Θ such that:
int_Q f(x) dx = p
and
f(x) > f(x*) for all x in Q and all x* not in Q.
For a continuous, unimodal density defined with respect to a single parameter (like the inverse-Gamma case considered here with parameters 0 < α < ∞, \,\, 0 < β < ∞), a HDR region Q of content p (with 0 < p < 1) is a unique, closed interval on the real half-line.
This function uses numerical methods to solve for the
boundaries of a HDR with content p for the
inverse-Gamma density, via repeated calls the functions
densigamma, pigamma and
qigamma. In particular, the function uniroot is used to
find points v and w such that
f(v) = f(w)
subject to the constraint
int_v^w f(x; alpha, beta) d theta = p.
densigamma gives the density, pigamma the
distribution function, qigamma the quantile function,
rigamma generates random samples, and igammaHDR gives
the lower and upper limits of the HDR, as defined above (NAs if
the optimization is not successful).
Simon Jackman simon.jackman@sydney.edu.au
alpha <- 4
beta <- 30
summary(rigamma(n=1000,alpha,beta))
xseq <- seq(.1,30,by=.1)
fx <- densigamma(xseq,alpha,beta)
plot(xseq,fx,type="n",
xlab="x",
ylab="f(x)",
ylim=c(0,1.01*max(fx)),
yaxs="i",
axes=FALSE)
axis(1)
title(substitute(list(alpha==a,beta==b),list(a=alpha,b=beta)))
q <- igammaHDR(alpha,beta,debug=TRUE)
xlo <- which.min(abs(q[1]-xseq))
xup <- which.min(abs(q[2]-xseq))
plotZero <- par()$usr[3]
polygon(x=xseq[c(xlo,xlo:xup,xup:xlo)],
y=c(plotZero,
fx[xlo:xup],
rep(plotZero,length(xlo:xup))),
border=FALSE,
col=gray(.45))
lines(xseq,fx,lwd=1.25)
## Not run:
alpha <- beta <- .1
xseq <- exp(seq(-7,30,length=1001))
fx <- densigamma(xseq,alpha,beta)
plot(xseq,fx,
log="xy",
type="l",
ylim=c(min(fx),1.01*max(fx)),
yaxs="i",
xlab="x, log scale",
ylab="f(x), log scale",
axes=FALSE)
axis(1)
title(substitute(list(alpha==a,beta==b),list(a=alpha,b=beta)))
q <- igammaHDR(alpha,beta,debug=TRUE)
xlo <- which.min(abs(q[1]-xseq))
xup <- which.min(abs(q[2]-xseq))
plotZero <- min(fx)
polygon(x=xseq[c(xlo,xlo:xup,xup:xlo)],
y=c(plotZero,
fx[xlo:xup],
rep(plotZero,length(xlo:xup))),
border=FALSE,
col=gray(.45))
lines(xseq,fx,lwd=1.25)
## End(Not run)Please choose more modern alternatives, such as Google Chrome or Mozilla Firefox.