Fix pdf/cdf domain handling in dist (StudentT negative cdf, Binomial overflow/infinite loop, off-support inputs)

src/statistics/dist.rs

@@ -306,7 +306,7 @@ use crate::special::function::*;
 use crate::traits::fp::FPVector;
 //use statistics::rand::ziggurat;
 use self::WeightedUniformError::*;
-use crate::statistics::{ops::C, stat::Statistics};
+use crate::statistics::stat::Statistics;
 use crate::structure::matrix::{matrix, Matrix, Row};
 use crate::util::non_macro::{linspace, seq};
 use crate::util::useful::{auto_zip, find_interval};
@@ -653,10 +653,13 @@ impl<T: PartialOrd + SampleUniform + Copy + Into<f64>> RNG for OPDist<T> {
     fn pdf<S: PartialOrd + SampleUniform + Copy + Into<f64>>(&self, x: S) -> f64 {
         match self {
             Bernoulli(prob) => {
-                if x.into() == 1f64 {
+                let x: f64 = x.into();
+                if x == 1f64 {
                     (*prob).into()
-                } else {
+                } else if x == 0f64 {
                     1f64 - (*prob).into()
+                } else {
+                    0f64
                 }
             }
             StudentT(nu) => {
@@ -690,7 +693,7 @@ impl<T: PartialOrd + SampleUniform + Copy + Into<f64>> RNG for OPDist<T> {
                     let x_t = nu / (x.powi(2) + nu);
                     1f64 - 0.5 * inc_beta(even_nu, 0.5, x_t)
                 } else if x < 0f64 {
-                    self.cdf(-x) - 0.5
+                    1f64 - self.cdf(-x)
                 } else {
                     0.5
                 }
@@ -814,8 +817,22 @@ impl<T: PartialOrd + SampleUniform + Copy + Into<f64>> RNG for TPDist<T> {
             Binomial(n, mu) => {
                 let n = *n;
                 let mu = (*mu).into();
-                let m = x.into() as usize;
-                (C(n, m) as f64) * mu.powi(m as i32) * (1f64 - mu).powi((n - m) as i32)
+                let x: f64 = x.into();
+                if x < 0f64 || x > n as f64 || x.fract() != 0f64 {
+                    return 0f64;
+                }
+                let m = x as usize;
+                // log-space: C(n, m) as an integer overflows from n = 63
+                let ln_c = ln_gamma(n as f64 + 1f64)
+                    - ln_gamma(m as f64 + 1f64)
+                    - ln_gamma((n - m) as f64 + 1f64);
+                let succ = if m == 0 { 0f64 } else { m as f64 * mu.ln() };
+                let fail = if m == n {
+                    0f64
+                } else {
+                    (n - m) as f64 * (1f64 - mu).ln()
+                };
+                (ln_c + succ + fail).exp()
             }
             Normal(m, s) => {
                 let mean = (*m).into();
@@ -869,8 +886,14 @@ impl<T: PartialOrd + SampleUniform + Copy + Into<f64>> RNG for TPDist<T> {
                 let n = *n;
                 let p = (*mu).into();
                 let q = 1f64 - p;
-                let k: f64 = x;
-                inc_beta(n as f64 - k, k + 1f64, q)
+                let k = x.floor();
+                if k < 0f64 {
+                    0f64
+                } else if k >= n as f64 {
+                    1f64
+                } else {
+                    inc_beta(n as f64 - k, k + 1f64, q)
+                }
             }
             Normal(m, s) => phi((x - (*m).into()) / (*s).into()),
             Beta(a, b) => {

src/statistics/ops.rs

@@ -11,10 +11,13 @@
 ///
 /// assert_eq!(factorial(5), 120);
 /// ```
+///
+/// # Panics
+/// Panics if the result overflows `usize` (`n > 20` on 64-bit).
 pub fn factorial(n: usize) -> usize {
     let mut p = 1usize;
     for i in 1..(n + 1) {
-        p *= i;
+        p = p.checked_mul(i).expect("factorial: overflow");
     }
     p
 }
@@ -36,32 +39,43 @@ pub fn double_factorial(n: usize) -> usize {
     let mut s = 1usize;
     let mut n = n;
     while n >= 2 {
-        s *= n;
+        s = s.checked_mul(n).expect("double_factorial: overflow");
         n -= 2;
     }
     s
 }
 
 /// Permutation
 ///
+/// Returns 0 when `r > n` (the falling factorial contains a zero factor).
+///
 /// # Usage
 ///
 /// ```
 /// extern crate peroxide;
 /// use peroxide::fuga::*;
 ///
 /// assert_eq!(P(5,3), 60);
+/// assert_eq!(P(3,5), 0);
 /// ```
+///
+/// # Panics
+/// Panics if the result overflows `usize`.
 #[allow(non_snake_case)]
 pub fn P(n: usize, r: usize) -> usize {
+    if r > n {
+        return 0;
+    }
     let mut p = 1usize;
     for i in 0..r {
-        p *= n - i;
+        p = p.checked_mul(n - i).expect("P: overflow");
     }
     p
 }
 
-/// Combination
+/// Combination (binomial coefficient)
+///
+/// Returns 0 when `r > n`.
 ///
 /// # Usage
 ///
@@ -70,14 +84,26 @@ pub fn P(n: usize, r: usize) -> usize {
 /// use peroxide::fuga::*;
 ///
 /// assert_eq!(C(10, 9), 10);
+/// assert_eq!(C(30, 15), 155117520);
+/// assert_eq!(C(5, 7), 0);
 /// ```
+///
+/// # Panics
+/// Panics if an intermediate product overflows `usize`. The multiplicative
+/// formula keeps intermediates no larger than `r * C(n, r)`, so every
+/// `C(n, r)` with `n <= 62` is exact on 64-bit.
 #[allow(non_snake_case)]
 pub fn C(n: usize, r: usize) -> usize {
-    if r > n / 2 {
-        return C(n, n - r);
+    if r > n {
+        return 0;
     }
-
-    P(n, r) / factorial(r)
+    let r = r.min(n - r);
+    let mut c = 1usize;
+    for i in 1..=r {
+        // c * (n - r + i) is divisible by i: c = C(n - r + i - 1, i - 1) here
+        c = c.checked_mul(n - r + i).expect("C: overflow") / i;
+    }
+    c
 }
 
 /// Combination with Repetition

tests/dist.rs

@@ -11,6 +11,65 @@ fn test_binomial() {
     assert!(nearly_eq(b.var(), 16f64));
 }
 
+// Reference values below cross-checked against scipy 1.17.1
+
+#[test]
+fn test_bernoulli_pdf_support() {
+    let b = Bernoulli(0.3);
+    assert!(nearly_eq(b.pdf(1.0), 0.3));
+    assert!(nearly_eq(b.pdf(0.0), 0.7));
+    // off-support values returned 1 - p before
+    assert_eq!(b.pdf(-1.0), 0f64);
+    assert_eq!(b.pdf(0.5), 0f64);
+    assert_eq!(b.pdf(2.0), 0f64);
+}
+
+#[test]
+fn test_studentt_cdf_negative_half() {
+    let t3 = StudentT(3f64);
+    // negative branch used cdf(-x) - 0.5 instead of 1 - cdf(-x)
+    assert!(nearly_eq(t3.cdf(-5.0), 0.007696219036651147));
+    assert!(nearly_eq(t3.cdf(-2.0), 0.06966298427942152));
+    assert!(nearly_eq(t3.cdf(-1.0), 0.19550110947788524));
+    assert!(nearly_eq(t3.cdf(-0.001), 0.49963244748473046));
+    assert!(nearly_eq(t3.cdf(0.0), 0.5));
+    assert!(nearly_eq(t3.cdf(1.0), 1.0 - t3.cdf(-1.0)));
+    // StudentT(1) is Cauchy: cdf(-1) = 1/4 exactly
+    assert!(nearly_eq(StudentT(1f64).cdf(-1.0), 0.25));
+}
+
+#[test]
+fn test_binomial_pdf() {
+    let b = Binomial(10, 0.3);
+    assert!(nearly_eq(b.pdf(0.0), 0.0282475249));
+    assert!(nearly_eq(b.pdf(3.0), 0.2668279320));
+    assert!(nearly_eq(b.pdf(10.0), 5.9049e-6));
+    // pdf(x) for x > n looped forever (C(n, n - r) usize underflow)
+    assert_eq!(b.pdf(12.0), 0f64);
+    assert_eq!(b.pdf(-1.0), 0f64);
+    assert_eq!(b.pdf(5.5), 0f64);
+    // C(30, 15) overflowed u64: pdf was 0.0131 in release, panic in debug
+    assert!(nearly_eq(Binomial(30, 0.5).pdf(15.0), 0.14446444809436781));
+    // far past u64 range
+    assert!(nearly_eq(Binomial(100, 0.5).pdf(50.0), 0.07958923738717871));
+    // degenerate p
+    assert!(nearly_eq(Binomial(5, 0.0).pdf(0.0), 1.0));
+    assert!(nearly_eq(Binomial(5, 1.0).pdf(5.0), 1.0));
+}
+
+#[test]
+fn test_binomial_cdf_domain() {
+    let b = Binomial(10, 0.3);
+    assert!(nearly_eq(b.cdf(3.0), 0.6496107184));
+    // k = n, k > n, k < 0 panicked in inc_beta before
+    assert!(nearly_eq(b.cdf(10.0), 1.0));
+    assert!(nearly_eq(b.cdf(11.0), 1.0));
+    assert_eq!(b.cdf(-1.0), 0f64);
+    // fractional k floors (was interpolated)
+    assert!(nearly_eq(b.cdf(5.5), 0.9526510126));
+    assert!(nearly_eq(Binomial(100, 0.5).cdf(50.0), 0.5397946186935891));
+}
+
 #[test]
 fn test_dirichlet() {
     let dir = MVDist::Dirichlet(vec![1.0, 2.0, 3.0]);

tests/ops.rs

@@ -0,0 +1,31 @@
+extern crate peroxide;
+use peroxide::fuga::*;
+
+#[test]
+fn test_combination() {
+    assert_eq!(C(10, 0), 1);
+    assert_eq!(C(10, 9), 10);
+    assert_eq!(C(10, 10), 1);
+    assert_eq!(C(52, 5), 2598960);
+    // overflowed u64 through P(30, 15) before (returned 14052247 in release)
+    assert_eq!(C(30, 15), 155117520);
+    assert_eq!(C(61, 30), 232714176627630544);
+}
+
+#[test]
+fn test_combination_r_greater_than_n() {
+    // C(n, n - r) with r > n underflowed usize -> infinite recursion in release
+    assert_eq!(C(10, 12), 0);
+    assert_eq!(C(0, 1), 0);
+    assert_eq!(C(5, 100), 0);
+}
+
+#[test]
+fn test_permutation() {
+    assert_eq!(P(5, 0), 1);
+    assert_eq!(P(5, 3), 60);
+    assert_eq!(P(5, 5), 120);
+    assert_eq!(P(30, 10), 109027350432000);
+    assert_eq!(P(20, 20), 2432902008176640000);
+    assert_eq!(P(3, 5), 0);
+}