# inverse-Gamma distribution Density, distribution function, quantile function, and highest density region calculation for the inverse-Gamma distribution with parameters `alpha` and `beta`. ```r 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)%(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%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 c("$p\n$", "$ \\in (0,1]$") is a region (or set of regions) c("$\\mathcal{Q}\n$", "$ \\subseteq \\Theta$") such that: $$ \int_\mathcal{Q} f(x) dx = p%int_Q f(x) dx = p $$ and $$ f(x) > f(x^*) \, \forall\ x \in \mathcal{Q},x^* \not\in \mathcal{Q}.%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 c("$0 <\n$", "$ \\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.%int_v^w f(x; alpha, beta) d theta = p. $$ ## Returns `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 (`NA`s 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 simon.jackman@sydney.edu.au ## See Also `gamma`, `dgamma`, `pgamma`, `qgamma`, `uniroot` ## Examples ```r 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) ```

igamma function