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

Issue:
Consider improving error handling and `StatsError` type
`enum StatsError` is supposed to be an

> Enumeration of possible errors thrown within the `statrs` library.

This indicates that it should not try to be a generic error type for any kind of statistics calculations, but instead only concern itself with the errors produced in `statrs`.

With that basic assumption, there are currently some inconsistencies and outdated practices though ([API guidelines](https://rust-lang.github.io/api-guidelines/interoperability.html#c-good-err) for reference):
- Add test that `StatsError` implements both `Sync` and `Send` (this is more of a formality, but is trivial to implement and good future-proofing)
- [`Error::description`](https://doc.rust-lang.org/std/error/trait.Error.html#method.description) is deprecated and should not be implemented
- There are a few unused variants (`ArgNotNegative`, `ArgIntervalExclMin`, etc.), these are leftovers from older versions and should be removed because they're no longer needed by `statrs`.
- There's at least one case of a function returning a `Result<T>` that cannot return an error: `Empirical::new`. There's no real reason for this infallible API to return a `Result<T, E>`. There might be others.

I realize that most of these are breaking changes, but seeing that the crate is pre-1.0, I don't think there's a big problem doing this.

Other things that could be improved:
`StatsError` is big: 40 bytes on a Linux x64 target. This is because there are variants which contain 2x `&str` (2 x 16 = 32 bytes plus discriminant and padding). Is it really necessary to have strings in the error type? The implementation could be replaced mostly 1:1 with some sort of `ArgName` enum, but there might be an even better solution that does not need this argument name wrapping at all.

All `new` functions seem to just return `StatsError::BadParams` whenever the params are outside of the defined/allowed range. Is there a good reason for these to be so vague when compared to the more specific errors returned by other functions? After all, the more specific errors already exist, why not return more exact error information? There might even be value in providing multiple error types, to have errors that are more specific to the exact API being used.


Discussion / hints:
I do see good reason for all of these and I'd be open to making changes for all but the infallible `new` not returning a result. 

All other public structs defined in the `distribution` module have a `new` method implemented that returns Result so it does provide consistency. Perhaps if `Empirical` were in a different module or it's `new` had a different name?
> All other public structs defined in the distribution module have a new method implemented that returns Result so it does provide consistency. Perhaps if Empirical were in a different module or it's new had a different name?

I see your point about consistency. I would personally value the expressiveness ("this call cannot fail") over the consistency ("all constructors return a Result and might need error handling"), but it doesn't matter much tbh.

Hmm I'm not sure about renaming the `new` function, it's a widespread naming convention in the Rust ecosystem and what most users would expect. Maybe a similar name like `new_<something>` so people can quickly find it in their IDEs.
Perhaps an `impl Default` over having `Empirical::new`? 

Regardless, the overall discussion you bring up on the error type is valid. I'd merge an effort that does any of
- scaling it down 
- making it adhere closer to API guidelines
- returning `StatsError::BadParams` less where possible.
During these changes, we should also work on removing support for degenerate distributions #102 and returning errors.
Tried a different approach, even wrote out an [example](https://github.com/YeungOnion/statrs/blob/feat/errors/examples/error_trials.rs) to test using it, but I didn't have a great motivating case, so I just made something up, but I wanted to see what it was like to actually handle those such errors.
Maybe we should figure out what the errors should be able to convey first. Since @FreezyLemon posted on #242 I've realized I've not done enough of this bigger picture stuff for decisions. I think it's fair to narrow the discussion to errors from distribution construction since that's at a module scope and most of where errors are handled.

The value in #258 is to make the errors flat, and delegates identifying what we call "bad parameters" to the caller. This is quite reasonable, as the conditions for errors are statically specified (up to decisions about what constitutes degeneracy), but there is no program state in which parameter conditions change. However it means that users of the library may repeat their work among projects where there's some input-validation cycle if they want useful error messages. 

The option I have in #247 is pretty dense, and it means all distributions' `new` have a different type for the Err variant. Would need quite a few imports just to do something that creates multiple distributions from different families. But it has the value of granularity.

---

A few lists to help organize my thoughts:

## Assumptions
- Users who are writing analyses will likely `unwrap` or pair `?` with `anyhow`.
- Some end users will benefit from specificity of erroneous input.
- Library authors using `statrs` wish for choice between dynamic and static dispatch.
- Library authors will not prefer to duplicate effort of input checking.
- Developers do not prefer being exhaustive with writing source code for error handling outside of tests. 

## Possible targeted behaviors adjacent to error API
- [ ] distribution family generic code
  - e.g. the application is fitting a distribution function to data where the end user chooses a distribution family
  - dynamic dispatch for creating distribution type is the only option we can support without changing beyond error API
- [ ] correcting while widening and narrowing of parameters
  - e.g. arguing again via fitting, multiparameter distributions could have jacobians that lead to invalid values for one (or multiple) parameters going out of bounds because types are not statically constraining onto intervals.
- [ ] matrix updating
   - #180 mentions possible memory performance benefit, a further benefit could be updates that have nice approximations as updates, for $M' = M + \epsilon A$ under appropriate constraints on $A$, perhaps symmetry, $A^T = A$ or single rank $A = x x$
- [ ] dimension checking/mismatch for dynamically sized `nalgebra` matrices
   - we should at least `panic` on things like this for now. Evaluating pdf for `Dyn` sizes may fail the [unwrapping here](https://github.com/statrs-dev/statrs/blob/039c54509c35444b7a785d7760ebf852240231c4/src/distribution/multivariate_normal.rs#L254)

I'm very open to modifying these lists.
I realize yet another option is to have the distribution error have too many variants for any one distribution, but use shared language `location`, `scale`, `shape`, etc. in the variants. Suppose this is called `ParameterError`. Then in each distribution's module, `impl From<super::ParameterError>` for the distribution's specific error or use some generics to impl Error with source converting to the right error.

This allows a shallow error, small step from unwrap for those writing analyses. And allows for deeper error which can rely more on handling programmatically. But I don't think this is a common idiom, the more I think about this, the more certain there's something easier I'm missing here and I don't think it's needed to replicate `std::io`'s pattern of error handling.
> the more I think about this, the more certain there's something easier I'm missing here and I don't think it's needed to replicate std::io's pattern of error handling.

Most usages of `StatsError` are currently for the return value of `SomeDistribution::new(...)`. Since this function is implemented per-distribution and not to satisfy some trait, the simplest solution here would be creating a concrete error type for each and every distribution in case the construction fails. For example:

Let's turn

```rs
// beta.rs

#[derive(Copy, Clone, PartialEq, Debug)]
pub struct Beta {
    shape_a: f64,
    shape_b: f64,
}

impl Beta {
    pub fn new(shape_a: f64, shape_b: f64) -> Result<Beta> {
        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 })
    }
}
```

into

```rs
#[derive(Copy, Clone, PartialEq, Debug)]
pub struct Beta {
    shape_a: f64,
    shape_b: f64,
}

#[derive(Copy, Clone, PartialEq, Eq, Debug)]
pub enum BetaError {
    ShapeANan,
    ShapeAInfinite,
    ShapeALessThanZero,
    ShapeBNan,
    ShapeBInfinite,
    ShapeBLessThanZero,
}

impl std::fmt::Display for BetaError {
    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
        match self {
            BetaError::ShapeANan => write!(f, "Shape A is NaN"),
            // ...
        }
    }
}

impl std::error::Error for BetaError {}

impl Beta {
    pub fn new(shape_a: f64, shape_b: f64) -> Result<Beta, BetaError> {
        if shape_a.is_nan() {
            Err(BetaError::ShapeANan)
        } else if shape_a.is_infinite() {
            Err(BetaError::ShapeAInfinite)
        } else if shape_a <= 0.0 {
            Err(BetaError::ShapeALessThanZero)
        } else if shape_b.is_nan() {
            Err(BetaError::ShapeBNan)
        } else if shape_b.is_infinite() {
            Err(BetaError::ShapeBInfinite)
        } else if shape_b <= 0.0 {
            Err(BetaError::ShapeBLessThanZero)
        } else {
            Ok(Beta { shape_a, shape_b })
        }
    }
}
```

Maybe this example is too verbose (i.e. maybe `BetaError::ShapeAInvalid` would be enough), but it gets the point across.

This will return very specific errors, and users that don't care about it can still `.unwrap()` or use pattern matching to handle whichever errors. This will not impact the distributions' traits (because `new` is not a part of any trait) and should allow for dynamic dispatch.

If new methods of creating a distribution are added later, the related error type could be extended or a separate error type, specific to the method of construction, could be added.

It will need a fair bit of boiler plate, especially now that it needs to be implemented all at once, but is otherwise very simple to implement. A crate like `thiserror` could cut down on boilerplate, but is honestly not needed. The creation errors are simple and `thiserror` is/has a proc-macro which means it has a negative impact on compilation times.

I might be missing something here, but this is the obvious KISS solution IMHO
Ah, it's true that the generic part can only come up once the type is correctly initialized.

I'm good with only specifying which parameter is invalid, so we can just have one variant per parameter, the one message denoting that one of its (potentially multiple) constraints was violated.
```rs
impl std::fmt::Display for BetaError {
    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
        match self {
            BetaError::ShapeAInvalid => write!(f, "Shape A must be positive finite"),
            BetaError::ShapeBInvalid => write!(f, "Shape B must be positive finite"),
        }
    }
}
```

It will be a little boilerplatey, but it need not be flexible, sounds like a text editor problem.
Alright, I'll prepare another draft pr to compare against the rest if you don't mind
Thanks for this, I'll keep an eye out. 

Fix patch:
diff --git a/src/distribution/multivariate_normal.rs b/src/distribution/multivariate_normal.rs
index b336d9db..33df4969 100644
--- a/src/distribution/multivariate_normal.rs
+++ b/src/distribution/multivariate_normal.rs
@@ -1,42 +1,37 @@
 use crate::distribution::Continuous;
 use crate::statistics::{Max, MeanN, Min, Mode, VarianceN};
-use crate::StatsError;
 use nalgebra::{Cholesky, Const, DMatrix, DVector, Dim, DimMin, Dyn, OMatrix, OVector};
 use std::f64;
 use std::f64::consts::{E, PI};
 
-/// computes both the normalization and exponential argument in the normal distribution
-/// # Errors
-/// will error on dimension mismatch
+/// Computes both the normalization and exponential argument in the normal
+/// distribution, returning `None` on dimension mismatch.
 pub(super) fn density_normalization_and_exponential<D>(
     mu: &OVector<f64, D>,
     cov: &OMatrix<f64, D, D>,
     precision: &OMatrix<f64, D, D>,
     x: &OVector<f64, D>,
-) -> std::result::Result<(f64, f64), StatsError>
+) -> Option<(f64, f64)>
 where
     D: DimMin<D, Output = D>,
     nalgebra::DefaultAllocator: nalgebra::allocator::Allocator<f64, D>
         + nalgebra::allocator::Allocator<f64, D, D>
         + nalgebra::allocator::Allocator<(usize, usize), D>,
 {
-    Ok((
+    Some((
         density_distribution_pdf_const(mu, cov)?,
         density_distribution_exponential(mu, precision, x)?,
     ))
 }
 
-/// computes the argument of the exponential term in the normal distribution
-/// ```text
-/// ```
-/// # Errors
-/// will error on dimension mismatch
+/// Computes the argument of the exponential term in the normal distribution,
+/// returning `None` on dimension mismatch.
 #[inline]
-pub(super) fn density_distribution_exponential<D>(
+fn density_distribution_exponential<D>(
     mu: &OVector<f64, D>,
     precision: &OMatrix<f64, D, D>,
     x: &OVector<f64, D>,
-) -> std::result::Result<f64, StatsError>
+) -> Option<f64>
 where
     D: Dim,
     nalgebra::DefaultAllocator:
@@ -46,22 +41,18 @@ where
         || x.shape_generic().0 != mu.shape_generic().0
         || !precision.is_square()
     {
-        return Err(StatsError::ContainersMustBeSameLength);
+        return None;
     }
+
     let dv = x - mu;
     let exp_term: f64 = -0.5 * (precision * &dv).dot(&dv);
-    Ok(exp_term)
-    // TODO update to dimension mismatch error
+    Some(exp_term)
 }
 
-/// computes the argument of the normalization term in the normal distribution
-/// # Errors
-/// will error on dimension mismatch
+/// Computes the argument of the normalization term in the normal distribution,
+/// returning `None` on dimension mismatch.
 #[inline]
-pub(super) fn density_distribution_pdf_const<D>(
-    mu: &OVector<f64, D>,
-    cov: &OMatrix<f64, D, D>,
-) -> std::result::Result<f64, StatsError>
+fn density_distribution_pdf_const<D>(mu: &OVector<f64, D>, cov: &OMatrix<f64, D, D>) -> Option<f64>
 where
     D: DimMin<D, Output = D>,
     nalgebra::DefaultAllocator: nalgebra::allocator::Allocator<f64, D>
@@ -69,12 +60,14 @@ where
         + nalgebra::allocator::Allocator<(usize, usize), D>,
 {
     if cov.shape_generic().0 != mu.shape_generic().0 || !cov.is_square() {
-        return Err(StatsError::ContainersMustBeSameLength);
+        return None;
     }
     let cov_det = cov.determinant();
-    Ok(((2. * PI).powi(mu.nrows() as i32) * cov_det.abs())
-        .recip()
-        .sqrt())
+    Some(
+        ((2. * PI).powi(mu.nrows() as i32) * cov_det.abs())
+            .recip()
+            .sqrt(),
+    )
 }
 
 /// Implements the [Multivariate Normal](https://en.wikipedia.org/wiki/Multivariate_normal_distribution)
diff --git a/src/distribution/multivariate_students_t.rs b/src/distribution/multivariate_students_t.rs
index 73e8f8f2..05532431 100644
--- a/src/distribution/multivariate_students_t.rs
+++ b/src/distribution/multivariate_students_t.rs
@@ -340,10 +340,13 @@ where
     /// [Γ(ν+p)/2] / [Γ(ν/2) ((ν * π)^p det(Σ))^(1 / 2)] * [1 + 1/ν (x - μ)ᵀ inv(Σ) (x - μ)]^(-(ν+p)/2)
     /// ```
     ///
-    /// where `ν` is the degrees of freedom,  `μ` is the mean, `Γ`
-    /// is the Gamma function, `inv(Σ)`
-    /// is the precision matrix, `det(Σ)` is the determinant
-    /// of the scale matrix, and `k` is the dimension of the distribution.
+    /// where
+    /// - `ν` is the degrees of freedom,
+    /// - `μ` is the mean,
+    /// - `Γ` is the Gamma function,
+    /// - `inv(Σ)` is the precision matrix,
+    /// - `det(Σ)` is the determinant of the scale matrix, and
+    /// - `k` is the dimension of the distribution.
     fn pdf(&self, x: &'a OVector<f64, D>) -> f64 {
         if self.freedom.is_infinite() {
             use super::multivariate_normal::density_normalization_and_exponential;
diff --git a/src/function/beta.rs b/src/function/beta.rs
index 794a8504..188c09fc 100644
--- a/src/function/beta.rs
+++ b/src/function/beta.rs
@@ -1,12 +1,37 @@
 //! Provides the [beta](https://en.wikipedia.org/wiki/Beta_function) and related
 //! function
 
-use crate::error::StatsError;
 use crate::function::gamma;
 use crate::prec;
-use crate::Result;
 use std::f64;
 
+/// Represents the errors that can occur when computing the natural logarithm
+/// of the beta function or the regularized lower incomplete beta function.
+#[derive(Copy, Clone, PartialEq, Eq, Debug, Hash)]
+#[non_exhaustive]
+pub enum BetaFuncError {
+    /// `a` is zero or less than zero.
+    ANotGreaterThanZero,
+
+    /// `b` is zero or less than zero.
+    BNotGreaterThanZero,
+
+    /// `x` is not in `[0, 1]`.
+    XOutOfRange,
+}
+
+impl std::fmt::Display for BetaFuncError {
+    fn fmt(&self, f: &mut std::fmt::Formatter) -> std::fmt::Result {
+        match self {
+            BetaFuncError::ANotGreaterThanZero => write!(f, "a is zero or less than zero"),
+            BetaFuncError::BNotGreaterThanZero => write!(f, "b is zero or less than zero"),
+            BetaFuncError::XOutOfRange => write!(f, "x is not in [0, 1]"),
+        }
+    }
+}
+
+impl std::error::Error for BetaFuncError {}
+
 /// Computes the natural logarithm
 /// of the beta function
 /// where `a` is the first beta parameter
@@ -29,11 +54,11 @@ pub fn ln_beta(a: f64, b: f64) -> f64 {
 /// # Errors
 ///
 /// if `a <= 0.0` or `b <= 0.0`
-pub fn checked_ln_beta(a: f64, b: f64) -> Result<f64> {
+pub fn checked_ln_beta(a: f64, b: f64) -> Result<f64, BetaFuncError> {
     if a <= 0.0 {
-        Err(StatsError::ArgMustBePositive("a"))
+        Err(BetaFuncError::ANotGreaterThanZero)
     } else if b <= 0.0 {
-        Err(StatsError::ArgMustBePositive("b"))
+        Err(BetaFuncError::BNotGreaterThanZero)
     } else {
         Ok(gamma::ln_gamma(a) + gamma::ln_gamma(b) - gamma::ln_gamma(a + b))
     }
@@ -59,7 +84,7 @@ pub fn beta(a: f64, b: f64) -> f64 {
 /// # Errors
 ///
 /// if `a <= 0.0` or `b <= 0.0`
-pub fn checked_beta(a: f64, b: f64) -> Result<f64> {
+pub fn checked_beta(a: f64, b: f64) -> Result<f64, BetaFuncError> {
     checked_ln_beta(a, b).map(|x| x.exp())
 }
 
@@ -83,7 +108,7 @@ pub fn beta_inc(a: f64, b: f64, x: f64) -> f64 {
 /// # Errors
 ///
 /// If `a <= 0.0`, `b <= 0.0`, `x < 0.0`, or `x > 1.0`
-pub fn checked_beta_inc(a: f64, b: f64, x: f64) -> Result<f64> {
+pub fn checked_beta_inc(a: f64, b: f64, x: f64) -> Result<f64, BetaFuncError> {
     checked_beta_reg(a, b, x).and_then(|x| checked_beta(a, b).map(|y| x * y))
 }
 
@@ -109,96 +134,100 @@ pub fn beta_reg(a: f64, b: f64, x: f64) -> f64 {
 /// # Errors
 ///
 /// if `a <= 0.0`, `b <= 0.0`, `x < 0.0`, or `x > 1.0`
-pub fn checked_beta_reg(a: f64, b: f64, x: f64) -> Result<f64> {
+pub fn checked_beta_reg(a: f64, b: f64, x: f64) -> Result<f64, BetaFuncError> {
     if a <= 0.0 {
-        Err(StatsError::ArgMustBePositive("a"))
-    } else if b <= 0.0 {
-        Err(StatsError::ArgMustBePositive("b"))
-    } else if !(0.0..=1.0).contains(&x) {
-        Err(StatsError::ArgIntervalIncl("x", 0.0, 1.0))
+        return Err(BetaFuncError::ANotGreaterThanZero);
+    }
+
+    if b <= 0.0 {
+        return Err(BetaFuncError::BNotGreaterThanZero);
+    }
+
+    if !(0.0..=1.0).contains(&x) {
+        return Err(BetaFuncError::XOutOfRange);
+    }
+
+    let bt = if x == 0.0 || ulps_eq!(x, 1.0) {
+        0.0
     } else {
-        let bt = if x == 0.0 || ulps_eq!(x, 1.0) {
-            0.0
-        } else {
-            (gamma::ln_gamma(a + b) - gamma::ln_gamma(a) - gamma::ln_gamma(b)
-                + a * x.ln()
-                + b * (1.0 - x).ln())
-            .exp()
-        };
-        let symm_transform = x >= (a + 1.0) / (a + b + 2.0);
-        let eps = prec::F64_PREC;
-        let fpmin = f64::MIN_POSITIVE / eps;
-
-        let mut a = a;
-        let mut b = b;
-        let mut x = x;
-        if symm_transform {
-            let swap = a;
-            x = 1.0 - x;
-            a = b;
-            b = swap;
-        }
+        (gamma::ln_gamma(a + b) - gamma::ln_gamma(a) - gamma::ln_gamma(b)
+            + a * x.ln()
+            + b * (1.0 - x).ln())
+        .exp()
+    };
+    let symm_transform = x >= (a + 1.0) / (a + b + 2.0);
+    let eps = prec::F64_PREC;
+    let fpmin = f64::MIN_POSITIVE / eps;
+
+    let mut a = a;
+    let mut b = b;
+    let mut x = x;
+    if symm_transform {
+        let swap = a;
+        x = 1.0 - x;
+        a = b;
+        b = swap;
+    }
+
+    let qab = a + b;
+    let qap = a + 1.0;
+    let qam = a - 1.0;
+    let mut c = 1.0;
+    let mut d = 1.0 - qab * x / qap;
+
+    if d.abs() < fpmin {
+        d = fpmin;
+    }
+    d = 1.0 / d;
+    let mut h = d;
 
-        let qab = a + b;
-        let qap = a + 1.0;
-        let qam = a - 1.0;
-        let mut c = 1.0;
-        let mut d = 1.0 - qab * x / qap;
+    for m in 1..141 {
+        let m = f64::from(m);
+        let m2 = m * 2.0;
+        let mut aa = m * (b - m) * x / ((qam + m2) * (a + m2));
+        d = 1.0 + aa * d;
 
         if d.abs() < fpmin {
             d = fpmin;
         }
-        d = 1.0 / d;
-        let mut h = d;
-
-        for m in 1..141 {
-            let m = f64::from(m);
-            let m2 = m * 2.0;
-            let mut aa = m * (b - m) * x / ((qam + m2) * (a + m2));
-            d = 1.0 + aa * d;
-
-            if d.abs() < fpmin {
-                d = fpmin;
-            }
 
-            c = 1.0 + aa / c;
-            if c.abs() < fpmin {
-                c = fpmin;
-            }
+        c = 1.0 + aa / c;
+        if c.abs() < fpmin {
+            c = fpmin;
+        }
 
-            d = 1.0 / d;
-            h = h * d * c;
-            aa = -(a + m) * (qab + m) * x / ((a + m2) * (qap + m2));
-            d = 1.0 + aa * d;
+        d = 1.0 / d;
+        h = h * d * c;
+        aa = -(a + m) * (qab + m) * x / ((a + m2) * (qap + m2));
+        d = 1.0 + aa * d;
 
-            if d.abs() < fpmin {
-                d = fpmin;
-            }
+        if d.abs() < fpmin {
+            d = fpmin;
+        }
 
-            c = 1.0 + aa / c;
+        c = 1.0 + aa / c;
 
-            if c.abs() < fpmin {
-                c = fpmin;
-            }
+        if c.abs() < fpmin {
+            c = fpmin;
+        }
 
-            d = 1.0 / d;
-            let del = d * c;
-            h *= del;
+        d = 1.0 / d;
+        let del = d * c;
+        h *= del;
 
-            if (del - 1.0).abs() <= eps {
-                return if symm_transform {
-                    Ok(1.0 - bt * h / a)
-                } else {
-                    Ok(bt * h / a)
-                };
-            }
+        if (del - 1.0).abs() <= eps {
+            return if symm_transform {
+                Ok(1.0 - bt * h / a)
+            } else {
+                Ok(bt * h / a)
+            };
         }
+    }
 
-        if symm_transform {
-            Ok(1.0 - bt * h / a)
-        } else {
-            Ok(bt * h / a)
-        }
+    if symm_transform {
+        Ok(1.0 - bt * h / a)
+    } else {
+        Ok(bt * h / a)
     }
 }
 
@@ -396,6 +425,8 @@ pub fn inv_beta_reg(mut a: f64, mut b: f64, mut x: f64) -> f64 {
 #[rustfmt::skip]
 #[cfg(test)]
 mod tests {
+    use super::*;
+
     #[test]
     fn test_ln_beta() {
         assert_almost_eq!(super::ln_beta(0.5, 0.5), 1.144729885849400174144, 1e-15);
@@ -597,4 +628,10 @@ mod tests {
     fn test_checked_beta_reg_x_gt_1() {
         assert!(super::checked_beta_reg(1.0, 1.0, 2.0).is_err());
     }
+
+    #[test]
+    fn test_error_is_sync_send() {
+        fn assert_sync_send<T: Sync + Send>() {}
+        assert_sync_send::<BetaFuncError>();
+    }
 }
diff --git a/src/function/exponential.rs b/src/function/exponential.rs
index 55280d7c..4d786f4d 100644
--- a/src/function/exponential.rs
+++ b/src/function/exponential.rs
@@ -1,14 +1,12 @@
 //! Provides functions related to exponential calculations
 
-use crate::{consts, Result, StatsError};
+use crate::consts;
 
 /// Computes the generalized Exponential Integral function
 /// where `x` is the argument and `n` is the integer power of the
 /// denominator term.
 ///
-/// # Errors
-///
-/// Returns an error if `x < 0.0` or the computation could not
+/// Returns `None` if `x < 0.0` or the computation could not
 /// converge after 100 iterations
 ///
 /// # Remarks
@@ -26,7 +24,7 @@ use crate::{consts, Result, StatsError};
 /// The continued fraction approach is used for `x > 1.0` while the taylor
 /// series expansions is used for `0.0 < x <= 1`.
 // TODO: Add examples
-pub fn integral(x: f64, n: u64) -> Result<f64> {
+pub fn integral(x: f64, n: u64) -> Option<f64> {
     let eps = 0.00000000000000001;
     let max_iter = 100;
     let nf64 = n as f64;
@@ -34,10 +32,10 @@ pub fn integral(x: f64, n: u64) -> Result<f64> {
 
     // special cases
     if n == 0 {
-        return Ok((-1.0 * x).exp() / x);
+        return Some((-1.0 * x).exp() / x);
     }
     if x == 0.0 {
-        return Ok(1.0 / (nf64 - 1.0));
+        return Some(1.0 / (nf64 - 1.0));
     }
 
     if x > 1.0 {
@@ -53,10 +51,10 @@ pub fn integral(x: f64, n: u64) -> Result<f64> {
             let del = c * d;
             h *= del;
             if (del - 1.0).abs() < eps {
-                return Ok(h * (-x).exp());
+                return Some(h * (-x).exp());
             }
         }
-        Err(StatsError::ComputationFailedToConverge)
+        None
     } else {
         let mut factorial = 1.0;
         let mut result = if n - 1 != 0 {
@@ -77,10 +75,10 @@ pub fn integral(x: f64, n: u64) -> Result<f64> {
             };
             result += del;
             if del.abs() < result.abs() * eps {
-                return Ok(result);
+                return Some(result);
             }
         }
-        Err(StatsError::ComputationFailedToConverge)
+        None
     }
 }
 
diff --git a/src/function/factorial.rs b/src/function/factorial.rs
index 77bcaa5b..08287aab 100644
--- a/src/function/factorial.rs
+++ b/src/function/factorial.rs
@@ -1,9 +1,7 @@
 //! Provides functions related to factorial calculations (e.g. binomial
 //! coefficient, factorial, multinomial)
 
-use crate::error::StatsError;
 use crate::function::gamma;
-use crate::Result;
 
 /// The maximum factorial representable
 /// by a 64-bit floating point without
@@ -74,17 +72,16 @@ pub fn multinomial(n: u64, ni: &[u64]) -> f64 {
 
 /// Computes the multinomial coefficient: `n choose n1, n2, n3, ...`
 ///
-/// # Errors
-///
-/// If the elements in `ni` do not sum to `n`
-pub fn checked_multinomial(n: u64, ni: &[u64]) -> Result<f64> {
+/// Returns `None` if the elements in `ni` do not sum to `n`.
+pub fn checked_multinomial(n: u64, ni: &[u64]) -> Option<f64> {
     let (sum, ret) = ni.iter().fold((0, ln_factorial(n)), |acc, &x| {
         (acc.0 + x, acc.1 - ln_factorial(x))
     });
-    if sum != n {
-        Err(StatsError::ContainerExpectedSumVar("ni", "n"))
+
+    if sum == n {
+        Some((0.5 + ret.exp()).floor())
     } else {
-        Ok((0.5 + ret.exp()).floor())
+        None
     }
 }
 
@@ -179,6 +176,6 @@ mod tests {
 
     #[test]
     fn test_checked_multinomial_bad_ni() {
-        assert!(checked_multinomial(1, &[1, 1]).is_err());
+        assert!(checked_multinomial(1, &[1, 1]).is_none());
     }
 }
diff --git a/src/function/gamma.rs b/src/function/gamma.rs
index ced63871..ab15fc09 100644
--- a/src/function/gamma.rs
+++ b/src/function/gamma.rs
@@ -2,11 +2,32 @@
 //! related functions
 
 use crate::consts;
-use crate::error::StatsError;
 use crate::prec;
-use crate::Result;
 use std::f64;
 
+/// Represents the errors that can occur when computing any of the incomplete
+/// gamma functions.
+#[derive(Copy, Clone, PartialEq, Eq, Debug, Hash)]
+#[non_exhaustive]
+pub enum GammaFuncError {
+    /// `a` is infinite, zero or less than zero.
+    AInvalid,
+
+    /// `x` is infinite, zero or less than zero.
+    XInvalid,
+}
+
+impl std::fmt::Display for GammaFuncError {
+    fn fmt(&self, f: &mut std::fmt::Formatter) -> std::fmt::Result {
+        match self {
+            GammaFuncError::AInvalid => write!(f, "a is infinite, zero or less than zero"),
+            GammaFuncError::XInvalid => write!(f, "x is infinite, zero or less than zero"),
+        }
+    }
+}
+
+impl std::error::Error for GammaFuncError {}
+
 /// Auxiliary variable when evaluating the `gamma_ln` function
 const GAMMA_R: f64 = 10.900511;
 
@@ -104,7 +125,7 @@ pub fn gamma_ui(a: f64, x: f64) -> f64 {
 /// # Errors
 ///
 /// if `a` or `x` are not in `(0, +inf)`
-pub fn checked_gamma_ui(a: f64, x: f64) -> Result<f64> {
+pub fn checked_gamma_ui(a: f64, x: f64) -> Result<f64, GammaFuncError> {
     checked_gamma_ur(a, x).map(|x| x * gamma(a))
 }
 
@@ -130,7 +151,7 @@ pub fn gamma_li(a: f64, x: f64) -> f64 {
 /// # Errors
 ///
 /// if `a` or `x` are not in `(0, +inf)`
-pub fn checked_gamma_li(a: f64, x: f64) -> Result<f64> {
+pub fn checked_gamma_li(a: f64, x: f64) -> Result<f64, GammaFuncError> {
     checked_gamma_lr(a, x).map(|x| x * gamma(a))
 }
 
@@ -162,15 +183,15 @@ pub fn gamma_ur(a: f64, x: f64) -> f64 {
 /// # Errors
 ///
 /// if `a` or `x` are not in `(0, +inf)`
-pub fn checked_gamma_ur(a: f64, x: f64) -> Result<f64> {
+pub fn checked_gamma_ur(a: f64, x: f64) -> Result<f64, GammaFuncError> {
     if a.is_nan() || x.is_nan() {
         return Ok(f64::NAN);
     }
     if a <= 0.0 || a == f64::INFINITY {
-        return Err(StatsError::ArgIntervalExcl("a", 0.0, f64::INFINITY));
+        return Err(GammaFuncError::AInvalid);
     }
     if x <= 0.0 || x == f64::INFINITY {
-        return Err(StatsError::ArgIntervalExcl("x", 0.0, f64::INFINITY));
+        return Err(GammaFuncError::XInvalid);
     }
 
     let eps = 0.000000000000001;
@@ -256,15 +277,15 @@ pub fn gamma_lr(a: f64, x: f64) -> f64 {
 /// # Errors
 ///
 /// if `a` or `x` are not in `(0, +inf)`
-pub fn checked_gamma_lr(a: f64, x: f64) -> Result<f64> {
+pub fn checked_gamma_lr(a: f64, x: f64) -> Result<f64, GammaFuncError> {
     if a.is_nan() || x.is_nan() {
         return Ok(f64::NAN);
     }
     if a <= 0.0 || a == f64::INFINITY {
-        return Err(StatsError::ArgIntervalExcl("a", 0.0, f64::INFINITY));
+        return Err(GammaFuncError::AInvalid);
     }
     if x <= 0.0 || x == f64::INFINITY {
-        return Err(StatsError::ArgIntervalExcl("x", 0.0, f64::INFINITY));
+        return Err(GammaFuncError::XInvalid);
     }
 
     let eps = 0.000000000000001;
@@ -423,7 +444,9 @@ fn signum(x: f64) -> f64 {
 #[rustfmt::skip]
 #[cfg(test)]
 mod tests {
-    use std::f64::{self, consts};
+    use super::*;
+
+    use std::f64::consts;
 
     #[test]
     fn test_gamma() {
@@ -804,4 +827,10 @@ mod tests {
         assert_almost_eq!(super::inv_digamma(1.6110931485817511237336268416044190359814435699427405), 5.5, 1e-14);
         assert_almost_eq!(super::inv_digamma(2.2622143570941481235561593642219403924532310597356171), 10.1, 1e-13);
     }
+
+    #[test]
+    fn test_error_is_sync_send() {
+        fn assert_sync_send<T: Sync + Send>() {}
+        assert_sync_send::<GammaFuncError>();
+    }
 }
diff --git a/src/function/logistic.rs b/src/function/logistic.rs
index 8adc8b79..57ccdd3e 100644
--- a/src/function/logistic.rs
+++ b/src/function/logistic.rs
@@ -1,9 +1,6 @@
 //! Provides the [logistic](http://en.wikipedia.org/wiki/Logistic_function) and
 //! related functions
 
-use crate::error::StatsError;
-use crate::Result;
-
 /// Computes the logistic function
 pub fn logistic(p: f64) -> f64 {
     1.0 / ((-p).exp() + 1.0)
@@ -18,16 +15,12 @@ pub fn logit(p: f64) -> f64 {
     checked_logit(p).unwrap()
 }
 
-/// Computes the logit function
-///
-/// # Errors
-///
-/// If `p < 0.0` or `p > 1.0`
-pub fn checked_logit(p: f64) -> Result<f64> {
-    if !(0.0..=1.0).contains(&p) {
-        Err(StatsError::ArgIntervalIncl("p", 0.0, 1.0))
+/// Computes the logit function, returning `None` if `p < 0.0` or `p > 1.0`.
+pub fn checked_logit(p: f64) -> Option<f64> {
+    if (0.0..=1.0).contains(&p) {
+        Some((p / (1.0 - p)).ln())
     } else {
-        Ok((p / (1.0 - p)).ln())
+        None
     }
 }
 
@@ -76,11 +69,11 @@ mod tests {
 
     #[test]
     fn test_checked_logit_p_lt_0() {
-        assert!(super::checked_logit(-1.0).is_err());
+        assert!(super::checked_logit(-1.0).is_none());
     }
 
     #[test]
     fn test_checked_logit_p_gt_1() {
-        assert!(super::checked_logit(2.0).is_err());
+        assert!(super::checked_logit(2.0).is_none());
     }
 }
diff --git a/src/lib.rs b/src/lib.rs
index 7ca3157c..634a2354 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -74,11 +74,3 @@ pub mod generate;
 pub mod prec;
 pub mod statistics;
 pub mod stats_tests;
-
-mod error;
-
-pub use crate::error::StatsError;
-
-/// Result type for the statrs library package that returns
-/// either a result type `T` or a `StatsError`
-pub type Result<T> = std::result::Result<T, StatsError>;
diff --git a/src/statistics/iter_statistics.rs b/src/statistics/iter_statistics.rs
index 97b7fd17..2dea1a9c 100644
--- a/src/statistics/iter_statistics.rs
+++ b/src/statistics/iter_statistics.rs
@@ -1,4 +1,3 @@
-use crate::error::StatsError;
 use crate::statistics::*;
 use std::borrow::Borrow;
 use std::f64;
@@ -175,7 +174,7 @@ where
         for x in self {
             let borrow = *x.borrow();
             let borrow2 = match iter.next() {
-                None => panic!("{}", StatsError::ContainersMustBeSameLength),
+                None => panic!("Iterators must have the same length"),
                 Some(x) => *x.borrow(),
             };
             let old_mean2 = mean2;
@@ -185,7 +184,7 @@ where
             comoment += (borrow - mean1) * (borrow2 - old_mean2);
         }
         if iter.next().is_some() {
-            panic!("{}", StatsError::ContainersMustBeSameLength);
+            panic!("Iterators must have the same length");
         }
 
         if n > 1.0 {
@@ -205,7 +204,7 @@ where
         for x in self {
             let borrow = *x.borrow();
             let borrow2 = match iter.next() {
-                None => panic!("{}", StatsError::ContainersMustBeSameLength),
+                None => panic!("Iterators must have the same length"),
                 Some(x) => *x.borrow(),
             };
             let old_mean2 = mean2;
@@ -215,7 +214,7 @@ where
             comoment += (borrow - mean1) * (borrow2 - old_mean2);
         }
         if iter.next().is_some() {
-            panic!("{}", StatsError::ContainersMustBeSameLength)
+            panic!("Iterators must have the same length")
         }
         if n > 0.0 {
             comoment / n


Test patch: