statmod: Probability Calculations for the Inverse Gaussian Distribution

TL;DR

This paper introduces the statmod package for R, providing fast, accurate probability functions for the inverse Gaussian distribution, including a novel method for computing quantiles with guaranteed convergence.

Contribution

It develops reliable, machine-precision probability functions for the inverse Gaussian distribution in R, including a new quantile computation method with proven convergence.

Findings

Probability functions achieve near full machine accuracy.

Newton's method for quantiles converges monotonically from the mode.

The package is available on CRAN for practical use.

Abstract

The inverse Gaussian distribution (IGD) is a well known and often used probability distribution for which fully reliable numerical algorithms have not been available. Our aim in this article is to develop software for this distribution for the R programming environment. We develop fast, reliable basic probability functions (dinvgauss, pinvgauss, qinvgauss and rinvgauss) that work for all possible parameter values and which achieve close to full machine accuracy. The most challenging task is to compute quantiles for given cumulative probabilities and we develop a simple but elegant mathematical solution to this problem. We show that Newton's method for finding the quantiles of a IGD always converges monotonically when started from the mode of the distribution. Simple Taylor series expansions are used to improve accuracy on the log-scale. The IGD probability functions provide the same options and obey the same conventions as do probability functions provided in the standard R stats package. The IGD functions are part of the statmod package available from the CRAN repository.