Repository: statrs-dev/statrs (language: rust)

Issue:
StudentT::inverse_cdf() ignores location and scale parameters
The current implementation hard-codes them as 0 and 1 respectively.  Obviously the right thing would be to use the actual parameters.  At minimum we should document the current behaviour.


Discussion / hints:

Fix patch:
diff --git a/Cargo.toml b/Cargo.toml
index 8278bc4c..a71228cf 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -22,6 +22,7 @@ rand = "0.7.3"
 nalgebra = "0.23.0"
 approx = "0.4.0"
 num-traits = "0.2.14"
+lazy_static = "1.4.0"
 
 [dev-dependencies]
 criterion = "0.3.3"
diff --git a/src/distribution/beta.rs b/src/distribution/beta.rs
index 6f1245fd..d384bc5b 100644
--- a/src/distribution/beta.rs
+++ b/src/distribution/beta.rs
@@ -3,8 +3,8 @@ use crate::function::{beta, gamma};
 use crate::is_zero;
 use crate::statistics::*;
 use crate::{Result, StatsError};
+use core::f64::INFINITY as INF;
 use rand::Rng;
-// use std::f64;
 
 /// Implements the [Beta](https://en.wikipedia.org/wiki/Beta_distribution)
 /// distribution
@@ -47,12 +47,15 @@ impl Beta {
     /// assert!(result.is_err());
     /// ```
     pub fn new(shape_a: f64, shape_b: f64) -> Result<Beta> {
-        let is_nan = shape_a.is_nan() || shape_b.is_nan();
-        match (shape_a, shape_b, is_nan) {
-            (_, _, true) => Err(StatsError::BadParams),
-            (_, _, false) if shape_a <= 0.0 || shape_b <= 0.0 => Err(StatsError::BadParams),
-            (_, _, false) => Ok(Beta { shape_a, shape_b }),
-        }
+        if shape_a.is_nan()
+            || shape_b.is_nan()
+            || shape_a.is_infinite() && shape_b.is_infinite()
+            || shape_a <= 0.0
+            || shape_b <= 0.0
+        {
+            return Err(StatsError::BadParams);
+        };
+        Ok(Beta { shape_a, shape_b })
     }
 
     /// Returns the shapeA (α) of the beta distribution
@@ -111,12 +114,6 @@ impl ContinuousCDF<f64, f64> for Beta {
             0.0
         } else if x >= 1.0 {
             1.0
-        } else if self.shape_a.is_infinite() && self.shape_b.is_infinite() {
-            if x < 0.5 {
-                0.0
-            } else {
-                1.0
-            }
         } else if self.shape_a.is_infinite() {
             if x < 1.0 {
                 0.0
@@ -174,9 +171,7 @@ impl Distribution<f64> for Beta {
     ///
     /// where `α` is shapeA and `β` is shapeB
     fn mean(&self) -> Option<f64> {
-        let mean = if self.shape_a.is_infinite() && self.shape_b.is_infinite() {
-            return None;
-        } else if self.shape_a.is_infinite() {
+        let mean = if self.shape_a.is_infinite() {
             1.0
         } else {
             self.shape_a / (self.shape_a + self.shape_b)
@@ -187,10 +182,6 @@ impl Distribution<f64> for Beta {
     ///
     /// # Remarks
     ///
-    /// Returns `None` if both `shape_a` and `shape_b` are
-    /// positive infinity, since this limit cannot be consistently
-    /// defined.
-    ///
     /// # Formula
     ///
     /// ```ignore
@@ -199,9 +190,7 @@ impl Distribution<f64> for Beta {
     ///
     /// where `α` is shapeA and `β` is shapeB
     fn variance(&self) -> Option<f64> {
-        let var = if self.shape_a.is_infinite() && self.shape_b.is_infinite() {
-            return None;
-        } else if self.shape_a.is_infinite() || self.shape_b.is_infinite() {
+        let var = if self.shape_a.is_infinite() || self.shape_b.is_infinite() {
             0.0
         } else {
             self.shape_a * self.shape_b
@@ -242,9 +231,7 @@ impl Distribution<f64> for Beta {
     ///
     /// where `α` is shapeA and `β` is shapeB
     fn skewness(&self) -> Option<f64> {
-        let skew = if self.shape_a.is_infinite() && self.shape_b.is_infinite() {
-            0.0
-        } else if self.shape_a.is_infinite() {
+        let skew = if self.shape_a.is_infinite() {
             -2.0
         } else if self.shape_b.is_infinite() {
             2.0
@@ -280,10 +267,7 @@ impl Mode<Option<f64>> for Beta {
     fn mode(&self) -> Option<f64> {
         // TODO: perhaps relax constraint in order to allow calculation
         // of 'anti-mode;
-        if self.shape_a <= 1.0
-            || self.shape_b <= 1.0
-            || self.shape_a.is_infinite() && self.shape_b.is_infinite()
-        {
+        if self.shape_a <= 1.0 || self.shape_b <= 1.0 {
             None
         } else if self.shape_a.is_infinite() {
             Some(1.0)
@@ -309,21 +293,15 @@ impl Continuous<f64, f64> for Beta {
     fn pdf(&self, x: f64) -> f64 {
         if x < 0.0 || x > 1.0 {
             0.0
-        } else if self.shape_a.is_infinite() && self.shape_b.is_infinite() {
-            if ulps_eq!(x, 0.5) {
-                f64::INFINITY
-            } else {
-                0.0
-            }
         } else if self.shape_a.is_infinite() {
             if ulps_eq!(x, 1.0) {
-                f64::INFINITY
+                INF
             } else {
                 0.0
             }
         } else if self.shape_b.is_infinite() {
             if is_zero(x) {
-                f64::INFINITY
+                INF
             } else {
                 0.0
             }
@@ -352,24 +330,18 @@ impl Continuous<f64, f64> for Beta {
     /// where `α` is shapeA, `β` is shapeB, and `Γ` is the gamma function
     fn ln_pdf(&self, x: f64) -> f64 {
         if x < 0.0 || x > 1.0 {
-            f64::NEG_INFINITY
-        } else if self.shape_a.is_infinite() && self.shape_b.is_infinite() {
-            if ulps_eq!(x, 0.5) {
-                f64::INFINITY
-            } else {
-                f64::NEG_INFINITY
-            }
+            -INF
         } else if self.shape_a.is_infinite() {
             if ulps_eq!(x, 1.0) {
-                f64::INFINITY
+                INF
             } else {
-                f64::NEG_INFINITY
+                -INF
             }
         } else if self.shape_b.is_infinite() {
             if is_zero(x) {
-                f64::INFINITY
+                INF
             } else {
-                f64::NEG_INFINITY
+                -INF
             }
         } else if ulps_eq!(self.shape_a, 1.0) && ulps_eq!(self.shape_b, 1.0) {
             0.0
@@ -380,14 +352,14 @@ impl Continuous<f64, f64> for Beta {
             let bb = if ulps_eq!(self.shape_a, 1.0) && is_zero(x) {
                 0.0
             } else if is_zero(x) {
-                f64::NEG_INFINITY
+                -INF
             } else {
                 (self.shape_a - 1.0) * x.ln()
             };
             let cc = if ulps_eq!(self.shape_b, 1.0) && ulps_eq!(x, 1.0) {
                 0.0
             } else if ulps_eq!(x, 1.0) {
-                f64::NEG_INFINITY
+                -INF
             } else {
                 (self.shape_b - 1.0) * (1.0 - x).ln()
             };
@@ -399,277 +371,244 @@ impl Continuous<f64, f64> for Beta {
 #[rustfmt::skip]
 #[cfg(test)]
 mod tests {
+    use super::*;
+    use crate::consts::ACC;
     use crate::distribution::internal::*;
-    use crate::distribution::{Beta, Continuous, ContinuousCDF};
     use crate::statistics::*;
-    use crate::consts::ACC;
+    use crate::testing_boiler;
 
-    fn try_create(shape_a: f64, shape_b: f64) -> Beta {
-        let n = Beta::new(shape_a, shape_b);
-        assert!(n.is_ok());
-        n.unwrap()
-    }
-
-    fn create_case(shape_a: f64, shape_b: f64) {
-        let n = try_create(shape_a, shape_b);
-        assert_eq!(n.shape_a(), shape_a);
-        assert_eq!(n.shape_b(), shape_b);
-    }
-
-    fn bad_create_case(shape_a: f64, shape_b: f64) {
-        let n = Beta::new(shape_a, shape_b);
-        assert!(n.is_err());
-    }
-
-    fn get_value<F,T>(shape_a: f64, shape_b: f64, eval: F) -> T
-    where
-        F: Fn(Beta) -> T,
-    {
-        let n = try_create(shape_a, shape_b);
-        eval(n)
-    }
-
-    fn test_case<F>(shape_a: f64, shape_b: f64, expected: f64, eval: F)
-    where F: Fn(Beta) -> f64
-    {
-        let x = get_value(shape_a, shape_b, eval);
-        assert_eq!(expected, x);
-    }
-
-    fn test_case_none<F>(shape_a: f64, shape_b: f64, eval: F)
-    where F: Fn(Beta) -> Option<f64>
-    {
-        let x = get_value(shape_a, shape_b, eval);
-        assert_eq!(None, x);
-    }
-
-    fn test_almost<F>(shape_a: f64, shape_b: f64, expected: f64, acc: f64, eval: F)
-    where F: Fn(Beta) -> f64
-    {
-        let x = get_value(shape_a, shape_b, eval);
-        assert_almost_eq!(expected, x, acc);
-    }
-
-    fn test_is_nan<F>(shape_a: f64, shape_b: f64, eval: F)
-    where
-        F: Fn(Beta) -> f64,
-    {
-        assert!(get_value(shape_a, shape_b, eval).is_nan())
-    }
+    testing_boiler!((f64, f64), Beta);
 
     #[test]
     fn test_create() {
-        create_case(1.0, 1.0);
-        create_case(9.0, 1.0);
-        create_case(5.0, 100.0);
-        create_case(1.0, f64::INFINITY);
-        create_case(f64::INFINITY, 1.0);
+        let valid = [(1.0, 1.0), (9.0, 1.0), (5.0, 100.0), (1.0, INF), (INF, 1.0)];
+        for &arg in valid.iter() {
+            try_create(arg);
+        }
     }
 
     #[test]
     fn test_bad_create() {
-        bad_create_case(0.0, 0.0);
-        bad_create_case(0.0, 0.1);
-        bad_create_case(1.0, 0.0);
-        bad_create_case(0.0, f64::INFINITY);
-        bad_create_case(f64::INFINITY, 0.0);
-        bad_create_case(f64::NAN, 1.0);
-        bad_create_case(1.0, f64::NAN);
-        bad_create_case(f64::NAN, f64::NAN);
-        bad_create_case(1.0, -1.0);
-        bad_create_case(-1.0, 1.0);
-        bad_create_case(-1.0, -1.0);
+        let invalid = [
+            (0.0, 0.0),
+            (0.0, 0.1),
+            (1.0, 0.0),
+            (0.0, INF),
+            (INF, 0.0),
+            (f64::NAN, 1.0),
+            (1.0, f64::NAN),
+            (f64::NAN, f64::NAN),
+            (1.0, -1.0),
+            (-1.0, 1.0),
+            (-1.0, -1.0),
+            (INF, INF),
+        ];
+        for &arg in invalid.iter() {
+            bad_create_case(arg);
+        }
     }
 
     #[test]
     fn test_mean() {
-        let mean = |x: Beta| x.mean().unwrap();
-        test_case(1.0, 1.0, 0.5, mean);
-        test_case(9.0, 1.0, 0.9, mean);
-        test_case(5.0, 100.0, 0.047619047619047619047616, mean);
-        test_case(1.0, f64::INFINITY, 0.0, mean);
-        test_case(f64::INFINITY, 1.0, 1.0, mean);
-        let mean = |x: Beta| x.mean();
-        test_case_none(f64::INFINITY, f64::INFINITY, mean);
+        let f = |x: Beta| x.mean().unwrap();
+        let test = [
+            ((1.0, 1.0), 0.5),
+            ((9.0, 1.0), 0.9),
+            ((5.0, 100.0), 0.047619047619047619047616),
+            ((1.0, INF), 0.0),
+            ((INF, 1.0), 1.0),
+        ];
+        for &(arg, res) in test.iter() {
+            test_case(arg, res, f);
+        }
     }
 
     #[test]
     fn test_variance() {
-        let variance = |x: Beta| x.variance().unwrap();
-        test_case(1.0, 1.0, 1.0 / 12.0, variance);
-        test_case(9.0, 1.0, 9.0 / 1100.0, variance);
-        test_case(5.0, 100.0, 500.0 / 1168650.0, variance);
-        test_case(1.0, f64::INFINITY, 0.0, variance);
-        test_case(f64::INFINITY, 1.0, 0.0, variance);
-        let variance = |x: Beta| x.variance();
-        test_case_none(f64::INFINITY, f64::INFINITY, variance);
+        let f = |x: Beta| x.variance().unwrap();
+        let test = [
+            ((1.0, 1.0), 1.0 / 12.0),
+            ((9.0, 1.0), 9.0 / 1100.0),
+            ((5.0, 100.0), 500.0 / 1168650.0),
+            ((1.0, INF), 0.0),
+            ((INF, 1.0), 0.0),
+        ];
+        for &(arg, res) in test.iter() {
+            test_case(arg, res, f);
+        }
     }
 
     #[test]
     fn test_entropy() {
-        let entropy = |x: Beta| x.entropy().unwrap();
-        test_almost(1.0, 1.0, 0.0, 1e-14, entropy);
-        test_almost(9.0, 1.0, -1.3083356884473304939016015849561625204060922267565917, 1e-13, entropy);
-        test_almost(5.0, 100.0, -2.5201623187602743679459255108827601222133603091753153, 1e-13, entropy);
+        let f = |x: Beta| x.entropy().unwrap();
+        let test = [
+            ((9.0, 1.0), -1.3083356884473304939016015),
+            ((5.0, 100.0), -2.52016231876027436794592),
+        ];
+        for &(arg, res) in test.iter() {
+            test_case(arg, res, f);
+        }
+        test_case_special((1.0, 1.0), 0.0, 1e-14, f);
         let entropy = |x: Beta| x.entropy();
-        test_case_none(1.0, f64::INFINITY, entropy);
-        test_case_none(f64::INFINITY, 1.0, entropy);
-        test_case_none(f64::INFINITY, f64::INFINITY, entropy);
+        test_none((1.0, INF), entropy);
+        test_none((INF, 1.0), entropy);
     }
 
     #[test]
     fn test_skewness() {
         let skewness = |x: Beta| x.skewness().unwrap();
-        test_case(1.0, 1.0, 0.0, skewness);
-        test_almost(9.0, 1.0, -1.4740554623801777107177478829647496373009282424841579, 1e-15, skewness);
-        test_almost(5.0, 100.0, 0.81759410927553430354583159143895018978562196953345572, 1e-15, skewness);
-        test_case(1.0, f64::INFINITY, 2.0, skewness);
-        test_case(f64::INFINITY, 1.0, -2.0, skewness);
-        test_case(f64::INFINITY, f64::INFINITY, 0.0, skewness);
+        test_case((1.0, 1.0), 0.0, skewness);
+        test_case((9.0, 1.0), -1.4740554623801777107177478829, skewness);
+        test_case((5.0, 100.0), 0.817594109275534303545831591, skewness);
+        test_case((1.0, INF), 2.0, skewness);
+        test_case((INF, 1.0), -2.0, skewness);
     }
 
     #[test]
     fn test_mode() {
         let mode = |x: Beta| x.mode().unwrap();
-        test_case(5.0, 100.0, 0.038834951456310676243255386452801758423447608947753906, mode);
-        test_case(2.0, f64::INFINITY, 0.0, mode);
-        test_case(f64::INFINITY, 2.0, 1.0, mode);
-        let mode = |x: Beta| x.mode();
-        test_case_none(f64::INFINITY, f64::INFINITY, mode);
+        test_case((5.0, 100.0), 0.038834951456310676243255386, mode);
+        test_case((92.0, INF), 0.0, mode);
+        test_case((INF, 2.0), 1.0, mode);
     }
 
     #[test]
     #[should_panic]
     fn test_mode_shape_a_lte_1() {
         let mode = |x: Beta| x.mode().unwrap();
-        get_value(1.0, 5.0, mode);
+        get_value((1.0, 5.0), mode);
     }
 
     #[test]
     #[should_panic]
     fn test_mode_shape_b_lte_1() {
         let mode = |x: Beta| x.mode().unwrap();
-        get_value(5.0, 1.0, mode);
+        get_value((5.0, 1.0), mode);
     }
 
     #[test]
     fn test_min_max() {
         let min = |x: Beta| x.min();
         let max = |x: Beta| x.max();
-        test_case(1.0, 1.0, 0.0, min);
-        test_case(1.0, 1.0, 1.0, max);
+        test_case((1.0, 1.0), 0.0, min);
+        test_case((1.0, 1.0), 1.0, max);
     }
 
-
     #[test]
     fn test_pdf() {
-        let pdf = |arg: f64| move |x: Beta| x.pdf(arg);
-        test_case(1.0, 1.0, 1.0, pdf(0.0));
-        test_case(1.0, 1.0, 1.0, pdf(0.5));
-        test_case(1.0, 1.0, 1.0, pdf(1.0));
-        test_case(9.0, 1.0, 0.0, pdf(0.0));
-        test_almost(9.0, 1.0, 0.03515625, 1e-15, pdf(0.5));
-        test_almost(9.0, 1.0, 9.0, 1e-13, pdf(1.0));
-        test_case(5.0, 100.0, 0.0, pdf(0.0));
-        test_almost(5.0, 100.0, 4.534102298350337661e-23, 1e-35, pdf(0.5));
-        test_case(5.0, 100.0, 0.0, pdf(1.0));
-        test_case(5.0, 100.0, 0.0, pdf(1.0));
-        test_case(1.0, f64::INFINITY, f64::INFINITY, pdf(0.0));
-        test_case(1.0, f64::INFINITY, 0.0, pdf(0.5));
-        test_case(1.0, f64::INFINITY, 0.0, pdf(1.0));
-        test_case(f64::INFINITY, 1.0, 0.0, pdf(0.0));
-        test_case(f64::INFINITY, 1.0, 0.0, pdf(0.5));
-        test_case(f64::INFINITY, 1.0, f64::INFINITY, pdf(1.0));
-        test_case(f64::INFINITY, f64::INFINITY, 0.0, pdf(0.0));
-        test_case(f64::INFINITY, f64::INFINITY, f64::INFINITY, pdf(0.5));
-        test_case(f64::INFINITY, f64::INFINITY, 0.0, pdf(1.0));
+        let f = |arg: f64| move |x: Beta| x.pdf(arg);
+        let test = [
+            ((1.0, 1.0), 0.0, 1.0),
+            ((1.0, 1.0), 0.5, 1.0),
+            ((1.0, 1.0), 1.0, 1.0),
+            ((9.0, 1.0), 0.0, 0.0),
+            ((9.0, 1.0), 0.5, 0.03515625),
+            ((9.0, 1.0), 1.0, 9.0),
+            ((5.0, 100.0), 0.0, 0.0),
+            ((5.0, 100.0), 0.5, 4.534102298350337661e-23),
+            ((5.0, 100.0), 1.0, 0.0),
+            ((5.0, 100.0), 1.0, 0.0),
+            ((1.0, INF), 0.0, INF),
+            ((1.0, INF), 0.5, 0.0),
+            ((1.0, INF), 1.0, 0.0),
+            ((INF, 1.0), 0.0, 0.0),
+            ((INF, 1.0), 0.5, 0.0),
+            ((INF, 1.0), 1.0, INF),
+        ];
+        for &(arg, x, expect) in test.iter() {
+            test_case(arg, expect, f(x));
+        }
     }
 
     #[test]
     fn test_pdf_input_lt_0() {
         let pdf = |arg: f64| move |x: Beta| x.pdf(arg);
-        test_case(1.0, 1.0, 0.0, pdf(-1.0));
+        test_case((1.0, 1.0), 0.0, pdf(-1.0));
     }
 
     #[test]
     fn test_pdf_input_gt_0() {
         let pdf = |arg: f64| move |x: Beta| x.pdf(arg);
-        test_case(1.0, 1.0, 0.0, pdf(2.0));
+        test_case((1.0, 1.0), 0.0, pdf(2.0));
     }
 
     #[test]
     fn test_ln_pdf() {
-        let ln_pdf = |arg: f64| move |x: Beta| x.ln_pdf(arg);
-        test_case(1.0, 1.0, 0.0, ln_pdf(0.0));
-        test_case(1.0, 1.0, 0.0, ln_pdf(0.5));
-        test_case(1.0, 1.0, 0.0, ln_pdf(1.0));
-        test_case(9.0, 1.0, f64::NEG_INFINITY, ln_pdf(0.0));
-        test_almost(9.0, 1.0, -3.3479528671433430925473664978203611353090199592365458, 1e-13, ln_pdf(0.5));
-        test_almost(9.0, 1.0, 2.1972245773362193827904904738450514092949811156454996, 1e-13, ln_pdf(1.0));
-        test_case(5.0, 100.0, f64::NEG_INFINITY, ln_pdf(0.0));
-        test_almost(5.0, 100.0, -51.447830024537682154565870837960406410586196074573801, 1e-12, ln_pdf(0.5));
-        test_case(5.0, 100.0, f64::NEG_INFINITY, ln_pdf(1.0));
-        test_case(1.0, f64::INFINITY, f64::INFINITY, ln_pdf(0.0));
-        test_case(1.0, f64::INFINITY, f64::NEG_INFINITY, ln_pdf(0.5));
-        test_case(1.0, f64::INFINITY, f64::NEG_INFINITY, ln_pdf(1.0));
-        test_case(f64::INFINITY, 1.0, f64::NEG_INFINITY, ln_pdf(0.0));
-        test_case(f64::INFINITY, 1.0, f64::NEG_INFINITY, ln_pdf(0.5));
-        test_case(f64::INFINITY, 1.0, f64::INFINITY, ln_pdf(1.0));
-        test_case(f64::INFINITY, f64::INFINITY, f64::NEG_INFINITY, ln_pdf(0.0));
-        test_case(f64::INFINITY, f64::INFINITY, f64::INFINITY, ln_pdf(0.5));
-        test_case(f64::INFINITY, f64::INFINITY, f64::NEG_INFINITY, ln_pdf(1.0));
+        let f = |arg: f64| move |x: Beta| x.ln_pdf(arg);
+        let test = [
+            ((1.0, 1.0), 0.0, 0.0),
+            ((1.0, 1.0), 0.5, 0.0),
+            ((1.0, 1.0), 1.0, 0.0),
+            ((9.0, 1.0), 0.0, -INF),
+            ((9.0, 1.0), 0.5, -3.347952867143343092547366497),
+            ((9.0, 1.0), 1.0, 2.1972245773362193827904904738),
+            ((5.0, 100.0), 0.0, -INF),
+            ((5.0, 100.0), 0.5, -51.447830024537682154565870),
+            ((5.0, 100.0), 1.0, -INF),
+            ((1.0, INF), 0.0, INF),
+            ((1.0, INF), 0.5, -INF),
+            ((1.0, INF), 1.0, -INF),
+            ((INF, 1.0), 0.0, -INF),
+            ((INF, 1.0), 0.5, -INF),
+            ((INF, 1.0), 1.0, INF),
+        ];
+        for &(arg, x, expect) in test.iter() {
+            test_case(arg, expect, f(x));
+        }
     }
 
     #[test]
     fn test_ln_pdf_input_lt_0() {
         let ln_pdf = |arg: f64| move |x: Beta| x.ln_pdf(arg);
-        test_case(1.0, 1.0, f64::NEG_INFINITY, ln_pdf(-1.0));
+        test_case((1.0, 1.0), -INF, ln_pdf(-1.0));
     }
 
     #[test]
     fn test_ln_pdf_input_gt_1() {
         let ln_pdf = |arg: f64| move |x: Beta| x.ln_pdf(arg);
-        test_case(1.0, 1.0, f64::NEG_INFINITY, ln_pdf(2.0));
+        test_case((1.0, 1.0), -INF, ln_pdf(2.0));
     }
 
     #[test]
     fn test_cdf() {
         let cdf = |arg: f64| move |x: Beta| x.cdf(arg);
-        test_case(1.0, 1.0, 0.0, cdf(0.0));
-        test_case(1.0, 1.0, 0.5, cdf(0.5));
-        test_case(1.0, 1.0, 1.0, cdf(1.0));
-        test_case(9.0, 1.0, 0.0, cdf(0.0));
-        test_almost(9.0, 1.0, 0.001953125, 1e-16, cdf(0.5));
-        test_case(9.0, 1.0, 1.0, cdf(1.0));
-        test_case(5.0, 100.0, 0.0, cdf(0.0));
-        test_case(5.0, 100.0, 1.0, cdf(0.5));
-        test_case(5.0, 100.0, 1.0, cdf(1.0));
-        test_case(1.0, f64::INFINITY, 1.0, cdf(0.0));
-        test_case(1.0, f64::INFINITY, 1.0, cdf(0.5));
-        test_case(1.0, f64::INFINITY, 1.0, cdf(1.0));
-        test_case(f64::INFINITY, 1.0, 0.0, cdf(0.0));
-        test_case(f64::INFINITY, 1.0, 0.0, cdf(0.5));
-        test_case(f64::INFINITY, 1.0, 1.0, cdf(1.0));
-        test_case(f64::INFINITY, f64::INFINITY, 0.0, cdf(0.0));
-        test_case(f64::INFINITY, f64::INFINITY, 1.0, cdf(0.5));
-        test_case(f64::INFINITY, f64::INFINITY, 1.0, cdf(1.0));
+        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.5, 0.001953125),
+            ((9.0, 1.0), 1.0, 1.0),
+            ((5.0, 100.0), 0.0, 0.0),
+            ((5.0, 100.0), 0.5, 1.0),
+            ((5.0, 100.0), 1.0, 1.0),
+            ((1.0, INF), 0.0, 1.0),
+            ((1.0, INF), 0.5, 1.0),
+            ((1.0, INF), 1.0, 1.0),
+            ((INF, 1.0), 0.0, 0.0),
+            ((INF, 1.0), 0.5, 0.0),
+            ((INF, 1.0), 1.0, 1.0),
+        ];
+        for &(arg, x, expect) in test.iter() {
+            test_case(arg, expect, cdf(x));
+        }
     }
 
     #[test]
     fn test_cdf_input_lt_0() {
         let cdf = |arg: f64| move |x: Beta| x.cdf(arg);
-        test_case(1.0, 1.0, 0.0, cdf(-1.0));
+        test_case((1.0, 1.0), 0.0, cdf(-1.0));
     }
 
     #[test]
     fn test_cdf_input_gt_1() {
         let cdf = |arg: f64| move |x: Beta| x.cdf(arg);
-        test_case(1.0, 1.0, 1.0, cdf(2.0));
+        test_case((1.0, 1.0), 1.0, cdf(2.0));
     }
 
     #[test]
     fn test_continuous() {
-        tests::check_continuous_distribution(&try_create(1.2, 3.4), 0.0, 1.0);
-        tests::check_continuous_distribution(&try_create(4.5, 6.7), 0.0, 1.0);
+        test::check_continuous_distribution(&try_create((1.2, 3.4)), 0.0, 1.0);
+        test::check_continuous_distribution(&try_create((4.5, 6.7)), 0.0, 1.0);
     }
 }
diff --git a/src/distribution/binomial.rs b/src/distribution/binomial.rs
index 35f2cf7c..285fb64d 100644
--- a/src/distribution/binomial.rs
+++ b/src/distribution/binomial.rs
@@ -512,7 +512,7 @@ mod tests {
 
     #[test]
     fn test_discrete() {
-        tests::check_discrete_distribution(&try_create(0.3, 5), 5);
-        tests::check_discrete_distribution(&try_create(0.7, 10), 10);
+        test::check_discrete_distribution(&try_create(0.3, 5), 5);
+        test::check_discrete_distribution(&try_create(0.7, 10), 10);
     }
 }
diff --git a/src/distribution/categorical.rs b/src/distribution/categorical.rs
index 46de9690..ea53eae1 100644
--- a/src/distribution/categorical.rs
+++ b/src/distribution/categorical.rs
@@ -488,7 +488,7 @@ mod tests {
 
     #[test]
     fn test_discrete() {
-        tests::check_discrete_distribution(&try_create(&[1.0, 2.0, 3.0, 4.0]), 4);
-        tests::check_discrete_distribution(&try_create(&[0.0, 1.0, 2.0, 3.0, 4.0]), 5);
+        test::check_discrete_distribution(&try_create(&[1.0, 2.0, 3.0, 4.0]), 4);
+        test::check_discrete_distribution(&try_create(&[0.0, 1.0, 2.0, 3.0, 4.0]), 5);
     }
 }
diff --git a/src/distribution/cauchy.rs b/src/distribution/cauchy.rs
index 61b0b931..dd56d28e 100644
--- a/src/distribution/cauchy.rs
+++ b/src/distribution/cauchy.rs
@@ -412,7 +412,7 @@ mod tests {
 
     #[test]
     fn test_continuous() {
-        tests::check_continuous_distribution(&try_create(-1.2, 3.4), -1500.0, 1500.0);
-        tests::check_continuous_distribution(&try_create(-4.5, 6.7), -5000.0, 5000.0);
+        test::check_continuous_distribution(&try_create(-1.2, 3.4), -1500.0, 1500.0);
+        test::check_continuous_distribution(&try_create(-4.5, 6.7), -5000.0, 5000.0);
     }
 }
diff --git a/src/distribution/chi.rs b/src/distribution/chi.rs
index 02d32995..75590026 100644
--- a/src/distribution/chi.rs
+++ b/src/distribution/chi.rs
@@ -546,8 +546,8 @@ mod tests {
 
     #[test]
     fn test_continuous() {
-        tests::check_continuous_distribution(&try_create(1.0), 0.0, 10.0);
-        tests::check_continuous_distribution(&try_create(2.0), 0.0, 10.0);
-        tests::check_continuous_distribution(&try_create(5.0), 0.0, 10.0);
+        test::check_continuous_distribution(&try_create(1.0), 0.0, 10.0);
+        test::check_continuous_distribution(&try_create(2.0), 0.0, 10.0);
+        test::check_continuous_distribution(&try_create(5.0), 0.0, 10.0);
     }
 }
diff --git a/src/distribution/chi_squared.rs b/src/distribution/chi_squared.rs
index 324ecb0c..2e1683f8 100644
--- a/src/distribution/chi_squared.rs
+++ b/src/distribution/chi_squared.rs
@@ -307,8 +307,8 @@ mod tests {
     #[test]
     fn test_continuous() {
         // TODO: figure out why this test fails:
-        //tests::check_continuous_distribution(&try_create(1.0), 0.0, 10.0);
-        tests::check_continuous_distribution(&try_create(2.0), 0.0, 10.0);
-        tests::check_continuous_distribution(&try_create(5.0), 0.0, 50.0);
+        //test::check_continuous_distribution(&try_create(1.0), 0.0, 10.0);
+        test::check_continuous_distribution(&try_create(2.0), 0.0, 10.0);
+        test::check_continuous_distribution(&try_create(5.0), 0.0, 50.0);
     }
 }
diff --git a/src/distribution/erlang.rs b/src/distribution/erlang.rs
index 38bd5272..9b1c80b0 100644
--- a/src/distribution/erlang.rs
+++ b/src/distribution/erlang.rs
@@ -293,8 +293,8 @@ mod tests {
 
     #[test]
     fn test_continuous() {
-        tests::check_continuous_distribution(&try_create(1, 2.5), 0.0, 20.0);
-        tests::check_continuous_distribution(&try_create(2, 1.5), 0.0, 20.0);
-        tests::check_continuous_distribution(&try_create(3, 0.5), 0.0, 20.0);
+        test::check_continuous_distribution(&try_create(1, 2.5), 0.0, 20.0);
+        test::check_continuous_distribution(&try_create(2, 1.5), 0.0, 20.0);
+        test::check_continuous_distribution(&try_create(3, 0.5), 0.0, 20.0);
     }
 }
diff --git a/src/distribution/exponential.rs b/src/distribution/exponential.rs
index 19e95320..0ffb64ef 100644
--- a/src/distribution/exponential.rs
+++ b/src/distribution/exponential.rs
@@ -447,8 +447,8 @@ mod tests {
 
     #[test]
     fn test_continuous() {
-        tests::check_continuous_distribution(&try_create(0.5), 0.0, 10.0);
-        tests::check_continuous_distribution(&try_create(1.5), 0.0, 20.0);
-        tests::check_continuous_distribution(&try_create(2.5), 0.0, 50.0);
+        test::check_continuous_distribution(&try_create(0.5), 0.0, 10.0);
+        test::check_continuous_distribution(&try_create(1.5), 0.0, 20.0);
+        test::check_continuous_distribution(&try_create(2.5), 0.0, 50.0);
     }
 }
diff --git a/src/distribution/fisher_snedecor.rs b/src/distribution/fisher_snedecor.rs
index e8e0c4c9..741b5772 100644
--- a/src/distribution/fisher_snedecor.rs
+++ b/src/distribution/fisher_snedecor.rs
@@ -538,6 +538,6 @@ mod tests {
 
     #[test]
     fn test_continuous() {
-        tests::check_continuous_distribution(&try_create(10.0, 10.0), 0.0, 10.0);
+        test::check_continuous_distribution(&try_create(10.0, 10.0), 0.0, 10.0);
     }
 }
diff --git a/src/distribution/gamma.rs b/src/distribution/gamma.rs
index 6c4389eb..cbdf6d92 100644
--- a/src/distribution/gamma.rs
+++ b/src/distribution/gamma.rs
@@ -2,8 +2,8 @@ use crate::distribution::{Continuous, ContinuousCDF};
 use crate::function::gamma;
 use crate::statistics::*;
 use crate::{Result, StatsError};
+use core::f64::INFINITY as INF;
 use rand::Rng;
-use std::f64;
 
 /// Implements the [Gamma](https://en.wikipedia.org/wiki/Gamma_distribution)
 /// distribution
@@ -31,7 +31,7 @@ impl Gamma {
     ///
     /// # Errors
     ///
-    /// Returns an error if `shape` or `rate` are `NaN`.
+    /// Returns an error if `shape` is 'NaN' or inf or `rate` is `NaN` or inf.
     /// Also returns an error if `shape <= 0.0` or `rate <= 0.0`
     ///
     /// # Examples
@@ -46,12 +46,15 @@ impl Gamma {
     /// assert!(result.is_err());
     /// ```
     pub fn new(shape: f64, rate: f64) -> Result<Gamma> {
-        let is_nan = shape.is_nan() || rate.is_nan();
-        match (shape, rate, is_nan) {
-            (_, _, true) => Err(StatsError::BadParams),
-            (_, _, false) if shape <= 0.0 || rate <= 0.0 => Err(StatsError::BadParams),
-            (_, _, false) => Ok(Gamma { shape, rate }),
+        if shape.is_nan()
+            || rate.is_nan()
+            || shape.is_infinite() && rate.is_infinite()
+            || shape <= 0.0
+            || rate <= 0.0
+        {
+            return Err(StatsError::BadParams);
         }
+        Ok(Gamma { shape, rate })
     }
 
     /// Returns the shape (α) of the gamma distribution
@@ -143,7 +146,7 @@ impl Max<f64> for Gamma {
     /// INF
     /// ```
     fn max(&self) -> f64 {
-        f64::INFINITY
+        INF
     }
 }
 
@@ -319,210 +322,221 @@ pub fn sample_unchecked<R: Rng + ?Sized>(rng: &mut R, shape: f64, rate: f64) ->
     }
 }
 
-#[rustfmt::skip]
 #[cfg(test)]
 mod tests {
-    use crate::statistics::*;
-    use crate::distribution::{ContinuousCDF, Continuous, Gamma};
-    use crate::distribution::internal::*;
+    use super::*;
     use crate::consts::ACC;
+    use crate::distribution::internal::*;
+    use crate::testing_boiler;
 
-    fn try_create(shape: f64, rate: f64) -> Gamma {
-        let n = Gamma::new(shape, rate);
-        assert!(n.is_ok());
-        n.unwrap()
-    }
-
-    fn create_case(shape: f64, rate: f64) {
-        let n = try_create(shape, rate);
-        assert_eq!(shape, n.shape());
-        assert_eq!(rate, n.rate());
-    }
-
-    fn bad_create_case(shape: f64, rate: f64) {
-        let n = Gamma::new(shape, rate);
-        assert!(n.is_err());
-    }
-
-    fn get_value<F>(shape: f64, rate: f64, eval: F) -> f64
-        where F: Fn(Gamma) -> f64
-    {
-        let n = try_create(shape, rate);
-        eval(n)
-    }
-
-    fn test_case<F>(shape: f64, rate: f64, expected: f64, eval: F)
-        where F: Fn(Gamma) -> f64
-    {
-        let x = get_value(shape, rate, eval);
-        assert_eq!(expected, x);
-    }
-
-    fn test_almost<F>(shape: f64, rate: f64, expected: f64, acc: f64, eval: F)
-        where F: Fn(Gamma) -> f64
-    {
-        let x = get_value(shape, rate, eval);
-        assert_almost_eq!(expected, x, acc);
-    }
-
-    fn test_is_nan<F>(shape: f64, rate: f64, eval: F)
-        where F: Fn(Gamma) -> f64
-    {
-        let x = get_value(shape, rate, eval);
-        assert!(x.is_nan());
-    }
+    testing_boiler!((f64, f64), Gamma);
 
     #[test]
     fn test_create() {
-        create_case(1.0, 0.1);
-        create_case(1.0, 1.0);
-        create_case(10.0, 10.0);
-        create_case(10.0, 1.0);
-        create_case(10.0, f64::INFINITY);
+        let valid = [
+            (1.0, 0.1),
+            (1.0, 1.0),
+            (10.0, 10.0),
+            (10.0, 1.0),
+            (10.0, INF),
+        ];
+
+        for &arg in valid.iter() {
+            try_create(arg);
+        }
     }
 
     #[test]
     fn test_bad_create() {
-        bad_create_case(0.0, 0.0);
-        bad_create_case(1.0, f64::NAN);
-        bad_create_case(1.0, -1.0);
-        bad_create_case(-1.0, 1.0);
-        bad_create_case(-1.0, -1.0);
-        bad_create_case(-1.0, f64::NAN);
+        let invalid = [
+            (0.0, 0.0),
+            (1.0, f64::NAN),
+            (1.0, -1.0),
+            (-1.0, 1.0),
+            (-1.0, -1.0),
+            (-1.0, f64::NAN),
+        ];
+        for &arg in invalid.iter() {
+            bad_create_case(arg);
+        }
     }
 
     #[test]
     fn test_mean() {
-        let mean = |x: Gamma| x.mean().unwrap();
-        test_case(1.0, 0.1, 10.0, mean);
-        test_case(1.0, 1.0, 1.0, mean);
-        test_case(10.0, 10.0, 1.0, mean);
-        test_case(10.0, 1.0, 10.0, mean);
-        test_case(10.0, f64::INFINITY, 0.0, mean);
+        let f = |x: Gamma| x.mean().unwrap();
+        let test = [
+            ((1.0, 0.1), 10.0),
+            ((1.0, 1.0), 1.0),
+            ((10.0, 10.0), 1.0),
+            ((10.0, 1.0), 10.0),
+            ((10.0, INF), 0.0),
+        ];
+        for &(arg, res) in test.iter() {
+            test_case(arg, res, f);
+        }
     }
 
     #[test]
     fn test_variance() {
-        let variance = |x: Gamma| x.variance().unwrap();
-        test_almost(1.0, 0.1, 100.0, 1e-13, variance);
-        test_case(1.0, 1.0, 1.0, variance);
-        test_case(10.0, 10.0, 0.1, variance);
-        test_case(10.0, 1.0, 10.0, variance);
-        test_case(10.0, f64::INFINITY, 0.0, variance);
-    }
-
-    #[test]
-    fn test_std_dev() {
-        let std_dev = |x: Gamma| x.std_dev().unwrap();
-        test_case(1.0, 0.1, 10.0, std_dev);
-        test_case(1.0, 1.0, 1.0, std_dev);
-        test_case(10.0, 10.0, 0.31622776601683794197697302588502426416723164097476643, std_dev);
-        test_case(10.0, 1.0, 3.1622776601683793319988935444327185337195551393252168, std_dev);
-        test_case(10.0, f64::INFINITY, 0.0, std_dev);
+        let f = |x: Gamma| x.variance().unwrap();
+        let test = [
+            ((1.0, 0.1), 100.0),
+            ((1.0, 1.0), 1.0),
+            ((10.0, 10.0), 0.1),
+            ((10.0, 1.0), 10.0),
+            ((10.0, INF), 0.0),
+        ];
+        for &(arg, res) in test.iter() {
+            test_case(arg, res, f);
+        }
     }
 
     #[test]
     fn test_entropy() {
-        let entropy = |x: Gamma| x.entropy().unwrap();
-        test_almost(1.0, 0.1, 3.3025850929940456285068402234265387271634735938763824, 1e-15, entropy);
-        test_almost(1.0, 1.0, 1.0, 1e-15, entropy);
-        test_almost(10.0, 10.0, 0.23346908548693395836262094490967812177376750477943892, 1e-13, entropy);
-        test_almost(10.0, 1.0, 2.5360541784809796423806123995940423293748689934081866, 1e-13, entropy);
-        test_case(10.0, f64::INFINITY, f64::NEG_INFINITY, entropy);
+        let f = |x: Gamma| x.entropy().unwrap();
+        let test = [
+            ((1.0, 0.1), 3.302585092994045628506840223),
+            ((1.0, 1.0), 1.0),
+            ((10.0, 10.0), 0.2334690854869339583626209),
+            ((10.0, 1.0), 2.53605417848097964238061239),
+            ((10.0, INF), f64::NEG_INFINITY),
+        ];
+        for &(arg, res) in test.iter() {
+            test_case(arg, res, f);
+        }
     }
 
     #[test]
     fn test_skewness() {
-        let skewness = |x: Gamma| x.skewness().unwrap();
-        test_case(1.0, 0.1, 2.0, skewness);
-        test_case(1.0, 1.0, 2.0, skewness);
-        test_case(10.0, 10.0, 0.63245553203367586639977870888654370674391102786504337, skewness);
-        test_case(10.0, 1.0, 0.63245553203367586639977870888654370674391102786504337, skewness);
-        test_case(10.0, f64::INFINITY, 0.63245553203367586639977870888654370674391102786504337, skewness);
+        let f = |x: Gamma| x.skewness().unwrap();
+        let test = [
+            ((1.0, 0.1), 2.0),
+            ((1.0, 1.0), 2.0),
+            ((10.0, 10.0), 0.6324555320336758663997787),
+            ((10.0, 1.0), 0.63245553203367586639977870),
+            ((10.0, INF), 0.6324555320336758),
+        ];
+        for &(arg, res) in test.iter() {
+            test_case(arg, res, f);
+        }
     }
 
     #[test]
     fn test_mode() {
-        let mode = |x: Gamma| x.mode().unwrap();
-        test_case(1.0, 0.1, 0.0, mode);
-        test_case(1.0, 1.0, 0.0, mode);
-        test_case(10.0, 10.0, 0.9, mode);
-        test_case(10.0, 1.0, 9.0, mode);
-        test_case(10.0, f64::INFINITY, 0.0, mode);
+        let f = |x: Gamma| x.mode().unwrap();
+        let test = [((1.0, 0.1), 0.0), ((1.0, 1.0), 0.0)];
+        for &(arg, res) in test.iter() {
+            test_case_special(arg, res, 10e-6, f);
+        }
+        let test = [((10.0, 10.0), 0.9), ((10.0, 1.0), 9.0), ((10.0, INF), 0.0)];
+        for &(arg, res) in test.iter() {
+            test_case(arg, res, f);
+        }
     }
 
     #[test]
     fn test_min_max() {
-        let min = |x: Gamma| x.min();
-        let max = |x: Gamma| x.max();
-        test_case(1.0, 0.1, 0.0, min);
-        test_case(1.0, 1.0, 0.0, min);
-        test_case(10.0, 10.0, 0.0, min);
-        test_case(10.0, 1.0, 0.0, min);
-        test_case(10.0, f64::INFINITY, 0.0, min);
-        test_case(1.0, 0.1, f64::INFINITY, max);
-        test_case(1.0, 1.0, f64::INFINITY, max);
-        test_case(10.0, 10.0, f64::INFINITY, max);
-        test_case(10.0, 1.0, f64::INFINITY, max);
-        test_case(10.0, f64::INFINITY, f64::INFINITY, max);
+        let f = |x: Gamma| x.min();
+        let test = [
+            ((1.0, 0.1), 0.0),
+            ((1.0, 1.0), 0.0),
+            ((10.0, 10.0), 0.0),
+            ((10.0, 1.0), 0.0),
+            ((10.0, INF), 0.0),
+        ];
+        for &(arg, res) in test.iter() {
+            test_case(arg, res, f);
+        }
+        let f = |x: Gamma| x.max();
+        let test = [
+            ((1.0, 0.1), INF),
+            ((1.0, 1.0), INF),
+            ((10.0, 10.0), INF),
+            ((10.0, 1.0), INF),
+            ((10.0, INF), INF),
+        ];
+        for &(arg, res) in test.iter() {
+            test_case(arg, res, f);
+        }
     }
 
     #[test]
     fn test_pdf() {
-        test_case(1.0, 0.1, 0.090483741803595961836995913651194571475319347018875963, |x| x.pdf(1.0));
-        test_case(1.0, 0.1, 0.036787944117144234201693506390001264039984687455876246, |x| x.pdf(10.0));
-        test_case(1.0, 1.0, 0.36787944117144232159552377016146086744581113103176804, |x| x.pdf(1.0));
-        test_case(1.0, 1.0, 0.000045399929762484851535591515560550610237918088866564953, |x| x.pdf(10.0));
-        test_almost(10.0, 10.0, 1.2511003572113329898476497894772544708420990097708588, 1e-14, |x| x.pdf(1.0));
-        test_almost(10.0, 10.0, 1.0251532120868705806216092933926141802686541811003037e-30, 1e-44, |x| x.pdf(10.0));
-        test_almost(10.0, 1.0, 0.0000010137771196302974029859010421116095333052555418644397, 1e-20, |x| x.pdf(1.0));
-        test_almost(10.0, 1.0, 0.12511003572113329898476497894772544708420990097708601, 1e-15, |x| x.pdf(10.0));
-        test_is_nan(10.0, f64::INFINITY, |x| x.pdf(1.0)); // is this really the behavior we want?
-        test_case(10.0, f64::INFINITY, 0.0, |x| x.pdf(f64::INFINITY));
+        let f = |arg: f64| move |x: Gamma| x.pdf(arg);
+        let test = [
+            ((1.0, 0.1), 1.0, 0.090483741803595961836995),
+            ((1.0, 0.1), 10.0, 0.036787944117144234201693),
+            ((1.0, 1.0), 1.0, 0.367879441171442321595523),
+            ((1.0, 1.0), 10.0, 0.000045399929762484851535),
+            ((10.0, 10.0), 1.0, 1.251100357211332989847649),
+            ((10.0, 10.0), 10.0, 1.025153212086870580621609e-30),
+            ((10.0, 1.0), 1.0, 0.000001013777119630297402),
+            ((10.0, 1.0), 10.0, 0.125110035721133298984764),
+        ];
+        for &(arg, x, res) in test.iter() {
+            test_case(arg, res, f(x));
+        }
+        //TODO: test special
+        // test_is_nan((10.0, INF), pdf(1.0)); // is this really the behavior we want?
+        //TODO: test special
+        // (10.0, INF, INF, 0.0, pdf(INF)),];
     }
 
     #[test]
     fn test_pdf_at_zero() {
-        test_almost(1.0, 0.1, 0.1, 1e-10, |x| x.pdf(0.0));
-        test_almost(1.0, 0.1, 0.1f64.ln(), 1e-10, |x| x.ln_pdf(0.0));
+        test_case((1.0, 0.1), 0.1, |x| x.pdf(0.0));
+        test_case((1.0, 0.1), 0.1f64.ln(), |x| x.ln_pdf(0.0));
     }
 
     #[test]
     fn test_ln_pdf() {
-        test_case(1.0, 0.1, -2.402585092994045634057955346552321429281631934330484, |x| x.ln_pdf(1.0));
-        test_case(1.0, 0.1, -3.3025850929940456285068402234265387271634735938763824, |x| x.ln_pdf(10.0));
-        test_case(1.0, 1.0, -1.0, |x| x.ln_pdf(1.0));
-        test_case(1.0, 1.0, -10.0, |x| x.ln_pdf(10.0));
-        test_almost(10.0, 10.0, 0.22402344985898722897219667227693591172986563062456522, 1e-15, |x| x.ln_pdf(1.0));
-        test_case(10.0, 10.0, -69.052710713194601614865880235563786219860220971716511, |x| x.ln_pdf(10.0));
-        test_almost(10.0, 1.0, -13.801827480081469611207717874566706164281149255663166, 1e-14, |x| x.ln_pdf(1.0));
-        test_almost(10.0, 1.0,  -2.0785616431350584550457947824074282958712358580042068, 1e-14, |x| x.ln_pdf(10.0));
-        test_is_nan(10.0, f64::INFINITY, |x| x.ln_pdf(1.0)); // is this really the behavior we want?
-        test_case(10.0, f64::INFINITY, f64::NEG_INFINITY, |x| x.ln_pdf(f64::INFINITY));
+        let f = |arg: f64| move |x: Gamma| x.ln_pdf(arg);
+        let test = [
+            ((1.0, 0.1), 1.0, -2.40258509299404563405795),
+            ((1.0, 0.1), 10.0, -3.30258509299404562850684),
+            ((1.0, 1.0), 1.0, -1.0),
+            ((1.0, 1.0), 10.0, -10.0),
+            ((10.0, 10.0), 1.0, 0.224023449858987228972196),
+            ((10.0, 10.0), 10.0, -69.0527107131946016148658),
+            ((10.0, 1.0), 1.0, -13.8018274800814696112077),
+            ((10.0, 1.0), 10.0, -2.07856164313505845504579),
+            ((10.0, INF), INF, f64::NEG_INFINITY),
+        ];
+        for &(arg, x, res) in test.iter() {
+            test_case(arg, res, f(x));
+        }
+        // TODO: test special
+        // test_is_nan((10.0, INF), f(1.0)); // is this really the behavior we want?
     }
 
     #[test]
     fn test_cdf() {
-        test_almost(1.0, 0.1, 0.095162581964040431858607615783064404690935346242622848, 1e-16, |x| x.cdf(1.0));
-        test_almost(1.0, 0.1, 0.63212055882855767840447622983853913255418886896823196, 1e-15, |x| x.cdf(10.0));
-        test_almost(1.0, 1.0, 0.63212055882855767840447622983853913255418886896823196, 1e-15, |x| x.cdf(1.0));
-        test_case(1.0, 1.0, 0.99995460007023751514846440848443944938976208191113396,|x| x.cdf(10.0));
-        test_almost(10.0, 10.0, 0.54207028552814779168583514294066541824736464003242184, 1e-15, |x| x.cdf(1.0));
-        test_case(10.0, 10.0, 0.99999999999999999999999999999988746526039157266114706, |x| x.cdf(10.0));
-        test_almost(10.0, 1.0, 0.00000011142547833872067735305068724025236288094949815466035, 1e-21, |x| x.cdf(1.0));
-        test_almost(10.0, 1.0, 0.54207028552814779168583514294066541824736464003242184, 1e-15, |x| x.cdf(10.0));
-        test_case(10.0, f64::INFINITY, 0.0, |x| x.cdf(1.0));
-        test_case(10.0, f64::INFINITY, 1.0, |x| x.cdf(10.0));
+        let f = |arg: f64| move |x: Gamma| x.cdf(arg);
+        let test = [
+            ((1.0, 0.1), 1.0, 0.095162581964040431858607),
+            ((1.0, 0.1), 10.0, 0.632120558828557678404476),
+            ((1.0, 1.0), 1.0, 0.632120558828557678404476),
+            ((1.0, 1.0), 10.0, 0.999954600070237515148464),
+            ((10.0, 10.0), 1.0, 0.542070285528147791685835),
+            ((10.0, 10.0), 10.0, 0.999999999999999999999999),
+            ((10.0, 1.0), 1.0, 0.000000111425478338720677),
+            ((10.0, 1.0), 10.0, 0.542070285528147791685835),
+            ((10.0, INF), 1.0, 0.0),
+            ((10.0, INF), 10.0, 1.0),
+        ];
+        for &(arg, x, res) in test.iter() {
+            test_case(arg, res, f(x));
+        }
     }
 
     #[test]
     fn test_cdf_at_zero() {
-        test_case(1.0, 0.1, 0.0, |x| x.cdf(0.0));
+        test_case((1.0, 0.1), 0.0, |x| x.cdf(0.0));
     }
 
     #[test]
     fn test_continuous() {
-        tests::check_continuous_distribution(&try_create(1.0, 0.5), 0.0, 20.0);
-        tests::check_continuous_distribution(&try_create(9.0, 2.0), 0.0, 20.0);
+        test::check_continuous_distribution(&try_create((1.0, 0.5)), 0.0, 20.0);
+        test::check_continuous_distribution(&try_create((9.0, 2.0)), 0.0, 20.0);
     }
 }
diff --git a/src/distribution/geometric.rs b/src/distribution/geometric.rs
index d2195e59..70321cdc 100644
--- a/src/distribution/geometric.rs
+++ b/src/distribution/geometric.rs
@@ -412,8 +412,8 @@ mod tests {
 
     #[test]
     fn test_discrete() {
-        tests::check_discrete_distribution(&try_create(0.3), 100);
-        tests::check_discrete_distribution(&try_create(0.6), 100);
-        tests::check_discrete_distribution(&try_create(1.0), 1);
+        test::check_discrete_distribution(&try_create(0.3), 100);
+        test::check_discrete_distribution(&try_create(0.6), 100);
+        test::check_discrete_distribution(&try_create(1.0), 1);
     }
 }
diff --git a/src/distribution/hypergeometric.rs b/src/distribution/hypergeometric.rs
index e15ee86c..a11469fc 100644
--- a/src/distribution/hypergeometric.rs
+++ b/src/distribution/hypergeometric.rs
@@ -527,7 +527,7 @@ mod tests {
 
     #[test]
     fn test_discrete() {
-        tests::check_discrete_distribution(&try_create(5, 4, 3), 4);
-        tests::check_discrete_distribution(&try_create(3, 2, 1), 2);
+        test::check_discrete_distribution(&try_create(5, 4, 3), 4);
+        test::check_discrete_distribution(&try_create(3, 2, 1), 2);
     }
 }
diff --git a/src/distribution/internal.rs b/src/distribution/internal.rs
index 86cb6710..70401449 100644
--- a/src/distribution/internal.rs
+++ b/src/distribution/internal.rs
@@ -13,11 +13,62 @@ pub fn is_valid_multinomial(arr: &[f64], incl_zero: bool) -> bool {
 }
 
 #[cfg(test)]
-pub mod tests {
+pub mod test {
     use super::is_valid_multinomial;
     use crate::consts::ACC;
     use crate::distribution::{Continuous, ContinuousCDF, Discrete, DiscreteCDF};
 
+    #[macro_export]
+    macro_rules! testing_boiler {
+        ($arg:ty, $dist:ty) => {
+            fn try_create(arg: $arg) -> $dist {
+                let n = <$dist>::new.call_once(arg);
+                assert!(n.is_ok());
+                n.unwrap()
+            }
+
+            fn bad_create_case(arg: $arg) {
+                let n = <$dist>::new.call(arg);
+                assert!(n.is_err());
+            }
+
+            fn get_value<F, T>(arg: $arg, eval: F) -> T
+            where
+                F: Fn($dist) -> T,
+            {
+                let n = try_create(arg);
+                eval(n)
+            }
+
+            fn test_case<F, T>(arg: $arg, expected: T, eval: F)
+            where
+                F: Fn($dist) -> T,
+                T: ::core::fmt::Debug + ::approx::RelativeEq<Epsilon = f64>,
+            {
+                let x = get_value(arg, eval);
+                assert_relative_eq!(expected, x, max_relative = ACC);
+            }
+
+            fn test_case_special<F, T>(arg: $arg, expected: T, acc: f64, eval: F)
+            where
+                F: Fn($dist) -> T,
+                T: ::core::fmt::Debug + ::approx::AbsDiffEq<Epsilon = f64>,
+            {
+                let x = get_value(arg, eval);
+                assert_abs_diff_eq!(expected, x, epsilon = acc);
+            }
+
+            fn test_none<F, T>(arg: $arg, eval: F)
+            where
+                F: Fn($dist) -> Option<T>,
+                T: ::core::cmp::PartialEq + ::core::fmt::Debug,
+            {
+                let x = get_value(arg, eval);
+                assert_eq!(None, x);
+            }
+        };
+    }
+
     /// cdf should be the integral of the pdf
     fn check_integrate_pdf_is_cdf<D: ContinuousCDF<f64, f64> + Continuous<f64, f64>>(
         dist: &D,
diff --git a/src/distribution/inverse_gamma.rs b/src/distribution/inverse_gamma.rs
index ae9eff0c..222a70dd 100644
--- a/src/distribution/inverse_gamma.rs
+++ b/src/distribution/inverse_gamma.rs
@@ -442,7 +442,7 @@ mod tests {
 
     #[test]
     fn test_continuous() {
-        tests::check_continuous_distribution(&try_create(1.0, 0.5), 0.0, 100.0);
-        tests::check_continuous_distribution(&try_create(9.0, 2.0), 0.0, 100.0);
+        test::check_continuous_distribution(&try_create(1.0, 0.5), 0.0, 100.0);
+        test::check_continuous_distribution(&try_create(9.0, 2.0), 0.0, 100.0);
     }
 }
diff --git a/src/distribution/laplace.rs b/src/distribution/laplace.rs
new file mode 100644
index 00000000..fc670450
--- /dev/null
+++ b/src/distribution/laplace.rs
@@ -0,0 +1,467 @@
+use crate::distribution::{Continuous, ContinuousCDF};
+use crate::statistics::*;
+use crate::{Result, StatsError};
+use rand::Rng;
+use std::f64;
+
+/// Implements the [Laplace](https://en.wikipedia.org/wiki/Laplace_distribution)
+/// distribution.
+///
+/// # Examples
+///
+/// ```
+/// use statrs::distribution::{Laplace, Continuous};
+/// use statrs::statistics::Mode;
+///
+/// let n = Laplace::new(0.0, 1.0).unwrap();
+/// assert_eq!(n.mode().unwrap(), 0.0);
+/// assert_eq!(n.pdf(1.0), 0.18393972058572117);
+/// ```
+#[derive(Debug, Copy, Clone, PartialEq)]
+pub struct Laplace {
+    location: f64,
+    scale: f64,
+}
+
+impl Laplace {
+    /// Constructs a new laplace distribution with the given
+    /// location and scale.
+    ///
+    /// # Errors
+    ///
+    /// Returns an error if location or scale are `NaN` or `scale <= 0.0`
+    ///
+    /// # Examples
+    ///
+    /// ```
+    /// use statrs::distribution::Laplace;
+    ///
+    /// let mut result = Laplace::new(0.0, 1.0);
+    /// assert!(result.is_ok());
+    ///
+    /// result = Laplace::new(0.0, -1.0);
+    /// assert!(result.is_err());
+    /// ```
+    pub fn new(location: f64, scale: f64) -> Result<Laplace> {
+        if location.is_nan() || scale.is_nan() || scale <= 0.0 {
+            Err(StatsError::BadParams)
+        } else {
+            Ok(Laplace { location, scale })
+        }
+    }
+
+    /// Returns the location of the laplace distribution
+    ///
+    /// # Examples
+    ///
+    /// ```
+    /// use statrs::distribution::Laplace;
+    ///
+    /// let n = Laplace::new(0.0, 1.0).unwrap();
+    /// assert_eq!(n.location(), 0.0);
+    /// ```
+    pub fn location(&self) -> f64 {
+        self.location
+    }
+
+    /// Returns the scale of the laplace distribution
+    ///
+    /// # Examples
+    ///
+    /// ```
+    /// use statrs::distribution::Laplace;
+    ///
+    /// let n = Laplace::new(0.0, 1.0).unwrap();
+    /// assert_eq!(n.scale(), 1.0);
+    /// ```
+    pub fn scale(&self) -> f64 {
+        self.scale
+    }
+}
+
+impl ::rand::distributions::Distribution<f64> for Laplace {
+    fn sample<R: Rng + ?Sized>(&self, rng: &mut R) -> f64 {
+        let x: f64 = rng.gen_range(-0.5, 0.5);
+        self.location - self.scale * x.signum() * (1. - 2. * x).ln()
+    }
+}
+
+impl ContinuousCDF<f64, f64> for Laplace {
+    /// Calculates the cumulative distribution function for the
+    /// laplace distribution at `x`
+    ///
+    /// # Formula
+    ///
+    /// ```ignore
+    /// (1 / 2) * (1 + signum(x - μ)) - signum(x - μ) * exp(-|x - μ| / b)
+    /// ```
+    ///
+    /// where `μ` is the location, `b` is the scale
+    fn cdf(&self, x: f64) -> f64 {
+        let y = (-(x - self.location).abs() / self.scale).exp() / 2.;
+        if x >= self.location {
+            1. - y
+        } else {
+            y
+        }
+    }
+    /// Calculates the inverse cumulative distribution function for the
+    /// laplace distribution at `p`
+    ///
+    /// # Formula
+    ///
+    /// if p <= 1/2
+    /// ```ignore
+    /// μ + b * ln(2p)
+    /// ```
+    /// if p >= 1/2
+    /// ```ignore
+    /// μ - b * ln(2 - 2p)
+    /// ```
+    ///
+    /// where `μ` is the location, `b` is the scale
+    fn inverse_cdf(&self, p: f64) -> f64 {
+        if p <= 0. || 1. <= p {
+            panic!("p must be in [0, 1]");
+        };
+        if p <= 0.5 {
+            self.location + self.scale * (2. * p)
+        } else {
+            self.location - self.scale * (2. - 2. * p)
+        }
+    }
+}
+
+impl Min<f64> for Laplace {
+    /// Returns the minimum value in the domain of the laplace
+    /// distribution representable by a double precision float
+    ///
+    /// # Formula
+    ///
+    /// ```ignore
+    /// NEG_INF
+    /// ```
+    fn min(&self) -> f64 {
+        f64::NEG_INFINITY
+    }
+}
+
+impl Max<f64> for Laplace {
+    /// Returns the maximum value in the domain of the laplace
+    /// distribution representable by a double precision float
+    ///
+    /// # Formula
+    ///
+    /// ```ignore
+    /// INF
+    /// ```
+    fn max(&self) -> f64 {
+        f64::INFINITY
+    }
+}
+
+impl Distribution<f64> for Laplace {
+    /// Returns the mode of the laplace distribution
+    ///
+    /// # Formula
+    ///
+    /// ```ignore
+    /// μ
+    /// ```
+    ///
+    /// where `μ` is the location
+    fn mean(&self) -> Option<f64> {
+        Some(self.location)
+    }
+    /// Returns the variance of the laplace distribution
+    ///
+    /// # Formula
+    ///
+    /// ```ignore
+    /// 2*b^2
+    /// ```
+    ///
+    /// where `b` is the scale
+    fn variance(&self) -> Option<f64> {
+        Some(2. * self.scale * self.scale)
+    }
+    /// Returns the entropy of the laplace distribution
+    ///
+    /// # Formula
+    ///
+    /// ```ignore
+    /// ln(2be)
+    /// ```
+    ///
+    /// where `b` is the scale
+    fn entropy(&self) -> Option<f64> {
+        Some((2. * self.scale).ln() + 1.)
+    }
+    /// Returns the skewness of the laplace distribution
+    ///
+    /// # Formula
+    ///
+    /// ```ignore
+    /// 0
+    /// ```
+    fn skewness(&self) -> Option<f64> {
+        Some(0.)
+    }
+}
+
+impl Median<f64> for Laplace {
+    /// Returns the median of the laplace distribution
+    ///
+    /// # Formula
+    ///
+    /// ```ignore
+    /// μ
+    /// ```
+    ///
+    /// where `μ` is the location
+    fn median(&self) -> f64 {
+        self.location
+    }
+}
+
+impl Mode<Option<f64>> for Laplace {
+    /// Returns the mode of the laplace distribution
+    ///
+    /// # Formula
+    ///
+    /// ```ignore
+    /// μ
+    /// ```
+    ///
+    /// where `μ` is the location
+    fn mode(&self) -> Option<f64> {
+        Some(self.location)
+    }
+}
+
+impl Continuous<f64, f64> for Laplace {
+    /// Calculates the probability density function for the laplace
+    /// distribution at `x`
+    ///
+    /// # Formula
+    ///
+    /// ```ignore
+    /// (1 / 2b) * exp(-|x - μ| / b)
+    /// ```
+    /// where `μ` is the location and `b` is the scale
+    fn pdf(&self, x: f64) -> f64 {
+        (-(x - self.location).abs() / self.scale).exp() / (2. * self.scale)
+    }
+
+    /// Calculates the log probability density function for the laplace
+    /// distribution at `x`
+    ///
+    /// # Formula
+    ///
+    /// ```ignore
+    /// ln((1 / 2b) * exp(-|x - μ| / b))
+    /// ```
+    ///
+    /// where `μ` is the location and `b` is the scale
+    fn ln_pdf(&self, x: f64) -> f64 {
+        ((-(x - self.location).abs() / self.scale).exp() / (2. * self.scale)).ln()
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use core::f64::INFINITY as INF;
+    use rand::thread_rng;
+
+    fn try_create(location: f64, scale: f64) -> Laplace {
+        let n = Laplace::new(location, scale);
+        assert!(n.is_ok());
+        n.unwrap()
+    }
+
+    fn create_case(location: f64, scale: f64) {
+        let n = try_create(location, scale);
+        assert_eq!(location, n.location);
+        assert_eq!(scale, n.scale);
+    }
+
+    fn bad_create_case(location: f64, scale: f64) {
+        let n = Laplace::new(location, scale);
+        assert!(n.is_err());
+    }
+
+    fn test_case<F>(location: f64, scale: f64, expected: f64, eval: F)
+    where
+        F: Fn(Laplace) -> f64,
+    {
+        let n = try_create(location, scale);
+        let x = eval(n);
+        assert_eq!(expected, x);
+    }
+
+    fn test_is_nan<F>(location: f64, scale: f64, eval: F)
+    where
+        F: Fn(Laplace) -> f64,
+    {
+        let n = try_create(location, scale);
+        let x = eval(n);
+        assert!(x.is_nan());
+    }
+
+    fn test_almost<F>(location: f64, scale: f64, expected: f64, acc: f64, eval: F)
+    where
+        F: Fn(Laplace) -> f64,
+    {
+        let n = try_create(location, scale);
+        let x = eval(n);
+        assert_almost_eq!(expected, x, acc);
+    }
+
+    #[test]
+    fn test_create() {
+        try_create(1.0, 2.0);
+        try_create(-INF, 0.1);
+        try_create(-5.0 - 1.0, 1.0);
+        try_create(0.0, 5.0);
+        try_create(1.0, 7.0);
+        try_create(5.0, 10.0);
+        try_create(INF, INF);
+    }
+
+    #[test]
+    fn test_bad_create() {
+        bad_create_case(2.0, -1.0);
+        bad_create_case(f64::NAN, 1.0);
+        bad_create_case(f64::NAN, -1.0);
+    }
+
+    #[test]
+    fn test_mean() {
+        let mean = |x: Laplace| x.mean().unwrap();
+        test_case(-INF, 0.1, -INF, mean);
+        test_case(-5.0 - 1.0, 1.0, -6.0, mean);
+        test_case(0.0, 5.0, 0.0, mean);
+        test_case(1.0, 10.0, 1.0, mean);
+        test_case(INF, INF, INF, mean);
+    }
+
+    #[test]
+    fn test_variance() {
+        let variance = |x: Laplace| x.variance().unwrap();
+        test_almost(-INF, 0.1, 0.02, 1E-12, variance);
+        test_almost(-5.0 - 1.0, 1.0, 2.0, 1E-12, variance);
+        test_almost(0.0, 5.0, 50.0, 1E-12, variance);
+        test_almost(1.0, 7.0, 98.0, 1E-12, variance);
+        test_almost(5.0, 10.0, 200.0, 1E-12, variance);
+        test_almost(INF, INF, INF, 1E-12, variance);
+    }
+    #[test]
+    fn test_entropy() {
+        let entropy = |x: Laplace| x.entropy().unwrap();
+        test_almost(-INF, 0.1, (2.0 * f64::consts::E * 0.1).ln(), 1E-12, entropy);
+        test_almost(-6.0, 1.0, (2.0 * f64::consts::E).ln(), 1E-12, entropy);
+        test_almost(1.0, 7.0, (2.0 * f64::consts::E * 7.0).ln(), 1E-12, entropy);
+        test_almost(5., 10., (2. * f64::consts::E * 10.).ln(), 1E-12, entropy);
+        test_almost(INF, INF, INF, 1E-12, entropy);
+    }
+
+    #[test]
+    fn test_skewness() {
+        let skewness = |x: Laplace| x.skewness().unwrap();
+        test_case(-INF, 0.1, 0.0, skewness);
+        test_case(-6.0, 1.0, 0.0, skewness);
+        test_case(1.0, 7.0, 0.0, skewness);
+        test_case(5.0, 10.0, 0.0, skewness);
+        test_case(INF, INF, 0.0, skewness);
+    }
+
+    #[test]
+    fn test_mode() {
+        let mode = |x: Laplace| x.mode().unwrap();
+        test_case(-INF, 0.1, -INF, mode);
+        test_case(-6.0, 1.0, -6.0, mode);
+        test_case(1.0, 7.0, 1.0, mode);
+        test_case(5.0, 10.0, 5.0, mode);
+        test_case(INF, INF, INF, mode);
+    }
+
+    #[test]
+    fn test_median() {
+        let median = |x: Laplace| x.median();
+        test_case(-INF, 0.1, -INF, median);
+        test_case(-6.0, 1.0, -6.0, median);
+        test_case(1.0, 7.0, 1.0, median);
+        test_case(5.0, 10.0, 5.0, median);
+        test_case(INF, INF, INF, median);
+    }
+
+    #[test]
+    fn test_min() {
+        test_case(0.0, 1.0, -INF, |l| l.min());
+    }
+
+    #[test]
+    fn test_max() {
+        test_case(0.0, 1.0, INF, |l| l.max());
+    }
+
+    #[test]
+    fn test_density() {
+        let pdf = |arg: f64| move |x: Laplace| x.pdf(arg);
+        test_almost(0.0, 0.1, 1.529511602509129e-06, 1E-12, pdf(1.5));
+        test_almost(1.0, 0.1, 7.614989872356341e-08, 1E-12, pdf(2.8));
+        test_almost(-1.0, 0.1, 3.8905661205668983e-19, 1E-12, pdf(-5.4));
+        test_almost(5.0, 0.1, 5.056107463052243e-43, 1E-12, pdf(-4.9));
+        test_almost(-5.0, 0.1, 1.9877248679543235e-30, 1E-12, pdf(2.0));
+        test_almost(INF, 0.1, 0.0, 1E-12, pdf(5.5));
+        test_almost(-INF, 0.1, 0.0, 1E-12, pdf(-0.0));
+        test_almost(0.0, 1.0, 0.0, 1E-12, pdf(INF));
+        test_almost(1.0, 1.0, 0.00915781944436709, 1E-12, pdf(5.0));
+        test_almost(-1.0, 1.0, 0.5, 1E-12, pdf(-1.0));
+        test_almost(5.0, 1.0, 0.0012393760883331792, 1E-12, pdf(-1.0));
+        test_almost(-5.0, 1.0, 0.0002765421850739168, 1E-12, pdf(2.5));
+        test_almost(INF, 0.1, 0.0, 1E-12, pdf(2.0));
+        test_almost(-INF, 0.1, 0.0, 1E-12, pdf(15.0));
+        test_almost(0.0, INF, 0.0, 1E-12, pdf(89.3));
+        test_almost(1.0, INF, 0.0, 1E-12, pdf(-0.1));
+        test_almost(-1.0, INF, 0.0, 1E-12, pdf(0.1));
+        test_almost(5.0, INF, 0.0, 1E-12, pdf(-6.1));
+        test_almost(-5.0, INF, 0.0, 1E-12, pdf(-10.0));
+        test_is_nan(INF, INF, pdf(2.0));
+        test_is_nan(-INF, INF, pdf(-5.1));
+    }
+
+    #[test]
+    fn test_ln_density() {
+        let ln_pdf = |arg: f64| move |x: Laplace| x.ln_pdf(arg);
+        test_almost(0.0, 0.1, -13.3905620875659, 1E-12, ln_pdf(1.5));
+        test_almost(1.0, 0.1, -16.390562087565897, 1E-12, ln_pdf(2.8));
+        test_almost(-1.0, 0.1, -42.39056208756591, 1E-12, ln_pdf(-5.4));
+        test_almost(5.0, 0.1, -97.3905620875659, 1E-12, ln_pdf(-4.9));
+        test_almost(-5.0, 0.1, -68.3905620875659, 1E-12, ln_pdf(2.0));
+        test_case(INF, 0.1, -INF, ln_pdf(5.5));
+        test_case(-INF, 0.1, -INF, ln_pdf(-0.0));
+        test_case(0.0, 1.0, -INF, ln_pdf(INF));
+        test_almost(1.0, 1.0, -4.693147180559945, 1E-12, ln_pdf(5.0));
+        test_almost(-1.0, 1.0, -f64::consts::LN_2, 1E-12, ln_pdf(-1.0));
+        test_almost(5.0, 1.0, -6.693147180559945, 1E-12, ln_pdf(-1.0));
+        test_almost(-5.0, 1.0, -8.193147180559945, 1E-12, ln_pdf(2.5));
+        test_case(INF, 0.1, -INF, ln_pdf(2.0));
+        test_case(-INF, 0.1, -INF, ln_pdf(15.0));
+        test_case(0.0, INF, -INF, ln_pdf(89.3));
+        test_case(1.0, INF, -INF, ln_pdf(-0.1));
+        test_case(-1.0, INF, -INF, ln_pdf(0.1));
+        test_case(5.0, INF, -INF, ln_pdf(-6.1));
+        test_case(-5.0, INF, -INF, ln_pdf(-10.0));
+        test_is_nan(INF, INF, ln_pdf(2.0));
+        test_is_nan(-INF, INF, ln_pdf(-5.1));
+    }
+
+    #[test]
+    fn test_sample() {
+        use ::rand::distributions::Distribution;
+        let l = try_create(0.1, 0.5);
+        l.sample(&mut thread_rng());
+    }
+}
diff --git a/src/distribution/log_normal.rs b/src/distribution/log_normal.rs
index 8926a83b..5f84b494 100644
--- a/src/distribution/log_normal.rs
+++ b/src/distribution/log_normal.rs
@@ -598,7 +598,7 @@ mod tests {
 
     #[test]
     fn test_continuous() {
-        tests::check_continuous_distribution(&try_create(0.0, 0.25), 0.0, 10.0);
-        tests::check_continuous_distribution(&try_create(0.0, 0.5), 0.0, 10.0);
+        test::check_continuous_distribution(&try_create(0.0, 0.25), 0.0, 10.0);
+        test::check_continuous_distribution(&try_create(0.0, 0.5), 0.0, 10.0);
     }
 }
diff --git a/src/distribution/mod.rs b/src/distribution/mod.rs
index 3943fb0f..ce924d33 100644
--- a/src/distribution/mod.rs
+++ b/src/distribution/mod.rs
@@ -22,6 +22,7 @@ pub use self::gamma::Gamma;
 pub use self::geometric::Geometric;
 pub use self::hypergeometric::Hypergeometric;
 pub use self::inverse_gamma::InverseGamma;
+pub use self::laplace::Laplace;
 pub use self::log_normal::LogNormal;
 pub use self::multinomial::Multinomial;
 pub use self::multivariate_normal::MultivariateNormal;
@@ -53,6 +54,7 @@ mod geometric;
 mod hypergeometric;
 mod internal;
 mod inverse_gamma;
+mod laplace;
 mod log_normal;
 mod multinomial;
 mod multivariate_normal;
diff --git a/src/distribution/negative_binomial.rs b/src/distribution/negative_binomial.rs
index 4c1f189e..494a7068 100644
--- a/src/distribution/negative_binomial.rs
+++ b/src/distribution/negative_binomial.rs
@@ -430,7 +430,7 @@ mod tests {
     //
     // #[test]
     // fn test_discrete() {
-    //     tests::check_discrete_distribution(&try_create(5.0, 0.3), 35);
-    //     tests::check_discrete_distribution(&try_create(10.0, 0.7), 21);
+    //     test::check_discrete_distribution(&try_create(5.0, 0.3), 35);
+    //     test::check_discrete_distribution(&try_create(10.0, 0.7), 21);
     // }
 }
diff --git a/src/distribution/normal.rs b/src/distribution/normal.rs
index 015831c3..2c8c7de5 100644
--- a/src/distribution/normal.rs
+++ b/src/distribution/normal.rs
@@ -446,8 +446,8 @@ mod tests {
 
     #[test]
     fn test_continuous() {
-        tests::check_continuous_distribution(&try_create(0.0, 1.0), -10.0, 10.0);
-        tests::check_continuous_distribution(&try_create(20.0, 0.5), 10.0, 30.0);
+        test::check_continuous_distribution(&try_create(0.0, 1.0), -10.0, 10.0);
+        test::check_continuous_distribution(&try_create(20.0, 0.5), 10.0, 30.0);
     }
 
     #[test]
diff --git a/src/distribution/pareto.rs b/src/distribution/pareto.rs
index 0298c2db..ff089e4e 100644
--- a/src/distribution/pareto.rs
+++ b/src/distribution/pareto.rs
@@ -489,7 +489,7 @@ mod tests {
 
     #[test]
     fn test_continuous() {
-        tests::check_continuous_distribution(&try_create(1.0, 10.0), 1.0, 10.0);
-        tests::check_continuous_distribution(&try_create(0.1, 2.0), 0.1, 100.0);
+        test::check_continuous_distribution(&try_create(1.0, 10.0), 1.0, 10.0);
+        test::check_continuous_distribution(&try_create(0.1, 2.0), 0.1, 100.0);
     }
 }
diff --git a/src/distribution/poisson.rs b/src/distribution/poisson.rs
index 37e82cdd..10ae475d 100644
--- a/src/distribution/poisson.rs
+++ b/src/distribution/poisson.rs
@@ -445,7 +445,7 @@ mod tests {
 
     #[test]
     fn test_discrete() {
-        tests::check_discrete_distribution(&try_create(0.3), 10);
-        tests::check_discrete_distribution(&try_create(4.5), 30);
+        test::check_discrete_distribution(&try_create(0.3), 10);
+        test::check_discrete_distribution(&try_create(4.5), 30);
     }
 }
diff --git a/src/distribution/students_t.rs b/src/distribution/students_t.rs
index d4fd0949..0c2d13e1 100644
--- a/src/distribution/students_t.rs
+++ b/src/distribution/students_t.rs
@@ -152,19 +152,23 @@ impl ContinuousCDF<f64, f64> for StudentsT {
     }
 
     /// Calculates the inverse cumulative distribution function for the
-    /// student's t-distribution at `x`
+    /// Student's T-distribution at `x`
     fn inverse_cdf(&self, x: f64) -> f64 {
+        // first calculate inverse_cdf for normal Student's T
         assert!(x >= 0.0 && x <= 1.0);
         let x = 2. * x.min(1. - x);
         let a = 0.5 * self.freedom;
         let b = 0.5;
-        let y = beta::inv_beta_reg(a, b, x);
-        let y = (self.freedom * (1. - y) / y).sqrt();
-        if x <= 0.5 {
-            y
-        } else {
-            -y
-        }
+        let mut y = beta::inv_beta_reg(a, b, x);
+        y = (self.freedom * (1. - y) / y).sqrt();
+        y = if x <= 0.5 { y } else { -y };
+        // generalised Student's T is related to normal Student's T by `Y = μ + σ X`
+        // where `X` is distributed as Student's T, so this result has to be scaled and shifted back
+        // formally: F_Y(t) = P(Y <= t) = P(X <= (t - μ) / σ) = F_X((t - μ) / σ)
+        // F_Y^{-1}(p) = inf { t' | F_Y(t') >= p } = inf { t' = μ + σ t | F_X((t' - μ) / σ) >= p }
+        // because scale is positive: loc + scale * t is strictly monotonic function
+        // = μ + σ inf { t | F_X(t) >= p } = μ + σ F_X^{-1}(p)
+        self.location + self.scale * y
     }
 }
 
@@ -606,9 +610,455 @@ mod tests {
 
     #[test]
     fn test_continuous() {
-        tests::check_continuous_distribution(&try_create(0.0, 1.0, 3.0), -30.0, 30.0);
-        tests::check_continuous_distribution(&try_create(0.0, 1.0, 10.0), -10.0, 10.0);
-        tests::check_continuous_distribution(&try_create(20.0, 0.5, 10.0), 10.0, 30.0);
+        test::check_continuous_distribution(&try_create(0.0, 1.0, 3.0), -30.0, 30.0);
+        test::check_continuous_distribution(&try_create(0.0, 1.0, 10.0), -10.0, 10.0);
+        test::check_continuous_distribution(&try_create(20.0, 0.5, 10.0), 10.0, 30.0);
+    }
+
+    #[test]
+    fn test_inv_cdf() {
+        let test = |x: f64, freedom: f64, expected: f64| {
+            use approx::*;
+            let d = StudentsT::new(0., 1., freedom).unwrap();
+            // Checks that left == right to 4 significant figures, unlike
+            // test_almost() which uses decimal places
+            assert_relative_eq!(d.inverse_cdf(x), expected, max_relative = 0.001);
+        };
+
+        // This test checks our implementation against the whole t-table
+        // copied from https://en.wikipedia.org/wiki/Student's_t-distribution
+
+        test(0.75, 1.0, 1.000);
+        test(0.8, 1.0, 1.376);
+        test(0.85, 1.0, 1.963);
+        test(0.9, 1.0, 3.078);
+        test(0.95, 1.0, 6.314);
+        test(0.975, 1.0, 12.71);
+        test(0.99, 1.0, 31.82);
+        test(0.995, 1.0, 63.66);
+        test(0.9975, 1.0, 127.3);
+        test(0.999, 1.0, 318.3);
+        test(0.9995, 1.0, 636.6);
+
+        test(0.75, 002.0, 0.816);
+        // test(0.8, 002.0, 1.080);  // We get 1.061 for some reason...
+        test(0.85, 002.0, 1.386);
+        test(0.9, 002.0, 1.886);
+        test(0.95, 002.0, 2.920);
+        test(0.975, 002.0, 4.303);
+        test(0.99, 002.0, 6.965);
+        test(0.995, 002.0, 9.925);
+        test(0.9975, 002.0, 14.09);
+        test(0.999, 002.0, 22.33);
+        test(0.9995, 002.0, 31.60);
+
+        test(0.75, 003.0, 0.765);
+        test(0.8, 003.0, 0.978);
+        test(0.85, 003.0, 1.250);
+        test(0.9, 003.0, 1.638);
+        test(0.95, 003.0, 2.353);
+        test(0.975, 003.0, 3.182);
+        test(0.99, 003.0, 4.541);
+        test(0.995, 003.0, 5.841);
+        test(0.9975, 003.0, 7.453);
+        test(0.999, 003.0, 10.21);
+        test(0.9995, 003.0, 12.92);
+
+        test(0.75, 004.0, 0.741);
+        test(0.8, 004.0, 0.941);
+        test(0.85, 004.0, 1.190);
+        test(0.9, 004.0, 1.533);
+        test(0.95, 004.0, 2.132);
+        test(0.975, 004.0, 2.776);
+        test(0.99, 004.0, 3.747);
+        test(0.995, 004.0, 4.604);
+        test(0.9975, 004.0, 5.598);
+        test(0.999, 004.0, 7.173);
+        test(0.9995, 004.0, 8.610);
+
+        test(0.75, 005.0, 0.727);
+        test(0.8, 005.0, 0.920);
+        test(0.85, 005.0, 1.156);
+        test(0.9, 005.0, 1.476);
+        test(0.95, 005.0, 2.015);
+        test(0.975, 005.0, 2.571);
+        test(0.99, 005.0, 3.365);
+        test(0.995, 005.0, 4.032);
+        test(0.9975, 005.0, 4.773);
+        test(0.999, 005.0, 5.893);
+        test(0.9995, 005.0, 6.869);
+
+        test(0.75, 006.0, 0.718);
+        test(0.8, 006.0, 0.906);
+        test(0.85, 006.0, 1.134);
+        test(0.9, 006.0, 1.440);
+        test(0.95, 006.0, 1.943);
+        test(0.975, 006.0, 2.447);
+        test(0.99, 006.0, 3.143);
+        test(0.995, 006.0, 3.707);
+        test(0.9975, 006.0, 4.317);
+        test(0.999, 006.0, 5.208);
+        test(0.9995, 006.0, 5.959);
+
+        test(0.75, 007.0, 0.711);
+        test(0.8, 007.0, 0.896);
+        test(0.85, 007.0, 1.119);
+        test(0.9, 007.0, 1.415);
+        test(0.95, 007.0, 1.895);
+        test(0.975, 007.0, 2.365);
+        test(0.99, 007.0, 2.998);
+        test(0.995, 007.0, 3.499);
+        test(0.9975, 007.0, 4.029);
+        test(0.999, 007.0, 4.785);
+        test(0.9995, 007.0, 5.408);
+
+        test(0.75, 008.0, 0.706);
+        test(0.8, 008.0, 0.889);
+        test(0.85, 008.0, 1.108);
+        test(0.9, 008.0, 1.397);
+        test(0.95, 008.0, 1.860);
+        test(0.975, 008.0, 2.306);
+        test(0.99, 008.0, 2.896);
+        test(0.995, 008.0, 3.355);
+        test(0.9975, 008.0, 3.833);
+        test(0.999, 008.0, 4.501);
+        test(0.9995, 008.0, 5.041);
+
+        test(0.75, 009.0, 0.703);
+        test(0.8, 009.0, 0.883);
+        test(0.85, 009.0, 1.100);
+        test(0.9, 009.0, 1.383);
+        test(0.95, 009.0, 1.833);
+        test(0.975, 009.0, 2.262);
+        test(0.99, 009.0, 2.821);
+        test(0.995, 009.0, 3.250);
+        test(0.9975, 009.0, 3.690);
+        test(0.999, 009.0, 4.297);
+        test(0.9995, 009.0, 4.781);
+
+        test(0.75, 010.0, 0.700);
+        test(0.8, 010.0, 0.879);
+        test(0.85, 010.0, 1.093);
+        test(0.9, 010.0, 1.372);
+        test(0.95, 010.0, 1.812);
+        test(0.975, 010.0, 2.228);
+        test(0.99, 010.0, 2.764);
+        test(0.995, 010.0, 3.169);
+        test(0.9975, 010.0, 3.581);
+        test(0.999, 010.0, 4.144);
+        test(0.9995, 010.0, 4.587);
+
+        test(0.75, 011.0, 0.697);
+        test(0.8, 011.0, 0.876);
+        test(0.85, 011.0, 1.088);
+        test(0.9, 011.0, 1.363);
+        test(0.95, 011.0, 1.796);
+        test(0.975, 011.0, 2.201);
+        test(0.99, 011.0, 2.718);
+        test(0.995, 011.0, 3.106);
+        test(0.9975, 011.0, 3.497);
+        test(0.999, 011.0, 4.025);
+        test(0.9995, 011.0, 4.437);
+
+        test(0.75, 012.0, 0.695);
+        test(0.8, 012.0, 0.873);
+        test(0.85, 012.0, 1.083);
+        test(0.9, 012.0, 1.356);
+        test(0.95, 012.0, 1.782);
+        test(0.975, 012.0, 2.179);
+        test(0.99, 012.0, 2.681);
+        test(0.995, 012.0, 3.055);
+        test(0.9975, 012.0, 3.428);
+        test(0.999, 012.0, 3.930);
+        test(0.9995, 012.0, 4.318);
+
+        test(0.75, 013.0, 0.694);
+        test(0.8, 013.0, 0.870);
+        test(0.85, 013.0, 1.079);
+        test(0.9, 013.0, 1.350);
+        test(0.95, 013.0, 1.771);
+        test(0.975, 013.0, 2.160);
+        test(0.99, 013.0, 2.650);
+        test(0.995, 013.0, 3.012);
+        test(0.9975, 013.0, 3.372);
+        test(0.999, 013.0, 3.852);
+        test(0.9995, 013.0, 4.221);
+
+        test(0.75, 014.0, 0.692);
+        test(0.8, 014.0, 0.868);
+        test(0.85, 014.0, 1.076);
+        test(0.9, 014.0, 1.345);
+        test(0.95, 014.0, 1.761);
+        test(0.975, 014.0, 2.145);
+        test(0.99, 014.0, 2.624);
+        test(0.995, 014.0, 2.977);
+        test(0.9975, 014.0, 3.326);
+        test(0.999, 014.0, 3.787);
+        test(0.9995, 014.0, 4.140);
+
+        test(0.75, 015.0, 0.691);
+        test(0.8, 015.0, 0.866);
+        test(0.85, 015.0, 1.074);
+        test(0.9, 015.0, 1.341);
+        test(0.95, 015.0, 1.753);
+        test(0.975, 015.0, 2.131);
+        test(0.99, 015.0, 2.602);
+        test(0.995, 015.0, 2.947);
+        test(0.9975, 015.0, 3.286);
+        test(0.999, 015.0, 3.733);
+        test(0.9995, 015.0, 4.073);
+
+        test(0.75, 016.0, 0.690);
+        test(0.8, 016.0, 0.865);
+        test(0.85, 016.0, 1.071);
+        test(0.9, 016.0, 1.337);
+        test(0.95, 016.0, 1.746);
+        test(0.975, 016.0, 2.120);
+        test(0.99, 016.0, 2.583);
+        test(0.995, 016.0, 2.921);
+        test(0.9975, 016.0, 3.252);
+        test(0.999, 016.0, 3.686);
+        test(0.9995, 016.0, 4.015);
+
+        test(0.75, 017.0, 0.689);
+        test(0.8, 017.0, 0.863);
+        test(0.85, 017.0, 1.069);
+        test(0.9, 017.0, 1.333);
+        test(0.95, 017.0, 1.740);
+        test(0.975, 017.0, 2.110);
+        test(0.99, 017.0, 2.567);
+        test(0.995, 017.0, 2.898);
+        test(0.9975, 017.0, 3.222);
+        test(0.999, 017.0, 3.646);
+        test(0.9995, 017.0, 3.965);
+
+        test(0.75, 018.0, 0.688);
+        test(0.8, 018.0, 0.862);
+        test(0.85, 018.0, 1.067);
+        test(0.9, 018.0, 1.330);
+        test(0.95, 018.0, 1.734);
+        test(0.975, 018.0, 2.101);
+        test(0.99, 018.0, 2.552);
+        test(0.995, 018.0, 2.878);
+        test(0.9975, 018.0, 3.197);
+        test(0.999, 018.0, 3.610);
+        test(0.9995, 018.0, 3.922);
+
+        test(0.75, 019.0, 0.688);
+        test(0.8, 019.0, 0.861);
+        test(0.85, 019.0, 1.066);
+        test(0.9, 019.0, 1.328);
+        test(0.95, 019.0, 1.729);
+        test(0.975, 019.0, 2.093);
+        test(0.99, 019.0, 2.539);
+        test(0.995, 019.0, 2.861);
+        test(0.9975, 019.0, 3.174);
+        test(0.999, 019.0, 3.579);
+        test(0.9995, 019.0, 3.883);
+
+        test(0.75, 020.0, 0.687);
+        test(0.8, 020.0, 0.860);
+        test(0.85, 020.0, 1.064);
+        test(0.9, 020.0, 1.325);
+        test(0.95, 020.0, 1.725);
+        test(0.975, 020.0, 2.086);
+        test(0.99, 020.0, 2.528);
+        test(0.995, 020.0, 2.845);
+        test(0.9975, 020.0, 3.153);
+        test(0.999, 020.0, 3.552);
+        test(0.9995, 020.0, 3.850);
+
+        test(0.75, 021.0, 0.686);
+        test(0.8, 021.0, 0.859);
+        test(0.85, 021.0, 1.063);
+        test(0.9, 021.0, 1.323);
+        test(0.95, 021.0, 1.721);
+        test(0.975, 021.0, 2.080);
+        test(0.99, 021.0, 2.518);
+        test(0.995, 021.0, 2.831);
+        test(0.9975, 021.0, 3.135);
+        test(0.999, 021.0, 3.527);
+        test(0.9995, 021.0, 3.819);
+
+        test(0.75, 022.0, 0.686);
+        test(0.8, 022.0, 0.858);
+        test(0.85, 022.0, 1.061);
+        test(0.9, 022.0, 1.321);
+        test(0.95, 022.0, 1.717);
+        test(0.975, 022.0, 2.074);
+        test(0.99, 022.0, 2.508);
+        test(0.995, 022.0, 2.819);
+        test(0.9975, 022.0, 3.119);
+        test(0.999, 022.0, 3.505);
+        test(0.9995, 022.0, 3.792);
+
+        test(0.75, 023.0, 0.685);
+        test(0.8, 023.0, 0.858);
+        test(0.85, 023.0, 1.060);
+        test(0.9, 023.0, 1.319);
+        test(0.95, 023.0, 1.714);
+        test(0.975, 023.0, 2.069);
+        test(0.99, 023.0, 2.500);
+        test(0.995, 023.0, 2.807);
+        test(0.9975, 023.0, 3.104);
+        test(0.999, 023.0, 3.485);
+        test(0.9995, 023.0, 3.767);
+
+        test(0.75, 024.0, 0.685);
+        test(0.8, 024.0, 0.857);
+        test(0.85, 024.0, 1.059);
+        test(0.9, 024.0, 1.318);
+        test(0.95, 024.0, 1.711);
+        test(0.975, 024.0, 2.064);
+        test(0.99, 024.0, 2.492);
+        test(0.995, 024.0, 2.797);
+        test(0.9975, 024.0, 3.091);
+        test(0.999, 024.0, 3.467);
+        test(0.9995, 024.0, 3.745);
+
+        test(0.75, 025.0, 0.684);
+        test(0.8, 025.0, 0.856);
+        test(0.85, 025.0, 1.058);
+        test(0.9, 025.0, 1.316);
+        test(0.95, 025.0, 1.708);
+        test(0.975, 025.0, 2.060);
+        test(0.99, 025.0, 2.485);
+        test(0.995, 025.0, 2.787);
+        test(0.9975, 025.0, 3.078);
+        test(0.999, 025.0, 3.450);
+        test(0.9995, 025.0, 3.725);
+
+        test(0.75, 026.0, 0.684);
+        test(0.8, 026.0, 0.856);
+        test(0.85, 026.0, 1.058);
+        test(0.9, 026.0, 1.315);
+        test(0.95, 026.0, 1.706);
+        test(0.975, 026.0, 2.056);
+        test(0.99, 026.0, 2.479);
+        test(0.995, 026.0, 2.779);
+        test(0.9975, 026.0, 3.067);
+        test(0.999, 026.0, 3.435);
+        test(0.9995, 026.0, 3.707);
+
+        test(0.75, 027.0, 0.684);
+        test(0.8, 027.0, 0.855);
+        test(0.85, 027.0, 1.057);
+        test(0.9, 027.0, 1.314);
+        test(0.95, 027.0, 1.703);
+        test(0.975, 027.0, 2.052);
+        test(0.99, 027.0, 2.473);
+        test(0.995, 027.0, 2.771);
+        test(0.9975, 027.0, 3.057);
+        test(0.999, 027.0, 3.421);
+        test(0.9995, 027.0, 3.690);
+
+        test(0.75, 028.0, 0.683);
+        test(0.8, 028.0, 0.855);
+        test(0.85, 028.0, 1.056);
+        test(0.9, 028.0, 1.313);
+        test(0.95, 028.0, 1.701);
+        test(0.975, 028.0, 2.048);
+        test(0.99, 028.0, 2.467);
+        test(0.995, 028.0, 2.763);
+        test(0.9975, 028.0, 3.047);
+        test(0.999, 028.0, 3.408);
+        test(0.9995, 028.0, 3.674);
+
+        test(0.75, 029.0, 0.683);
+        test(0.8, 029.0, 0.854);
+        test(0.85, 029.0, 1.055);
+        test(0.9, 029.0, 1.311);
+        test(0.95, 029.0, 1.699);
+        test(0.975, 029.0, 2.045);
+        test(0.99, 029.0, 2.462);
+        test(0.995, 029.0, 2.756);
+        test(0.9975, 029.0, 3.038);
+        test(0.999, 029.0, 3.396);
+        test(0.9995, 029.0, 3.659);
+
+        test(0.75, 030.0, 0.683);
+        test(0.8, 030.0, 0.854);
+        test(0.85, 030.0, 1.055);
+        test(0.9, 030.0, 1.310);
+        test(0.95, 030.0, 1.697);
+        test(0.975, 030.0, 2.042);
+        test(0.99, 030.0, 2.457);
+        test(0.995, 030.0, 2.750);
+        test(0.9975, 030.0, 3.030);
+        test(0.999, 030.0, 3.385);
+        test(0.9995, 030.0, 3.646);
+
+        test(0.75, 040.0, 0.681);
+        test(0.8, 040.0, 0.851);
+        test(0.85, 040.0, 1.050);
+        test(0.9, 040.0, 1.303);
+        test(0.95, 040.0, 1.684);
+        test(0.975, 040.0, 2.021);
+        test(0.99, 040.0, 2.423);
+        test(0.995, 040.0, 2.704);
+        test(0.9975, 040.0, 2.971);
+        test(0.999, 040.0, 3.307);
+        test(0.9995, 040.0, 3.551);
+
+        test(0.75, 050.0, 0.679);
+        test(0.8, 050.0, 0.849);
+        test(0.85, 050.0, 1.047);
+        test(0.9, 050.0, 1.299);
+        test(0.95, 050.0, 1.676);
+        test(0.975, 050.0, 2.009);
+        test(0.99, 050.0, 2.403);
+        test(0.995, 050.0, 2.678);
+        test(0.9975, 050.0, 2.937);
+        test(0.999, 050.0, 3.261);
+        test(0.9995, 050.0, 3.496);
+
+        test(0.75, 060.0, 0.679);
+        test(0.8, 060.0, 0.848);
+        test(0.85, 060.0, 1.045);
+        test(0.9, 060.0, 1.296);
+        test(0.95, 060.0, 1.671);
+        test(0.975, 060.0, 2.000);
+        test(0.99, 060.0, 2.390);
+        test(0.995, 060.0, 2.660);
+        test(0.9975, 060.0, 2.915);
+        test(0.999, 060.0, 3.232);
+        test(0.9995, 060.0, 3.460);
+
+        test(0.75, 080.0, 0.678);
+        test(0.8, 080.0, 0.846);
+        test(0.85, 080.0, 1.043);
+        test(0.9, 080.0, 1.292);
+        test(0.95, 080.0, 1.664);
+        test(0.975, 080.0, 1.990);
+        test(0.99, 080.0, 2.374);
+        test(0.995, 080.0, 2.639);
+        test(0.9975, 080.0, 2.887);
+        test(0.999, 080.0, 3.195);
+        test(0.9995, 080.0, 3.416);
+
+        test(0.75, 100.0, 0.677);
+        test(0.8, 100.0, 0.845);
+        test(0.85, 100.0, 1.042);
+        test(0.9, 100.0, 1.290);
+        test(0.95, 100.0, 1.660);
+        test(0.975, 100.0, 1.984);
+        test(0.99, 100.0, 2.364);
+        test(0.995, 100.0, 2.626);
+        test(0.9975, 100.0, 2.871);
+        test(0.999, 100.0, 3.174);
+        test(0.9995, 100.0, 3.390);
+
+        test(0.75, 120.0, 0.677);
+        test(0.8, 120.0, 0.845);
+        test(0.85, 120.0, 1.041);
+        test(0.9, 120.0, 1.289);
+        test(0.95, 120.0, 1.658);
+        test(0.975, 120.0, 1.980);
+        test(0.99, 120.0, 2.358);
+        test(0.995, 120.0, 2.617);
+        test(0.9975, 120.0, 2.860);
+        test(0.999, 120.0, 3.160);
+        test(0.9995, 120.0, 3.373);
     }
 
     #[test]
diff --git a/src/distribution/triangular.rs b/src/distribution/triangular.rs
index 6149c176..99067cf1 100644
--- a/src/distribution/triangular.rs
+++ b/src/distribution/triangular.rs
@@ -474,7 +474,7 @@ mod tests {
 
     #[test]
     fn test_continuous() {
-        tests::check_continuous_distribution(&try_create(-5.0, 5.0, 0.0), -5.0, 5.0);
-        tests::check_continuous_distribution(&try_create(-15.0, -2.0, -3.0), -15.0, -2.0);
+        test::check_continuous_distribution(&try_create(-5.0, 5.0, 0.0), -5.0, 5.0);
+        test::check_continuous_distribution(&try_create(-15.0, -2.0, -3.0), -15.0, -2.0);
     }
 }
diff --git a/src/distribution/uniform.rs b/src/distribution/uniform.rs
index 9e604662..c909a45a 100644
--- a/src/distribution/uniform.rs
+++ b/src/distribution/uniform.rs
@@ -409,8 +409,8 @@ mod tests {
 
     #[test]
     fn test_continuous() {
-        tests::check_continuous_distribution(&try_create(0.0, 10.0), 0.0, 10.0);
-        tests::check_continuous_distribution(&try_create(-2.0, 15.0), -2.0, 15.0);
+        test::check_continuous_distribution(&try_create(0.0, 10.0), 0.0, 10.0);
+        test::check_continuous_distribution(&try_create(-2.0, 15.0), -2.0, 15.0);
     }
 
     #[test]
diff --git a/src/distribution/weibull.rs b/src/distribution/weibull.rs
index eb617ae1..1399233b 100644
--- a/src/distribution/weibull.rs
+++ b/src/distribution/weibull.rs
@@ -485,6 +485,6 @@ mod tests {
 
     #[test]
     fn test_continuous() {
-        tests::check_continuous_distribution(&try_create(1.0, 0.2), 0.0, 10.0);
+        test::check_continuous_distribution(&try_create(1.0, 0.2), 0.0, 10.0);
     }
 }
diff --git a/src/function/factorial.rs b/src/function/factorial.rs
index 07bc88d2..eac59bcf 100644
--- a/src/function/factorial.rs
+++ b/src/function/factorial.rs
@@ -4,13 +4,12 @@
 use crate::error::StatsError;
 use crate::function::gamma;
 use crate::Result;
-use std::f64;
-use std::sync::Once;
+use core::f64::INFINITY as INF;
 
 /// The maximum factorial representable
 /// by a 64-bit floating point without
 /// overflowing
-pub const MAX_ARG: u64 = 170;
+pub const MAX_FACTORIAL: usize = 170;
 
 /// Computes the factorial function `x -> x!` for
 /// `170 >= x >= 0`. All factorials larger than `170!`
@@ -20,11 +19,8 @@ pub const MAX_ARG: u64 = 170;
 ///
 /// Returns `f64::INFINITY` if `x > 170`
 pub fn factorial(x: u64) -> f64 {
-    if x > MAX_ARG {
-        f64::INFINITY
-    } else {
-        get_fcache()[x as usize]
-    }
+    let x = x as usize;
+    FCACHE.get(x).map_or(INF, |&fac| fac)
 }
 
 /// Computes the logarithmic factorial function `x -> ln(x!)`
@@ -34,13 +30,10 @@ pub fn factorial(x: u64) -> f64 {
 ///
 /// Returns `0.0` if `x <= 1`
 pub fn ln_factorial(x: u64) -> f64 {
-    if x <= 1 {
-        0.0
-    } else if x > MAX_ARG {
-        gamma::ln_gamma(x as f64 + 1.0)
-    } else {
-        get_fcache()[x as usize].ln()
-    }
+    let x = x as usize;
+    FCACHE
+        .get(x)
+        .map_or_else(|| gamma::ln_gamma(x as f64 + 1.0), |&fac| fac.ln())
 }
 
 /// Computes the binomial coefficient `n choose k`
@@ -98,89 +91,87 @@ pub fn checked_multinomial(n: u64, ni: &[u64]) -> Result<f64> {
 
 // Initialization for pre-computed cache of 171 factorial
 // values 0!...170!
-const CACHE_SIZE: usize = 171;
-
-static mut FCACHE: [f64; CACHE_SIZE] = [1.0; CACHE_SIZE];
-static START: Once = Once::new();
-
-fn get_fcache() -> [f64; CACHE_SIZE] {
-    START.call_once(|| {
-        (1..CACHE_SIZE).fold(1.0, |acc, i| {
-            let fac = acc * i as f64;
-            unsafe {
-                FCACHE[i] = fac;
-            }
-            fac
-        });
-    });
-    unsafe { FCACHE }
+lazy_static! {
+    static ref FCACHE: [f64; MAX_FACTORIAL + 1] = {
+        let mut fcache = [1.0; MAX_FACTORIAL + 1];
+        fcache
+            .iter_mut()
+            .enumerate()
+            .skip(1)
+            .fold(1.0, |acc, (i, elt)| {
+                let fac = acc * i as f64;
+                *elt = fac;
+                fac
+            });
+        fcache
+    };
 }
 
-#[rustfmt::skip]
 #[cfg(test)]
 mod tests {
-    use std::{f64, u64};
+    use super::*;
 
     #[test]
     fn test_factorial_and_ln_factorial() {
-        let mut factorial = 1.0;
+        let mut fac = 1.0;
+        assert_eq!(factorial(0), fac);
         for i in 1..171 {
-            factorial *= i as f64;
-            assert_eq!(super::factorial(i), factorial);
-            assert_eq!(super::ln_factorial(i), factorial.ln());
+            fac *= i as f64;
+            assert_eq!(factorial(i), fac);
+            assert_eq!(ln_factorial(i), fac.ln());
         }
     }
 
     #[test]
     fn test_factorial_overflow() {
-        assert_eq!(super::factorial(172), f64::INFINITY);
-        assert_eq!(super::factorial(u64::MAX), f64::INFINITY);
+        assert_eq!(factorial(172), INF);
+        assert_eq!(factorial(u64::MAX), INF);
     }
 
     #[test]
     fn test_ln_factorial_does_not_overflow() {
-        assert_eq!(super::ln_factorial(1 << 10), 6078.2118847500501140);
-        assert_almost_eq!(super::ln_factorial(1 << 12), 29978.648060844048236, 1e-11);
-        assert_eq!(super::ln_factorial(1 << 15), 307933.81973375485425);
-        assert_eq!(super::ln_factorial(1 << 17), 1413421.9939462073242);
+        assert_eq!(ln_factorial(1 << 10), 6078.2118847500501140);
+        assert_almost_eq!(ln_factorial(1 << 12), 29978.648060844048236, 1e-11);
+        assert_eq!(ln_factorial(1 << 15), 307933.81973375485425);
+        assert_eq!(ln_factorial(1 << 17), 1413421.9939462073242);
     }
 
     #[test]
     fn test_binomial() {
-        assert_eq!(super::binomial(1, 1), 1.0);
-        assert_eq!(super::binomial(5, 2), 10.0);
-        assert_eq!(super::binomial(7, 3), 35.0);
-        assert_eq!(super::binomial(1, 0), 1.0);
-        assert_eq!(super::binomial(0, 1), 0.0);
-        assert_eq!(super::binomial(5, 7), 0.0);
+        assert_eq!(binomial(1, 1), 1.0);
+        assert_eq!(binomial(5, 2), 10.0);
+        assert_eq!(binomial(7, 3), 35.0);
+        assert_eq!(binomial(1, 0), 1.0);
+        assert_eq!(binomial(0, 1), 0.0);
+        assert_eq!(binomial(5, 7), 0.0);
     }
 
     #[test]
     fn test_ln_binomial() {
-        assert_eq!(super::ln_binomial(1, 1), 1f64.ln());
-        assert_almost_eq!(super::ln_binomial(5, 2), 10f64.ln(), 1e-14);
-        assert_almost_eq!(super::ln_binomial(7, 3), 35f64.ln(), 1e-14);
-        assert_eq!(super::ln_binomial(1, 0), 1f64.ln());
-        assert_eq!(super::ln_binomial(0, 1), 0f64.ln());
-        assert_eq!(super::ln_binomial(5, 7), 0f64.ln());
+        assert_eq!(ln_binomial(1, 1), 1f64.ln());
+        assert_almost_eq!(ln_binomial(5, 2), 10f64.ln(), 1e-14);
+        assert_almost_eq!(ln_binomial(7, 3), 35f64.ln(), 1e-14);
+        assert_eq!(ln_binomial(1, 0), 1f64.ln());
+        assert_eq!(ln_binomial(0, 1), 0f64.ln());
+        assert_eq!(ln_binomial(5, 7), 0f64.ln());
     }
 
     #[test]
     fn test_multinomial() {
-        assert_eq!(1.0, super::multinomial(1, &[1, 0]));
-        assert_eq!(10.0, super::multinomial(5, &[3, 2]));
-        assert_eq!(10.0, super::multinomial(5, &[2, 3]));
-        assert_eq!(35.0, super::multinomial(7, &[3, 4]));
+        assert_eq!(1.0, multinomial(1, &[1, 0]));
+        assert_eq!(10.0, multinomial(5, &[3, 2]));
+        assert_eq!(10.0, multinomial(5, &[2, 3]));
+        assert_eq!(35.0, multinomial(7, &[3, 4]));
     }
 
     #[test]
     #[should_panic]
     fn test_multinomial_bad_ni() {
-        super::multinomial(1, &[1, 1]);
+        multinomial(1, &[1, 1]);
     }
 
     #[test]
     fn test_checked_multinomial_bad_ni() {
-        assert!(super::checked_multinomial(1, &[1, 1]).is_err());
+        assert!(checked_multinomial(1, &[1, 1]).is_err());
     }
 }
diff --git a/src/lib.rs b/src/lib.rs
index 83d2eee7..7481b662 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -27,12 +27,17 @@
 #![crate_name = "statrs"]
 #![allow(clippy::excessive_precision)]
 #![allow(clippy::many_single_char_names)]
-#![allow(dead_code)]
 #![allow(unused_imports)]
+#![forbid(unsafe_code)]
+#![cfg_attr(test, feature(unboxed_closures))]
+#![cfg_attr(test, feature(fn_traits))]
 
 #[macro_use]
 extern crate approx;
 
+#[macro_use]
+extern crate lazy_static;
+
 #[macro_export]
 macro_rules! assert_almost_eq {
     ($a:expr, $b:expr, $prec:expr) => {
diff --git a/src/statistics/mod.rs b/src/statistics/mod.rs
index 0b11a7f3..b156c6ec 100644
--- a/src/statistics/mod.rs
+++ b/src/statistics/mod.rs
@@ -2,12 +2,13 @@
 
 pub use self::iter_statistics::*;
 pub use self::order_statistics::*;
+pub use self::slice_statistics::*;
 pub use self::statistics::*;
 pub use self::traits::*;
 
 mod iter_statistics;
 mod order_statistics;
 // TODO: fix later
-// mod slice_statistics;
+mod slice_statistics;
 mod statistics;
 mod traits;
diff --git a/src/statistics/order_statistics.rs b/src/statistics/order_statistics.rs
index 9c231eb2..1655e681 100644
--- a/src/statistics/order_statistics.rs
+++ b/src/statistics/order_statistics.rs
@@ -17,15 +17,18 @@ pub trait OrderStatistics<T> {
     ///
     /// ```
     /// use statrs::statistics::OrderStatistics;
+    /// use statrs::statistics::Data;
     ///
-    /// let mut x = [];
+    /// let x = [];
+    /// let mut x = Data::new(x);
     /// assert!(x.order_statistic(1).is_nan());
     ///
-    /// let mut y = [0.0, 3.0, -2.0];
+    /// let y = [0.0, 3.0, -2.0];
+    /// let mut y = Data::new(y);
     /// assert!(y.order_statistic(0).is_nan());
     /// assert!(y.order_statistic(4).is_nan());
     /// assert_eq!(y.order_statistic(2), 0.0);
-    /// assert!(y != [0.0, 3.0, -2.0]);
+    /// assert!(y != Data::new([0.0, 3.0, -2.0]));
     /// ```
     fn order_statistic(&mut self, order: usize) -> T;
 
@@ -39,13 +42,16 @@ pub trait OrderStatistics<T> {
     ///
     /// ```
     /// use statrs::statistics::OrderStatistics;
+    /// use statrs::statistics::Data;
     ///
-    /// let mut x = [];
+    /// let x = [];
+    /// let mut x = Data::new(x);
     /// assert!(x.median().is_nan());
     ///
-    /// let mut y = [0.0, 3.0, -2.0];
+    /// let y = [0.0, 3.0, -2.0];
+    /// let mut y = Data::new(y);
     /// assert_eq!(y.median(), 0.0);
-    /// assert!(y != [0.0, 3.0, -2.0]);
+    /// assert!(y != Data::new([0.0, 3.0, -2.0]));
     fn median(&mut self) -> T;
 
     /// Estimates the tau-th quantile from the data. The tau-th quantile
@@ -62,15 +68,18 @@ pub trait OrderStatistics<T> {
     ///
     /// ```
     /// use statrs::statistics::OrderStatistics;
+    /// use statrs::statistics::Data;
     ///
-    /// let mut x = [];
+    /// let x = [];
+    /// let mut x = Data::new(x);
     /// assert!(x.quantile(0.5).is_nan());
     ///
-    /// let mut y = [0.0, 3.0, -2.0];
+    /// let y = [0.0, 3.0, -2.0];
+    /// let mut y = Data::new(y);
     /// assert!(y.quantile(-1.0).is_nan());
     /// assert!(y.quantile(2.0).is_nan());
     /// assert_eq!(y.quantile(0.5), 0.0);
-    /// assert!(y != [0.0, 3.0, -2.0]);
+    /// assert!(y != Data::new([0.0, 3.0, -2.0]));
     /// ```
     fn quantile(&mut self, tau: f64) -> T;
 
@@ -87,16 +96,19 @@ pub trait OrderStatistics<T> {
     ///
     /// ```
     /// use statrs::statistics::OrderStatistics;
+    /// use statrs::statistics::Data;
     ///
-    /// let mut x = [];
+    /// let x = [];
+    /// let mut x = Data::new(x);
     /// assert!(x.percentile(0).is_nan());
     ///
-    /// let mut y = [1.0, 5.0, 3.0, 4.0, 10.0, 9.0, 6.0, 7.0, 8.0, 2.0];
+    /// let y = [1.0, 5.0, 3.0, 4.0, 10.0, 9.0, 6.0, 7.0, 8.0, 2.0];
+    /// let mut y = Data::new(y);
     /// assert_eq!(y.percentile(0), 1.0);
     /// assert_eq!(y.percentile(50), 5.5);
     /// assert_eq!(y.percentile(100), 10.0);
     /// assert!(y.percentile(105).is_nan());
-    /// assert!(y != [1.0, 5.0, 3.0, 4.0, 10.0, 9.0, 6.0, 7.0, 8.0, 2.0]);
+    /// assert!(y != Data::new([1.0, 5.0, 3.0, 4.0, 10.0, 9.0, 6.0, 7.0, 8.0, 2.0]));
     /// ```
     fn percentile(&mut self, p: usize) -> T;
 
@@ -113,14 +125,17 @@ pub trait OrderStatistics<T> {
     /// extern crate statrs;
     ///
     /// use statrs::statistics::OrderStatistics;
+    /// use statrs::statistics::Data;
     ///
     /// # fn main() {
-    /// let mut x = [];
+    /// let x = [];
+    /// let mut x = Data::new(x);
     /// assert!(x.lower_quartile().is_nan());
     ///
-    /// let mut y = [2.0, 1.0, 3.0, 4.0];
+    /// let y = [2.0, 1.0, 3.0, 4.0];
+    /// let mut y = Data::new(y);
     /// assert_almost_eq!(y.lower_quartile(), 1.416666666666666, 1e-15);
-    /// assert!(y != [2.0, 1.0, 3.0, 4.0]);
+    /// assert!(y != Data::new([2.0, 1.0, 3.0, 4.0]));
     /// # }
     /// ```
     fn lower_quartile(&mut self) -> T;
@@ -138,14 +153,17 @@ pub trait OrderStatistics<T> {
     /// extern crate statrs;
     ///
     /// use statrs::statistics::OrderStatistics;
+    /// use statrs::statistics::Data;
     ///
     /// # fn main() {
-    /// let mut x = [];
+    /// let x = [];
+    /// let mut x = Data::new(x);
     /// assert!(x.upper_quartile().is_nan());
     ///
-    /// let mut y = [2.0, 1.0, 3.0, 4.0];
+    /// let y = [2.0, 1.0, 3.0, 4.0];
+    /// let mut y = Data::new(y);
     /// assert_almost_eq!(y.upper_quartile(), 3.5833333333333333, 1e-15);
-    /// assert!(y != [2.0, 1.0, 3.0, 4.0]);
+    /// assert!(y != Data::new([2.0, 1.0, 3.0, 4.0]));
     /// # }
     /// ```
     fn upper_quartile(&mut self) -> T;
@@ -162,15 +180,18 @@ pub trait OrderStatistics<T> {
     /// #[macro_use]
     /// extern crate statrs;
     ///
+    /// use statrs::statistics::Data;
     /// use statrs::statistics::OrderStatistics;
     ///
     /// # fn main() {
-    /// let mut x = [];
+    /// let x = [];
+    /// let mut x = Data::new(x);
     /// assert!(x.interquartile_range().is_nan());
     ///
-    /// let mut y = [2.0, 1.0, 3.0, 4.0];
+    /// let y = [2.0, 1.0, 3.0, 4.0];
+    /// let mut y = Data::new(y);
     /// assert_almost_eq!(y.interquartile_range(), 2.166666666666667, 1e-15);
-    /// assert!(y != [2.0, 1.0, 3.0, 4.0]);
+    /// assert!(y != Data::new([2.0, 1.0, 3.0, 4.0]));
     /// # }
     /// ```
     fn interquartile_range(&mut self) -> T;
@@ -181,14 +202,17 @@ pub trait OrderStatistics<T> {
     ///
     /// ```
     /// use statrs::statistics::{OrderStatistics, RankTieBreaker};
+    /// use statrs::statistics::Data;
     ///
-    /// let mut x = [];
+    /// let x = [];
+    /// let mut x = Data::new(x);
     /// assert_eq!(x.ranks(RankTieBreaker::Average).len(), 0);
     ///
     /// let y = [1.0, 3.0, 2.0, 2.0];
-    /// assert_eq!((&mut y.clone()).ranks(RankTieBreaker::Average), [1.0, 4.0,
+    /// let mut y = Data::new([1.0, 3.0, 2.0, 2.0]);
+    /// assert_eq!(y.clone().ranks(RankTieBreaker::Average), [1.0, 4.0,
     /// 2.5, 2.5]);
-    /// assert_eq!((&mut y.clone()).ranks(RankTieBreaker::Min), [1.0, 4.0, 2.0,
+    /// assert_eq!(y.clone().ranks(RankTieBreaker::Min), [1.0, 4.0, 2.0,
     /// 2.0]);
     /// ```
     fn ranks(&mut self, tie_breaker: RankTieBreaker) -> Vec<T>;
diff --git a/src/statistics/slice_statistics.rs b/src/statistics/slice_statistics.rs
index d39cea30..abb6cae7 100644
--- a/src/statistics/slice_statistics.rs
+++ b/src/statistics/slice_statistics.rs
@@ -1,23 +1,130 @@
 use crate::statistics::*;
-use std::f64;
+use core::ops::{Index, IndexMut};
+use rand::prelude::SliceRandom;
 
-impl OrderStatistics<f64> for [f64] {
+#[derive(Clone, Debug, PartialEq, Eq)]
+pub struct Data<D>(D);
+
+impl<D: AsRef<[f64]>> Index<usize> for Data<D> {
+    type Output = f64;
+    fn index(&self, i: usize) -> &f64 {
+        &self.0.as_ref()[i]
+    }
+}
+
+impl<D: AsMut<[f64]> + AsRef<[f64]>> IndexMut<usize> for Data<D> {
+    fn index_mut(&mut self, i: usize) -> &mut f64 {
+        &mut self.0.as_mut()[i]
+    }
+}
+
+impl<D: AsMut<[f64]> + AsRef<[f64]>> Data<D> {
+    pub fn new(data: D) -> Self {
+        Data(data)
+    }
+    pub fn swap(&mut self, i: usize, j: usize) {
+        self.0.as_mut().swap(i, j)
+    }
+    pub fn len(&self) -> usize {
+        self.0.as_ref().len()
+    }
+    pub fn is_empty(&self) -> bool {
+        self.0.as_ref().len() == 0
+    }
+    pub fn iter(&self) -> core::slice::Iter<'_, f64> {
+        self.0.as_ref().iter()
+    }
+    // Selection algorithm from Numerical Recipes
+    // See: https://en.wikipedia.org/wiki/Selection_algorithm
+    fn select_inplace(&mut self, rank: usize) -> f64 {
+        if rank == 0 {
+            return self.min();
+        }
+        if rank > self.len() - 1 {
+            return self.max();
+        }
+
+        let mut low = 0;
+        let mut high = self.len() - 1;
+        loop {
+            if high <= low + 1 {
+                if high == low + 1 && self[high] < self[low] {
+                    self.swap(low, high)
+                }
+                return self[rank];
+            }
+
+            let middle = (low + high) / 2;
+            self.swap(middle, low + 1);
+
+            if self[low] > self[high] {
+                self.swap(low, high);
+            }
+            if self[low + 1] > self[high] {
+                self.swap(low + 1, high);
+            }
+            if self[low] > self[low + 1] {
+                self.swap(low, low + 1);
+            }
+
+            let mut begin = low + 1;
+            let mut end = high;
+            let pivot = self[begin];
+            loop {
+                loop {
+                    begin += 1;
+                    if self[begin] >= pivot {
+                        break;
+                    }
+                }
+                loop {
+                    end -= 1;
+                    if self[end] <= pivot {
+                        break;
+                    }
+                }
+                if end < begin {
+                    break;
+                }
+                self.swap(begin, end);
+            }
+
+            self[low + 1] = self[end];
+            self[end] = pivot;
+
+            if end >= rank {
+                high = end - 1;
+            }
+            if end <= rank {
+                low = begin;
+            }
+        }
+    }
+}
+
+impl<D: AsRef<[f64]>> ::rand::distributions::Distribution<f64> for Data<D> {
+    fn sample<R: ::rand::Rng + ?Sized>(&self, rng: &mut R) -> f64 {
+        *self.0.as_ref().choose(rng).unwrap()
+    }
+}
+
+impl<D: AsMut<[f64]> + AsRef<[f64]>> OrderStatistics<f64> for Data<D> {
     fn order_statistic(&mut self, order: usize) -> f64 {
         let n = self.len();
         match order {
             1 => self.min(),
             _ if order == n => self.max(),
             _ if order < 1 || order > n => f64::NAN,
-            _ => select_inplace(self, order - 1),
+            _ => self.select_inplace(order - 1),
         }
     }
 
     fn median(&mut self) -> f64 {
         let k = self.len() / 2;
         if self.len() % 2 != 0 {
-            select_inplace(self, k)
+            self.select_inplace(k)
         } else {
-            (select_inplace(self, k.saturating_sub(1)) + select_inplace(self, k)) / 2.0
+            (self.select_inplace(k.saturating_sub(1)) + self.select_inplace(k)) / 2.0
         }
     }
 
@@ -36,8 +143,8 @@ impl OrderStatistics<f64> for [f64] {
             return self.max();
         }
 
-        let a = select_inplace(self, (hf as usize).saturating_sub(1));
-        let b = select_inplace(self, hf as usize);
+        let a = self.select_inplace((hf as usize).saturating_sub(1));
+        let b = self.select_inplace(hf as usize);
         a + (h - hf as f64) * (b - a)
     }
 
@@ -98,7 +205,7 @@ impl OrderStatistics<f64> for [f64] {
     }
 }
 
-impl Min<f64> for [f64] {
+impl<D: AsMut<[f64]> + AsRef<[f64]>> Min<f64> for Data<D> {
     /// Returns the minimum value in the data
     ///
     /// # Remarks
@@ -108,24 +215,27 @@ impl Min<f64> for [f64] {
     /// # Examples
     ///
     /// ```
-    /// use std::f64;
     /// use statrs::statistics::Min;
+    /// use statrs::statistics::Data;
     ///
-    /// let x: [f64; 0] = [];
+    /// let x = [];
+    /// let x = Data::new(x);
     /// assert!(x.min().is_nan());
     ///
     /// let y = [0.0, f64::NAN, 3.0, -2.0];
+    /// let y = Data::new(y);
     /// assert!(y.min().is_nan());
     ///
     /// let z = [0.0, 3.0, -2.0];
+    /// let z = Data::new(z);
     /// assert_eq!(z.min(), -2.0);
     /// ```
     fn min(&self) -> f64 {
-        Statistics::min(self)
+        Statistics::min(self.iter())
     }
 }
 
-impl Max<f64> for [f64] {
+impl<D: AsMut<[f64]> + AsRef<[f64]>> Max<f64> for Data<D> {
     /// Returns the maximum value in the data
     ///
     /// # Remarks
@@ -135,24 +245,27 @@ impl Max<f64> for [f64] {
     /// # Examples
     ///
     /// ```
-    /// use std::f64;
     /// use statrs::statistics::Max;
+    /// use statrs::statistics::Data;
     ///
-    /// let x: [f64; 0] = [];
+    /// let x = [];
+    /// let x = Data::new(x);
     /// assert!(x.max().is_nan());
     ///
     /// let y = [0.0, f64::NAN, 3.0, -2.0];
+    /// let y = Data::new(y);
     /// assert!(y.max().is_nan());
     ///
     /// let z = [0.0, 3.0, -2.0];
+    /// let z = Data::new(z);
     /// assert_eq!(z.max(), 3.0);
     /// ```
     fn max(&self) -> f64 {
-        Statistics::max(self)
+        Statistics::max(self.iter())
     }
 }
 
-impl Distribution<f64> for [f64] {
+impl<D: AsMut<[f64]> + AsRef<[f64]>> Distribution<f64> for Data<D> {
     /// Evaluates the sample mean, an estimate of the population
     /// mean.
     ///
@@ -166,22 +279,25 @@ impl Distribution<f64> for [f64] {
     /// #[macro_use]
     /// extern crate statrs;
     ///
-    /// use std::f64;
     /// use statrs::statistics::Distribution;
+    /// use statrs::statistics::Data;
     ///
     /// # fn main() {
     /// let x = [];
+    /// let x = Data::new(x);
     /// assert!(x.mean().unwrap().is_nan());
     ///
     /// let y = [0.0, f64::NAN, 3.0, -2.0];
+    /// let y = Data::new(y);
     /// assert!(y.mean().unwrap().is_nan());
     ///
     /// let z = [0.0, 3.0, -2.0];
+    /// let z = Data::new(z);
     /// assert_almost_eq!(z.mean().unwrap(), 1.0 / 3.0, 1e-15);
     /// # }
     /// ```
     fn mean(&self) -> Option<f64> {
-        Some(Statistics::mean(self))
+        Some(Statistics::mean(self.iter()))
     }
     /// Estimates the unbiased population variance from the provided samples
     ///
@@ -196,24 +312,27 @@ impl Distribution<f64> for [f64] {
     /// # Examples
     ///
     /// ```
-    /// use std::f64;
     /// use statrs::statistics::Distribution;
+    /// use statrs::statistics::Data;
     ///
     /// let x = [];
+    /// let x = Data::new(x);
     /// assert!(x.variance().unwrap().is_nan());
     ///
     /// let y = [0.0, f64::NAN, 3.0, -2.0];
+    /// let y = Data::new(y);
     /// assert!(y.variance().unwrap().is_nan());
     ///
     /// let z = [0.0, 3.0, -2.0];
+    /// let z = Data::new(z);
     /// assert_eq!(z.variance().unwrap(), 19.0 / 3.0);
     /// ```
     fn variance(&self) -> Option<f64> {
-        Some(Statistics::variance(self))
+        Some(Statistics::variance(self.iter()))
     }
 }
 
-impl Median<f64> for [f64] {
+impl<D: AsMut<[f64]> + AsRef<[f64]> + Clone> Median<f64> for Data<D> {
     /// Returns the median value from the data
     ///
     /// # Remarks
@@ -224,15 +343,18 @@ impl Median<f64> for [f64] {
     ///
     /// ```
     /// use statrs::statistics::Median;
+    /// use statrs::statistics::Data;
     ///
     /// let x = [];
+    /// let x = Data::new(x);
     /// assert!(x.median().is_nan());
     ///
     /// let y = [0.0, 3.0, -2.0];
+    /// let y = Data::new(y);
     /// assert_eq!(y.median(), 0.0);
     fn median(&self) -> f64 {
-        let mut copy = self.to_vec();
-        OrderStatistics::median(&mut *copy)
+        let mut v = self.clone();
+        OrderStatistics::median(&mut v)
     }
 }
 
@@ -255,82 +377,15 @@ fn handle_rank_ties(
     }
 }
 
-// Selection algorithm from Numerical Recipes
-// See: https://en.wikipedia.org/wiki/Selection_algorithm
-fn select_inplace(arr: &mut [f64], rank: usize) -> f64 {
-    if rank == 0 {
-        return arr.min();
-    }
-    if rank > arr.len() - 1 {
-        return arr.max();
-    }
-
-    let mut low = 0;
-    let mut high = arr.len() - 1;
-    loop {
-        if high <= low + 1 {
-            if high == low + 1 && arr[high] < arr[low] {
-                arr.swap(low, high)
-            }
-            return arr[rank];
-        }
-
-        let middle = (low + high) / 2;
-        arr.swap(middle, low + 1);
-
-        if arr[low] > arr[high] {
-            arr.swap(low, high);
-        }
-        if arr[low + 1] > arr[high] {
-            arr.swap(low + 1, high);
-        }
-        if arr[low] > arr[low + 1] {
-            arr.swap(low, low + 1);
-        }
-
-        let mut begin = low + 1;
-        let mut end = high;
-        let pivot = arr[begin];
-        loop {
-            loop {
-                begin += 1;
-                if arr[begin] >= pivot {
-                    break;
-                }
-            }
-            loop {
-                end -= 1;
-                if arr[end] <= pivot {
-                    break;
-                }
-            }
-            if end < begin {
-                break;
-            }
-            arr.swap(begin, end);
-        }
-
-        arr[low + 1] = arr[end];
-        arr[end] = pivot;
-
-        if end >= rank {
-            high = end - 1;
-        }
-        if end <= rank {
-            low = begin;
-        }
-    }
-}
-
-#[rustfmt::skip]
 #[cfg(test)]
 mod tests {
-    use std::f64;
+    use super::*;
     use crate::statistics::*;
 
     #[test]
     fn test_order_statistic_short() {
-        let mut data = [-1.0, 5.0, 0.0, -3.0, 10.0, -0.5, 4.0, 1.0, 6.0];
+        let data = [-1.0, 5.0, 0.0, -3.0, 10.0, -0.5, 4.0, 1.0, 6.0];
+        let mut data = Data::new(data);
         assert!(data.order_statistic(0).is_nan());
         assert_eq!(data.order_statistic(1), -3.0);
         assert_eq!(data.order_statistic(2), -1.0);
@@ -343,7 +398,8 @@ mod tests {
 
     #[test]
     fn test_quantile_short() {
-        let mut data = [-1.0, 5.0, 0.0, -3.0, 10.0, -0.5, 4.0, 0.2, 1.0, 6.0];
+        let data = [-1.0, 5.0, 0.0, -3.0, 10.0, -0.5, 4.0, 0.2, 1.0, 6.0];
+        let mut data = Data::new(data);
         assert_eq!(data.quantile(0.0), -3.0);
         assert_eq!(data.quantile(1.0), 10.0);
         assert_almost_eq!(data.quantile(0.5), 3.0 / 5.0, 1e-15);
@@ -357,44 +413,100 @@ mod tests {
 
     #[test]
     fn test_ranks() {
-        let mut sorted_distinct = [1.0, 2.0, 4.0, 7.0, 8.0, 9.0, 10.0, 12.0];
-        let mut sorted_ties = [1.0, 2.0, 2.0, 7.0, 9.0, 9.0, 10.0, 12.0];
-        assert_eq!(sorted_distinct.ranks(RankTieBreaker::Average), [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0]);
-        assert_eq!(sorted_ties.ranks(RankTieBreaker::Average), [1.0, 2.5, 2.5, 4.0, 5.5, 5.5, 7.0, 8.0]);
-        assert_eq!(sorted_distinct.ranks(RankTieBreaker::Min), [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0]);
-        assert_eq!(sorted_ties.ranks(RankTieBreaker::Min), [1.0, 2.0, 2.0, 4.0, 5.0, 5.0, 7.0, 8.0]);
-        assert_eq!(sorted_distinct.ranks(RankTieBreaker::Max), [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0]);
-        assert_eq!(sorted_ties.ranks(RankTieBreaker::Max), [1.0, 3.0, 3.0, 4.0, 6.0, 6.0, 7.0, 8.0]);
-        assert_eq!(sorted_distinct.ranks(RankTieBreaker::First), [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0]);
-        assert_eq!(sorted_ties.ranks(RankTieBreaker::First), [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0]);
+        let sorted_distinct = [1.0, 2.0, 4.0, 7.0, 8.0, 9.0, 10.0, 12.0];
+        let mut sorted_distinct = Data::new(sorted_distinct);
+        let sorted_ties = [1.0, 2.0, 2.0, 7.0, 9.0, 9.0, 10.0, 12.0];
+        let mut sorted_ties = Data::new(sorted_ties);
+        assert_eq!(
+            sorted_distinct.ranks(RankTieBreaker::Average),
+            [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0]
+        );
+        assert_eq!(
+            sorted_ties.ranks(RankTieBreaker::Average),
+            [1.0, 2.5, 2.5, 4.0, 5.5, 5.5, 7.0, 8.0]
+        );
+        assert_eq!(
+            sorted_distinct.ranks(RankTieBreaker::Min),
+            [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0]
+        );
+        assert_eq!(
+            sorted_ties.ranks(RankTieBreaker::Min),
+            [1.0, 2.0, 2.0, 4.0, 5.0, 5.0, 7.0, 8.0]
+        );
+        assert_eq!(
+            sorted_distinct.ranks(RankTieBreaker::Max),
+            [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0]
+        );
+        assert_eq!(
+            sorted_ties.ranks(RankTieBreaker::Max),
+            [1.0, 3.0, 3.0, 4.0, 6.0, 6.0, 7.0, 8.0]
+        );
+        assert_eq!(
+            sorted_distinct.ranks(RankTieBreaker::First),
+            [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0]
+        );
+        assert_eq!(
+            sorted_ties.ranks(RankTieBreaker::First),
+            [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0]
+        );
 
         let distinct = [1.0, 8.0, 12.0, 7.0, 2.0, 9.0, 10.0, 4.0];
+        let distinct = Data::new(distinct);
         let ties = [1.0, 9.0, 12.0, 7.0, 2.0, 9.0, 10.0, 2.0];
-        assert_eq!(distinct.clone().ranks(RankTieBreaker::Average), [1.0, 5.0, 8.0, 4.0, 2.0, 6.0, 7.0, 3.0]);
-        assert_eq!(ties.clone().ranks(RankTieBreaker::Average), [1.0, 5.5, 8.0, 4.0, 2.5, 5.5, 7.0, 2.5]);
-        assert_eq!(distinct.clone().ranks(RankTieBreaker::Min), [1.0, 5.0, 8.0, 4.0, 2.0, 6.0, 7.0, 3.0]);
-        assert_eq!(ties.clone().ranks(RankTieBreaker::Min), [1.0, 5.0, 8.0, 4.0, 2.0, 5.0, 7.0, 2.0]);
-        assert_eq!(distinct.clone().ranks(RankTieBreaker::Max), [1.0, 5.0, 8.0, 4.0, 2.0, 6.0, 7.0, 3.0]);
-        assert_eq!(ties.clone().ranks(RankTieBreaker::Max), [1.0, 6.0, 8.0, 4.0, 3.0, 6.0, 7.0, 3.0]);
-        assert_eq!(distinct.clone().ranks(RankTieBreaker::First), [1.0, 5.0, 8.0, 4.0, 2.0, 6.0, 7.0, 3.0]);
-        assert_eq!(ties.clone().ranks(RankTieBreaker::First), [1.0, 5.0, 8.0, 4.0, 2.0, 6.0, 7.0, 3.0]);
+        let ties = Data::new(ties);
+        assert_eq!(
+            distinct.clone().ranks(RankTieBreaker::Average),
+            [1.0, 5.0, 8.0, 4.0, 2.0, 6.0, 7.0, 3.0]
+        );
+        assert_eq!(
+            ties.clone().ranks(RankTieBreaker::Average),
+            [1.0, 5.5, 8.0, 4.0, 2.5, 5.5, 7.0, 2.5]
+        );
+        assert_eq!(
+            distinct.clone().ranks(RankTieBreaker::Min),
+            [1.0, 5.0, 8.0, 4.0, 2.0, 6.0, 7.0, 3.0]
+        );
+        assert_eq!(
+            ties.clone().ranks(RankTieBreaker::Min),
+            [1.0, 5.0, 8.0, 4.0, 2.0, 5.0, 7.0, 2.0]
+        );
+        assert_eq!(
+            distinct.clone().ranks(RankTieBreaker::Max),
+            [1.0, 5.0, 8.0, 4.0, 2.0, 6.0, 7.0, 3.0]
+        );
+        assert_eq!(
+            ties.clone().ranks(RankTieBreaker::Max),
+            [1.0, 6.0, 8.0, 4.0, 3.0, 6.0, 7.0, 3.0]
+        );
+        assert_eq!(
+            distinct.clone().ranks(RankTieBreaker::First),
+            [1.0, 5.0, 8.0, 4.0, 2.0, 6.0, 7.0, 3.0]
+        );
+        assert_eq!(
+            ties.clone().ranks(RankTieBreaker::First),
+            [1.0, 5.0, 8.0, 4.0, 2.0, 6.0, 7.0, 3.0]
+        );
     }
 
     #[test]
     fn test_median_short() {
         let even = [-1.0, 5.0, 0.0, -3.0, 10.0, -0.5, 4.0, 0.2, 1.0, 6.0];
+        let even = Data::new(even);
         assert_eq!(even.median(), 0.6);
 
         let odd = [-1.0, 5.0, 0.0, -3.0, 10.0, -0.5, 4.0, 0.2, 1.0];
+        let odd = Data::new(odd);
         assert_eq!(odd.median(), 0.2);
     }
 
     #[test]
     fn test_median_long_constant_seq() {
         let even = vec![2.0; 100000];
+        let even = Data::new(even);
         assert_eq!(2.0, even.median());
 
         let odd = vec![2.0; 100001];
+        let odd = Data::new(odd);
         assert_eq!(2.0, odd.median());
     }
 
@@ -402,20 +514,30 @@ mod tests {
 
     #[test]
     fn test_median_robust_on_infinities() {
-        let mut data3 = [2.0, f64::NEG_INFINITY, f64::INFINITY];
+        let data3 = [2.0, f64::NEG_INFINITY, f64::INFINITY];
+        let data3 = Data::new(data3);
         assert_eq!(data3.median(), 2.0);
         assert_eq!(data3.median(), 2.0);
 
-        data3 = [f64::NEG_INFINITY, 2.0, f64::INFINITY];
+        let data3 = [f64::NEG_INFINITY, 2.0, f64::INFINITY];
+        let data3 = Data::new(data3);
         assert_eq!(data3.median(), 2.0);
         assert_eq!(data3.median(), 2.0);
 
-        data3 = [f64::NEG_INFINITY, f64::INFINITY, 2.0];
+        let data3 = [f64::NEG_INFINITY, f64::INFINITY, 2.0];
+        let data3 = Data::new(data3);
         assert_eq!(data3.median(), 2.0);
         assert_eq!(data3.median(), 2.0);
 
         let data4 = [f64::NEG_INFINITY, 2.0, 3.0, f64::INFINITY];
+        let data4 = Data::new(data4);
         assert_eq!(data4.median(), 2.5);
         assert_eq!(data4.median(), 2.5);
     }
+    #[test]
+    fn test_foo() {
+        let arr = [0.0, 1.0, 2.0, 3.0];
+        let mut arr = Data::new(arr);
+        arr.order_statistic(2);
+    }
 }


Test patch: