invgamma: the inverse gamma distribution in R

Introduction

R is the programming lingua franca of academic statisticians and one of the most popular statistics and data science programming languages (R Core Team, 2016; Muenchen, 2023). In its base distribution it contains core functionality for all sorts of operations essential to a proper platform for statistical computing: reading and writing files, data manipulation, graphical capabilities, scientific computing data structures and algorithms (e.g., arithmetic, matrix algebra, sorting, and special mathematical functions), modeling capabilities (e.g., regression), and much more. R also has a thriving developer base creating contributed packages, which number in the tens of thousands at the time of this writing, hosted on R’s primary package repository, the Comprehensive R Archive Network (CRAN, https://cran.r-project.org).

One of the key components of R is a simple and efficient code base for working with probability distributions. R has dozens of built-in functions related to common probability distributions such as probability density functions (PDFs), cumulative distribution functions (CDFs), quantile functions (QFs), and random number generators (RNGs), typically implemented as compiled C executables. For a given family of probability distributions, these four functions follow a memorable naming pattern: [dpqr]<dist>(), where the letters d, p, q, and r signify the PDF, CDF, QF, and RNG of the distribution (respectively), and <dist> is the short-hand form of the family name that R uses. Thus, dnorm() is the PDF of the normal distribution; qbinom() is the quantile function of the binomial distribution; and rt() is the RNG for the $t$ distribution. We refer to these as the dpqr functions, and for brevity when referring to the collection of four for a particular family write ∗<dist>(), e.g., ∗gamma(). The specific parameterizations of these families are communicated in their corresponding documentation accessible with ?, e.g., ?dnorm. These functions are all provided by the stats package, one of R’s base packages that is loaded by default when an R session is initialized.

Conspicuously missing from the list, however, is the inverse gamma distribution and its related subfamilies: the inverse ${\chi }^2$ and inverse exponential distributions. In this short article we discuss invgamma, a low-dependency, efficient, and low-maintenance package available on CRAN that has provided the functions to fill this void for the past decade, with downloads totaling well over three hundred thousand (Csárdi, 2019). Although previously only a software implementation, we, the authors of invgamma, write this article now to reflect on the software, present the core mathematics behind it (some of which is sparsely documented elsewhere), and clarify a few design choices.

Getting invgamma

There are two ways to obtain invgamma. First, one can obtain it from CRAN using R’s built-in packaging tool install.packages():

Alternatively, one can install the most recent stable development version of the package using remotes (Csárdi et al., 2024):

The package is then loaded and attached to the current R session using the standard mechanism:

Mathematical underpinnings

Illustrative examples

In this section we briefly show how the functions can be used and demonstrate their correctness. We begin by visualizing the $Y\sim \mathrm{I}\mathrm{n}\mathrm{v}\text{-}\mathrm{G}\mathrm{a}\mathrm{m}\mathrm{m}\mathrm{a}\left(7,10\right)$ distribution (Wickham, 2009; Wickham et al., 2019). The corresponding PDF is included in Fig. 1.

This also illustrates that the function dinvgamma() is vectorized across its first argument, as R users expect. This is a really nice feature for likelihood estimation, where the likelihood function is evaluated at several points and multiplied. Typically computed on the log-scale, we can pass in the argument log = TRUE and apply sum() to affect a similar, but more numerically useful, result.

The consistency of pinvgamma() with dinvgamma() can be checked with the integrate() function, which applies numerical quadrature via adaptive Gauss–Kronrod quadrature. Here, we compute $F_Y(2)=P\left[Y\le 2\right]$ using both

integrate(dinvgamma()) and pgamma():

Since these values are equivalent, dinvgamma() and pinvgamma() are satisfying the expected relationship consistently. Checking qinvgamma()’s consistence is a simple step from here: we simply take the inverse of $P\left[Y\le 2\right]$, expecting 2.

The above illustrates that dinvgamma(), pinvgamma(), and qinvgamma() are consistent in the sense that they satisfy the probabilistic relationships demanded among their respective functions: they refer to the same probability distribution. To demonstrate that they are sampling from the correct distribution, we compute the same quantity as before, $F_Y(2)$, using Monte Carlo simulation via rinvgamma(). Since rinvgamma() simply uses 1/rgamma(), it samples from the correct distribution by definition; its numerics are addressed shortly.

Since 0.76209 is within the Monte Carlo error of 0.7621835 (the true value computed previously) when using 100,000 draws, it is clear that the implementations are correct. As more complex demonstrations, in Fig. 2 we superimpose a kernel density estimate of the 100,000 sampled values with the density given by dinvgamma() and present a QQ-plot of the draws, both of which strongly confirm the correctness of the sampler (Pedersen, 2024). Figure 2 is generated with the following code:

As a final note, like their ∗gamma() counterparts, ∗invgamma() functions accept both the rate (default) and scale parameters. Their consistency can be seen in the following example with dinvgamma():

The Inverse ${\mathit{\chi }}^2$ and Inverse Exponential Distributions

The mathematical arguments made in the previous section work for any continuous distribution with a smooth density defined on the positive real numbers. Particularly notable among this class of distributions are two commonly used families subsumed by the gamma distribution: the ${\chi }^2$ (chi-squared) distribution and the exponential distribution. Since these two are so closely related to the gamma distribution, we have implemented the entire dpqr line of functions for the inverses of both of these distributions in invgamma as well.

The chi-squared distribution with $\nu$ degrees of freedom, denoted ${\chi }_{\nu }^2$, is equal to the $\mathrm{G}\mathrm{a}\mathrm{m}\mathrm{m}\mathrm{a}\left(\alpha =\frac{\nu }{2},\lambda =\frac{1}{2}\right)$ distribution; it is the distribution of the sum of $\nu$ independent squared standard normal variates. In invgamma, the ∗invchisq() functions for the inverse ${\chi }^2$ distribution could be implemented through the ∗invgamma() functions, but we have chosen to model them off of the ∗chisq() functions in the same way that the ∗invgamma() functions are implemented using the ∗gamma() functions, for two reasons. First, the ∗chisq() functions allow for non-central ${\chi }^2$ distributions, which are not special cases of the inverse gamma but nevertheless can be obtained by logic similar to that of the previous section. Second, using the ∗chisq() functions insures that the ∗invchisq() functions are synchronized with base R’s implementation of the ${\chi }^2$ distribution, which is more natural than syncing them with the ∗gamma() implementations.

The exponential distribution, like the gamma distribution, also manifests in a rate and a scale parameterization; both are obtained by setting the gamma distribution’s shape parameter $\alpha$ to 1. However, only the rate parameterization is implemented in R’s ∗exp() functions. invgamma stays true to this convention in its ∗invexp() functions, which have a rate argument defaulted to 1. And, like the ∗invchisq() functions, the ∗invexp() functions wrap the ∗exp() functions, not the ∗invgamma() functions. As in the inverse ${\chi }^2$ case, this allows the functions to naturally sync with the ∗exp() functions.

A note on parameterization and argument names

The point of invgamma is to provide a simple, low-maintenance implementation of the d, p, q, and r functions of the inverse gamma, chi-squared, and exponential distributions. To do this, it uses the transformation theorem along with the implementations of the corresponding functions for the gamma distribution, which are already implemented in stats as dgamma() and so on. The argument names were intended to reflect this fact. Thus, rate in dinvgamma() refers to the rate parameter of dgamma(), the rate parameter of the gamma distribution.

This forms a kind of linguistic problem: the parameter rate of dgamma() is called rate because it is a rate parameter, which has a technical definition outside of the gamma/inverse gamma context. The inverse gamma distribution has a rate parameter in this other sense as well, but that rate parameter is not the same as the rate parameter of the gamma distribution. Worse: the rate parameter of the inverse gamma is the inverse of the rate parameter of the gamma distribution, which is the scale parameter. Thus, while (for example) dinvgamma() admits a parameter called rate, it is referring to the rate of the rate of the associated gamma distribution, not the rate parameter of the inverse gamma distribution, and this has naturally led to some confusion.

An analogy of the naming convention oddities may be helpful here to explain how these clashes occur. The lognormal distribution $LogN(\mu ,{\sigma }^2)$ is, by definition, the distribution of a random variable whose log is a normal $N(\mu ,{\sigma }^2)$ distribution, i.e., $Y=e^X$ with $X\sim N(\mu ,{\sigma }^2)$: a lognormal is the exponential of a normal. Throughout statistics $\mu$ conventionally denotes the mean of a random variable or distribution and ${\sigma }^2$ its variance. However, although the first parameter of the lognormal is conventionally denoted $\mu$, the mean of a lognormal random variable is not this $\mu$, and similarly ${\sigma }^2$ is not its variance. Users simply have to know that they are the mean and variance of the “distribution behind the distribution.” This is a bit confusing and goes against the convention that $\mu$ ( ${\sigma }^2$) is the mean (variance) of the distribution. It is a clash of conventions that the community deals with and accepts as a fact of life, albeit in some cases begrudgingly.

Base R implements the lognormal distribution as dlnorm() with formal arguments meanlog and sdlog to try to clarify this distinction, but this was not done when initially writing invgamma, and this has been a source of confusion. When writing invgamma, our original intention was to write other similar packages that do the same thing: provide tools for other common distributions not in R that are simple transformations of distributions that are. The convention was to maintain the parameters of the “distributions behind the distributions.” This makes sense, but it does make this part messy.

Unfortunately it is not something that can be changed at this point, as the package has been in use for nearly a decade with several hundred thousand downloads. Similar to the base distribution of R, we have adopted a long-term mindset lending significant deference to software interface continuity, and thus have opted to keep the naming the way it is. To avoid mishaps, invgamma issues a package startup message alerting users to this issue. It is possible that if a responsible way can be determined to gradually transition to a more transparent naming convention, the code base will shift, alerting users to the change.

Related works

As the inverse gamma distribution has important applications all over statistics, similar implementations of it have arisen in other user-contributed R packages. The two most prominent of these are the MCMCpack and actuar packages (Martin, Quinn & Park, 2011; Dutang, Goulet & Pigeon, 2008; Dutang, 2016). MCMCpack is designed for Markov chain Monte Carlo algorithms in Bayesian statistics, and actuar is designed for actuarial science in R. Both are large, domain-specific packages that are strange to list as dependencies in packages with applications in unrelated domains, e.g., pharmaceutical statistics. Moreover, both of them have little quirks that make them less than ideal: MCMCpack only implements dinvgamma() and rinvgamma(), with restricted argument lists and a hard-coded version of the PDF; and actuar uses custom C-level implementations of the functions.

While not dramatically different than these, invgamma is preferable to these two packages as a general purpose implementation of the inverse gamma distribution for three reasons. First, invgamma is a small package, domain neutral, and has a name that reflects what it is, which is ideal for loading in other packages and for finding it. Second, invgamma leverages mathematical relationships between the gamma distribution and the inverse gamma distribution to allow the ∗invgamma() functions to wrap the ∗gamma() functions. Consequently, the ∗invgamma() functions are both efficient and automatically keep pace with changes to the gamma distribution functions made by the core team of R maintainers. Third, invgamma follows the same naming conventions and parameterizations as base R and has the same arguments as the ∗gamma() functions, so that using the ∗invgamma() functions are natural to R users. We note that actuar’s implementation satisfies (3) as well. Additionally, invgamma includes ∗invchisq() and ∗invexp() functions for the closely related inverse ${\chi }^2$ (chi-squared) and exponential distributions.

Accuracy considerations and limitations

One additional point that warrants a note is numerics concerning the transformations. invgamma was intended to be a lightweight, low-maintenance package implementing the inverse gamma, inverse chi-square, and inverse exponential distributions. It uses the transformation theorem in all cases. One of the challenges to using this implementation strategy is that the numerics are not optimized for the particular situation at hand. Arithmetic on a computer is not the same as arithmetic in theory, and as a consequence the best computer implementations of mathematical facts/algorithms are tailored to the specific cases at hand. Since these are not optimized in this way, questions may arise to the extent of the validity of the functions. The biggest questions surround the random number generators.

One way to test the values produced by the RNGs is via goodness-of-fit (GoF) testing of their draws. From an inferential perspective this is analogous to the graphics presented in Fig. 2, but operationally it is more mechanical as it does not require judgment calls, at least not directly. The basic idea is to take a large number of draws $n$ from the presumed distribution and assess those with a GoF test where $H_0:F_X=F_0$, where $F_X$ represents the distribution of the draws and $F_0$ represents the claimed distribution. If we conduct such a procedure at points ranging across the $(\alpha ,\lambda )$ parameter space, we can get a good feel for which ranges of the parameter space the RNGs break down.

To code such a simulation, we begin with a function that performs such a GoF test, the Kolmogorov-Smirnov (KS) test, for a specific shape/rate combination. The following code accepts an $(\alpha ,\lambda )$ shape/rate parameterization of the inverse gamma and uses the built-in function ks.test() to generate a $p$-value corresponding to the KS test for the inverse gamma distribution:

We then create a grid across a wide range of the parameter space, using a logarithmic spacing for better representation. We conduct the simulation at each of the points on the grid. Under the null hypothesis that the sampler is working correctly, at the 5% level roughly one out of every 20 tests will reject by design, so if we see regions that seem to have a higher rejection rate, those regions should be considered suspect. Figure 3 contains the results based off of the following code.

If the RNG were correct, roughly 1 out of every 20 points should be purple. As this is violated for small shapes, the RNG is not to be trusted for small $\alpha$. Figure 3 clearly suggests that this is the case when $\alpha$ is less than about 0.01, regardless of what the value of $\lambda$ is. As a consequence, rinvgamma() issues a warning in these situations:

The reason for this failure is quite evident: when the shape is low, the gamma distribution piles its mass near 0. When the draws are inverted, they are enormous, frequently overflowing the floating point number system to return infinity. This problem is inherited from R’s RNG for the gamma distribution, for which the corresponding values are identically zero:

Similar R problems exist with duplications, even when the sampler is accurate:

Both of these illustrate that the problems that rinvgamma() has are limitations inherited from the R’s rgamma() itself and are quite mild. Because the functions of invgamma are tied to these implementations, they will be corrected automatically when rgamma() is.

Similar tests can be conducted for the inverse chi-squared and inverse exponential distributions. The inverse chi-squared distribution is similar to the inverse gamma. Figure 4 illustrates that the inverse chi-squared RNG fails when $\nu \le 0.01$ and $ncp\le 10$ or so, and thus the RNG issues a warning in this setting.

The inverse exponential is the easiest as it only admits a single parameter, rate. Unlike the others, the simulations illustrated in Fig. 5 suggest that the inverse exponential RNG is robust across a wide range of parameter values.

Conclusions

In this article we have surveyed the invgamma R package, along with documenting the mathematics motivating the functions provided and tests of their accuracy. The package functions ∗invgamma(), ∗invchisq(), and ∗invexp() are largely self-maintained through their connections to the ∗gamma(), ∗chisq(), and ∗exp() functions provided in the base distribution of R and efficient for the same reason. The package is light weight and ideal for including in other R packages.

Supplemental Information