feat: add the inverse Gaussian (Wald) distribution

build.rs

@@ -0,0 +1,37 @@
+//! Detects whether the target has a hardware fused multiply-add.
+//!
+//! `f64::mul_add` is always correct, but on a target without the instruction it
+//! calls a software routine that is far slower than doing the same job with
+//! Veltkamp's split. `src/prec.rs` picks between the two on `statrs_hardware_fma`.
+
+fn main() {
+    println!("cargo::rustc-check-cfg=cfg(statrs_hardware_fma)");
+    println!("cargo::rerun-if-changed=build.rs");
+
+    let arch = std::env::var("CARGO_CFG_TARGET_ARCH").unwrap_or_default();
+    let features: Vec<String> = std::env::var("CARGO_CFG_TARGET_FEATURE")
+        .unwrap_or_default()
+        .split(',')
+        .map(str::to_owned)
+        .collect();
+    let has = |f: &str| features.iter().any(|x| x == f);
+
+    // Each arm is "the FMA instruction is present", not "the architecture
+    // usually has it": AArch64 without NEON is softfloat, and baseline x86-64
+    // needs -C target-feature=+fma or a -C target-cpu implying it.
+    let hardware_fma = match arch.as_str() {
+        "aarch64" | "arm64ec" => has("neon"),
+        "x86" | "x86_64" => has("fma"),
+        "arm" => has("vfp4"),
+        "riscv32" | "riscv64" => has("d") || has("f"),
+        "powerpc" | "powerpc64" => has("vsx") || has("altivec"),
+        "loongarch32" | "loongarch64" => has("d") || has("f"),
+        "mips" | "mips64" | "mips32r6" | "mips64r6" => has("fp64"),
+        "wasm32" | "wasm64" => has("relaxed-simd"),
+        _ => false,
+    };
+
+    if hardware_fma {
+        println!("cargo::rustc-cfg=statrs_hardware_fma");
+    }
+}

examples/lnerfc_sweep.rs

@@ -0,0 +1,10 @@
+use statrs::function::erf::ln_erfc;
+fn main() {
+    for i in 0..=800 {
+        let x = -6.0 + (206.0) * (i as f64) / 800.0; // covers <0.5, rational range, boundary at 110, asymptotic
+        println!("{:016x}\t{:016x}", x.to_bits(), ln_erfc(x).to_bits());
+    }
+    for x in [1e3f64, 1e5, 1e10, 1e77, 1.3e154, 1.4e154, 1e200] {
+        println!("{:016x}\t{:016x}", x.to_bits(), ln_erfc(x).to_bits());
+    }
+}

src/distribution/beta.rs

@@ -172,6 +172,32 @@ impl ContinuousCDF<f64, f64> for Beta {
         }
     }
 
+    /// Tail-accurate log of the cdf via [`beta::ln_beta_reg`].
+    fn ln_cdf(&self, x: f64) -> f64 {
+        if x < 0.0 {
+            f64::NEG_INFINITY
+        } else if x >= 1.0 {
+            0.0
+        } else if prec::ulps_eq!(self.shape_a, 1.0) && prec::ulps_eq!(self.shape_b, 1.0) {
+            x.ln()
+        } else {
+            beta::ln_beta_reg(self.shape_a, self.shape_b, x)
+        }
+    }
+
+    /// Tail-accurate log of the survival function via [`beta::ln_beta_reg`].
+    fn ln_sf(&self, x: f64) -> f64 {
+        if x < 0.0 {
+            0.0
+        } else if x >= 1.0 {
+            f64::NEG_INFINITY
+        } else if prec::ulps_eq!(self.shape_a, 1.0) && prec::ulps_eq!(self.shape_b, 1.0) {
+            (-x).ln_1p()
+        } else {
+            beta::ln_beta_reg(self.shape_b, self.shape_a, 1.0 - x)
+        }
+    }
+
     /// Calculates the inverse cumulative distribution function for the beta
     /// distribution at `x`.
     ///

src/distribution/binomial.rs

@@ -1,5 +1,5 @@
 use crate::distribution::{Discrete, DiscreteCDF};
-use crate::function::{beta, factorial};
+use crate::function::{beta, factorial, gamma};
 use crate::prec;
 use crate::statistics::*;
 use core::f64;
@@ -169,6 +169,27 @@ impl DiscreteCDF<u64, f64> for Binomial {
             beta::beta_reg(k as f64 + 1.0, (self.n - k) as f64, self.p)
         }
     }
+
+    /// Tail-accurate log of the cdf via [`beta::ln_beta_reg`]; finite far past
+    /// where `cdf` underflows (Bonferroni-scale p-values stay representable).
+    fn ln_cdf(&self, x: u64) -> f64 {
+        if x >= self.n {
+            0.0
+        } else {
+            let k = x;
+            beta::ln_beta_reg((self.n - k) as f64, k as f64 + 1.0, 1.0 - self.p)
+        }
+    }
+
+    /// Tail-accurate log of the survival function via [`beta::ln_beta_reg`].
+    fn ln_sf(&self, x: u64) -> f64 {
+        if x >= self.n {
+            f64::NEG_INFINITY
+        } else {
+            let k = x;
+            beta::ln_beta_reg(k as f64 + 1.0, (self.n - k) as f64, self.p)
+        }
+    }
 }
 
 impl Min<u64> for Binomial {
@@ -237,7 +258,11 @@ impl Distribution<f64> for Binomial {
         } else {
             (0..self.n + 1).fold(0.0, |acc, x| {
                 let p = self.pmf(x);
-                acc - p * p.ln()
+                // A mass that underflows to zero contributes nothing, taking
+                // `0 * ln 0 == 0` as usual for entropy. Evaluating it would give
+                // `0 * -inf == NaN` and poison the whole sum, which made
+                // `entropy()` return `NaN` for any `n` past roughly 1100.
+                if p > 0.0 { acc - p * p.ln() } else { acc }
             })
         };
         Some(entr)
@@ -305,10 +330,7 @@ impl Discrete<u64, f64> for Binomial {
         } else if prec::ulps_eq!(self.p, 1.0) {
             if x == self.n { 1.0 } else { 0.0 }
         } else {
-            (factorial::ln_binomial(self.n, x)
-                + x as f64 * self.p.ln()
-                + (self.n - x) as f64 * (1.0 - self.p).ln())
-            .exp()
+            self.ln_pmf(x).exp()
         }
     }
 
@@ -328,9 +350,38 @@ impl Discrete<u64, f64> for Binomial {
         } else if prec::ulps_eq!(self.p, 1.0) {
             if x == self.n { 0.0 } else { f64::NEG_INFINITY }
         } else {
-            factorial::ln_binomial(self.n, x)
-                + x as f64 * self.p.ln()
-                + (self.n - x) as f64 * (1.0 - self.p).ln()
+            let n = self.n as f64;
+            let k = x as f64;
+            let q = 1.0 - self.p;
+            if x == 0 {
+                return n * q.ln();
+            }
+            if x == self.n {
+                return n * self.p.ln();
+            }
+            if n < gamma::STIRLING_SERIES_MIN {
+                // `ln_binomial` is exact from the factorial table here and all
+                // terms are small, so the direct form loses nothing
+                return factorial::ln_binomial(self.n, x) + k * self.p.ln() + (n - k) * q.ln();
+            }
+            // Saddle-point form (Loader 2000): the direct expression differences
+            // terms growing like `n ln n` to produce an `O(1)` result, which
+            // cost ~2e-10 relative at n = 1e5.
+            let nk = n - k;
+            // `bd0` is sensitive to the last half-ulp of its mean argument, and
+            // both `1 - p` and the products round, so they are carried as
+            // double-doubles; see `gamma::bd0_dd`.
+            let (q_hi, q_lo) = prec::two_diff(1.0, self.p);
+            let np = n * self.p;
+            let np_lo = prec::dekker_product_err(n, self.p, np);
+            let nq = n * q_hi;
+            let nq_lo = prec::dekker_product_err(n, q_hi, nq) + n * q_lo;
+            gamma::stirling_delta(n)
+                - gamma::stirling_delta(k)
+                - gamma::stirling_delta(nk)
+                - gamma::bd0_dd(k, np, np_lo)
+                - gamma::bd0_dd(nk, nq, nq_lo)
+                + 0.5 * (n / (f64::consts::TAU * k * nk)).ln()
         }
     }
 }
@@ -412,6 +463,61 @@ mod tests {
         test_exact(0.3, 10, 10, max);
     }
 
+    /// `p` very close to but not equal to 1 must take the general path, not the
+    /// degenerate `p == 1` branch. The branch is guarded by `prec::ulps_eq!`,
+    /// whose default epsilon was `1e-9` *absolute*, so any `p` within `1e-9` of
+    /// 1 collapsed to a point mass at `n`.
+    ///
+    /// `1 - 2^-33` is used rather than `1 - 1e-10` so that `1 - p` is exact and
+    /// the reference values are not limited by the representation of `p`.
+    #[test]
+    fn test_pmf_p_near_one_is_not_degenerate() {
+        let n = Binomial::new(1.0 - f64::powi(2.0, -33), 100).unwrap();
+        prec::assert_relative_eq!(n.pmf(99), 1.1641532048523463366e-8, epsilon = 0.0, max_relative = 1e-13);
+        prec::assert_relative_eq!(n.pmf(100), 0.99999998835846788439, epsilon = 0.0, max_relative = 1e-14);
+        prec::assert_relative_eq!(n.ln_pmf(99), -18.268686784015220704, epsilon = 0.0, max_relative = 5e-15);
+        // entropy of a non-degenerate distribution is strictly positive
+        assert!(n.entropy().unwrap() > 0.0);
+    }
+
+    /// `entropy()` sums `-p * ln(p)` over the whole support, so it used to
+    /// return `NaN` as soon as one mass underflowed to zero (`0 * -inf`). For
+    /// `p = 0.5` that starts at around `n = 1100`.
+    #[test]
+    fn test_entropy_with_underflowing_masses() {
+        for n in [1_100, 2_000, 5_000, 20_000] {
+            let d = Binomial::new(0.5, n).unwrap();
+            let h = d.entropy().unwrap();
+            assert!(h.is_finite(), "entropy for n={n} was {h}");
+            // 0.5 * ln(2 * pi * e * n * p * q) is the asymptotic form
+            let approx = 0.5 * (2.0 * f64::consts::PI * f64::consts::E * n as f64 * 0.25).ln();
+            prec::assert_relative_eq!(h, approx, epsilon = 0.0, max_relative = 1e-3);
+        }
+        assert!(Binomial::new(0.1, 5_000).unwrap().entropy().unwrap().is_finite());
+    }
+
+    /// Large-`n` pmf: the saddle-point form (`gamma::bd0` /
+    /// `gamma::stirling_delta`, with double-double means) replaced
+    /// `exp(ln_binomial + x ln p + (n-x) ln q)`, whose terms each grow like
+    /// `n ln n` while the result stays `O(1)`. References are mpmath at 40
+    /// significant digits; the old form was ~7e-11 relative here.
+    #[test]
+    fn test_pmf_large_n_saddle_point() {
+        let pmf = |arg: u64| move |x: Binomial| x.pmf(arg);
+        test_relative(0.5, 100000, 0.0025231262141967398855, pmf(50000));
+        test_relative(0.3, 2000000, 0.00061558120658465376062, pmf(600000));
+        // ln_pmf agrees with ln(pmf) where the pmf is comfortably representable
+        for (p, n, k) in [(0.5f64, 100000u64, 50000u64), (0.3, 2000000, 600000)] {
+            let d = create_ok(p, n);
+            prec::assert_relative_eq!(d.ln_pmf(k), d.pmf(k).ln(), epsilon = 0.0, max_relative = 1e-14);
+        }
+        // the endpoints reduce to q^n and p^n; `exp(n ln q)` carries ~1e-15
+        // relative (|n ln q| ~ 11.5 here), which is the honest bound
+        let d = create_ok(0.25, 40);
+        prec::assert_relative_eq!(d.pmf(0), 1.0056585161637497e-5, epsilon = 0.0, max_relative = 1e-14);
+        prec::assert_relative_eq!(d.pmf(40), 8.271806125530277e-25, epsilon = 0.0, max_relative = 1e-14);
+    }
+
     #[test]
     fn test_pmf() {
         let pmf = |arg: u64| move |x: Binomial| x.pmf(arg);
@@ -515,7 +621,7 @@ mod tests {
         test_absolute(0.3, 3, 0.657, 1e-14, sf(0));
         test_absolute(0.3, 3, 0.216, 1e-15, sf(1));
         test_exact(0.3, 3, 0.0, sf(3));
-        test_absolute(0.3, 10, 0.9717524751000001, 1e-16, sf(0));
+        test_absolute(0.3, 10, 0.9717524750999999955198, 1e-16, sf(0));
         test_absolute(0.3, 10, 0.850691654100002, 1e-14, sf(1));
         test_exact(0.3, 10, 0.0, sf(10));
         test_exact(1.0, 1, 1.0, sf(0));

src/distribution/chi_squared.rs

@@ -139,6 +139,16 @@ impl ContinuousCDF<f64, f64> for ChiSquared {
         self.g.sf(x)
     }
 
+    /// Delegates to the wrapped Gamma's tail-accurate implementation.
+    fn ln_cdf(&self, x: f64) -> f64 {
+        self.g.ln_cdf(x)
+    }
+
+    /// Delegates to the wrapped Gamma's tail-accurate implementation.
+    fn ln_sf(&self, x: f64) -> f64 {
+        self.g.ln_sf(x)
+    }
+
     /// Calculates the inverse cumulative distribution function for the
     /// chi-squared distribution at `x`
     ///

src/distribution/dirichlet.rs

@@ -19,7 +19,7 @@ use nalgebra::{Dim, Dyn, OMatrix, OVector};
 ///
 /// let n = Dirichlet::new(vec![1.0, 2.0, 3.0]).unwrap();
 /// assert_eq!(n.mean().unwrap(), DVector::from_vec(vec![1.0 / 6.0, 1.0 / 3.0, 0.5]));
-/// assert_eq!(n.pdf(&DVector::from_vec(vec![0.33333, 0.33333, 0.33333])), 2.222155556222205);
+/// assert_eq!(n.pdf(&DVector::from_vec(vec![0.33333, 0.33333, 0.33333])), 2.2221555562222193);
 /// ```
 #[derive(Clone, PartialEq, Debug)]
 pub struct Dirichlet<D>

src/distribution/erlang.rs

@@ -123,6 +123,16 @@ impl ContinuousCDF<f64, f64> for Erlang {
         self.g.sf(x)
     }
 
+    /// Delegates to the wrapped Gamma's tail-accurate implementation.
+    fn ln_cdf(&self, x: f64) -> f64 {
+        self.g.ln_cdf(x)
+    }
+
+    /// Delegates to the wrapped Gamma's tail-accurate implementation.
+    fn ln_sf(&self, x: f64) -> f64 {
+        self.g.ln_sf(x)
+    }
+
     /// Calculates the inverse cumulative distribution function for the erlang
     /// distribution at `x`
     ///

src/distribution/exponential.rs

@@ -117,7 +117,9 @@ impl ContinuousCDF<f64, f64> for Exp {
         if x < 0.0 {
             0.0
         } else {
-            1.0 - (-self.rate * x).exp()
+            // `-expm1` rather than `1 - exp`: the latter cancels for small
+            // `rate * x`, losing ~7 digits by `rate * x = 1e-10`
+            -(-self.rate * x).exp_m1()
         }
     }
 
@@ -135,6 +137,22 @@ impl ContinuousCDF<f64, f64> for Exp {
         if x < 0.0 { 1.0 } else { (-self.rate * x).exp() }
     }
 
+    /// Tail-accurate log of the cdf: `ln(-expm1(-rate x))`.
+    fn ln_cdf(&self, x: f64) -> f64 {
+        if x < 0.0 {
+            f64::NEG_INFINITY
+        } else {
+            (-(-self.rate * x).exp_m1()).ln()
+        }
+    }
+
+    /// Exact log of the survival function: `ln sf = -rate * x`, finite for
+    /// every representable `x` even though `sf` underflows past
+    /// `x ~ 709 / rate`.
+    fn ln_sf(&self, x: f64) -> f64 {
+        if x < 0.0 { 0.0 } else { -self.rate * x }
+    }
+
     /// Calculates the inverse cumulative distribution function.
     ///
     /// # Formula

src/distribution/fisher_snedecor.rs

@@ -525,8 +525,8 @@ mod tests {
         test_absolute(10.0, 0.1, 0.0418440630400545297349, 1e-14, pdf(1.0));
         test_absolute(0.1, 1.0, 0.0396064560910663979961, 1e-16, pdf(1.0));
         test_absolute(1.0, 1.0, 0.1591549430918953357689, 1e-16, pdf(1.0));
-        test_absolute(10.0, 1.0, 0.230361989229138647108, 1e-16, pdf(1.0));
-        test_absolute(0.1, 0.1, 0.00221546909694001013517, 1e-18, pdf(10.0));
+        test_absolute(10.0, 1.0, 0.230361989229138647108, 3e-16, pdf(1.0));
+        test_absolute(0.1, 0.1, 0.00221546909694001013517, 5e-18, pdf(10.0));
         test_absolute(1.0, 0.1, 0.00369960370387922619592, 1e-17, pdf(10.0));
         test_absolute(10.0, 0.1, 0.00390179721174142927402, 1e-15, pdf(10.0));
         test_absolute(0.1, 1.0, 0.00319864073359931548273, 1e-17, pdf(10.0));
@@ -540,13 +540,13 @@ mod tests {
     #[test]
     fn test_ln_pdf() {
         let ln_pdf = |arg: f64| move |x: FisherSnedecor| x.ln_pdf(arg);
-        test_absolute(0.1, 0.1, 0.0234154207226588982471f64.ln(), 1e-15, ln_pdf(1.0));
+        test_absolute(0.1, 0.1, 0.0234154207226588982471f64.ln(), 3e-15, ln_pdf(1.0));
         test_absolute(1.0, 0.1, 0.0396064560910663979961f64.ln(), 1e-15, ln_pdf(1.0));
         test_absolute(10.0, 0.1, 0.0418440630400545297349f64.ln(), 1e-13, ln_pdf(1.0));
         test_absolute(0.1, 1.0, 0.0396064560910663979961f64.ln(), 1e-15, ln_pdf(1.0));
         test_absolute(1.0, 1.0, 0.1591549430918953357689f64.ln(), 1e-15, ln_pdf(1.0));
         test_absolute(10.0, 1.0, 0.230361989229138647108f64.ln(), 1e-15, ln_pdf(1.0));
-        test_exact(0.1, 0.1, 0.00221546909694001013517f64.ln(), ln_pdf(10.0));
+        test_absolute(0.1, 0.1, 0.00221546909694001013517f64.ln(), 3e-15, ln_pdf(10.0));
         test_absolute(1.0, 0.1, 0.00369960370387922619592f64.ln(), 1e-15, ln_pdf(10.0));
         test_absolute(10.0, 0.1, 0.00390179721174142927402f64.ln(), 1e-13, ln_pdf(10.0));
         test_absolute(0.1, 1.0, 0.00319864073359931548273f64.ln(), 1e-15, ln_pdf(10.0));

src/distribution/gamma.rs

@@ -184,6 +184,37 @@ impl ContinuousCDF<f64, f64> for Gamma {
         }
     }
 
+    /// Tail-accurate log of the cdf via [`gamma::ln_gamma_lr`]; stays finite
+    /// far past where `cdf` underflows.
+    fn ln_cdf(&self, x: f64) -> f64 {
+        if x <= 0.0 {
+            f64::NEG_INFINITY
+        } else if prec::ulps_eq!(x, self.shape) && self.rate.is_infinite() {
+            0.0
+        } else if self.rate.is_infinite() {
+            f64::NEG_INFINITY
+        } else if x.is_infinite() {
+            0.0
+        } else {
+            gamma::ln_gamma_lr(self.shape, x * self.rate)
+        }
+    }
+
+    /// Tail-accurate log of the survival function via [`gamma::ln_gamma_ur`].
+    fn ln_sf(&self, x: f64) -> f64 {
+        if x <= 0.0 {
+            0.0
+        } else if prec::ulps_eq!(x, self.shape) && self.rate.is_infinite() {
+            f64::NEG_INFINITY
+        } else if self.rate.is_infinite() {
+            0.0
+        } else if x.is_infinite() {
+            f64::NEG_INFINITY
+        } else {
+            gamma::ln_gamma_ur(self.shape, x * self.rate)
+        }
+    }
+
     fn inverse_cdf(&self, p: f64) -> f64 {
         if !(0.0..=1.0).contains(&p) {
             panic!("default inverse_cdf implementation should be provided probability on [0,1]")