statmod: Probability Calculations for the Inverse Gaussian Distribution

• The IGD is used in reliability and survival analysis and for modeling non-negative, positively skewed data.Source: §1 · p. 1
• The density is easy to compute and random variates can be generated from chi-square and binomial variables, but cdf log-scale accuracy needs care and quantiles require iteration.Source: §1 · p. 2
• Earlier quantile solvers either converge slowly with limited precision or can diverge unpredictably, even from a close initial approximation.Source: §1 · p. 2
• Starting Newton’s method at the mode guarantees monotonic convergence for IGD quantiles, eliminates backtracking, and converges rapidly even in extreme tails.Convergence eventually becomes quadratic.Source: §1 · p. 2
• The statmod functions dinvgauss, pinvgauss, qinvgauss, and rinvgauss work for all parameter values at close to machine accuracy and test faster, more accurate, and more reliable than existing IGD functions.Source: §1 · pp. 2–3

• For x > 0, µ > 0, and φ > 0, IG(µ,φ) has mean µ and variance φµ^3; λ = 1/φ is the alternative shape parameter.Source: §2 · p. 3
• Scaling X by µ gives an inverse Gaussian with mean 1 and dispersion φµ, which also determines skewness and equals the squared coefficient of variation.Source: §2 · p. 3
• The density is unimodal with mode strictly between 0 and µ; when φµ is large, the mode becomes dependent on φ alone.Table 1 gives pdf values for special parameter cases, including theoretical limits for infinite parameters.Source: §2 · p. 4, Table 1 · p. 4
• When φ is very large, the density develops a spike at 0 regardless of µ.Source: §2 · p. 4
• statmod handles examples with missing or infinite parameters and out-of-support x values, whereas the existing functions error on those calls.Source: §2 · p. 5

3 Cumulative distribution function

• Underflow, overflow, and subtractive cancellation threaten the cdf’s unlogged sum, so statmod combines terms on the log scale to compute far-tail values reliably.The calculation removes a common factor before combining terms; log1p() avoids cancellation for small arguments.Source: §3 · pp. 5–6
• The right-tail probability has an asymptotic approximation that is very accurate when φ−1/2 m(qm−1) > 105, but for more modest values it gives only 2–3 significant figures correctly.Source: §3 · pp. 5–6
• -7146.914: statmod returns this log probability for q = 0.0001, whereas other packages cannot compute log probabilities below about -700.Existing functions also return zero for some very small left- and right-tail probabilities.Source: §3 · pp. 6–7
• The cdf handles limiting cases: infinite mean corresponds to an inverse-chi-square case, infinite dispersion to a spike at zero, and extreme x values yield 0 or 1.The example with missing mean and dispersion returns 0 for x = -1 and 1 for x = Inf.Source: §3 · p. 7
• 15 significant figures: statmod agrees with the chi-square reference at q1 = 0.01, while STAR reaches 3 and rmutil and SuppDists have relative error near 100%.For q1 = 0.1, statmod also agrees to 15 significant figures, versus 12 for the other packages.Source: §3 · p. 8

4 Inverting the cdf

• Starting Newton’s method at the mode guarantees monotonic convergence because the cdf is convex to its left and concave to its right.The iteration therefore needs no backtracking and is eventually quadratically convergent.Source: §4 · p. 9, Figure 2 · p. 10
• For tail probabilities below 10^-5, the method uses normal- or gamma-based starting values that remain within the domain of monotonic convergence.The mode remains the starting value otherwise; tail-specific starting values are used only when their iteration savings outweigh their computation cost.Source: §4 · pp. 9, 11
• When the target and current cdf probabilities are close, a Taylor expansion on the log scale avoids cancellation in the Newton-step probability difference.The expansion is used when |δ| < 10^-5, with δ defined as the difference between their log probabilities.Source: §4 · pp. 9, 11
• 16 significant figures: statmod::qinvgauss outperforms the other packages’ 6–8 figures, with probability error never above 2e-16 and relative quantile error never above 5e-16.The reported probability and relative quantile errors are separate measures.Source: §4 · pp. 11–12
• 40 times faster than rmutil or STAR and about 3 times faster than SuppDists, statmod::qinvgauss processed one million probabilities in the reported benchmark.The timings were measured on a Windows laptop running 64-bit R-devel.Source: §4 · p. 12
• statmod::qinvgauss runs successfully for all parameter values, while SuppDists can fail to converge and rmutil and STAR do not support right-tail probabilities.It also computes a right-tail quantile when the corresponding left-tail probability is indistinguishable from 1 on the unlogged scale.Source: §4 · pp. 12–13
• The statmod quantile function accepts either dispersion or shape and follows R probability-function conventions for tails and log probabilities.The functions also detect boundary or invalid probabilities and preserve input attributes under the stated length condition.Source: §4 · pp. 13–14

5 Random deviates

• statmod, SuppDists, and STAR generate chi-square deviates and select between the two corresponding IGD values using probabilities worked out by.The shared method generates chi-square deviates corresponding to (X − µ)^2/(φXµ^2).Source: §5 · p. 14
• rmutil::rinvgauss generates deviates by applying qinvgauss to uniform draws, making it slower and less accurate than the other functions.statmod generates one million deviates in about a quarter-second; STAR takes about half a second, while SuppDists is faster because it is implemented in C.Source: §5 · p. 14

6 Discussion

• The statmod functions are the first IGD probability functions reported to work for all parameter values and return close to full machine accuracy.They compute probabilities on the log-scale where possible and handle special limiting cases, including some not previously described.Source: §6 · p. 14
• Starting Newton iteration from the IGD mode gives monotonic convergence, avoiding backtracking and allowing a change in step sign to signal floating-point limits.The iteration stops on a sign change if its convergence criterion has not yet been met.Source: §6 · pp. 14–15
• qinvgauss is 40 times faster than rmutil and STAR quantile functions while returning 16 rather than 6 figures of accuracy.It is also 3 times faster than SuppDists::qinvGauss, despite that function being written in C.Source: §6 · p. 15
• The statmod functions remain in pure R; reimplementation in C could make them faster, while the current versions are easier to understand and maintain.They are only slightly slower than comparable functions such as qchisq and qt.Source: §6 · p. 15
• The quantile strategy may apply to any continuous unimodal distribution or one that can be transformed to be unimodal.The paper also notes that the current implementation prioritizes the benefits of understandability and easy maintenance over possible speed gains from C.Source: §6 · p. 15

Appendix: asymptotic right tail probabilities

• For large q, the appendix derives an asymptotic expression for the right-tail probability, assuming µ = 1 without loss of generality.The derivation begins by dropping the 1/x term in the density exponent for large x.Source: Appendix: asymptotic right tail… · p. 18
• The asymptotic derivation uses the large-x approximation to the density before transforming and approximating the integral for large q.The supplied passage describes these steps but does not state the resulting expression.Source: Appendix: asymptotic right tail… · p. 18
• After simplifying the density exponent, the derivation integrates the density to express the right-tail probability for large q.The appendix then transforms the variable of integration and approximates the resulting integral.Source: Appendix: asymptotic right tail… · p. 18