diff --git a/src/distribution/multinomial.rs b/src/distribution/multinomial.rs
index 0794d352..c128105e 100644
--- a/src/distribution/multinomial.rs
+++ b/src/distribution/multinomial.rs
@@ -1,7 +1,7 @@
 use crate::distribution::Discrete;
 use crate::function::factorial;
 use crate::statistics::*;
-use nalgebra::{DVector, Dim, Dyn, OMatrix, OVector};
+use nalgebra::{Cholesky, DVector, Dim, DimMin, Dyn, OMatrix, OVector};
 
 /// Implements the
 /// [Multinomial](https://en.wikipedia.org/wiki/Multinomial_distribution)
@@ -260,6 +260,46 @@ where
     }
 }
 
+impl<D> CentralMoment<f64> for Multinomial<D>
+where
+    D: DimMin<D>,
+    nalgebra::DefaultAllocator:
+        nalgebra::allocator::Allocator<f64, D> + nalgebra::allocator::Allocator<f64, D, D>,
+{
+    type Mu = OVector<f64, D>;
+    type Var = Cholesky<f64, D>;
+    type Kurt = ();
+    type Skew = ();
+
+    fn mean(&self) -> Self::Mu {
+        self.p.scale(self.n as f64)
+    }
+
+    fn variance(&self) -> Self::Var {
+        let mut cov = OMatrix::from_diagonal(&self.p.map(|x| x * (1.0 - x)));
+        let mut offdiag = |x: usize, y: usize| {
+            let elt = -self.p[x] * self.p[y];
+            // cov[(x, y)] = elt;
+            cov[(y, x)] = elt;
+        };
+
+        for j in 0..self.p.len() {
+            for i in 0..j {
+                offdiag(i, j);
+            }
+        }
+        Cholesky::new_unchecked(cov.scale(self.n as f64))
+    }
+
+    fn excess_kurtosis(&self) -> Self::Kurt {
+        unimplemented!()
+    }
+
+    fn skewness(&self) -> Self::Skew {
+        unimplemented!()
+    }
+}
+
 // impl Skewness<Vec<f64>> for Multinomial {
 //     /// Returns the skewness of the multinomial distribution
 //     ///
diff --git a/src/distribution/multivariate_normal.rs b/src/distribution/multivariate_normal.rs
index 6eb3b777..658cbb71 100644
--- a/src/distribution/multivariate_normal.rs
+++ b/src/distribution/multivariate_normal.rs
@@ -1,5 +1,5 @@
 use crate::distribution::Continuous;
-use crate::statistics::{Max, MeanN, Min, Mode, VarianceN};
+use crate::statistics::{CentralMoment, Max, Min, Mode};
 use nalgebra::{Cholesky, Const, DMatrix, DVector, Dim, DimMin, Dyn, OMatrix, OVector};
 use std::f64;
 use std::f64::consts::{E, PI};
@@ -216,14 +216,11 @@ where
     ///
     /// where `Σ` is the covariance matrix and `det` is the determinant
     pub fn entropy(&self) -> Option<f64> {
-        Some(
-            0.5 * self
-                .variance()
-                .unwrap()
-                .scale(2. * PI * E)
-                .determinant()
-                .ln(),
-        )
+        Some(0.5 * self.variance().scale(2. * PI * E).determinant().ln())
+    }
+
+    pub fn cov(&self) -> &OMatrix<f64, D, D> {
+        &self.cov
     }
 }
 
@@ -264,6 +261,34 @@ where
     }
 }
 
+impl<D> CentralMoment<f64> for MultivariateNormal<D>
+where
+    D: DimMin<D>,
+    nalgebra::DefaultAllocator:
+        nalgebra::allocator::Allocator<f64, D> + nalgebra::allocator::Allocator<f64, D, D>,
+{
+    type Mu = OVector<f64, D>;
+    type Var = OMatrix<f64, D, D>;
+    type Kurt = ();
+    type Skew = ();
+
+    fn mean(&self) -> Self::Mu {
+        self.mu.clone()
+    }
+
+    fn variance(&self) -> Self::Var {
+        self.cov.clone()
+    }
+
+    fn excess_kurtosis(&self) -> Self::Kurt {
+        unimplemented!()
+    }
+
+    fn skewness(&self) -> Self::Skew {
+        unimplemented!()
+    }
+}
+
 impl<D> Min<OVector<f64, D>> for MultivariateNormal<D>
 where
     D: Dim,
@@ -290,34 +315,6 @@ where
     }
 }
 
-impl<D> MeanN<OVector<f64, D>> for MultivariateNormal<D>
-where
-    D: Dim,
-    nalgebra::DefaultAllocator:
-        nalgebra::allocator::Allocator<f64, D> + nalgebra::allocator::Allocator<f64, D, D>,
-{
-    /// Returns the mean of the normal distribution
-    ///
-    /// # Remarks
-    ///
-    /// This is the same mean used to construct the distribution
-    fn mean(&self) -> Option<OVector<f64, D>> {
-        Some(self.mu.clone())
-    }
-}
-
-impl<D> VarianceN<OMatrix<f64, D, D>> for MultivariateNormal<D>
-where
-    D: Dim,
-    nalgebra::DefaultAllocator:
-        nalgebra::allocator::Allocator<f64, D> + nalgebra::allocator::Allocator<f64, D, D>,
-{
-    /// Returns the covariance matrix of the multivariate normal distribution
-    fn variance(&self) -> Option<OMatrix<f64, D, D>> {
-        Some(self.cov.clone())
-    }
-}
-
 impl<D> Mode<OVector<f64, D>> for MultivariateNormal<D>
 where
     D: Dim,
@@ -379,7 +376,7 @@ mod tests  {
 
     use crate::{
         distribution::{Continuous, MultivariateNormal},
-        statistics::{Max, MeanN, Min, Mode, VarianceN},
+        statistics::{Max, Min, Mode, CentralMoment},
     };
 
     use super::MultivariateNormalError;
@@ -404,8 +401,8 @@ mod tests  {
             + nalgebra::allocator::Allocator<(usize, usize), D>,
     {
         let mvn = try_create(mean.clone(), covariance.clone());
-        assert_eq!(mean, mvn.mean().unwrap());
-        assert_eq!(covariance, mvn.variance().unwrap());
+        assert_eq!(mean, mvn.mean());
+        assert_eq!(covariance, mvn.variance());
     }
 
     fn bad_create_case<D>(mean: OVector<f64, D>, covariance: OMatrix<f64, D, D>)
@@ -477,7 +474,7 @@ mod tests  {
 
     #[test]
     fn test_variance() {
-        let variance = |x: MultivariateNormal<_>| x.variance().unwrap();
+        let variance = |x: MultivariateNormal<_>| x.variance();
         test_case(
             vector![0., 0.],
             matrix![1., 0.; 0., 1.],
diff --git a/src/distribution/multivariate_students_t.rs b/src/distribution/multivariate_students_t.rs
index 919e3ebf..bb7e5a1e 100644
--- a/src/distribution/multivariate_students_t.rs
+++ b/src/distribution/multivariate_students_t.rs
@@ -1,6 +1,6 @@
 use crate::distribution::Continuous;
 use crate::function::gamma;
-use crate::statistics::{Max, MeanN, Min, Mode, VarianceN};
+use crate::statistics::{CentralMoment, Max, Min, Mode};
 use nalgebra::{Cholesky, Const, DMatrix, Dim, DimMin, Dyn, OMatrix, OVector};
 use std::f64::consts::PI;
 
@@ -259,52 +259,41 @@ where
     }
 }
 
-impl<D> MeanN<OVector<f64, D>> for MultivariateStudent<D>
+impl<D> CentralMoment<f64> for MultivariateStudent<D>
 where
-    D: Dim,
+    D: DimMin<D>,
     nalgebra::DefaultAllocator:
         nalgebra::allocator::Allocator<f64, D> + nalgebra::allocator::Allocator<f64, D, D>,
 {
-    /// Returns the mean of the student distribution.
-    ///
-    /// # Remarks
-    ///
-    /// This is the same mean used to construct the distribution if
-    /// the degrees of freedom is larger than 1.
-    fn mean(&self) -> Option<OVector<f64, D>> {
+    type Mu = Option<OVector<f64, D>>;
+    type Var = Option<OMatrix<f64, D, D>>;
+    type Kurt = ();
+    type Skew = ();
+
+    fn mean(&self) -> Self::Mu {
         if self.freedom > 1. {
             Some(self.location.clone())
         } else {
             None
         }
     }
-}
 
-impl<D> VarianceN<OMatrix<f64, D, D>> for MultivariateStudent<D>
-where
-    D: Dim,
-    nalgebra::DefaultAllocator:
-        nalgebra::allocator::Allocator<f64, D> + nalgebra::allocator::Allocator<f64, D, D>,
-{
-    /// Returns the covariance matrix of the multivariate student distribution.
-    ///
-    /// # Formula
-    ///
-    /// ```math
-    /// Σ ⋅ ν / (ν - 2)
-    /// ```
-    ///
-    /// where `Σ` is the scale matrix and `ν` is the degrees of freedom.
-    /// Only defined if freedom is larger than 2.
-    fn variance(&self) -> Option<OMatrix<f64, D, D>> {
+    fn variance(&self) -> Self::Var {
         if self.freedom > 2. {
             Some(self.scale.clone() * self.freedom / (self.freedom - 2.))
         } else {
             None
         }
     }
-}
 
+    fn excess_kurtosis(&self) -> Self::Kurt {
+        unimplemented!()
+    }
+
+    fn skewness(&self) -> Self::Skew {
+        unimplemented!()
+    }
+}
 impl<D> Mode<OVector<f64, D>> for MultivariateStudent<D>
 where
     D: Dim,
@@ -399,7 +388,7 @@ mod tests  {
 
     use crate::{
         distribution::{Continuous, MultivariateStudent, MultivariateNormal},
-        statistics::{Max, MeanN, Min, Mode, VarianceN},
+        statistics::{Max,  Min, Mode, CentralMoment},
     };
 
     use super::MultivariateStudentError;