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.
Highlights
The inverse Gaussian distribution (IGD) (Tweedie, 1957; Johnson and Kotz, 1970) is widely used in a variety of application areas including reliability and survival analysis (Whitmore, 1975; Chhikara and Folks, 1977; Bardsley, 1980; Chhikara, 1989; Wang and Xu, 2010; Balakrishna and Rahul, 2014)
We have found that none of these IGD functions work for all parameter values or return results to full machine accuracy
> pinvgauss(q[1], mean = mu, disp = phi) + + pinvgauss(q[2], mean = mu, disp = phi, lower.tail = FALSE) [1] 1.6427313604456183e-32 > rmutil::pinvgauss(q[1], m = mu, s = phi) + + 1 - rmutil::pinvgauss(q[2], m = mu, s = phi) [1] 0 > SuppDists::pinvGauss(q[1], nu = mu, lambda = 1/phi) + + SuppDists::pinvGauss(q[2], nu = mu, lambda = 1/phi, lower.tail = FALSE) [1] 8.2136568022278466e-33 > STAR::pinvgauss(q[1], mu = mu, sigma2 = phi) + + STAR::pinvgauss(q[2], mu = mu, sigma2 = phi, lower.tail = FALSE) [1] 1.6319986233795599e-32. It can be seen from the above that rmutil and SuppDists do not agree with pchisq to any significant figures, meaning that the relative error is close to 100%, while STAR manages 3 significant figures. statmod on the other hand continues to agree with pchisq to 15 significant figures
Summary
The inverse Gaussian distribution (IGD) (Tweedie, 1957; Johnson and Kotz, 1970) is widely used in a variety of application areas including reliability and survival analysis (Whitmore, 1975; Chhikara and Folks, 1977; Bardsley, 1980; Chhikara, 1989; Wang and Xu, 2010; Balakrishna and Rahul, 2014). One is to solve for the quantile using a general-purpose equation solver such as the uniroot function in R This is the approach taken by the qinvgauss functions in the rmutil and STAR packages. The other approach is to use Newton’s method to solve the equation after applying an initial approximation (Kallioras and Koutrouvelis, 2014) This approach was taken by one of the current authors when developing inverse Gaussian code for S-PLUS (Smyth, 1998). It is the approach taken by the qinvGauss function in the SuppDists package. Newton’s method cannot in general be guaranteed to converge, even when the initial approximation is close to the required value, and the parameter values for which divergence occurs are hard to predict. The starting value may be far from the required solution, the rapid convergence
Sign-in/Register to view highlights & summary Sign-in/Register to access full text options 
Free) 
Talk to us
Join us for a 30 min session where you can share your feedback and ask us any queries you have
Schedule a call
