igamma {pscl}R Documentation

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 α>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 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_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 < α < infty, ,, 0 < β < 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 jackman@stanford.edu

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)

[Package pscl version 1.03 Index]