diff --git a/CHANGELOG.md b/CHANGELOG.md
index 9a39cea0..8d1293e7 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,3 +1,16 @@
+v0.14.0
+- upgrade `rand` dependency to `0.8`
+- fix inaccurate sampling of `Gamma`
+- Implemented Empirical distribution
+- Implemented Laplace distribution
+- Removed Checked* traits
+- Almost clippy-clean
+- Almost fully enabled rustfmt
+- Begin applying consistent numeric relative-accuracy targets with the approx crate
+- Introduce macro to generate testing boilerplate, yet not all tests use this yet
+- Moved to dynamic vectors in the MultivariateNormal distribution
+- Reduced a number of distribution-specific traits into the Distribution and DiscreteDistribution traits
+
 v0.13.0
 - Implemented `MultivariateNormal` distribution (depends on `nalgebra 0.19`)
 - Implemented `Dirac` distribution
diff --git a/src/distribution/beta.rs b/src/distribution/beta.rs
index a477ff4d..17da2a04 100644
--- a/src/distribution/beta.rs
+++ b/src/distribution/beta.rs
@@ -291,7 +291,7 @@ impl Continuous<f64, f64> for Beta {
     ///
     /// where `α` is shapeA, `β` is shapeB, and `Γ` is the gamma function
     fn pdf(&self, x: f64) -> f64 {
-        if x < 0.0 || x > 1.0 {
+        if !(0.0..=1.0).contains(&x) {
             0.0
         } else if self.shape_a.is_infinite() {
             if ulps_eq!(x, 1.0) {
@@ -329,7 +329,7 @@ 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 {
+        if !(0.0..=1.0).contains(&x) {
             -INF
         } else if self.shape_a.is_infinite() {
             if ulps_eq!(x, 1.0) {
diff --git a/src/distribution/categorical.rs b/src/distribution/categorical.rs
index ea53eae1..488773ed 100644
--- a/src/distribution/categorical.rs
+++ b/src/distribution/categorical.rs
@@ -282,7 +282,7 @@ pub fn prob_mass_to_cdf(prob_mass: &[f64]) -> Vec<f64> {
 fn binary_index(search: &[f64], val: f64) -> usize {
     use std::cmp;
 
-    let mut low = 0 as isize;
+    let mut low = 0_isize;
     let mut high = search.len() as isize - 1;
     while low <= high {
         let mid = low + ((high - low) / 2);
diff --git a/src/distribution/chi.rs b/src/distribution/chi.rs
index fc674f76..3fa7f1a3 100644
--- a/src/distribution/chi.rs
+++ b/src/distribution/chi.rs
@@ -145,7 +145,7 @@ impl Distribution<f64> for Chi {
     /// where `k` is degrees of freedom and `Γ` is the gamma function
     fn mean(&self) -> Option<f64> {
         if self.freedom.is_infinite() {
-            return None;
+            None
         } else if self.freedom > 300.0 {
             // Large n approximation based on the Stirling series approximation to the Gamma function
             // This avoids call the Gamma function with large arguments and returning NaN
diff --git a/src/distribution/negative_binomial.rs b/src/distribution/negative_binomial.rs
index 494a7068..14fffde6 100644
--- a/src/distribution/negative_binomial.rs
+++ b/src/distribution/negative_binomial.rs
@@ -51,7 +51,7 @@ impl NegativeBinomial {
         if p.is_nan() || p < 0.0 || p > 1.0 || r.is_nan() || r < 0.0 {
             Err(StatsError::BadParams)
         } else {
-            Ok(NegativeBinomial { p, r })
+            Ok(NegativeBinomial { r, p })
         }
     }
 
diff --git a/src/distribution/normal.rs b/src/distribution/normal.rs
index 2c8c7de5..d6f64c1a 100644
--- a/src/distribution/normal.rs
+++ b/src/distribution/normal.rs
@@ -90,7 +90,7 @@ impl ContinuousCDF<f64, f64> for Normal {
     /// where `μ` is the mean, `σ` is the standard deviation and `erfc_inv` is
     /// the inverse of the complementary error function
     fn inverse_cdf(&self, x: f64) -> f64 {
-        if x < 0.0 || x > 1.0 {
+        if !(0.0..=1.0).contains(&x) {
             panic!("x must be in [0, 1]");
         } else {
             self.mean - (self.std_dev * f64::consts::SQRT_2 * erf::erfc_inv(2.0 * x))
diff --git a/src/distribution/students_t.rs b/src/distribution/students_t.rs
index 23036285..8e50fb2a 100644
--- a/src/distribution/students_t.rs
+++ b/src/distribution/students_t.rs
@@ -155,7 +155,7 @@ impl ContinuousCDF<f64, f64> for StudentsT {
     /// 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);
+        assert!((0.0..=1.0).contains(&x));
         let x = 2. * x.min(1. - x);
         let a = 0.5 * self.freedom;
         let b = 0.5;
diff --git a/src/function/beta.rs b/src/function/beta.rs
index 616e1847..23fb3430 100644
--- a/src/function/beta.rs
+++ b/src/function/beta.rs
@@ -115,7 +115,7 @@ pub fn checked_beta_reg(a: f64, b: f64, x: f64) -> Result<f64> {
         Err(StatsError::ArgMustBePositive("a"))
     } else if b <= 0.0 {
         Err(StatsError::ArgMustBePositive("b"))
-    } else if x < 0.0 || x > 1.0 {
+    } else if !(0.0..=1.0).contains(&x) {
         Err(StatsError::ArgIntervalIncl("x", 0.0, 1.0))
     } else {
         let bt = if is_zero(x) || ulps_eq!(x, 1.0) {
@@ -275,7 +275,7 @@ pub fn inv_beta_reg(mut a: f64, mut b: f64, mut x: f64) -> f64 {
     const SAE: i32 = -30;
     const FPU: f64 = 1e-30; // 10^SAE
 
-    debug_assert!(x >= 0.0 && x <= 1.0 && a > 0.0 && b > 0.0);
+    debug_assert!((0.0..=1.0).contains(&x) && a > 0.0 && b > 0.0);
 
     if x == 0.0 {
         return 0.0;
diff --git a/src/function/logistic.rs b/src/function/logistic.rs
index 27cbe9cc..8adc8b79 100644
--- a/src/function/logistic.rs
+++ b/src/function/logistic.rs
@@ -24,7 +24,7 @@ pub fn logit(p: f64) -> f64 {
 ///
 /// If `p < 0.0` or `p > 1.0`
 pub fn checked_logit(p: f64) -> Result<f64> {
-    if p < 0.0 || p > 1.0 {
+    if !(0.0..=1.0).contains(&p) {
         Err(StatsError::ArgIntervalIncl("p", 0.0, 1.0))
     } else {
         Ok((p / (1.0 - p)).ln())
diff --git a/src/statistics/slice_statistics.rs b/src/statistics/slice_statistics.rs
index abb6cae7..1d1b79cc 100644
--- a/src/statistics/slice_statistics.rs
+++ b/src/statistics/slice_statistics.rs
@@ -129,7 +129,7 @@ impl<D: AsMut<[f64]> + AsRef<[f64]>> OrderStatistics<f64> for Data<D> {
     }
 
     fn quantile(&mut self, tau: f64) -> f64 {
-        if tau < 0.0 || tau > 1.0 || self.is_empty() {
+        if !(0.0..=1.0).contains(&tau) || self.is_empty() {
             return f64::NAN;
         }
 
diff --git a/src/statistics/statistics.rs b/src/statistics/statistics.rs
index 7ab6836b..40b91c72 100644
--- a/src/statistics/statistics.rs
+++ b/src/statistics/statistics.rs
@@ -28,14 +28,14 @@ pub trait Statistics<T> {
     /// use std::f64;
     /// use statrs::statistics::Statistics;
     ///
-    /// let x: [f64; 0] = [];
-    /// assert!(x.min().is_nan());
+    /// let x = &[];
+    /// assert!(Statistics::min(x).is_nan());
     ///
-    /// let y = [0.0, f64::NAN, 3.0, -2.0];
-    /// assert!(y.min().is_nan());
+    /// let y = &[0.0, f64::NAN, 3.0, -2.0];
+    /// assert!(Statistics::min(y).is_nan());
     ///
-    /// let z = [0.0, 3.0, -2.0];
-    /// assert_eq!(z.min(), -2.0);
+    /// let z = &[0.0, 3.0, -2.0];
+    /// assert_eq!(Statistics::min(z), -2.0);
     /// ```
     fn min(self) -> T;
 
@@ -51,14 +51,14 @@ pub trait Statistics<T> {
     /// use std::f64;
     /// use statrs::statistics::Statistics;
     ///
-    /// let x: [f64; 0] = [];
-    /// assert!(x.max().is_nan());
+    /// let x = &[];
+    /// assert!(Statistics::max(x).is_nan());
     ///
-    /// let y = [0.0, f64::NAN, 3.0, -2.0];
-    /// assert!(y.max().is_nan());
+    /// let y = &[0.0, f64::NAN, 3.0, -2.0];
+    /// assert!(Statistics::max(y).is_nan());
     ///
-    /// let z = [0.0, 3.0, -2.0];
-    /// assert_eq!(z.max(), 3.0);
+    /// let z = &[0.0, 3.0, -2.0];
+    /// assert_eq!(Statistics::max(z), 3.0);
     /// ```
     fn max(self) -> T;
 
@@ -74,13 +74,13 @@ pub trait Statistics<T> {
     /// use std::f64;
     /// use statrs::statistics::Statistics;
     ///
-    /// let x = [];
+    /// let x = &[];
     /// assert!(x.abs_min().is_nan());
     ///
-    /// let y = [0.0, f64::NAN, 3.0, -2.0];
+    /// let y = &[0.0, f64::NAN, 3.0, -2.0];
     /// assert!(y.abs_min().is_nan());
     ///
-    /// let z = [0.0, 3.0, -2.0];
+    /// let z = &[0.0, 3.0, -2.0];
     /// assert_eq!(z.abs_min(), 0.0);
     /// ```
     fn abs_min(self) -> T;
@@ -97,13 +97,13 @@ pub trait Statistics<T> {
     /// use std::f64;
     /// use statrs::statistics::Statistics;
     ///
-    /// let x = [];
+    /// let x = &[];
     /// assert!(x.abs_max().is_nan());
     ///
-    /// let y = [0.0, f64::NAN, 3.0, -2.0];
+    /// let y = &[0.0, f64::NAN, 3.0, -2.0];
     /// assert!(y.abs_max().is_nan());
     ///
-    /// let z = [0.0, 3.0, -2.0, -8.0];
+    /// let z = &[0.0, 3.0, -2.0, -8.0];
     /// assert_eq!(z.abs_max(), 8.0);
     /// ```
     fn abs_max(self) -> T;
@@ -125,13 +125,13 @@ pub trait Statistics<T> {
     /// use statrs::statistics::Statistics;
     ///
     /// # fn main() {
-    /// let x = [];
+    /// let x = &[];
     /// assert!(x.mean().is_nan());
     ///
-    /// let y = [0.0, f64::NAN, 3.0, -2.0];
+    /// let y = &[0.0, f64::NAN, 3.0, -2.0];
     /// assert!(y.mean().is_nan());
     ///
-    /// let z = [0.0, 3.0, -2.0];
+    /// let z = &[0.0, 3.0, -2.0];
     /// assert_almost_eq!(z.mean(), 1.0 / 3.0, 1e-15);
     /// # }
     /// ```
@@ -155,19 +155,19 @@ pub trait Statistics<T> {
     /// use statrs::statistics::Statistics;
     ///
     /// # fn main() {
-    /// let x = [];
+    /// let x = &[];
     /// assert!(x.geometric_mean().is_nan());
     ///
-    /// let y = [0.0, f64::NAN, 3.0, -2.0];
+    /// let y = &[0.0, f64::NAN, 3.0, -2.0];
     /// assert!(y.geometric_mean().is_nan());
     ///
-    /// let mut z = [0.0, 3.0, -2.0];
+    /// let mut z = &[0.0, 3.0, -2.0];
     /// assert!(z.geometric_mean().is_nan());
     ///
-    /// z = [0.0, 3.0, 2.0];
+    /// z = &[0.0, 3.0, 2.0];
     /// assert_eq!(z.geometric_mean(), 0.0);
     ///
-    /// z = [1.0, 2.0, 3.0];
+    /// z = &[1.0, 2.0, 3.0];
     /// // test value from online calculator, could be more accurate
     /// assert_almost_eq!(z.geometric_mean(), 1.81712, 1e-5);
     /// # }
@@ -194,19 +194,19 @@ pub trait Statistics<T> {
     /// use statrs::statistics::Statistics;
     ///
     /// # fn main() {
-    /// let x = [];
+    /// let x = &[];
     /// assert!(x.harmonic_mean().is_nan());
     ///
-    /// let y = [0.0, f64::NAN, 3.0, -2.0];
+    /// let y = &[0.0, f64::NAN, 3.0, -2.0];
     /// assert!(y.harmonic_mean().is_nan());
     ///
-    /// let mut z = [0.0, 3.0, -2.0];
+    /// let mut z = &[0.0, 3.0, -2.0];
     /// assert!(z.harmonic_mean().is_nan());
     ///
-    /// z = [0.0, 3.0, 2.0];
+    /// z = &[0.0, 3.0, 2.0];
     /// assert_eq!(z.harmonic_mean(), 0.0);
     ///
-    /// z = [1.0, 2.0, 3.0];
+    /// z = &[1.0, 2.0, 3.0];
     /// // test value from online calculator, could be more accurate
     /// assert_almost_eq!(z.harmonic_mean(), 1.63636, 1e-5);
     /// # }
@@ -229,13 +229,13 @@ pub trait Statistics<T> {
     /// use std::f64;
     /// use statrs::statistics::Statistics;
     ///
-    /// let x = [];
+    /// let x = &[];
     /// assert!(x.variance().is_nan());
     ///
-    /// let y = [0.0, f64::NAN, 3.0, -2.0];
+    /// let y = &[0.0, f64::NAN, 3.0, -2.0];
     /// assert!(y.variance().is_nan());
     ///
-    /// let z = [0.0, 3.0, -2.0];
+    /// let z = &[0.0, 3.0, -2.0];
     /// assert_eq!(z.variance(), 19.0 / 3.0);
     /// ```
     fn variance(self) -> T;
@@ -257,13 +257,13 @@ pub trait Statistics<T> {
     /// use std::f64;
     /// use statrs::statistics::Statistics;
     ///
-    /// let x = [];
+    /// let x = &[];
     /// assert!(x.std_dev().is_nan());
     ///
-    /// let y = [0.0, f64::NAN, 3.0, -2.0];
+    /// let y = &[0.0, f64::NAN, 3.0, -2.0];
     /// assert!(y.std_dev().is_nan());
     ///
-    /// let z = [0.0, 3.0, -2.0];
+    /// let z = &[0.0, 3.0, -2.0];
     /// assert_eq!(z.std_dev(), (19f64 / 3.0).sqrt());
     /// ```
     fn std_dev(self) -> T;
@@ -283,13 +283,13 @@ pub trait Statistics<T> {
     /// use std::f64;
     /// use statrs::statistics::Statistics;
     ///
-    /// let x = [];
+    /// let x = &[];
     /// assert!(x.population_variance().is_nan());
     ///
-    /// let y = [0.0, f64::NAN, 3.0, -2.0];
+    /// let y = &[0.0, f64::NAN, 3.0, -2.0];
     /// assert!(y.population_variance().is_nan());
     ///
-    /// let z = [0.0, 3.0, -2.0];
+    /// let z = &[0.0, 3.0, -2.0];
     /// assert_eq!(z.population_variance(), 38.0 / 9.0);
     /// ```
     fn population_variance(self) -> T;
@@ -309,13 +309,13 @@ pub trait Statistics<T> {
     /// use std::f64;
     /// use statrs::statistics::Statistics;
     ///
-    /// let x = [];
+    /// let x = &[];
     /// assert!(x.population_std_dev().is_nan());
     ///
-    /// let y = [0.0, f64::NAN, 3.0, -2.0];
+    /// let y = &[0.0, f64::NAN, 3.0, -2.0];
     /// assert!(y.population_std_dev().is_nan());
     ///
-    /// let z = [0.0, 3.0, -2.0];
+    /// let z = &[0.0, 3.0, -2.0];
     /// assert_eq!(z.population_std_dev(), (38f64 / 9.0).sqrt());
     /// ```
     fn population_std_dev(self) -> T;
@@ -345,16 +345,16 @@ pub trait Statistics<T> {
     /// use statrs::statistics::Statistics;
     ///
     /// # fn main() {
-    /// let x = [];
+    /// let x = &[];
     /// assert!(x.covariance(&[]).is_nan());
     ///
-    /// let y1 = [0.0, f64::NAN, 3.0, -2.0];
-    /// let y2 = [-5.0, 4.0, 10.0, f64::NAN];
-    /// assert!(y1.covariance(&y2).is_nan());
+    /// let y1 = &[0.0, f64::NAN, 3.0, -2.0];
+    /// let y2 = &[-5.0, 4.0, 10.0, f64::NAN];
+    /// assert!(y1.covariance(y2).is_nan());
     ///
-    /// let z1 = [0.0, 3.0, -2.0];
-    /// let z2 = [-5.0, 4.0, 10.0];
-    /// assert_almost_eq!(z1.covariance(&z2), -5.5, 1e-14);
+    /// let z1 = &[0.0, 3.0, -2.0];
+    /// let z2 = &[-5.0, 4.0, 10.0];
+    /// assert_almost_eq!(z1.covariance(z2), -5.5, 1e-14);
     /// # }
     /// ```
     fn covariance(self, other: Self) -> T;
@@ -382,16 +382,16 @@ pub trait Statistics<T> {
     /// use statrs::statistics::Statistics;
     ///
     /// # fn main() {
-    /// let x = [];
+    /// let x = &[];
     /// assert!(x.population_covariance(&[]).is_nan());
     ///
-    /// let y1 = [0.0, f64::NAN, 3.0, -2.0];
-    /// let y2 = [-5.0, 4.0, 10.0, f64::NAN];
-    /// assert!(y1.population_covariance(&y2).is_nan());
+    /// let y1 = &[0.0, f64::NAN, 3.0, -2.0];
+    /// let y2 = &[-5.0, 4.0, 10.0, f64::NAN];
+    /// assert!(y1.population_covariance(y2).is_nan());
     ///
-    /// let z1 = [0.0, 3.0, -2.0];
-    /// let z2 = [-5.0, 4.0, 10.0];
-    /// assert_almost_eq!(z1.population_covariance(&z2), -11.0 / 3.0, 1e-14);
+    /// let z1 = &[0.0, 3.0, -2.0];
+    /// let z2 = &[-5.0, 4.0, 10.0];
+    /// assert_almost_eq!(z1.population_covariance(z2), -11.0 / 3.0, 1e-14);
     /// # }
     /// ```
     fn population_covariance(self, other: Self) -> T;
@@ -412,13 +412,13 @@ pub trait Statistics<T> {
     /// use statrs::statistics::Statistics;
     ///
     /// # fn main() {
-    /// let x = [];
+    /// let x = &[];
     /// assert!(x.quadratic_mean().is_nan());
     ///
-    /// let y = [0.0, f64::NAN, 3.0, -2.0];
+    /// let y = &[0.0, f64::NAN, 3.0, -2.0];
     /// assert!(y.quadratic_mean().is_nan());
     ///
-    /// let z = [0.0, 3.0, -2.0];
+    /// let z = &[0.0, 3.0, -2.0];
     /// // test value from online calculator, could be more accurate
     /// assert_almost_eq!(z.quadratic_mean(), 2.08167, 1e-5);
     /// # }