README.md
In ginormal: Generalized Inverse Normal Distribution Density and Generation

ginormal: Generalized Inverse Normal distribution

The ginormal package provides the density function and random variable generation from the generalized inverse normal (GIN) distribution introduced by Robert (1991). The GIN distribution is a way to generalize the distribution of the reciprocal of a normal random variable. That is, the distribution generalizes the distribution of the random variable $Z = 1/X$ where $X \sim \text{Normal}(\mu, \sigma^2)$. This distribution is different from the generalized inverse Gaussian (GIG) distribution (Jørgensen, 2012) despite the similarities in naming (see below).

The GIN distribution is supported on the entire real line $z \in (-\infty, \infty)$ and takes three parameters:

$\alpha > 1$, a degrees-of-freedom parameter,
$\mu \in (-\infty, \infty)$, similar to a location parameter, it shifts the density of the distribution left and right,
$\tau > 0$, similar to a scale parameter, it spreads the density of the distribution.

This package is the first to provide an efficient sampling algorithm for drawing from the GIN distribution. We provide similar routines for the GIN distribution truncated to the positive or negative reals. Further details of the distribution, theoretical guarantees and pseudo-code for the sampling algorithms, as well as an application to Bayesian estimation of network formation models can be found in the working paper Ding, Estrada and Montoya-Blandón (2023).

To install, type from within R:

install.packages("ginormal")

You could also install directly from the GitHub repository. To do so from within R, first install the devtools package and then type:

# install.packages("devtools") # Uncomment if devtools not installed
library(devtools)
install_github(repo = "smonto2/ginormal", subdir = "R_package", ref = "main")

Examples of how to use the ginormal package routines are available in the GitHub repository.

Provided with the package are four main routines:

dgin(z, alpha, mu, tau, log = TRUE, quasi = FALSE)
dtgin(z, alpha, mu, tau, sign, log = TRUE, quasi = FALSE)
rgin(size, alpha, mu, tau, algo)
rtgin(size, alpha, mu, tau, sign, algo)

The first two compute the densities and the last two are used for random number generation. Density routines take in the quantile z; parameters alpha, mu and tau; and two optional logical arguments:

log, should the logarithm of the density be returned? Defaults to TRUE.
quasi, should the value of the kernel (or quasi-density) be returned? Defaults to FALSE.

Generation routines take the same parameters but require a size argument determining the amount of random variates to generate. These routines only admit a parameter alpha larger than 2. They take an additional argument algo, which can be either "hormann" or "leydold", and defaults to "hormann" as our prefered method. See below for details on both points.

Those routines including "t" in their name work for the truncated variants. They take an additional logical argument sign, where sign = TRUE implies truncation to positive numbers $(z > 0)$ and sign = FALSE to negative numbers $(z < 0)$.

Let $Z \sim \text{GIN}(\alpha, \mu, \tau)$. The GIN density function is given by $$f_Z(z) = \frac{1}{C(\alpha, \mu, \tau)} |z|^{-\alpha}\exp\left[-\frac{1}{2\tau^2} \left( \frac{1}{z} - \mu \right)^2 \right] \equiv \frac{g(z; \alpha, \mu, \tau)}{C(\alpha, \mu, \tau)}$$ where $g(z; \alpha, \mu, \tau)$ is the kernel or quasi-density and the proportionality constant can be written in closed form as $$ C(\alpha, \mu, \tau) = (\sqrt{2} \tau)^{\alpha-1} \exp\left(- \frac{\mu^2}{2\tau^2} \right) \Gamma\left(\frac{\alpha-1}{2}\right) {}_1F_1 \left(\frac{\alpha-1}{2}; \frac{1}{2}; \frac{\mu^2}{\tau^2}\right) $$

where $\Gamma(x)$ is the Gamma function and ${}1F_1(a, b; x)$ is the confluent hypergeometric function. In addition to the density and generation routines for the GIN distribution, we provide similar routines for the GIN distribution truncated to positive or negative numbers. These are denoted by $\text{GIN}^{+}$ when truncated to $(0, \infty)$ and by $\text{GIN}^{-}$ when truncated to $(-\infty, 0)$. Let $Z^{+} \sim \text{GIN}^{+}(\alpha, \mu, \tau)$ and $Z^{-} \sim \text{GIN}^{-}(\alpha, \mu, \tau)$. Their densities are given by $$f{Z^{+}}(z) = \frac{g(z; \alpha, \mu, \tau)}{C^{+}(\alpha, \mu, \tau)} \mathbb{I}(z > 0)$$ $$f_{Z^{-}}(z) = \frac{g(z; \alpha, \mu, \tau)}{C^{-}(\alpha, \mu, \tau)} \mathbb{I}(z < 0)$$ with proportionality constants $$C^{+}(\alpha, \mu) = e^{-\frac{\mu^2}{4}} \Gamma(\alpha - 1) D_{-(\alpha-1)}(-\mu)$$ $$C^{-}(\alpha, \mu) = e^{-\frac{\mu^2}{4}} \Gamma(\alpha - 1) D_{-(\alpha-1)}(\mu)$$ where $\mathbb{I}(\cdot)$ is the indicator function that is 1 when its argument is true and 0 otherwise, and $D_\nu(x)$ is the parabolic cylinder function. [^1]

Ding, Estrada and Montoya-Blandón (2023) provide an efficient sampling algorithm for the GIN distribution and its truncated variants for the case of $\alpha > 2$. This restriction is not of concern if the goal is the perform Bayesian estimation using this distribution (see below for more details and Remark 2 in the paper). Generation is done using the ratio-of-uniforms method with mode shift (Kinderman and Monahan, 1977), which requires the computation of the minimal bounding rectangle. We implement two alternatives found in the literature:

Leydold (2001) that requires information on the proportionality constants.
Hörmann and Leydold (2014) that requires solving a cubic equation. This is our prefered method and the default in the package.

While the kernels — and therefore the sampling techniques — for the GIN and GIG distribution are similar, these two distribution share some important differences. The main is their conceptualization, as they both attempt to generalize the idea of an inverse normal distribution in different ways. The GIG distribution does so by choosing cumulants that are inverses to those of the normal distribution. The GIN distribution does so by directly using the density of the reciprocal after a change of variables. Another important difference comes from their use as conjugate priors in Bayesian analysis:

$\theta \sim \text{GIN}(\alpha, \mu, \tau)$ is the conjugate prior if observations are random samples from $Y \sim \text{Normal}(\theta, \theta^2)$
$\theta \sim \text{GIG}(\alpha, \mu, \tau)$ is the conjugate prior if observations are random samples from $Y \sim \text{Normal}(\theta, \theta)$

These are both mixture models with a similar structure but carry different interpretations and thus require different posterior sampling algorithms. This interpretation also shows why the restriction of $\alpha \geq 2$ is not binding if the goal is to perform Bayesian analysis. A prior $\theta \sim \text{GIN}(\alpha_0, \mu_0, \tau_0)$ with $\alpha_0 = 1 + \varepsilon$ is non-informative when $\varepsilon > 0$ is arbitrarily small. However, the posterior distribution will have degrees-of-freedom parameter $\alpha_N = N + 1 + \varepsilon$ where $N$ is the sample size. As $N \geq 1$ implies $\alpha_N > 2$, for a conjugate Bayesian analysis we are always drawing from the GIN distribution with $\alpha > 2$.

[^1]: In R, package BAS contains the confluent hypergeometric function. For the parabolic cylinder function, we use a Fortran subroutine provided in the SPECFUN library (Zhang and Jin, 1996) and our own R translation of this function.