feat: add specialised inverse cdf implementation for geometric distribution

src/distribution/geometric.rs

@@ -2,6 +2,7 @@ use crate::distribution::{Discrete, DiscreteCDF};
 use crate::prec;
 use crate::statistics::*;
 use core::f64;
+use num_traits::{One, Zero};
 
 /// Implements the
 /// [Geometric](https://en.wikipedia.org/wiki/Geometric_distribution)
@@ -153,6 +154,37 @@ impl DiscreteCDF<u64, f64> for Geometric {
             ((-self.p).ln_1p() * (x as f64)).exp()
         }
     }
+
+    /// Calculates the inverse cumulative distribution function for the
+    /// geometric distribution at `x`.
+    /// In other languages, such as R, this is known as the the quantile function.
+    ///
+    /// # Formula
+    ///
+    /// ```text
+    /// p = ceil(log(1-x)/log(1-p))
+    /// ```
+    ///
+    /// # Panics
+    /// panics if provided `x` not on interval [0.0, 1.0]
+    fn inverse_cdf(&self, x: f64) -> u64 {
+        // ceil(log(1-x)/log(1-self.p)) = ceil(log1p(-x)/log1p(-self.p))
+        //                              = ceil((-x).ln_1p()/(-self.p).ln_1p())
+        //                              = ((-x).ln_1p()/(-self.p).ln_1p()).ceil()
+        if x <= self.cdf(self.min()) {
+            // if p = 1 this branch will always be taken
+            // no matter the value of x
+            return self.min();
+        }
+        if x == <f64>::one() {
+            return self.max();
+        }
+        if !(<f64>::zero()..=<f64>::one()).contains(&x) {
+            core::panic!("p must be on [0, 1]")
+        }
+
+        ((-x).ln_1p() / (-self.p).ln_1p()).ceil() as u64
+    }
 }
 
 impl Min<u64> for Geometric {
@@ -520,5 +552,14 @@ mod tests {
         test_exact(1., 1, invcdf(1.));
         test_exact(0.2, 1, invcdf(0.2));
         test_exact(0.004, 173, invcdf(0.5));
+        test_exact(0.5, u64::MAX, invcdf(1.));
+        test_exact(0.5, 2, invcdf(0.75));
+    }
+
+    #[test]
+    #[should_panic]
+    fn test_inverse_cdf_panic() {
+        let invcdf = |arg: f64| move |x: Geometric| x.inverse_cdf(arg);
+        test_exact(1., 1, invcdf(2.));
     }
 }

src/function/kernel.rs

@@ -28,7 +28,7 @@
 //! assert!((e.evaluate(0.0) - 0.75).abs() < 1e-12);
 //! ```
 
-use std::f64::consts::{FRAC_PI_2, PI};
+use core::f64::consts::{FRAC_PI_2, PI};
 
 /// Common interface for kernel functions used in KDE and smoothing.
 pub trait Kernel {
@@ -66,11 +66,7 @@ pub struct Epanechnikov;
 impl Kernel for Epanechnikov {
     fn evaluate(&self, x: f64) -> f64 {
         let a = x.abs();
-        if a <= 1.0 {
-            0.75 * (1.0 - a * a)
-        } else {
-            0.0
-        }
+        if a <= 1.0 { 0.75 * (1.0 - a * a) } else { 0.0 }
     }
 
     fn support(&self) -> Option<(f64, f64)> {
@@ -85,11 +81,7 @@ pub struct Triangular;
 impl Kernel for Triangular {
     fn evaluate(&self, x: f64) -> f64 {
         let a = x.abs();
-        if a <= 1.0 {
-            1.0 - a
-        } else {
-            0.0
-        }
+        if a <= 1.0 { 1.0 - a } else { 0.0 }
     }
 
     fn support(&self) -> Option<(f64, f64)> {
@@ -143,11 +135,7 @@ pub struct Uniform;
 
 impl Kernel for Uniform {
     fn evaluate(&self, x: f64) -> f64 {
-        if x.abs() <= 1.0 {
-            0.5
-        } else {
-            0.0
-        }
+        if x.abs() <= 1.0 { 0.5 } else { 0.0 }
     }
 
     fn support(&self) -> Option<(f64, f64)> {