diff --git a/src/distribution/bernoulli.rs b/src/distribution/bernoulli.rs
index 94af5450..12466f42 100644
--- a/src/distribution/bernoulli.rs
+++ b/src/distribution/bernoulli.rs
@@ -10,13 +10,14 @@ use crate::statistics::*;
 /// # Examples
 ///
 /// ```
-/// use statrs::distribution::{Bernoulli, Discrete};
-/// use statrs::statistics::Distribution;
+/// use statrs::distribution::{Bernoulli, BinomialError, Discrete};
+/// use statrs::statistics::*;
 ///
-/// let n = Bernoulli::new(0.5).unwrap();
-/// assert_eq!(n.mean().unwrap(), 0.5);
+/// let n = Bernoulli::new(0.5)?;
+/// assert_eq!(n.mean(), 0.5);
 /// assert_eq!(n.pmf(0), 0.5);
 /// assert_eq!(n.pmf(1), 0.5);
+/// # Ok::<(), BinomialError>(())
 /// ```
 #[derive(Copy, Clone, PartialEq, Debug)]
 pub struct Bernoulli {
@@ -53,10 +54,11 @@ impl Bernoulli {
     /// # Examples
     ///
     /// ```
-    /// use statrs::distribution::Bernoulli;
+    /// use statrs::distribution::{Bernoulli, BinomialError};
     ///
-    /// let n = Bernoulli::new(0.5).unwrap();
+    /// let n = Bernoulli::new(0.5)?;
     /// assert_eq!(n.p(), 0.5);
+    /// # Ok::<(), BinomialError>(())
     /// ```
     pub fn p(&self) -> f64 {
         self.b.p()
@@ -68,10 +70,11 @@ impl Bernoulli {
     /// # Examples
     ///
     /// ```
-    /// use statrs::distribution::Bernoulli;
+    /// use statrs::distribution::{Bernoulli, BinomialError};
     ///
-    /// let n = Bernoulli::new(0.5).unwrap();
+    /// let n = Bernoulli::new(0.5)?;
     /// assert_eq!(n.n(), 1);
+    /// # Ok::<(), BinomialError>(())
     /// ```
     pub fn n(&self) -> u64 {
         1
diff --git a/src/distribution/beta.rs b/src/distribution/beta.rs
index 2e56ffee..5d030100 100644
--- a/src/distribution/beta.rs
+++ b/src/distribution/beta.rs
@@ -8,13 +8,14 @@ use crate::statistics::*;
 /// # Examples
 ///
 /// ```
-/// use statrs::distribution::{Beta, Continuous};
+/// use statrs::distribution::{Beta, Continuous, BetaError};
 /// use statrs::statistics::*;
 /// use statrs::prec;
 ///
-/// let n = Beta::new(2.0, 2.0).unwrap();
-/// assert_eq!(n.mean().unwrap(), 0.5);
+/// let n = Beta::new(2.0, 2.0)?;
+/// assert_eq!(n.mean(), 0.5);
 /// assert!(prec::almost_eq(n.pdf(0.5), 1.5, 1e-14));
+/// # Ok::<(), BetaError>(())
 /// ```
 #[derive(Copy, Clone, PartialEq, Debug)]
 pub struct Beta {
@@ -82,10 +83,11 @@ impl Beta {
     /// # Examples
     ///
     /// ```
-    /// use statrs::distribution::Beta;
+    /// use statrs::distribution::{Beta, BetaError};
     ///
-    /// let n = Beta::new(1.0, 2.0).unwrap();
+    /// let n = Beta::new(1.0, 2.0)?;
     /// assert_eq!(n.shape_a(), 1.0);
+    /// # Ok::<(), BetaError>(())
     /// ```
     pub fn shape_a(&self) -> f64 {
         self.shape_a
@@ -96,10 +98,11 @@ impl Beta {
     /// # Examples
     ///
     /// ```
-    /// use statrs::distribution::Beta;
+    /// use statrs::distribution::{Beta, BetaError};
     ///
-    /// let n = Beta::new(1.0, 2.0).unwrap();
+    /// let n = Beta::new(1.0, 2.0)?;
     /// assert_eq!(n.shape_b(), 2.0);
+    /// # Ok::<(), BetaError>(())
     /// ```
     pub fn shape_b(&self) -> f64 {
         self.shape_b
diff --git a/src/distribution/binomial.rs b/src/distribution/binomial.rs
index b7abca78..7d0553bf 100644
--- a/src/distribution/binomial.rs
+++ b/src/distribution/binomial.rs
@@ -10,13 +10,14 @@ use std::f64;
 /// # Examples
 ///
 /// ```
-/// use statrs::distribution::{Binomial, Discrete};
-/// use statrs::statistics::Distribution;
+/// use statrs::distribution::{Binomial, BinomialError, Discrete};
+/// use statrs::statistics::*;
 ///
-/// let n = Binomial::new(0.5, 5).unwrap();
-/// assert_eq!(n.mean().unwrap(), 2.5);
+/// let n = Binomial::new(0.5, 5)?;
+/// assert_eq!(n.mean(), 2.5);
 /// assert_eq!(n.pmf(0), 0.03125);
 /// assert_eq!(n.pmf(3), 0.3125);
+/// # Ok::<(), BinomialError>(())
 /// ```
 #[derive(Copy, Clone, PartialEq, Debug)]
 pub struct Binomial {
@@ -78,10 +79,11 @@ impl Binomial {
     /// # Examples
     ///
     /// ```
-    /// use statrs::distribution::Binomial;
+    /// use statrs::distribution::{Binomial, BinomialError};
     ///
-    /// let n = Binomial::new(0.5, 5).unwrap();
+    /// let n = Binomial::new(0.5, 5)?;
     /// assert_eq!(n.p(), 0.5);
+    /// # Ok::<(), BinomialError>(())
     /// ```
     pub fn p(&self) -> f64 {
         self.p
@@ -93,10 +95,11 @@ impl Binomial {
     /// # Examples
     ///
     /// ```
-    /// use statrs::distribution::Binomial;
+    /// use statrs::distribution::{Binomial, BinomialError};
     ///
-    /// let n = Binomial::new(0.5, 5).unwrap();
+    /// let n = Binomial::new(0.5, 5)?;
     /// assert_eq!(n.n(), 5);
+    /// # Ok::<(), BinomialError>(())
     /// ```
     pub fn n(&self) -> u64 {
         self.n
diff --git a/src/distribution/categorical.rs b/src/distribution/categorical.rs
index dcec087b..44371699 100644
--- a/src/distribution/categorical.rs
+++ b/src/distribution/categorical.rs
@@ -10,13 +10,14 @@ use std::f64;
 /// # Examples
 ///
 /// ```
-/// use statrs::distribution::{Categorical, Discrete};
-/// use statrs::statistics::Distribution;
+/// use statrs::distribution::{Categorical, CategoricalError, Discrete};
+/// use statrs::statistics::*;
 /// use statrs::prec;
 ///
-/// let n = Categorical::new(&[0.0, 1.0, 2.0]).unwrap();
-/// assert!(prec::almost_eq(n.mean().unwrap(), 5.0 / 3.0, 1e-15));
+/// let n = Categorical::new(&[0.0, 1.0, 2.0])?;
+/// assert!(prec::almost_eq(n.mean(), 5.0 / 3.0, 1e-15));
 /// assert_eq!(n.pmf(1), 1.0 / 3.0);
+/// # Ok::<(), CategoricalError>(())
 /// ```
 #[derive(Clone, PartialEq, Debug)]
 pub struct Categorical {
diff --git a/src/distribution/cauchy.rs b/src/distribution/cauchy.rs
index 6f1b7c7c..fbde1dd5 100644
--- a/src/distribution/cauchy.rs
+++ b/src/distribution/cauchy.rs
@@ -8,12 +8,13 @@ use std::f64;
 /// # Examples
 ///
 /// ```
-/// use statrs::distribution::{Cauchy, Continuous};
-/// use statrs::statistics::Mode;
+/// use statrs::distribution::{Cauchy, Continuous, CauchyError};
+/// use statrs::statistics::*;
 ///
-/// let n = Cauchy::new(0.0, 1.0).unwrap();
+/// let n = Cauchy::new(0.0, 1.0)?;
 /// assert_eq!(n.mode().unwrap(), 0.0);
 /// assert_eq!(n.pdf(1.0), 0.1591549430918953357689);
+/// # Ok::<(), CauchyError>(())
 /// ```
 #[derive(Copy, Clone, PartialEq, Debug)]
 pub struct Cauchy {
@@ -80,10 +81,11 @@ impl Cauchy {
     /// # Examples
     ///
     /// ```
-    /// use statrs::distribution::Cauchy;
+    /// use statrs::distribution::{Cauchy, CauchyError};
     ///
-    /// let n = Cauchy::new(0.0, 1.0).unwrap();
+    /// let n = Cauchy::new(0.0, 1.0)?;
     /// assert_eq!(n.location(), 0.0);
+    /// # Ok::<(), CauchyError>(())
     /// ```
     pub fn location(&self) -> f64 {
         self.location
@@ -94,10 +96,11 @@ impl Cauchy {
     /// # Examples
     ///
     /// ```
-    /// use statrs::distribution::Cauchy;
+    /// use statrs::distribution::{Cauchy, CauchyError};
     ///
-    /// let n = Cauchy::new(0.0, 1.0).unwrap();
+    /// let n = Cauchy::new(0.0, 1.0)?;
     /// assert_eq!(n.scale(), 1.0);
+    /// # Ok::<(), CauchyError>(())
     /// ```
     pub fn scale(&self) -> f64 {
         self.scale
diff --git a/src/distribution/chi.rs b/src/distribution/chi.rs
index 8b5a9576..11b6a2c5 100644
--- a/src/distribution/chi.rs
+++ b/src/distribution/chi.rs
@@ -9,13 +9,14 @@ use std::f64;
 /// # Examples
 ///
 /// ```
-/// use statrs::distribution::{Chi, Continuous};
-/// use statrs::statistics::Distribution;
+/// use statrs::distribution::{Chi, Continuous, ChiError};
+/// use statrs::statistics::*;
 /// use statrs::prec;
 ///
-/// let n = Chi::new(2.0).unwrap();
+/// let n = Chi::new(2.0)?;
 /// assert!(prec::almost_eq(n.mean().unwrap(), 1.25331413731550025121, 1e-14));
 /// assert!(prec::almost_eq(n.pdf(1.0), 0.60653065971263342360, 1e-15));
+/// # Ok::<(), ChiError>(())
 /// ```
 #[derive(Copy, Clone, PartialEq, Debug)]
 pub struct Chi {
@@ -77,10 +78,12 @@ impl Chi {
     /// # Examples
     ///
     /// ```
-    /// use statrs::distribution::Chi;
+    /// use statrs::distribution::{Chi, ChiError};
+    ///
     ///
-    /// let n = Chi::new(2.0).unwrap();
+    /// let n = Chi::new(2.0)?;
     /// assert_eq!(n.freedom(), 2.0);
+    /// # Ok::<(), ChiError>(())
     /// ```
     pub fn freedom(&self) -> f64 {
         self.freedom
diff --git a/src/distribution/chi_squared.rs b/src/distribution/chi_squared.rs
index de7b6534..8cfb3275 100644
--- a/src/distribution/chi_squared.rs
+++ b/src/distribution/chi_squared.rs
@@ -11,13 +11,14 @@ use std::f64;
 /// # Examples
 ///
 /// ```
-/// use statrs::distribution::{ChiSquared, Continuous};
-/// use statrs::statistics::Distribution;
+/// use statrs::distribution::{ChiSquared, Continuous, GammaError};
+/// use statrs::statistics::*;
 /// use statrs::prec;
 ///
-/// let n = ChiSquared::new(3.0).unwrap();
-/// assert_eq!(n.mean().unwrap(), 3.0);
+/// let n = ChiSquared::new(3.0)?;
+/// assert_eq!(n.mean(), 3.0);
 /// assert!(prec::almost_eq(n.pdf(4.0), 0.107981933026376103901, 1e-15));
+/// # Ok::<(), GammaError>(())
 /// ```
 #[derive(Copy, Clone, PartialEq, Debug)]
 pub struct ChiSquared {
@@ -56,10 +57,11 @@ impl ChiSquared {
     /// # Examples
     ///
     /// ```
-    /// use statrs::distribution::ChiSquared;
+    /// use statrs::distribution::{ChiSquared, GammaError};
     ///
-    /// let n = ChiSquared::new(3.0).unwrap();
+    /// let n = ChiSquared::new(3.0)?;
     /// assert_eq!(n.freedom(), 3.0);
+    /// # Ok::<(), GammaError>(())
     /// ```
     pub fn freedom(&self) -> f64 {
         self.freedom
@@ -70,10 +72,11 @@ impl ChiSquared {
     /// # Examples
     ///
     /// ```
-    /// use statrs::distribution::ChiSquared;
+    /// use statrs::distribution::{ChiSquared, GammaError};
     ///
-    /// let n = ChiSquared::new(3.0).unwrap();
+    /// let n = ChiSquared::new(3.0)?;
     /// assert_eq!(n.shape(), 3.0 / 2.0);
+    /// # Ok::<(), GammaError>(())
     /// ```
     pub fn shape(&self) -> f64 {
         self.g.shape()
@@ -84,10 +87,11 @@ impl ChiSquared {
     /// # Examples
     ///
     /// ```
-    /// use statrs::distribution::ChiSquared;
+    /// use statrs::distribution::{ChiSquared, GammaError};
     ///
-    /// let n = ChiSquared::new(3.0).unwrap();
+    /// let n = ChiSquared::new(3.0)?;
     /// assert_eq!(n.rate(), 0.5);
+    /// # Ok::<(), GammaError>(())
     /// ```
     pub fn rate(&self) -> f64 {
         self.g.rate()
diff --git a/src/distribution/dirac.rs b/src/distribution/dirac.rs
index 2b4b5401..acd5c0ba 100644
--- a/src/distribution/dirac.rs
+++ b/src/distribution/dirac.rs
@@ -7,11 +7,12 @@ use crate::statistics::*;
 /// # Examples
 ///
 /// ```
-/// use statrs::distribution::{Dirac, Continuous};
-/// use statrs::statistics::Distribution;
+/// use statrs::distribution::{Dirac, Continuous, DiracError};
+/// use statrs::statistics::*;
 ///
-/// let n = Dirac::new(3.0).unwrap();
-/// assert_eq!(n.mean().unwrap(), 3.0);
+/// let n = Dirac::new(3.0)?;
+/// assert_eq!(n.mean(), 3.0);
+/// # Ok::<(), DiracError>(())
 /// ```
 #[derive(Debug, Copy, Clone, PartialEq)]
 pub struct Dirac(f64);
diff --git a/src/distribution/dirichlet.rs b/src/distribution/dirichlet.rs
index 5e08d0a4..7fe40282 100644
--- a/src/distribution/dirichlet.rs
+++ b/src/distribution/dirichlet.rs
@@ -12,14 +12,14 @@ use std::f64;
 /// # Examples
 ///
 /// ```
-/// use statrs::distribution::{Dirichlet, Continuous};
-/// use statrs::statistics::Distribution;
+/// use statrs::distribution::{Dirichlet, Continuous, DirichletError};
+/// use statrs::statistics::*;
 /// use nalgebra::DVector;
-/// use statrs::statistics::MeanN;
 ///
-/// let n = Dirichlet::new(vec![1.0, 2.0, 3.0]).unwrap();
+/// let n = Dirichlet::new(vec![1.0, 2.0, 3.0])?;
 /// assert_eq!(n.mean().unwrap(), DVector::from_vec(vec![1.0 / 6.0, 1.0 / 3.0, 0.5]));
 /// assert_eq!(n.pdf(&DVector::from_vec(vec![0.33333, 0.33333, 0.33333])), 2.222155556222205);
+/// # Ok::<(), DirichletError>(())
 /// ```
 #[derive(Clone, PartialEq, Debug)]
 pub struct Dirichlet<D>
@@ -138,11 +138,12 @@ where
     /// # Examples
     ///
     /// ```
-    /// use statrs::distribution::Dirichlet;
+    /// use statrs::distribution::{Dirichlet, DirichletError};
     /// use nalgebra::DVector;
     ///
-    /// let n = Dirichlet::new(vec![1.0, 2.0, 3.0]).unwrap();
+    /// let n = Dirichlet::new(vec![1.0, 2.0, 3.0])?;
     /// assert_eq!(n.alpha(), &DVector::from_vec(vec![1.0, 2.0, 3.0]));
+    /// # Ok::<(), DirichletError>(())
     /// ```
     pub fn alpha(&self) -> &nalgebra::OVector<f64, D> {
         &self.alpha
diff --git a/src/distribution/discrete_uniform.rs b/src/distribution/discrete_uniform.rs
index 62c66a95..9b650da9 100644
--- a/src/distribution/discrete_uniform.rs
+++ b/src/distribution/discrete_uniform.rs
@@ -8,12 +8,13 @@ use crate::statistics::*;
 /// # Examples
 ///
 /// ```
-/// use statrs::distribution::{DiscreteUniform, Discrete};
-/// use statrs::statistics::Distribution;
+/// use statrs::distribution::{DiscreteUniform, Discrete, DiscreteUniformError};
+/// use statrs::statistics::*;
 ///
-/// let n = DiscreteUniform::new(0, 5).unwrap();
-/// assert_eq!(n.mean().unwrap(), 2.5);
+/// let n = DiscreteUniform::new(0, 5)?;
+/// assert_eq!(n.mean(), 2.5);
 /// assert_eq!(n.pmf(3), 1.0 / 6.0);
+/// # Ok::<(), DiscreteUniformError>(())
 /// ```
 #[derive(Copy, Clone, PartialEq, Debug)]
 pub struct DiscreteUniform {
diff --git a/src/distribution/empirical.rs b/src/distribution/empirical.rs
index 4bb1cbba..695c8999 100644
--- a/src/distribution/empirical.rs
+++ b/src/distribution/empirical.rs
@@ -27,7 +27,7 @@ impl<T: PartialOrd> Ord for NonNan<T> {
 ///
 /// ```
 /// use statrs::distribution::{Continuous, Empirical};
-/// use statrs::statistics::Distribution;
+/// use statrs::statistics::*;
 ///
 /// let samples = vec![0.0, 5.0, 10.0];
 ///
diff --git a/src/distribution/erlang.rs b/src/distribution/erlang.rs
index f14720dd..4b79f7c8 100644
--- a/src/distribution/erlang.rs
+++ b/src/distribution/erlang.rs
@@ -10,13 +10,14 @@ use crate::statistics::*;
 /// # Examples
 ///
 /// ```
-/// use statrs::distribution::{Erlang, Continuous};
-/// use statrs::statistics::Distribution;
+/// use statrs::distribution::{Erlang, Continuous, GammaError};
+/// use statrs::statistics::*;
 /// use statrs::prec;
 ///
-/// let n = Erlang::new(3, 1.0).unwrap();
-/// assert_eq!(n.mean().unwrap(), 3.0);
+/// let n = Erlang::new(3, 1.0)?;
+/// assert_eq!(n.mean(), 3.0);
 /// assert!(prec::almost_eq(n.pdf(2.0), 0.270670566473225383788, 1e-15));
+/// # Ok::<(), GammaError>(())
 /// ```
 #[derive(Copy, Clone, PartialEq, Debug)]
 pub struct Erlang {
@@ -52,10 +53,11 @@ impl Erlang {
     /// # Examples
     ///
     /// ```
-    /// use statrs::distribution::Erlang;
+    /// use statrs::distribution::{Erlang, GammaError};
     ///
-    /// let n = Erlang::new(3, 1.0).unwrap();
+    /// let n = Erlang::new(3, 1.0)?;
     /// assert_eq!(n.shape(), 3);
+    /// # Ok::<(), GammaError>(())
     /// ```
     pub fn shape(&self) -> u64 {
         self.g.shape() as u64
@@ -66,10 +68,11 @@ impl Erlang {
     /// # Examples
     ///
     /// ```
-    /// use statrs::distribution::Erlang;
+    /// use statrs::distribution::{Erlang, GammaError};
     ///
     /// let n = Erlang::new(3, 1.0).unwrap();
     /// assert_eq!(n.rate(), 1.0);
+    /// # Ok::<(), GammaError>(())
     /// ```
     pub fn rate(&self) -> f64 {
         self.g.rate()
diff --git a/src/distribution/exponential.rs b/src/distribution/exponential.rs
index bc40025d..7198ce9e 100644
--- a/src/distribution/exponential.rs
+++ b/src/distribution/exponential.rs
@@ -11,12 +11,13 @@ use std::f64;
 /// # Examples
 ///
 /// ```
-/// use statrs::distribution::{Exp, Continuous};
-/// use statrs::statistics::Distribution;
+/// use statrs::distribution::{Exp, Continuous, ExpError};
+/// use statrs::statistics::*;
 ///
-/// let n = Exp::new(1.0).unwrap();
-/// assert_eq!(n.mean().unwrap(), 1.0);
+/// let n = Exp::new(1.0)?;
+/// assert_eq!(n.mean(), 1.0);
 /// assert_eq!(n.pdf(1.0), 0.3678794411714423215955);
+/// # Ok::<(), ExpError>(())
 /// ```
 #[derive(Copy, Clone, PartialEq, Debug)]
 pub struct Exp {
@@ -74,10 +75,11 @@ impl Exp {
     /// # Examples
     ///
     /// ```
-    /// use statrs::distribution::Exp;
+    /// use statrs::distribution::{Exp, ExpError};
     ///
-    /// let n = Exp::new(1.0).unwrap();
+    /// let n = Exp::new(1.0)?;
     /// assert_eq!(n.rate(), 1.0);
+    /// # Ok::<(), ExpError>(())
     /// ```
     pub fn rate(&self) -> f64 {
         self.rate
diff --git a/src/distribution/fisher_snedecor.rs b/src/distribution/fisher_snedecor.rs
index e4c19256..b5f7d01d 100644
--- a/src/distribution/fisher_snedecor.rs
+++ b/src/distribution/fisher_snedecor.rs
@@ -10,13 +10,14 @@ use std::f64;
 /// # Examples
 ///
 /// ```
-/// use statrs::distribution::{FisherSnedecor, Continuous};
-/// use statrs::statistics::Distribution;
+/// use statrs::distribution::{FisherSnedecor, Continuous, FisherSnedecorError};
+/// use statrs::statistics::*;
 /// use statrs::prec;
 ///
-/// let n = FisherSnedecor::new(3.0, 3.0).unwrap();
+/// let n = FisherSnedecor::new(3.0, 3.0)?;
 /// assert_eq!(n.mean().unwrap(), 3.0);
 /// assert!(prec::almost_eq(n.pdf(1.0), 0.318309886183790671538, 1e-15));
+/// # Ok::<(), FisherSnedecorError>(())
 /// ```
 #[derive(Debug, Copy, Clone, PartialEq)]
 pub struct FisherSnedecor {
@@ -92,10 +93,11 @@ impl FisherSnedecor {
     /// # Examples
     ///
     /// ```
-    /// use statrs::distribution::FisherSnedecor;
+    /// use statrs::distribution::{FisherSnedecor, FisherSnedecorError};
     ///
-    /// let n = FisherSnedecor::new(2.0, 3.0).unwrap();
+    /// let n = FisherSnedecor::new(2.0, 3.0)?;
     /// assert_eq!(n.freedom_1(), 2.0);
+    /// # Ok::<(), FisherSnedecorError>(())
     /// ```
     pub fn freedom_1(&self) -> f64 {
         self.freedom_1
@@ -107,10 +109,11 @@ impl FisherSnedecor {
     /// # Examples
     ///
     /// ```
-    /// use statrs::distribution::FisherSnedecor;
+    /// use statrs::distribution::{FisherSnedecor, FisherSnedecorError};
     ///
     /// let n = FisherSnedecor::new(2.0, 3.0).unwrap();
     /// assert_eq!(n.freedom_2(), 3.0);
+    /// # Ok::<(), FisherSnedecorError>(())
     /// ```
     pub fn freedom_2(&self) -> f64 {
         self.freedom_2
diff --git a/src/distribution/gamma.rs b/src/distribution/gamma.rs
index 83122202..cf2a6c38 100644
--- a/src/distribution/gamma.rs
+++ b/src/distribution/gamma.rs
@@ -9,13 +9,14 @@ use crate::statistics::*;
 /// # Examples
 ///
 /// ```
-/// use statrs::distribution::{Gamma, Continuous};
-/// use statrs::statistics::Distribution;
+/// use statrs::distribution::{Gamma, Continuous, GammaError};
+/// use statrs::statistics::*;
 /// use statrs::prec;
 ///
-/// let n = Gamma::new(3.0, 1.0).unwrap();
-/// assert_eq!(n.mean().unwrap(), 3.0);
+/// let n = Gamma::new(3.0, 1.0)?;
+/// assert_eq!(n.mean(), 3.0);
 /// assert!(prec::almost_eq(n.pdf(2.0), 0.270670566473225383788, 1e-15));
+/// # Ok::<(), GammaError>(())
 /// ```
 #[derive(Copy, Clone, PartialEq, Debug)]
 pub struct Gamma {
@@ -91,10 +92,11 @@ impl Gamma {
     /// # Examples
     ///
     /// ```
-    /// use statrs::distribution::Gamma;
+    /// use statrs::distribution::{Gamma, GammaError};
     ///
-    /// let n = Gamma::new(3.0, 1.0).unwrap();
+    /// let n = Gamma::new(3.0, 1.0)?;
     /// assert_eq!(n.shape(), 3.0);
+    /// # Ok::<(), GammaError>(())
     /// ```
     pub fn shape(&self) -> f64 {
         self.shape
@@ -105,10 +107,11 @@ impl Gamma {
     /// # Examples
     ///
     /// ```
-    /// use statrs::distribution::Gamma;
+    /// use statrs::distribution::{Gamma, GammaError};
     ///
     /// let n = Gamma::new(3.0, 1.0).unwrap();
     /// assert_eq!(n.rate(), 1.0);
+    /// # Ok::<(), GammaError>(())
     /// ```
     pub fn rate(&self) -> f64 {
         self.rate
diff --git a/src/distribution/geometric.rs b/src/distribution/geometric.rs
index 10bb2a2b..d86cf942 100644
--- a/src/distribution/geometric.rs
+++ b/src/distribution/geometric.rs
@@ -9,13 +9,14 @@ use std::f64;
 /// # Examples
 ///
 /// ```
-/// use statrs::distribution::{Geometric, Discrete};
-/// use statrs::statistics::Distribution;
+/// use statrs::distribution::{Geometric, GeometricError, Discrete};
+/// use statrs::statistics::*;
 ///
-/// let n = Geometric::new(0.3).unwrap();
-/// assert_eq!(n.mean().unwrap(), 1.0 / 0.3);
+/// let n = Geometric::new(0.3)?;
+/// assert_eq!(n.mean(), 1.0 / 0.3);
 /// assert_eq!(n.pmf(1), 0.3);
 /// assert_eq!(n.pmf(2), 0.21);
+/// # Ok::<(), GeometricError>(())
 /// ```
 #[derive(Copy, Clone, PartialEq, Debug)]
 pub struct Geometric {
@@ -75,9 +76,11 @@ impl Geometric {
     ///
     /// ```
     /// use statrs::distribution::Geometric;
+    /// # use statrs::distribution::GeometricError;
     ///
-    /// let n = Geometric::new(0.5).unwrap();
+    /// let n = Geometric::new(0.5)?;
     /// assert_eq!(n.p(), 0.5);
+    /// # Ok::<(), GeometricError>(())
     /// ```
     pub fn p(&self) -> f64 {
         self.p
diff --git a/src/distribution/gumbel.rs b/src/distribution/gumbel.rs
index 1be77a45..3465159f 100644
--- a/src/distribution/gumbel.rs
+++ b/src/distribution/gumbel.rs
@@ -9,13 +9,14 @@ use std::f64::{self, consts::PI};
 /// # Examples
 ///
 /// ```
-/// use statrs::distribution::{Gumbel, Continuous};
-/// use statrs::{consts::EULER_MASCHERONI, statistics::Distribution};
+/// use statrs::distribution::{Gumbel, Continuous, GumbelError};
+/// use statrs::{consts::EULER_MASCHERONI, statistics::*};
 ///
-/// let n = Gumbel::new(0.0, 1.0).unwrap();
+/// let n = Gumbel::new(0.0, 1.0)?;
 /// assert_eq!(n.location(), 0.0);
-/// assert_eq!(n.skewness().unwrap(), 1.13955);
+/// assert_eq!(n.mean(), 0.5772156649015329);
 /// assert_eq!(n.pdf(0.0), 0.36787944117144233);
+/// # Ok::<(), GumbelError>(())
 /// ```
 #[derive(Copy, Clone, PartialEq, Debug)]
 pub struct Gumbel {
@@ -82,10 +83,11 @@ impl Gumbel {
     /// # Examples
     ///
     /// ```
-    /// use statrs::distribution::Gumbel;
+    /// use statrs::distribution::{Gumbel, GumbelError};
     ///
-    /// let n = Gumbel::new(0.0, 1.0).unwrap();
+    /// let n = Gumbel::new(0.0, 1.0)?;
     /// assert_eq!(n.location(), 0.0);
+    /// # Ok::<(), GumbelError>(())
     /// ```
     pub fn location(&self) -> f64 {
         self.location
@@ -96,10 +98,11 @@ impl Gumbel {
     /// # Examples
     ///
     /// ```
-    /// use statrs::distribution::Gumbel;
+    /// use statrs::distribution::{Gumbel, GumbelError};
     ///
-    /// let n = Gumbel::new(0.0, 1.0).unwrap();
+    /// let n = Gumbel::new(0.0, 1.0)?;
     /// assert_eq!(n.scale(), 1.0);
+    /// # Ok::<(), GumbelError>(())
     /// ```
     pub fn scale(&self) -> f64 {
         self.scale
diff --git a/src/distribution/hypergeometric.rs b/src/distribution/hypergeometric.rs
index b165c2ea..0d6b7bac 100644
--- a/src/distribution/hypergeometric.rs
+++ b/src/distribution/hypergeometric.rs
@@ -11,14 +11,15 @@ use std::f64;
 /// # Examples
 ///
 /// ```
-/// use statrs::distribution::{Hypergeometric, Discrete};
-/// use statrs::statistics::Distribution;
+/// use statrs::distribution::{Hypergeometric, Discrete, HypergeometricError};
+/// use statrs::statistics::*;
 /// use statrs::prec;
 ///
-/// let n = Hypergeometric::new(500, 50, 100).unwrap();
+/// let n = Hypergeometric::new(500, 50, 100)?;
 /// assert_eq!(n.mean().unwrap(), 10.);
 /// assert!(prec::almost_eq(n.pmf(10), 0.14736784, 1e-8));
 /// assert!(prec::almost_eq(n.pmf(25), 3.537e-7, 1e-10));
+/// # Ok::<(), HypergeometricError>(())
 /// ```
 #[derive(Copy, Clone, PartialEq, Eq, Debug)]
 pub struct Hypergeometric {
@@ -97,10 +98,11 @@ impl Hypergeometric {
     /// # Examples
     ///
     /// ```
-    /// use statrs::distribution::Hypergeometric;
+    /// use statrs::distribution::{Hypergeometric, HypergeometricError};
     ///
-    /// let n = Hypergeometric::new(10, 5, 3).unwrap();
+    /// let n = Hypergeometric::new(10, 5, 3)?;
     /// assert_eq!(n.population(), 10);
+    /// # Ok::<(), HypergeometricError>(())
     /// ```
     pub fn population(&self) -> u64 {
         self.population
@@ -112,10 +114,11 @@ impl Hypergeometric {
     /// # Examples
     ///
     /// ```
-    /// use statrs::distribution::Hypergeometric;
+    /// use statrs::distribution::{Hypergeometric, HypergeometricError};
     ///
-    /// let n = Hypergeometric::new(10, 5, 3).unwrap();
+    /// let n = Hypergeometric::new(10, 5, 3)?;
     /// assert_eq!(n.successes(), 5);
+    /// # Ok::<(), HypergeometricError>(())
     /// ```
     pub fn successes(&self) -> u64 {
         self.successes
@@ -127,10 +130,11 @@ impl Hypergeometric {
     /// # Examples
     ///
     /// ```
-    /// use statrs::distribution::Hypergeometric;
+    /// use statrs::distribution::{Hypergeometric, HypergeometricError};
     ///
-    /// let n = Hypergeometric::new(10, 5, 3).unwrap();
+    /// let n = Hypergeometric::new(10, 5, 3)?;
     /// assert_eq!(n.draws(), 3);
+    /// # Ok::<(), HypergeometricError>(())
     /// ```
     pub fn draws(&self) -> u64 {
         self.draws
diff --git a/src/distribution/inverse_gamma.rs b/src/distribution/inverse_gamma.rs
index 9a8c8e58..e9cdfe21 100644
--- a/src/distribution/inverse_gamma.rs
+++ b/src/distribution/inverse_gamma.rs
@@ -10,13 +10,14 @@ use std::f64;
 /// # Examples
 ///
 /// ```
-/// use statrs::distribution::{InverseGamma, Continuous};
-/// use statrs::statistics::Distribution;
+/// use statrs::distribution::{InverseGamma, InverseGammaError, Continuous};
+/// use statrs::statistics::*;
 /// use statrs::prec;
 ///
-/// let n = InverseGamma::new(1.1, 0.1).unwrap();
+/// let n = InverseGamma::new(1.1, 0.1)?;
 /// assert!(prec::almost_eq(n.mean().unwrap(), 1.0, 1e-14));
 /// assert_eq!(n.pdf(1.0), 0.07554920138253064);
+/// # Ok::<(), InverseGammaError>(())
 /// ```
 #[derive(Copy, Clone, PartialEq, Debug)]
 pub struct InverseGamma {
@@ -88,10 +89,11 @@ impl InverseGamma {
     /// # Examples
     ///
     /// ```
-    /// use statrs::distribution::InverseGamma;
+    /// use statrs::distribution::{InverseGamma, InverseGammaError};
     ///
-    /// let n = InverseGamma::new(3.0, 1.0).unwrap();
+    /// let n = InverseGamma::new(3.0, 1.0)?;
     /// assert_eq!(n.shape(), 3.0);
+    /// # Ok::<(), InverseGammaError>(())
     /// ```
     pub fn shape(&self) -> f64 {
         self.shape
@@ -102,10 +104,11 @@ impl InverseGamma {
     /// # Examples
     ///
     /// ```
-    /// use statrs::distribution::InverseGamma;
+    /// use statrs::distribution::{InverseGamma, InverseGammaError};
     ///
-    /// let n = InverseGamma::new(3.0, 1.0).unwrap();
+    /// let n = InverseGamma::new(3.0, 1.0)?;
     /// assert_eq!(n.rate(), 1.0);
+    /// # Ok::<(), InverseGammaError>(())
     /// ```
     pub fn rate(&self) -> f64 {
         self.rate
diff --git a/src/distribution/laplace.rs b/src/distribution/laplace.rs
index e29db1fa..caadb279 100644
--- a/src/distribution/laplace.rs
+++ b/src/distribution/laplace.rs
@@ -8,12 +8,13 @@ use std::f64;
 /// # Examples
 ///
 /// ```
-/// use statrs::distribution::{Laplace, Continuous};
-/// use statrs::statistics::Mode;
+/// use statrs::distribution::{Laplace, Continuous, LaplaceError};
+/// use statrs::statistics::*;
 ///
-/// let n = Laplace::new(0.0, 1.0).unwrap();
+/// let n = Laplace::new(0.0, 1.0)?;
 /// assert_eq!(n.mode().unwrap(), 0.0);
 /// assert_eq!(n.pdf(1.0), 0.18393972058572117);
+/// # Ok::<(), LaplaceError>(())
 /// ```
 #[derive(Copy, Clone, PartialEq, Debug)]
 pub struct Laplace {
@@ -80,10 +81,11 @@ impl Laplace {
     /// # Examples
     ///
     /// ```
-    /// use statrs::distribution::Laplace;
+    /// use statrs::distribution::{Laplace, LaplaceError};
     ///
-    /// let n = Laplace::new(0.0, 1.0).unwrap();
+    /// let n = Laplace::new(0.0, 1.0)?;
     /// assert_eq!(n.location(), 0.0);
+    /// # Ok::<(), LaplaceError>(())
     /// ```
     pub fn location(&self) -> f64 {
         self.location
@@ -94,10 +96,11 @@ impl Laplace {
     /// # Examples
     ///
     /// ```
-    /// use statrs::distribution::Laplace;
+    /// use statrs::distribution::{Laplace, LaplaceError};
     ///
     /// let n = Laplace::new(0.0, 1.0).unwrap();
     /// assert_eq!(n.scale(), 1.0);
+    /// # Ok::<(), LaplaceError>(())
     /// ```
     pub fn scale(&self) -> f64 {
         self.scale
diff --git a/src/distribution/log_normal.rs b/src/distribution/log_normal.rs
index 0b10f5c2..12ce7413 100644
--- a/src/distribution/log_normal.rs
+++ b/src/distribution/log_normal.rs
@@ -11,13 +11,14 @@ use std::f64;
 /// # Examples
 ///
 /// ```
-/// use statrs::distribution::{LogNormal, Continuous};
-/// use statrs::statistics::Distribution;
+/// use statrs::distribution::{LogNormal, LogNormalError, Continuous};
+/// use statrs::statistics::*;
 /// use statrs::prec;
 ///
-/// let n = LogNormal::new(0.0, 1.0).unwrap();
-/// assert_eq!(n.mean().unwrap(), (0.5f64).exp());
+/// let n = LogNormal::new(0.0, 1.0)?;
+/// assert_eq!(n.mean(), (0.5f64).exp());
 /// assert!(prec::almost_eq(n.pdf(1.0), 0.3989422804014326779399, 1e-16));
+/// # Ok::<(), LogNormalError>(())
 /// ```
 #[derive(Copy, Clone, PartialEq, Debug)]
 pub struct LogNormal {
diff --git a/src/distribution/mod.rs b/src/distribution/mod.rs
index 5aa5ff5d..56b33173 100644
--- a/src/distribution/mod.rs
+++ b/src/distribution/mod.rs
@@ -2,8 +2,8 @@
 //! and provides
 //! concrete implementations for a variety of distributions.
 use super::statistics::{Max, Min};
-use ::num_traits::{Float, Num};
 use num_traits::NumAssignOps;
+use num_traits::{Float, Num};
 
 pub use self::bernoulli::Bernoulli;
 pub use self::beta::{Beta, BetaError};
@@ -101,10 +101,11 @@ pub trait ContinuousCDF<K: Float, T: Float>: Min<K> + Max<K> {
     /// # Examples
     ///
     /// ```
-    /// use statrs::distribution::{ContinuousCDF, Uniform};
+    /// use statrs::distribution::{ContinuousCDF, Uniform, UniformError};
     ///
-    /// let n = Uniform::new(0.0, 1.0).unwrap();
+    /// let n = Uniform::new(0.0, 1.0)?;
     /// assert_eq!(0.5, n.cdf(0.5));
+    /// # Ok::<(), UniformError>(())
     /// ```
     fn cdf(&self, x: K) -> T;
 
@@ -115,10 +116,11 @@ pub trait ContinuousCDF<K: Float, T: Float>: Min<K> + Max<K> {
     /// # Examples
     ///
     /// ```
-    /// use statrs::distribution::{ContinuousCDF, Uniform};
+    /// use statrs::distribution::{ContinuousCDF, Uniform, UniformError};
     ///
-    /// let n = Uniform::new(0.0, 1.0).unwrap();
+    /// let n = Uniform::new(0.0, 1.0)?;
     /// assert_eq!(0.5, n.sf(0.5));
+    /// # Ok::<(), UniformError>(())
     /// ```
     fn sf(&self, x: K) -> T {
         T::one() - self.cdf(x)
@@ -175,9 +177,11 @@ pub trait DiscreteCDF<K: Sized + Num + Ord + Clone + NumAssignOps, T: Float>:
     ///
     /// ```
     /// use statrs::distribution::{DiscreteCDF, DiscreteUniform};
+    /// # use statrs::distribution::DiscreteUniformError;
     ///
-    /// let n = DiscreteUniform::new(1, 10).unwrap();
+    /// let n = DiscreteUniform::new(1, 10)?;
     /// assert_eq!(0.6, n.cdf(6));
+    /// # Ok::<(), DiscreteUniformError>(())
     /// ```
     fn cdf(&self, x: K) -> T;
 
@@ -188,9 +192,11 @@ pub trait DiscreteCDF<K: Sized + Num + Ord + Clone + NumAssignOps, T: Float>:
     ///
     /// ```
     /// use statrs::distribution::{DiscreteCDF, DiscreteUniform};
+    /// # use statrs::distribution::DiscreteUniformError;
     ///
-    /// let n = DiscreteUniform::new(1, 10).unwrap();
+    /// let n = DiscreteUniform::new(1, 10)?;
     /// assert_eq!(0.4, n.sf(6));
+    /// # Ok::<(), DiscreteUniformError>(())
     /// ```
     fn sf(&self, x: K) -> T {
         T::one() - self.cdf(x)
@@ -238,9 +244,11 @@ pub trait Continuous<K, T> {
     ///
     /// ```
     /// use statrs::distribution::{Continuous, Uniform};
+    /// # use statrs::distribution::UniformError;
     ///
-    /// let n = Uniform::new(0.0, 1.0).unwrap();
+    /// let n = Uniform::new(0.0, 1.0)?;
     /// assert_eq!(1.0, n.pdf(0.5));
+    /// # Ok::<(), UniformError>(())
     /// ```
     fn pdf(&self, x: K) -> T;
 
@@ -252,9 +260,11 @@ pub trait Continuous<K, T> {
     ///
     /// ```
     /// use statrs::distribution::{Continuous, Uniform};
+    /// # use statrs::distribution::UniformError;
     ///
-    /// let n = Uniform::new(0.0, 1.0).unwrap();
+    /// let n = Uniform::new(0.0, 1.0)?;
     /// assert_eq!(0.0, n.ln_pdf(0.5));
+    /// # Ok::<(), UniformError>(())
     /// ```
     fn ln_pdf(&self, x: K) -> T;
 }
@@ -265,8 +275,7 @@ pub trait Continuous<K, T> {
 /// # Remarks
 ///
 /// All methods provided by the `Discrete` trait are unchecked, meaning
-/// they can panic if in an invalid state or encountering invalid input
-/// depending on the implementing distribution.
+/// they can panic if in an invalid state or encountering invalid input depending on the implementing distribution.
 pub trait Discrete<K, T> {
     /// Returns the probability mass function calculated at `x` for a given
     /// distribution.
@@ -277,9 +286,11 @@ pub trait Discrete<K, T> {
     /// ```
     /// use statrs::distribution::{Discrete, Binomial};
     /// use statrs::prec;
+    /// # use statrs::distribution::BinomialError;
     ///
-    /// let n = Binomial::new(0.5, 10).unwrap();
+    /// let n = Binomial::new(0.5, 10)?;
     /// assert!(prec::almost_eq(n.pmf(5), 0.24609375, 1e-15));
+    /// # Ok::<(), BinomialError>(())
     /// ```
     fn pmf(&self, x: K) -> T;
 
@@ -293,8 +304,11 @@ pub trait Discrete<K, T> {
     /// use statrs::distribution::{Discrete, Binomial};
     /// use statrs::prec;
     ///
-    /// let n = Binomial::new(0.5, 10).unwrap();
+    /// # use statrs::distribution::BinomialError;
+    ///
+    /// let n = Binomial::new(0.5, 10)?;
     /// assert!(prec::almost_eq(n.ln_pmf(5), (0.24609375f64).ln(), 1e-15));
+    /// # Ok::<(), BinomialError>(())
     /// ```
     fn ln_pmf(&self, x: K) -> T;
 }
diff --git a/src/distribution/multinomial.rs b/src/distribution/multinomial.rs
index bed93ef5..243ead4b 100644
--- a/src/distribution/multinomial.rs
+++ b/src/distribution/multinomial.rs
@@ -12,12 +12,13 @@ use nalgebra::{Cholesky, DVector, Dim, DimMin, Dyn, OMatrix, OVector};
 /// # Examples
 ///
 /// ```
-/// use statrs::distribution::Multinomial;
-/// use statrs::statistics::MeanN;
+/// use statrs::distribution::{Multinomial, MultinomialError};
+/// use statrs::statistics::StandardizedMoment;
 /// use nalgebra::vector;
 ///
-/// let n = Multinomial::new_from_nalgebra(vector![0.3, 0.7], 5).unwrap();
-/// assert_eq!(n.mean().unwrap(), (vector![1.5, 3.5]));
+/// let n = Multinomial::new_from_nalgebra(vector![0.3, 0.7], 5)?;
+/// assert_eq!(n.mean(), (vector![1.5, 3.5]));
+/// # Ok::<(), MultinomialError>(())
 /// ```
 #[derive(Debug, Clone, PartialEq)]
 pub struct Multinomial<D>
@@ -123,11 +124,12 @@ where
     /// # Examples
     ///
     /// ```
-    /// use statrs::distribution::Multinomial;
+    /// use statrs::distribution::{Multinomial, MultinomialError};
     /// use nalgebra::dvector;
     ///
-    /// let n = Multinomial::new(vec![0.0, 1.0, 2.0], 3).unwrap();
+    /// let n = Multinomial::new(vec![0.0, 1.0, 2.0], 3)?;
     /// assert_eq!(*n.p(), dvector![0.0, 1.0/3.0, 2.0/3.0]);
+    /// # Ok::<(), MultinomialError>(())
     /// ```
     pub fn p(&self) -> &OVector<f64, D> {
         &self.p
@@ -139,10 +141,11 @@ where
     /// # Examples
     ///
     /// ```
-    /// use statrs::distribution::Multinomial;
+    /// use statrs::distribution::{Multinomial, MultinomialError};
     ///
-    /// let n = Multinomial::new(vec![0.0, 1.0, 2.0], 3).unwrap();
+    /// let n = Multinomial::new(vec![0.0, 1.0, 2.0], 3)?;
     /// assert_eq!(n.n(), 3);
+    /// # Ok::<(), MultinomialError>(())
     /// ```
     pub fn n(&self) -> u64 {
         self.n
diff --git a/src/distribution/multivariate_normal.rs b/src/distribution/multivariate_normal.rs
index 4fa2ff4d..949e15cc 100644
--- a/src/distribution/multivariate_normal.rs
+++ b/src/distribution/multivariate_normal.rs
@@ -76,14 +76,16 @@ where
 /// # Examples
 ///
 /// ```
-/// use statrs::distribution::{MultivariateNormal, Continuous};
+/// use statrs::distribution::{MultivariateNormal, MultivariateNormalError, Continuous};
 /// use nalgebra::{matrix, vector};
-/// use statrs::statistics::{MeanN, VarianceN};
+/// use statrs::statistics::*;
 ///
-/// let mvn = MultivariateNormal::new_from_nalgebra(vector![0., 0.], matrix![1., 0.; 0., 1.]).unwrap();
-/// assert_eq!(mvn.mean().unwrap(), vector![0., 0.]);
-/// assert_eq!(mvn.variance().unwrap(), matrix![1., 0.; 0., 1.]);
+/// let mvn = MultivariateNormal::new_from_nalgebra(vector![0., 0.], matrix![1., 0.; 0.,
+/// 1.])?;
+/// assert_eq!(mvn.mean(), vector![0., 0.]);
+/// assert_eq!(mvn.variance(), matrix![1., 0.; 0., 1.]);
 /// assert_eq!(mvn.pdf(&vector![1.,  1.]), 0.05854983152431917);
+/// # Ok::<(), MultivariateNormalError>(())
 /// ```
 #[derive(Clone, PartialEq, Debug)]
 pub struct MultivariateNormal<D>
diff --git a/src/distribution/multivariate_students_t.rs b/src/distribution/multivariate_students_t.rs
index 734ac568..081adc37 100644
--- a/src/distribution/multivariate_students_t.rs
+++ b/src/distribution/multivariate_students_t.rs
@@ -12,14 +12,15 @@ use std::f64::consts::PI;
 /// # Examples
 ///
 /// ```
-/// use statrs::distribution::{MultivariateStudent, Continuous};
+/// use statrs::distribution::{MultivariateStudent, Continuous, MultivariateStudentError};
 /// use nalgebra::{DVector, DMatrix};
-/// use statrs::statistics::{MeanN, VarianceN};
+/// use statrs::statistics::*;
 ///
-/// let mvs = MultivariateStudent::new(vec![0., 0.], vec![1., 0., 0., 1.], 4.).unwrap();
+/// let mvs = MultivariateStudent::new(vec![0., 0.], vec![1., 0., 0., 1.], 4.)?;
 /// assert_eq!(mvs.mean().unwrap(), DVector::from_vec(vec![0., 0.]));
 /// assert_eq!(mvs.variance().unwrap(), DMatrix::from_vec(2, 2, vec![2., 0., 0., 2.]));
 /// assert_eq!(mvs.pdf(&DVector::from_vec(vec![1.,  1.])), 0.04715702017537655);
+/// # Ok::<(), MultivariateStudentError>(())
 /// ```
 #[derive(Debug, Clone, PartialEq)]
 pub struct MultivariateStudent<D>
diff --git a/src/distribution/negative_binomial.rs b/src/distribution/negative_binomial.rs
index 8949f2b6..8dbf7cff 100644
--- a/src/distribution/negative_binomial.rs
+++ b/src/distribution/negative_binomial.rs
@@ -24,14 +24,15 @@ use std::f64;
 /// # Examples
 ///
 /// ```
-/// use statrs::distribution::{NegativeBinomial, Discrete};
-/// use statrs::statistics::DiscreteDistribution;
+/// use statrs::distribution::{NegativeBinomial, NegativeBinomialError, Discrete};
+/// use statrs::statistics::*;
 /// use statrs::prec::almost_eq;
 ///
-/// let r = NegativeBinomial::new(4.0, 0.5).unwrap();
-/// assert_eq!(r.mean().unwrap(), 4.0);
+/// let r = NegativeBinomial::new(4.0, 0.5)?;
+/// assert_eq!(r.mean(), 4.0);
 /// assert!(almost_eq(r.pmf(0), 0.0625, 1e-8));
 /// assert!(almost_eq(r.pmf(3), 0.15625, 1e-8));
+/// # Ok::<(), NegativeBinomialError>(())
 /// ```
 #[derive(Copy, Clone, PartialEq, Debug)]
 pub struct NegativeBinomial {
@@ -104,10 +105,11 @@ impl NegativeBinomial {
     /// # Examples
     ///
     /// ```
-    /// use statrs::distribution::NegativeBinomial;
+    /// use statrs::distribution::{NegativeBinomial, NegativeBinomialError};
     ///
-    /// let r = NegativeBinomial::new(5.0, 0.5).unwrap();
+    /// let r = NegativeBinomial::new(5.0, 0.5)?;
     /// assert_eq!(r.p(), 0.5);
+    /// # Ok::<(), NegativeBinomialError>(())
     /// ```
     pub fn p(&self) -> f64 {
         self.p
@@ -119,10 +121,11 @@ impl NegativeBinomial {
     /// # Examples
     ///
     /// ```
-    /// use statrs::distribution::NegativeBinomial;
+    /// use statrs::distribution::{NegativeBinomial, NegativeBinomialError};
     ///
-    /// let r = NegativeBinomial::new(5.0, 0.5).unwrap();
+    /// let r = NegativeBinomial::new(5.0, 0.5)?;
     /// assert_eq!(r.r(), 5.0);
+    /// # Ok::<(), NegativeBinomialError>(())
     /// ```
     pub fn r(&self) -> f64 {
         self.r
diff --git a/src/distribution/normal.rs b/src/distribution/normal.rs
index 1ecc2a3b..009acf4c 100644
--- a/src/distribution/normal.rs
+++ b/src/distribution/normal.rs
@@ -10,12 +10,13 @@ use std::f64;
 /// # Examples
 ///
 /// ```
-/// use statrs::distribution::{Normal, Continuous};
-/// use statrs::statistics::Distribution;
+/// use statrs::distribution::{Normal, Continuous, NormalError};
+/// use statrs::statistics::*;
 ///
-/// let n = Normal::new(0.0, 1.0).unwrap();
-/// assert_eq!(n.mean().unwrap(), 0.0);
+/// let n = Normal::new(0.0, 1.0)?;
+/// assert_eq!(n.mean(), 0.0);
 /// assert_eq!(n.pdf(1.0), 0.2419707245191433497978);
+/// # Ok::<(), NormalError>(())
 /// ```
 #[derive(Copy, Clone, PartialEq, Debug)]
 pub struct Normal {
diff --git a/src/distribution/pareto.rs b/src/distribution/pareto.rs
index 66f9dd56..65b49075 100644
--- a/src/distribution/pareto.rs
+++ b/src/distribution/pareto.rs
@@ -8,13 +8,14 @@ use std::f64;
 /// # Examples
 ///
 /// ```
-/// use statrs::distribution::{Pareto, Continuous};
-/// use statrs::statistics::Distribution;
+/// use statrs::distribution::{Pareto, Continuous, ParetoError};
+/// use statrs::statistics::*;
 /// use statrs::prec;
 ///
-/// let p = Pareto::new(1.0, 2.0).unwrap();
+/// let p = Pareto::new(1.0, 2.0)?;
 /// assert_eq!(p.mean().unwrap(), 2.0);
 /// assert!(prec::almost_eq(p.pdf(2.0), 0.25, 1e-15));
+/// # Ok::<(), ParetoError>(())
 /// ```
 #[derive(Copy, Clone, PartialEq, Debug)]
 pub struct Pareto {
@@ -82,10 +83,11 @@ impl Pareto {
     /// # Examples
     ///
     /// ```
-    /// use statrs::distribution::Pareto;
+    /// use statrs::distribution::{Pareto, ParetoError};
     ///
-    /// let n = Pareto::new(1.0, 2.0).unwrap();
+    /// let n = Pareto::new(1.0, 2.0)?;
     /// assert_eq!(n.scale(), 1.0);
+    /// # Ok::<(), ParetoError>(())
     /// ```
     pub fn scale(&self) -> f64 {
         self.scale
@@ -96,10 +98,11 @@ impl Pareto {
     /// # Examples
     ///
     /// ```
-    /// use statrs::distribution::Pareto;
+    /// use statrs::distribution::{Pareto, ParetoError};
     ///
-    /// let n = Pareto::new(1.0, 2.0).unwrap();
+    /// let n = Pareto::new(1.0, 2.0)?;
     /// assert_eq!(n.shape(), 2.0);
+    /// # Ok::<(), ParetoError>(())
     /// ```
     pub fn shape(&self) -> f64 {
         self.shape
diff --git a/src/distribution/poisson.rs b/src/distribution/poisson.rs
index 2146413a..901534ab 100644
--- a/src/distribution/poisson.rs
+++ b/src/distribution/poisson.rs
@@ -9,13 +9,14 @@ use std::f64;
 /// # Examples
 ///
 /// ```
-/// use statrs::distribution::{Poisson, Discrete};
-/// use statrs::statistics::Distribution;
+/// use statrs::distribution::{Poisson, PoissonError, Discrete};
+/// use statrs::statistics::*;
 /// use statrs::prec;
 ///
-/// let n = Poisson::new(1.0).unwrap();
-/// assert_eq!(n.mean().unwrap(), 1.0);
+/// let n = Poisson::new(1.0)?;
+/// assert_eq!(n.mean(), 1.0);
 /// assert!(prec::almost_eq(n.pmf(1), 0.367879441171442, 1e-15));
+/// # Ok::<(), PoissonError>(())
 /// ```
 #[derive(Copy, Clone, PartialEq, Debug)]
 pub struct Poisson {
@@ -73,10 +74,11 @@ impl Poisson {
     /// # Examples
     ///
     /// ```
-    /// use statrs::distribution::Poisson;
+    /// use statrs::distribution::{Poisson, PoissonError};
     ///
-    /// let n = Poisson::new(1.0).unwrap();
+    /// let n = Poisson::new(1.0)?;
     /// assert_eq!(n.lambda(), 1.0);
+    /// # Ok::<(), PoissonError>(())
     /// ```
     pub fn lambda(&self) -> f64 {
         self.lambda
diff --git a/src/distribution/students_t.rs b/src/distribution/students_t.rs
index 1225c2cb..d33c08fc 100644
--- a/src/distribution/students_t.rs
+++ b/src/distribution/students_t.rs
@@ -9,13 +9,14 @@ use std::f64;
 /// # Examples
 ///
 /// ```
-/// use statrs::distribution::{StudentsT, Continuous};
-/// use statrs::statistics::Distribution;
+/// use statrs::distribution::{StudentsT, Continuous, StudentsTError};
+/// use statrs::statistics::*;
 /// use statrs::prec;
 ///
-/// let n = StudentsT::new(0.0, 1.0, 2.0).unwrap();
+/// let n = StudentsT::new(0.0, 1.0, 2.0)?;
 /// assert_eq!(n.mean().unwrap(), 0.0);
 /// assert!(prec::almost_eq(n.pdf(0.0), 0.353553390593274, 1e-15));
+/// # Ok::<(), StudentsTError>(())
 /// ```
 #[derive(Copy, Clone, PartialEq, Debug)]
 pub struct StudentsT {
@@ -98,10 +99,11 @@ impl StudentsT {
     /// # Examples
     ///
     /// ```
-    /// use statrs::distribution::StudentsT;
+    /// use statrs::distribution::{StudentsT, StudentsTError};
     ///
-    /// let n = StudentsT::new(0.0, 1.0, 2.0).unwrap();
+    /// let n = StudentsT::new(0.0, 1.0, 2.0)?;
     /// assert_eq!(n.location(), 0.0);
+    /// # Ok::<(), StudentsTError>(())
     /// ```
     pub fn location(&self) -> f64 {
         self.location
@@ -112,10 +114,11 @@ impl StudentsT {
     /// # Examples
     ///
     /// ```
-    /// use statrs::distribution::StudentsT;
+    /// use statrs::distribution::{StudentsT, StudentsTError};
     ///
     /// let n = StudentsT::new(0.0, 1.0, 2.0).unwrap();
     /// assert_eq!(n.scale(), 1.0);
+    /// # Ok::<(), StudentsTError>(())
     /// ```
     pub fn scale(&self) -> f64 {
         self.scale
@@ -126,10 +129,11 @@ impl StudentsT {
     /// # Examples
     ///
     /// ```
-    /// use statrs::distribution::StudentsT;
+    /// use statrs::distribution::{StudentsT, StudentsTError};
     ///
     /// let n = StudentsT::new(0.0, 1.0, 2.0).unwrap();
     /// assert_eq!(n.freedom(), 2.0);
+    /// # Ok::<(), StudentsTError>(())
     /// ```
     pub fn freedom(&self) -> f64 {
         self.freedom
diff --git a/src/distribution/triangular.rs b/src/distribution/triangular.rs
index 363620fa..9fb146d1 100644
--- a/src/distribution/triangular.rs
+++ b/src/distribution/triangular.rs
@@ -9,12 +9,13 @@ use std::f64;
 /// # Examples
 ///
 /// ```
-/// use statrs::distribution::{Triangular, Continuous};
-/// use statrs::statistics::Distribution;
+/// use statrs::distribution::{Triangular, Continuous, TriangularError};
+/// use statrs::statistics::*;
 ///
-/// let n = Triangular::new(0.0, 5.0, 2.5).unwrap();
-/// assert_eq!(n.mean().unwrap(), 7.5 / 3.0);
+/// let n = Triangular::new(0.0, 5.0, 2.5)?;
+/// assert_eq!(n.mean(), 7.5 / 3.0);
 /// assert_eq!(n.pdf(2.5), 5.0 / 12.5);
+/// # Ok::<(), TriangularError>(())
 /// ```
 #[derive(Copy, Clone, PartialEq, Debug)]
 pub struct Triangular {
diff --git a/src/distribution/uniform.rs b/src/distribution/uniform.rs
index c13ba454..021f8b68 100644
--- a/src/distribution/uniform.rs
+++ b/src/distribution/uniform.rs
@@ -10,12 +10,13 @@ use std::fmt::Debug;
 /// # Examples
 ///
 /// ```
-/// use statrs::distribution::{Uniform, Continuous};
-/// use statrs::statistics::Distribution;
+/// use statrs::distribution::{Uniform, Continuous, UniformError};
+/// use statrs::statistics::*;
 ///
-/// let n = Uniform::new(0.0, 1.0).unwrap();
-/// assert_eq!(n.mean().unwrap(), 0.5);
+/// let n = Uniform::new(0.0, 1.0)?;
+/// assert_eq!(n.mean(), 0.5);
 /// assert_eq!(n.pdf(0.5), 1.0);
+/// # Ok::<(), UniformError>(())
 /// ```
 #[derive(Debug, Copy, Clone, PartialEq)]
 pub struct Uniform {
diff --git a/src/distribution/weibull.rs b/src/distribution/weibull.rs
index ea5054a3..ffe2f4e8 100644
--- a/src/distribution/weibull.rs
+++ b/src/distribution/weibull.rs
@@ -10,14 +10,15 @@ use std::f64;
 /// # Examples
 ///
 /// ```
-/// use statrs::distribution::{Weibull, Continuous};
-/// use statrs::statistics::Distribution;
+/// use statrs::distribution::{Weibull, Continuous, WeibullError};
+/// use statrs::statistics::*;
 /// use statrs::prec;
 ///
-/// let n = Weibull::new(10.0, 1.0).unwrap();
-/// assert!(prec::almost_eq(n.mean().unwrap(),
+/// let n = Weibull::new(10.0, 1.0)?;
+/// assert!(prec::almost_eq(n.mean(),
 /// 0.95135076986687318362924871772654021925505786260884, 1e-15));
 /// assert_eq!(n.pdf(1.0), 3.6787944117144232159552377016146086744581113103177);
+/// # Ok::<(), WeibullError>(())
 /// ```
 #[derive(Copy, Clone, PartialEq, Debug)]
 pub struct Weibull {
@@ -90,10 +91,11 @@ impl Weibull {
     /// # Examples
     ///
     /// ```
-    /// use statrs::distribution::Weibull;
+    /// use statrs::distribution::{Weibull, WeibullError};
     ///
-    /// let n = Weibull::new(10.0, 1.0).unwrap();
+    /// let n = Weibull::new(10.0, 1.0)?;
     /// assert_eq!(n.shape(), 10.0);
+    /// # Ok::<(), WeibullError>(())
     /// ```
     pub fn shape(&self) -> f64 {
         self.shape
@@ -104,10 +106,11 @@ impl Weibull {
     /// # Examples
     ///
     /// ```
-    /// use statrs::distribution::Weibull;
+    /// use statrs::distribution::{Weibull, WeibullError};
     ///
-    /// let n = Weibull::new(10.0, 1.0).unwrap();
+    /// let n = Weibull::new(10.0, 1.0)?;
     /// assert_eq!(n.scale(), 1.0);
+    /// # Ok::<(), WeibullError>(())
     /// ```
     pub fn scale(&self) -> f64 {
         self.scale
diff --git a/src/generate.rs b/src/generate.rs
index c353f45f..0f585012 100644
--- a/src/generate.rs
+++ b/src/generate.rs
@@ -265,20 +265,15 @@ impl InfiniteTriangle {
     /// # Examples
     ///
     /// ```
-    /// #[macro_use]
-    /// extern crate statrs;
-    ///
+    /// use approx::assert_relative_eq;
     /// use statrs::generate::InfiniteTriangle;
-    ///
-    /// # fn main() {
     /// let x = InfiniteTriangle::new(4, 7, 1.0, -1.0,
     /// 1).take(12).collect::<Vec<f64>>();
     /// let expected: [f64; 12] = [-0.714, -1.0, -0.5, 0.0, 0.5, 1.0, 0.714,
     /// 0.429, 0.143, -0.143, -0.429, -0.714];
     /// for (&left, &right) in x.iter().zip(expected.iter()) {
-    ///     assert_almost_eq!(left, right, 1e-3);
+    ///     assert_relative_eq!(left, right, epsilon=1e-3);
     /// }
-    /// # }
     /// ```
     pub fn new(
         raise_duration: i64,
diff --git a/src/lib.rs b/src/lib.rs
index d4d80ebc..8938f626 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -19,27 +19,29 @@
 //! # Introspecting distributions
 //! Statrs also comes with a number of useful utility traits for more detailed introspection of distributions.
 //! ```
-//! use statrs::distribution::{Exp, Continuous, ContinuousCDF}; // `cdf` and `pdf`
-//! use statrs::statistics::Distribution; // statistical moments and entropy
+//! use statrs::distribution::{Exp, Continuous, ContinuousCDF, ExpError}; // `cdf` and `pdf`
+//! use statrs::statistics::*;
 //!
-//! let n = Exp::new(1.0).unwrap();
-//! assert_eq!(n.mean(), Some(1.0));
-//! assert_eq!(n.variance(), Some(1.0));
-//! assert_eq!(n.entropy(), Some(1.0));
-//! assert_eq!(n.skewness(), Some(2.0));
+//! let n = Exp::new(1.0)?;
+//! assert_eq!(n.mean(), 1.0);
+//! assert_eq!(n.variance(), 1.0);
+//! assert_eq!(n.entropy(), 1.0);
+//! assert_eq!(n.skewness(), 2.0);
 //! assert_eq!(n.cdf(1.0), 0.6321205588285576784045);
 //! assert_eq!(n.pdf(1.0), 0.3678794411714423215955);
+//! # Ok::<(), ExpError>(())
 //! ```
 //!
 //! # Utility functions
 //! as well as utility functions including `erf`, `gamma`, `ln_gamma`, `beta`, etc.
 //!
 //! ```
-//! use statrs::distribution::FisherSnedecor;
-//! use statrs::statistics::Distribution;
+//! use statrs::distribution::{FisherSnedecor, FisherSnedecorError};
+//! use statrs::statistics::*;
 //!
-//! let n = FisherSnedecor::new(1.0, 1.0).unwrap();
+//! let n = FisherSnedecor::new(1.0, 1.0)?;
 //! assert!(n.variance().is_none());
+//! # Ok::<(), FisherSnedecorError>(())
 //! ```
 //! ## Distributions implemented
 //! Statrs comes with a number of commonly used distributions including Normal, Gamma, Student's T, Exponential, Weibull, etc. view all implemented in `distributions` module.
diff --git a/src/statistics/slice_statistics.rs b/src/statistics/slice_statistics.rs
index 1f81e2e7..3b1c811a 100644
--- a/src/statistics/slice_statistics.rs
+++ b/src/statistics/slice_statistics.rs
@@ -317,13 +317,13 @@ impl<D: AsMut<[f64]> + AsRef<[f64]>> StandardizedMoment<f64> for Data<D> {
     /// #[macro_use]
     /// extern crate statrs;
     ///
-    /// use statrs::statistics::Distribution;
+    /// use statrs::statistics::StandardizedMoment;
     /// use statrs::statistics::Data;
     ///
     /// # fn main() {
     /// let x = [];
     /// let x = Data::new(x);
-    /// assert!(x.mean().unwrap().is_nan());
+    /// assert!(x.mean().is_none());
     ///
     /// let y = [0.0, f64::NAN, 3.0, -2.0];
     /// let y = Data::new(y);
@@ -362,12 +362,12 @@ impl<D: AsMut<[f64]> + AsRef<[f64]>> StandardizedMoment<f64> for Data<D> {
     /// # Examples
     ///
     /// ```
-    /// use statrs::statistics::Distribution;
+    /// use statrs::statistics::StandardizedMoment;
     /// use statrs::statistics::Data;
     ///
     /// let x = [];
     /// let x = Data::new(x);
-    /// assert!(x.variance().unwrap().is_nan());
+    /// assert!(x.variance().is_none());
     ///
     /// let y = [0.0, f64::NAN, 3.0, -2.0];
     /// let y = Data::new(y);
@@ -378,19 +378,18 @@ impl<D: AsMut<[f64]> + AsRef<[f64]>> StandardizedMoment<f64> for Data<D> {
     /// assert_eq!(z.variance().unwrap(), 19.0 / 3.0);
     /// ```
     fn variance(&self) -> Option<f64> {
-        if self.0.as_ref().is_empty() {
-            return None;
-        }
+        // TODO: rewrite iterator statistics and rely on them instead
+        let mut iter = self.0.as_ref().iter().enumerate();
 
-        let mut sum = 0.0;
+        let mut sum = *iter.next()?.1;
         let mut variance = 0.0;
 
-        for (i, &x) in self.0.as_ref().iter().enumerate() {
+        for (i, &x) in iter {
             if x.is_nan() {
                 return Some(f64::NAN);
             }
-            sum += x;
             let diff = (i + 1) as f64 * x - sum;
+            sum += x;
             variance += diff * diff / (i * (i + 1)) as f64
         }
         Some(variance / (self.0.as_ref().len() - 1) as f64)
diff --git a/src/statistics/traits.rs b/src/statistics/traits.rs
index deb8948c..9c9428f7 100644
--- a/src/statistics/traits.rs
+++ b/src/statistics/traits.rs
@@ -9,10 +9,11 @@ pub trait Min<T> {
     ///
     /// ```
     /// use statrs::statistics::Min;
-    /// use statrs::distribution::Uniform;
+    /// use statrs::distribution::{Uniform, UniformError};
     ///
-    /// let n = Uniform::new(0.0, 1.0).unwrap();
+    /// let n = Uniform::new(0.0, 1.0)?;
     /// assert_eq!(0.0, n.min());
+    /// # Ok::<(), UniformError>(())
     /// ```
     fn min(&self) -> T;
 }
@@ -26,10 +27,11 @@ pub trait Max<T> {
     ///
     /// ```
     /// use statrs::statistics::Max;
-    /// use statrs::distribution::Uniform;
+    /// use statrs::distribution::{Uniform, UniformError};
     ///
-    /// let n = Uniform::new(0.0, 1.0).unwrap();
+    /// let n = Uniform::new(0.0, 1.0)?;
     /// assert_eq!(1.0, n.max());
+    /// # Ok::<(), UniformError>(())
     /// ```
     fn max(&self) -> T;
 }
@@ -238,10 +240,11 @@ pub trait Median<T> {
     ///
     /// ```
     /// use statrs::statistics::Median;
-    /// use statrs::distribution::Uniform;
+    /// use statrs::distribution::{Uniform, UniformError};
     ///
-    /// let n = Uniform::new(0.0, 1.0).unwrap();
+    /// let n = Uniform::new(0.0, 1.0)?;
     /// assert_eq!(0.5, n.median());
+    /// # Ok::<(), UniformError>(())
     /// ```
     fn median(&self) -> T;
 }
@@ -255,10 +258,11 @@ pub trait Mode<T> {
     ///
     /// ```
     /// use statrs::statistics::Mode;
-    /// use statrs::distribution::Uniform;
+    /// use statrs::distribution::{Uniform, UniformError};
     ///
-    /// let n = Uniform::new(0.0, 1.0).unwrap();
+    /// let n = Uniform::new(0.0, 1.0)?;
     /// assert_eq!(Some(0.5), n.mode());
+    /// # Ok::<(), UniformError>(())
     /// ```
     fn mode(&self) -> T;
 }