Add some more inverse CDF functions

src/distribution/beta.rs

@@ -164,6 +164,26 @@
             beta::beta_reg(self.shape_b, self.shape_a, 1.0 - x)
         }
     }
+
+    /// Calculates the inverse cumulative distribution function for the beta
+    /// distribution
+    /// at `x`
+    ///
+    /// # Formula
+    ///
+    /// ```text
+    /// I^{-1}_x(α, β)
+    /// ```
+    ///
+    /// where `α` is shapeA, `β` is shapeB, and `I_x` is the inverse of the
+    /// regularized lower incomplete beta function
+    fn inverse_cdf(&self, x: f64) -> f64 {
+        if !(0.0..=1.0).contains(&x) {
+            panic!("x must be in [0, 1]");
+        } else {
+            beta::inv_beta_reg(self.shape_a, self.shape_b, x)
+        }
+    }
 }
 
 impl Min<f64> for Beta {
@@ -657,6 +677,27 @@
         }
     }
 
+    #[test]
+    fn test_inverse_cdf() {
+        // let inverse_cdf = |arg: f64| move |x: Beta| x.inverse_cdf(arg);
+        let func = |arg: f64| move |x: Beta| x.inverse_cdf(x.cdf(arg));
+        let test = [
+            ((1.0, 1.0), 0.0, 0.0),
+            ((1.0, 1.0), 0.5, 0.5),
+            ((1.0, 1.0), 1.0, 1.0),
+            ((9.0, 1.0), 0.0, 0.0),
+            ((9.0, 1.0), 0.001953125, 0.001953125),
+            ((9.0, 1.0), 0.5, 0.5),
+            ((9.0, 1.0), 1.0, 1.0),
+            ((5.0, 100.0), 0.0, 0.0),
+            ((5.0, 100.0), 0.01, 0.01),
+            ((5.0, 100.0), 1.0, 1.0),
+        ];
+        for ((a, b), x, expect) in test {
+           test_case(a, b, expect, func(x));
+        };
+    }
+
     #[test]
     fn test_cdf_input_lt_0() {
         let cdf = |arg: f64| move |x: Beta| x.cdf(arg);

src/distribution/cauchy.rs

@@ -121,6 +121,24 @@
     fn sf(&self, x: f64) -> f64 {
         (1.0 / f64::consts::PI) * ((self.location - x) / self.scale).atan() + 0.5
     }
+
+    /// Calculates the inverse cumulative distribution function for the
+    /// cauchy distribution at `x`
+    ///
+    /// # Formula
+    ///
+    /// ```text
+    /// x_0 + γ tan((x - 0.5) π)
+    /// ```
+    ///
+    /// where `x_0` is the location and `γ` is the scale
+    fn inverse_cdf(&self, x: f64) -> f64 {
+        if !(0.0..=1.0).contains(&x) {
+            panic!("x must be in [0, 1]");
+        } else {
+            self.location + self.scale * (f64::consts::PI * (x - 0.5)).tan()
+        }
+    }
 }
 
 impl Min<f64> for Cauchy {
@@ -466,6 +484,31 @@
         test_case(f64::INFINITY, 1.0, 1.0, sf(5.0));
     }
 
+    #[test]
+    fn test_inverse_cdf() {
+        let func = |arg: f64| move |x: Cauchy| x.inverse_cdf(x.cdf(arg));
+        test_almost(0.0, 0.1, -5.0, 1e-10, func(-5.0));
+        test_almost(0.0, 0.1, -1.0, 1e-14, func(-1.0));
+        test_case(0.0, 0.1, 0.0, func(0.0));
+        test_almost(0.0, 0.1, 1.0, 1e-14, func(1.0));
+        test_almost(0.0, 0.1, 5.0, 1e-10, func(5.0));
+        test_almost(0.0, 1.0, -5.0, 1e-14, func(-5.0));
+        test_almost(0.0, 1.0, -1.0, 1e-15, func(-1.0));
+        test_case(0.0, 1.0, 0.0, func(0.0));
+        test_almost(0.0, 1.0, 1.0, 1e-15, func(1.0));
+        test_almost(0.0, 1.0, 5.0, 1e-14, func(5.0));
+        test_almost(0.0, 10.0, -5.0, 1e-14, func(-5.0));
+        test_almost(0.0, 10.0, -1.0, 1e-14, func(-1.0));
+        test_case(0.0, 10.0, 0.0, func(0.0));
+        test_almost(0.0, 10.0, 1.0, 1e-14, func(1.0));
+        test_almost(0.0, 10.0, 5.0, 1e-14, func(5.0));
+        test_case(-5.0, 100.0, -5.0, func(-5.0));
+        test_almost(-5.0, 100.0, -1.0, 1e-10, func(-1.0));
+        test_almost(-5.0, 100.0, 0.0, 1e-14, func(0.0));
+        test_almost(-5.0, 100.0, 1.0, 1e-14, func(1.0));
+        test_almost(-5.0, 100.0, 5.0, 1e-10, func(5.0));
+    }
+
     #[test]
     fn test_continuous() {
         test::check_continuous_distribution(&try_create(-1.2, 3.4), -1500.0, 1500.0);

src/distribution/chi_squared.rs

@@ -138,6 +138,21 @@
     fn sf(&self, x: f64) -> f64 {
         self.g.sf(x)
     }
+
+    /// Calculates the inverse cumulative distribution function for the
+    /// chi-squared distribution at `x`
+    ///
+    /// # Formula
+    ///
+    /// ```text
+    /// γ^{-1}(k / 2, x * Γ(k / 2) / 2)
+    /// ```
+    ///
+    /// where `k` is the degrees of freedom, `Γ` is the gamma function,
+    /// and `γ` is the lower incomplete gamma function
+    fn inverse_cdf(&self, p: f64) -> f64 {
+        self.g.inverse_cdf(p)
+    }
 }
 
 impl Min<f64> for ChiSquared {

src/distribution/erlang.rs

@@ -122,6 +122,21 @@
     fn sf(&self, x: f64) -> f64 {
         self.g.sf(x)
     }
+
+    /// Calculates the inverse cumulative distribution function for the erlang
+    /// distribution at `x`
+    ///
+    /// # Formula
+    ///
+    /// ```text
+    /// γ^{-1}(k, (k - 1)! x) / λ
+    /// ```
+    ///
+    /// where `k` is the shape, `λ` is the rate, and `γ` is the upper
+    /// incomplete gamma function
+    fn inverse_cdf(&self, p: f64) -> f64 {
+        self.g.inverse_cdf(p)
+    }
 }
 
 impl Min<f64> for Erlang {

src/distribution/fisher_snedecor.rs

@@ -155,6 +155,28 @@
             )
         }
     }
+
+    /// Calculates the inverse cumulative distribution function for the
+    /// fisher-snedecor distribution at `x`
+    ///
+    /// # Formula
+    ///
+    /// ```text
+    /// z = I^{-1}_(x)(d1 / 2, d2 / 2)
+    /// d2 / (d1 (1 / z - 1))
+    /// ```
+    ///
+    /// where `d1` is the first degree of freedom, `d2` is
+    /// the second degree of freedom, and `I` is the regularized incomplete
+    /// beta function
+    fn inverse_cdf(&self, x: f64) -> f64 {
+        if !(0.0..=1.0).contains(&x) {
+            panic!("x must be in [0, 1]");
+        } else {
+            let z = beta::inv_beta_reg(self.freedom_1 / 2.0, self.freedom_2 / 2.0, x);
+            self.freedom_2 / (self.freedom_1 * (1.0 / z - 1.0))
+        }
+    }
 }
 
 impl Min<f64> for FisherSnedecor {
@@ -586,6 +608,23 @@
         test_almost(10.0, 1.0, 0.65910686769794, 1e-12, sf(1.0));
     }
 
+    #[test]
+    fn test_inverse_cdf() {
+        let func = |arg: f64| move |x: FisherSnedecor| x.inverse_cdf(x.cdf(arg));
+        test_almost(0.1, 0.1, 0.1, 1e-12, func(0.1));
+        test_almost(1.0, 0.1, 0.1, 1e-12, func(0.1));
+        test_almost(10.0, 0.1, 0.1, 1e-12, func(0.1));
+        test_almost(0.1, 1.0, 0.1, 1e-12, func(0.1));
+        test_almost(1.0, 1.0, 0.1, 1e-12, func(0.1));
+        test_almost(10.0, 1.0, 0.1, 1e-12, func(0.1));
+        test_almost(0.1, 0.1, 1.0, 1e-13, func(1.0));
+        test_almost(1.0, 0.1, 1.0, 1e-12, func(1.0));
+        test_almost(10.0, 0.1, 1.0, 1e-12, func(1.0));
+        test_almost(0.1, 1.0, 1.0, 1e-12, func(1.0));
+        test_almost(1.0, 1.0, 1.0, 1e-12, func(1.0));
+        test_almost(10.0, 1.0, 1.0, 1e-12, func(1.0));
+    }
+
     #[test]
     fn test_sf_lower_bound() {
         let sf = |arg: f64| move |x: FisherSnedecor| x.sf(arg);

src/distribution/pareto.rs

@@ -141,6 +141,24 @@
             (self.scale / x).powf(self.shape)
         }
     }
+
+    /// Calculates the inverse cumulative distribution function for the Pareto
+    /// distribution at `x`
+    ///
+    /// # Formula
+    ///
+    /// ```text
+    /// x_m / (1 - x)^(1 / α)
+    /// ```
+    ///
+    /// where `x_m` is the scale and `α` is the shape
+    fn inverse_cdf(&self, p: f64) -> f64 {
+        if !(0.0..=1.0).contains(&p) {
+            panic!("x must be in [0, 1]");
+        } else {
+            self.scale * (1.0 - p).powf(-1.0 / self.shape)
+        }
+    }
 }
 
 impl Min<f64> for Pareto {
@@ -524,12 +542,23 @@
         test_case(1.0, 1.0, 1.0, sf(1.0));
         test_case(5.0, 5.0, 1.0, sf(2.0));
         test_almost(7.0, 7.0, 0.08235429999999999, 1e-14, sf(10.0));
-        test_almost(10.0, 10.0, 0.16150558288984573, 1e14, sf(12.0));
+        test_almost(10.0, 10.0, 0.16150558288984573, 1e-14, sf(12.0));
         test_case(5.0, 1.0, 0.5, sf(10.0));
         test_almost(3.0, 10.0, 0.0009765625, 1e-14, sf(6.0));
         test_case(1.0, 1.0, 0.0, sf(f64::INFINITY));
     }
 
+    #[test]
+    fn test_inverse_cdf() {
+        let func = |arg: f64| move |x: Pareto| x.inverse_cdf(x.cdf(arg));
+        test_case(0.1, 0.1, 0.1, func(0.1));
+        test_case(1.0, 1.0, 1.0, func(1.0));
+        test_case(7.0, 7.0, 10.0, func(10.0));
+        test_case(10.0, 10.0, 12.0, func(12.0));
+        test_case(5.0, 1.0, 10.0, func(10.0));
+        test_case(3.0, 10.0, 6.0, func(6.0));
+    }
+
     #[test]
     fn test_continuous() {
         test::check_continuous_distribution(&try_create(1.0, 10.0), 1.0, 10.0);

src/distribution/triangular.rs

@@ -84,7 +84,7 @@
     /// } if min < x <= mode {
     ///     (x - min)^2 / ((max - min) * (mode - min))
     /// } else if mode < x < max {
-    ///     1 - (max - min)^2 / ((max - min) * (max - mode))
+    ///     1 - (max - x)^2 / ((max - min) * (max - mode))
     /// } else {
     ///     1
     /// }
@@ -134,6 +134,34 @@
             0.0
         }
     }
+
+    /// Calculates the inverse cumulative distribution function for the triangular
+    /// distribution
+    /// at `x`
+    ///
+    /// # Formula
+    ///
+    /// ```text
+    /// if x < (mode - min) / (max - min) {
+    ///     min + ((max - min) * (mode - min) * x)^(1 / 2)
+    /// } else {
+    ///     max - ((max - min) * (max - mode) * (1 - x))^(1 / 2)
+    /// }
+    /// ```
+    fn inverse_cdf(&self, p: f64) -> f64 {
+        let a = self.min;
+        let b = self.max;
+        let c = self.mode;
+        if !(0.0..=1.0).contains(&p) {
+            panic!("x must be in [0, 1]");
+        }
+
+        if p < (c - a) / (b - a) {
+            a + ((c - a) * (b - a) * p).sqrt()
+        } else {
+            b - ((b - a) * (b - c) * (1.0 - p)).sqrt()
+        }
+    }
 }
 
 impl Min<f64> for Triangular {
@@ -538,6 +566,20 @@
         test_case(0.0, 3.0, 1.5, 0.0, sf(5.0));
     }
 
+    #[test]
+    fn test_inverse_cdf() {
+        let func = |arg: f64| move |x: Triangular| x.inverse_cdf(x.cdf(arg));
+        test_almost(0.0, 1.0, 0.5, 0.25, 1e-15, func(0.25));
+        test_almost(0.0, 1.0, 0.5, 0.5, 1e-15, func(0.5));
+        test_almost(0.0, 1.0, 0.5, 0.75, 1e-15, func(0.75));
+        test_almost(-5.0, 8.0, -3.5, -4.0, 1e-15, func(-4.0));
+        test_almost(-5.0, 8.0, -3.5, -3.5, 1e-15, func(-3.5));
+        test_almost(-5.0, 8.0, -3.5, 4.0, 1e-15, func(4.0));
+        test_almost(-5.0, -3.0, -4.0, -4.5, 1e-15, func(-4.5));
+        test_almost(-5.0, -3.0, -4.0, -4.0, 1e-15, func(-4.0));
+        test_almost(-5.0, -3.0, -4.0, -3.5, 1e-15, func(-3.5));
+    }
+
     #[test]
     fn test_continuous() {
         test::check_continuous_distribution(&try_create(-5.0, 5.0, 0.0), -5.0, 5.0);

src/distribution/weibull.rs

@@ -139,6 +139,24 @@
             (-x.powf(self.shape) * self.scale_pow_shape_inv).exp()
         }
     }
+
+    /// Calculates the inverse cumulative distribution function for the weibull
+    /// distribution at `x`
+    ///
+    /// # Formula
+    ///
+    /// ```text
+    /// λ (-ln(1 - x))^(1 / k)
+    /// ```
+    ///
+    /// where `k` is the shape and `λ` is the scale
+    fn inverse_cdf(&self, p: f64) -> f64 {
+        if !(0.0..=1.0).contains(&p) {
+            panic!("x must be in [0, 1]");
+        }
+
+        (-((-p).ln_1p() / self.scale_pow_shape_inv)).powf(1.0 / self.shape)
+    }
 }
 
 impl Min<f64> for Weibull {
@@ -526,6 +544,21 @@
         test_case(10.0, 1.0, 0.0, sf(10.0));
     }
 
+    #[test]
+    fn test_inverse_cdf() {
+        let func = |arg: f64| move |x: Weibull| x.inverse_cdf(x.cdf(arg));
+        test_case(1.0, 0.1, 0.0, func(0.0));
+        test_almost(1.0, 0.1, 1.0, 1e-13, func(1.0));
+        test_case(1.0, 1.0, 0.0, func(0.0));
+        test_case(1.0, 1.0, 1.0, func(1.0));
+        test_almost(1.0, 1.0, 10.0, 1e-10, func(10.0));
+        test_case(10.0, 10.0, 0.0, func(0.0));
+        test_almost(10.0, 10.0, 1.0, 1e-5, func(1.0));
+        test_almost(10.0, 10.0, 10.0, 1e-10, func(10.0));
+        test_case(10.0, 1.0, 0.0, func(0.0));
+        test_case(10.0, 1.0, 1.0, func(1.0));
+    }
+
     #[test]
     fn test_continuous() {
         test::check_continuous_distribution(&try_create(1.0, 0.2), 0.0, 10.0);