igamma: inverse-Gamma distribution

Description

Density, distribution function, quantile function, and highest density region calculation for the inverse-Gamma distribution with parameters alpha and beta.

Usage

densigamma(x,alpha,beta)
  	pigamma(q,alpha,beta)
	qigamma(p,alpha,beta)
	rigamma(n,alpha,beta)
	igammaHDR(alpha,beta,content=.95,debug=FALSE)

Value

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).

Arguments

x,q: vector of quantiles
p: vector of probabilities
n: number of random samples in rigamma
alpha,beta: rate and shape parameters of the inverse-Gamma density, both positive
content: scalar, 0 < content < 1, volume of highest density region
debug: logical; if TRUE, debugging information from the search for the HDR is printed

Author

Simon Jackman simon.jackman@sydney.edu.au

Details

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 $\alpha>0$ and $\beta>0$ is $$f(x) = \frac{\beta^\alpha}{\Gamma(\alpha)} x^{-\alpha-1} \exp(-\beta/x)$$ where $\Gamma(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 \sim IG(\alpha,\beta)$ then $1/x \sim G(\alpha,\beta)$.

Highest Density Regions. In general, suppose $x$ has a density $f(x)$, where $x \in \Theta$. Then a highest density region (HDR) for $x$ with content $p \in (0,1]$ is a region (or set of regions) $\mathcal{Q} \subseteq \Theta$ such that: $$\int_\mathcal{Q} f(x) dx = p$$ and $$f(x) > f(x^*) \, \forall\ x \in \mathcal{Q}, x^* \not\in \mathcal{Q}.$$ For a continuous, unimodal density defined with respect to a single parameter (like the inverse-Gamma case considered here with parameters $0 < \alpha < \infty, \,\, 0 < \beta < \infty$), 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.$$

Examples

Run this code

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)


if (FALSE) {
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)
}

Run the code above in your browser using DataCamp Workspace