igamma {pscl}

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)

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

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.

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

Note

The densigamma is named so as not to conflict with the digamma function in the R base package (the derivative of the gamma function).

Author(s)

Simon Jackman [e-mail]

See Also

gamma, dgamma, pgamma, qgamma, uniroot

Examples

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)