diff --git a/Cargo.toml b/Cargo.toml
index a8b6e4ca..15bd9d64 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -9,7 +9,7 @@ keywords = ["probability", "statistics", "stats", "distribution", "math"]
 categories = ["science"]
 homepage = "https://github.com/statrs-dev/statrs"
 repository = "https://github.com/statrs-dev/statrs"
-edition = "2018"
+edition = "2021"
 
 include = ["CHANGELOG.md", "LICENSE.md", "src/"]
 
@@ -22,12 +22,13 @@ nightly = []
 
 [dependencies]
 rand = "0.8"
-nalgebra = { version = "0.32", features = ["rand"] }
+nalgebra = { version = "0.32", default_features = false, features = ["rand", "std"] }
 approx = "0.5.0"
 num-traits = "0.2.14"
 
 [dev-dependencies]
 criterion = "0.3.3"
+nalgebra = { version = "0.32", default_features = false, features = ["macros"] }
 
 [[bench]]
 name = "order_statistics"
diff --git a/src/distribution/multivariate_normal.rs b/src/distribution/multivariate_normal.rs
index d6f400ec..c75765df 100644
--- a/src/distribution/multivariate_normal.rs
+++ b/src/distribution/multivariate_normal.rs
@@ -3,10 +3,9 @@ use crate::distribution::Normal;
 use crate::statistics::{Max, MeanN, Min, Mode, VarianceN};
 use crate::{Result, StatsError};
 use nalgebra::{
-    base::allocator::Allocator, base::dimension::DimName, Cholesky, DefaultAllocator, Dim, DimMin,
-    LU, U1,
+    base::allocator::Allocator, Cholesky, Const, DMatrix, DVector, DefaultAllocator, Dim, DimMin,
+    Dyn, OMatrix, OVector,
 };
-use nalgebra::{DMatrix, DVector};
 use rand::Rng;
 use std::f64;
 use std::f64::consts::{E, PI};
@@ -18,26 +17,30 @@ use std::f64::consts::{E, PI};
 ///
 /// ```
 /// use statrs::distribution::{MultivariateNormal, Continuous};
-/// use nalgebra::{DVector, DMatrix};
+/// use nalgebra::{matrix, vector};
 /// use statrs::statistics::{MeanN, VarianceN};
 ///
-/// let mvn = MultivariateNormal::new(vec![0., 0.], vec![1., 0., 0., 1.]).unwrap();
-/// assert_eq!(mvn.mean().unwrap(), DVector::from_vec(vec![0., 0.]));
-/// assert_eq!(mvn.variance().unwrap(), DMatrix::from_vec(2, 2, vec![1., 0., 0., 1.]));
-/// assert_eq!(mvn.pdf(&DVector::from_vec(vec![1.,  1.])), 0.05854983152431917);
+/// 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.]);
+/// assert_eq!(mvn.pdf(&vector![1.,  1.]), 0.05854983152431917);
 /// ```
 #[derive(Debug, Clone, PartialEq)]
-pub struct MultivariateNormal {
-    dim: usize,
-    cov_chol_decomp: DMatrix<f64>,
-    mu: DVector<f64>,
-    cov: DMatrix<f64>,
-    precision: DMatrix<f64>,
+pub struct MultivariateNormal<D>
+where
+    D: Dim,
+    nalgebra::DefaultAllocator:
+        nalgebra::allocator::Allocator<f64, D> + nalgebra::allocator::Allocator<f64, D, D>,
+{
+    cov_chol_decomp: OMatrix<f64, D, D>,
+    mu: OVector<f64, D>,
+    cov: OMatrix<f64, D, D>,
+    precision: OMatrix<f64, D, D>,
     pdf_const: f64,
 }
 
-impl MultivariateNormal {
-    /// Constructs a new multivariate normal distribution with a mean of `mean`
+impl MultivariateNormal<Dyn> {
+    ///  Constructs a new multivariate normal distribution with a mean of `mean`
     /// and covariance matrix `cov`
     ///
     /// # Errors
@@ -47,19 +50,26 @@ impl MultivariateNormal {
     pub fn new(mean: Vec<f64>, cov: Vec<f64>) -> Result<Self> {
         let mean = DVector::from_vec(mean);
         let cov = DMatrix::from_vec(mean.len(), mean.len(), cov);
-        return MultivariateNormal::new_from_nalgebra(mean, cov)
+        MultivariateNormal::new_from_nalgebra(mean, cov)
     }
+}
 
+impl<D> MultivariateNormal<D>
+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>,
+{
     /// Constructs a new multivariate normal distribution with a mean of `mean`
-    /// and covariance matrix `cov`, but with explicitly using nalgebras
-    /// DVector and DMatrix instead of Vec<f64>
+    /// and covariance matrix `cov` using `nalgebra` `OVector` and `OMatrix`
+    /// instead of `Vec<f64>`
     ///
     /// # Errors
     ///
     /// Returns an error if the given covariance matrix is not
     /// symmetric or positive-definite
-    pub fn new_from_nalgebra(mean: DVector<f64>, cov: DMatrix<f64>) -> Result<Self> {
-        let dim = mean.len();
+    pub fn new_from_nalgebra(mean: OVector<f64, D>, cov: OMatrix<f64, D, D>) -> Result<Self> {
         // Check that the provided covariance matrix is symmetric
         if cov.lower_triangle() != cov.upper_triangle().transpose()
         // Check that mean and covariance do not contain NaN
@@ -81,7 +91,6 @@ impl MultivariateNormal {
             Some(cholesky_decomp) => {
                 let precision = cholesky_decomp.inverse();
                 Ok(MultivariateNormal {
-                    dim,
                     cov_chol_decomp: cholesky_decomp.unpack(),
                     mu: mean,
                     cov,
@@ -113,7 +122,12 @@ impl MultivariateNormal {
     }
 }
 
-impl ::rand::distributions::Distribution<DVector<f64>> for MultivariateNormal {
+impl<D> ::rand::distributions::Distribution<OVector<f64, D>> for MultivariateNormal<D>
+where
+    D: Dim,
+    nalgebra::DefaultAllocator:
+        nalgebra::allocator::Allocator<f64, D> + nalgebra::allocator::Allocator<f64, D, D>,
+{
     /// Samples from the multivariate normal distribution
     ///
     /// # Formula
@@ -125,52 +139,73 @@ impl ::rand::distributions::Distribution<DVector<f64>> for MultivariateNormal {
     /// `Z` is a vector of normally distributed random variables, and
     /// `μ` is the mean vector
 
-    fn sample<R: Rng + ?Sized>(&self, rng: &mut R) -> DVector<f64> {
+    fn sample<R: Rng + ?Sized>(&self, rng: &mut R) -> OVector<f64, D> {
         let d = Normal::new(0., 1.).unwrap();
-        let z = DVector::<f64>::from_distribution(self.dim, &d, rng);
+        let z = OVector::from_distribution_generic(self.mu.shape_generic().0, Const::<1>, &d, rng);
         (&self.cov_chol_decomp * z) + &self.mu
     }
 }
 
-impl Min<DVector<f64>> for MultivariateNormal {
+impl<D> Min<OVector<f64, D>> for MultivariateNormal<D>
+where
+    D: Dim,
+    nalgebra::DefaultAllocator:
+        nalgebra::allocator::Allocator<f64, D> + nalgebra::allocator::Allocator<f64, D, D>,
+{
     /// Returns the minimum value in the domain of the
     /// multivariate normal distribution represented by a real vector
-    fn min(&self) -> DVector<f64> {
-        DVector::from_vec(vec![f64::NEG_INFINITY; self.dim])
+    fn min(&self) -> OVector<f64, D> {
+        OMatrix::repeat_generic(self.mu.shape_generic().0, Const::<1>, f64::NEG_INFINITY)
     }
 }
 
-impl Max<DVector<f64>> for MultivariateNormal {
+impl<D> Max<OVector<f64, D>> for MultivariateNormal<D>
+where
+    D: Dim,
+    nalgebra::DefaultAllocator:
+        nalgebra::allocator::Allocator<f64, D> + nalgebra::allocator::Allocator<f64, D, D>,
+{
     /// Returns the maximum value in the domain of the
     /// multivariate normal distribution represented by a real vector
-    fn max(&self) -> DVector<f64> {
-        DVector::from_vec(vec![f64::INFINITY; self.dim])
+    fn max(&self) -> OVector<f64, D> {
+        OMatrix::repeat_generic(self.mu.shape_generic().0, Const::<1>, f64::INFINITY)
     }
 }
 
-impl MeanN<DVector<f64>> for MultivariateNormal {
+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<DVector<f64>> {
-        let mut vec = vec![];
-        for elt in self.mu.clone().into_iter() {
-            vec.push(*elt);
-        }
-        Some(DVector::from_vec(vec))
+    fn mean(&self) -> Option<OVector<f64, D>> {
+        Some(self.mu.clone())
     }
 }
 
-impl VarianceN<DMatrix<f64>> for MultivariateNormal {
+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<DMatrix<f64>> {
+    fn variance(&self) -> Option<OMatrix<f64, D, D>> {
         Some(self.cov.clone())
     }
 }
 
-impl Mode<DVector<f64>> for MultivariateNormal {
+impl<D> Mode<OVector<f64, D>> for MultivariateNormal<D>
+where
+    D: Dim,
+    nalgebra::DefaultAllocator:
+        nalgebra::allocator::Allocator<f64, D> + nalgebra::allocator::Allocator<f64, D, D>,
+{
     /// Returns the mode of the multivariate normal distribution
     ///
     /// # Formula
@@ -180,12 +215,18 @@ impl Mode<DVector<f64>> for MultivariateNormal {
     /// ```
     ///
     /// where `μ` is the mean
-    fn mode(&self) -> DVector<f64> {
+    fn mode(&self) -> OVector<f64, D> {
         self.mu.clone()
     }
 }
 
-impl<'a> Continuous<&'a DVector<f64>, f64> for MultivariateNormal {
+impl<'a, D> Continuous<&'a OVector<f64, D>, f64> for MultivariateNormal<D>
+where
+    D: Dim,
+    nalgebra::DefaultAllocator: nalgebra::allocator::Allocator<f64, D>
+        + nalgebra::allocator::Allocator<f64, D, D>
+        + nalgebra::allocator::Allocator<f64, nalgebra::Const<1>, D>,
+{
     /// Calculates the probability density function for the multivariate
     /// normal distribution at `x`
     ///
@@ -197,7 +238,7 @@ impl<'a> Continuous<&'a DVector<f64>, f64> for MultivariateNormal {
     ///
     /// where `μ` is the mean, `inv(Σ)` is the precision matrix, `det(Σ)` is the determinant
     /// of the covariance matrix, and `k` is the dimension of the distribution
-    fn pdf(&self, x: &'a DVector<f64>) -> f64 {
+    fn pdf(&self, x: &'a OVector<f64, D>) -> f64 {
         let dv = x - &self.mu;
         let exp_term = -0.5
             * *(&dv.transpose() * &self.precision * &dv)
@@ -205,9 +246,10 @@ impl<'a> Continuous<&'a DVector<f64>, f64> for MultivariateNormal {
                 .unwrap();
         self.pdf_const * exp_term.exp()
     }
+
     /// Calculates the log probability density function for the multivariate
     /// normal distribution at `x`. Equivalent to pdf(x).ln().
-    fn ln_pdf(&self, x: &'a DVector<f64>) -> f64 {
+    fn ln_pdf(&self, x: &'a OVector<f64, D>) -> f64 {
         let dv = x - &self.mu;
         let exp_term = -0.5
             * *(&dv.transpose() * &self.precision * &dv)
@@ -217,7 +259,7 @@ impl<'a> Continuous<&'a DVector<f64>, f64> for MultivariateNormal {
     }
 }
 
-impl Continuous<Vec<f64>, f64> for MultivariateNormal {
+impl Continuous<Vec<f64>, f64> for MultivariateNormal<Dyn> {
     /// Calculates the probability density function for the multivariate
     /// normal distribution at `x`
     ///
@@ -232,6 +274,7 @@ impl Continuous<Vec<f64>, f64> for MultivariateNormal {
     fn pdf(&self, x: Vec<f64>) -> f64 {
         self.pdf(&DVector::from(x))
     }
+
     /// Calculates the log probability density function for the multivariate
     /// normal distribution at `x`. Equivalent to pdf(x).ln().
     fn ln_pdf(&self, x: Vec<f64>) -> f64 {
@@ -239,155 +282,322 @@ impl Continuous<Vec<f64>, f64> for MultivariateNormal {
     }
 }
 
-#[rustfmt::skip]
 #[cfg(all(test, feature = "nightly"))]
-mod tests  {
-    use crate::distribution::{Continuous, MultivariateNormal};
-    use crate::statistics::*;
-    use crate::consts::ACC;
+mod tests {
     use core::fmt::Debug;
-    use nalgebra::base::allocator::Allocator;
-    use nalgebra::{
-        DefaultAllocator, Dim, DimMin, DimName, Matrix2, Matrix3, Vector2, Vector3,
-        U1, U2,
+
+    use nalgebra::{dmatrix, dvector, matrix, vector, DimMin, OMatrix, OVector};
+
+    use crate::{
+        distribution::{Continuous, MultivariateNormal},
+        statistics::{Max, MeanN, Min, Mode, VarianceN},
     };
 
-    fn try_create(mean: Vec<f64>, covariance: Vec<f64>) -> MultivariateNormal
+    fn try_create<D>(mean: OVector<f64, D>, covariance: OMatrix<f64, D, D>) -> MultivariateNormal<D>
+    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>,
     {
-        let mvn = MultivariateNormal::new(mean, covariance);
+        let mvn = MultivariateNormal::new_from_nalgebra(mean, covariance);
         assert!(mvn.is_ok());
         mvn.unwrap()
     }
 
-    fn create_case(mean: Vec<f64>, covariance: Vec<f64>)
+    fn create_case<D>(mean: OVector<f64, D>, covariance: OMatrix<f64, D, D>)
+    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>,
     {
         let mvn = try_create(mean.clone(), covariance.clone());
-        assert_eq!(DVector::from_vec(mean.clone()), mvn.mean().unwrap());
-        assert_eq!(DMatrix::from_vec(mean.len(), mean.len(), covariance), mvn.variance().unwrap());
+        assert_eq!(mean, mvn.mean().unwrap());
+        assert_eq!(covariance, mvn.variance().unwrap());
     }
 
-    fn bad_create_case(mean: Vec<f64>, covariance: Vec<f64>)
+    fn bad_create_case<D>(mean: OVector<f64, D>, covariance: OMatrix<f64, D, D>)
+    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>,
     {
-        let mvn = MultivariateNormal::new(mean, covariance);
+        let mvn = MultivariateNormal::new_from_nalgebra(mean, covariance);
         assert!(mvn.is_err());
     }
 
-    fn test_case<T, F>(mean: Vec<f64>, covariance: Vec<f64>, expected: T, eval: F)
-    where
+    fn test_case<T, F, D>(
+        mean: OVector<f64, D>, covariance: OMatrix<f64, D, D>, expected: T, eval: F,
+    ) where
         T: Debug + PartialEq,
-        F: FnOnce(MultivariateNormal) -> T,
+        F: FnOnce(MultivariateNormal<D>) -> T,
+        D: DimMin<D, Output = D>,
+        nalgebra::DefaultAllocator: nalgebra::allocator::Allocator<f64, D>
+            + nalgebra::allocator::Allocator<f64, D, D>
+            + nalgebra::allocator::Allocator<(usize, usize), D>,
     {
         let mvn = try_create(mean, covariance);
         let x = eval(mvn);
         assert_eq!(expected, x);
     }
 
-    fn test_almost<F>(
-        mean: Vec<f64>,
-        covariance: Vec<f64>,
-        expected: f64,
-        acc: f64,
-        eval: F,
+    fn test_almost<F, D>(
+        mean: OVector<f64, D>, covariance: OMatrix<f64, D, D>, expected: f64, acc: f64, eval: F,
     ) where
-        F: FnOnce(MultivariateNormal) -> f64,
+        F: FnOnce(MultivariateNormal<D>) -> f64,
+        D: DimMin<D, Output = D>,
+        nalgebra::DefaultAllocator: nalgebra::allocator::Allocator<f64, D>
+            + nalgebra::allocator::Allocator<f64, D, D>
+            + nalgebra::allocator::Allocator<(usize, usize), D>,
     {
         let mvn = try_create(mean, covariance);
         let x = eval(mvn);
         assert_almost_eq!(expected, x, acc);
     }
 
-    use super::*;
-
-    macro_rules! dvec {
-        ($($x:expr),*) => (DVector::from_vec(vec![$($x),*]));
-    }
-
-    macro_rules! mat2 {
-        ($x11:expr, $x12:expr, $x21:expr, $x22:expr) => (DMatrix::from_vec(2,2,vec![$x11, $x12, $x21, $x22]));
-    }
-
-    // macro_rules! mat3 {
-    //     ($x11:expr, $x12:expr, $x13:expr, $x21:expr, $x22:expr, $x23:expr, $x31:expr, $x32:expr, $x33:expr) => (DMatrix::from_vec(3,3,vec![$x11, $x12, $x13, $x21, $x22, $x23, $x31, $x32, $x33]));
-    // }
-
     #[test]
     fn test_create() {
-        create_case(vec![0., 0.], vec![1., 0., 0., 1.]);
-        create_case(vec![10.,  5.], vec![2., 1., 1., 2.]);
-        create_case(vec![4., 5., 6.], vec![2., 1., 0., 1., 2., 1., 0., 1., 2.]);
-        create_case(vec![0., f64::INFINITY], vec![1., 0., 0., 1.]);
-        create_case(vec![0., 0.], vec![f64::INFINITY, 0., 0., f64::INFINITY]);
+        create_case(vector![0., 0.], matrix![1., 0.; 0., 1.]);
+        create_case(vector![10., 5.], matrix![2., 1.; 1., 2.]);
+        create_case(
+            vector![4., 5., 6.],
+            matrix![2., 1., 0.; 1., 2., 1.; 0., 1., 2.],
+        );
+        create_case(dvector![0., f64::INFINITY], dmatrix![1., 0.; 0., 1.]);
+        create_case(
+            dvector![0., 0.],
+            dmatrix![f64::INFINITY, 0.; 0., f64::INFINITY],
+        );
     }
 
     #[test]
     fn test_bad_create() {
         // Covariance not symmetric
-        bad_create_case(vec![0., 0.], vec![1., 1., 0., 1.]);
+        bad_create_case(vector![0., 0.], matrix![1., 1.; 0., 1.]);
         // Covariance not positive-definite
-        bad_create_case(vec![0., 0.], vec![1., 2., 2., 1.]);
+        bad_create_case(vector![0., 0.], matrix![1., 2.; 2., 1.]);
         // NaN in mean
-        bad_create_case(vec![0., f64::NAN], vec![1., 0., 0., 1.]);
+        bad_create_case(dvector![0., f64::NAN], dmatrix![1., 0.; 0., 1.]);
         // NaN in Covariance Matrix
-        bad_create_case(vec![0., 0.], vec![1., 0., 0., f64::NAN]);
+        bad_create_case(dvector![0., 0.], dmatrix![1., 0.; 0., f64::NAN]);
     }
 
     #[test]
     fn test_variance() {
-        let variance = |x: MultivariateNormal| x.variance().unwrap();
-        test_case(vec![0., 0.], vec![1., 0., 0., 1.], mat2![1., 0., 0., 1.], variance);
-        test_case(vec![0., 0.], vec![f64::INFINITY, 0., 0., f64::INFINITY], mat2![f64::INFINITY, 0., 0., f64::INFINITY], variance);
+        let variance = |x: MultivariateNormal<_>| x.variance().unwrap();
+        test_case(
+            vector![0., 0.],
+            matrix![1., 0.; 0., 1.],
+            matrix![1., 0.; 0., 1.],
+            variance,
+        );
+        test_case(
+            vector![0., 0.],
+            matrix![f64::INFINITY, 0.; 0., f64::INFINITY],
+            matrix![f64::INFINITY, 0.; 0., f64::INFINITY],
+            variance,
+        );
     }
 
     #[test]
     fn test_entropy() {
-        let entropy = |x: MultivariateNormal| x.entropy().unwrap();
-        test_case(vec![0., 0.], vec![1., 0., 0., 1.], 2.8378770664093453, entropy);
-        test_case(vec![0., 0.], vec![1., 0.5, 0.5, 1.], 2.694036030183455, entropy);
-        test_case(vec![0., 0.], vec![f64::INFINITY, 0., 0., f64::INFINITY], f64::INFINITY, entropy);
+        let entropy = |x: MultivariateNormal<_>| x.entropy().unwrap();
+        test_case(
+            dvector![0., 0.],
+            dmatrix![1., 0.; 0., 1.],
+            2.8378770664093453,
+            entropy,
+        );
+        test_case(
+            dvector![0., 0.],
+            dmatrix![1., 0.5; 0.5, 1.],
+            2.694036030183455,
+            entropy,
+        );
+        test_case(
+            dvector![0., 0.],
+            dmatrix![f64::INFINITY, 0.; 0., f64::INFINITY],
+            f64::INFINITY,
+            entropy,
+        );
     }
 
     #[test]
     fn test_mode() {
-        let mode = |x: MultivariateNormal| x.mode();
-        test_case(vec![0., 0.], vec![1., 0., 0., 1.], dvec![0.,  0.], mode);
-        test_case(vec![f64::INFINITY, f64::INFINITY], vec![1., 0., 0., 1.], dvec![f64::INFINITY,  f64::INFINITY], mode);
+        let mode = |x: MultivariateNormal<_>| x.mode();
+        test_case(
+            vector![0., 0.],
+            matrix![1., 0.; 0., 1.],
+            vector![0., 0.],
+            mode,
+        );
+        test_case(
+            vector![f64::INFINITY, f64::INFINITY],
+            matrix![1., 0.; 0., 1.],
+            vector![f64::INFINITY, f64::INFINITY],
+            mode,
+        );
     }
 
     #[test]
     fn test_min_max() {
-        let min = |x: MultivariateNormal| x.min();
-        let max = |x: MultivariateNormal| x.max();
-        test_case(vec![0., 0.], vec![1., 0., 0., 1.], dvec![f64::NEG_INFINITY, f64::NEG_INFINITY], min);
-        test_case(vec![0., 0.], vec![1., 0., 0., 1.], dvec![f64::INFINITY, f64::INFINITY], max);
-        test_case(vec![10., 1.], vec![1., 0., 0., 1.], dvec![f64::NEG_INFINITY, f64::NEG_INFINITY], min);
-        test_case(vec![-3., 5.], vec![1., 0., 0., 1.], dvec![f64::INFINITY, f64::INFINITY], max);
+        let min = |x: MultivariateNormal<_>| x.min();
+        let max = |x: MultivariateNormal<_>| x.max();
+        test_case(
+            dvector![0., 0.],
+            dmatrix![1., 0.; 0., 1.],
+            dvector![f64::NEG_INFINITY, f64::NEG_INFINITY],
+            min,
+        );
+        test_case(
+            dvector![0., 0.],
+            dmatrix![1., 0.; 0., 1.],
+            dvector![f64::INFINITY, f64::INFINITY],
+            max,
+        );
+        test_case(
+            dvector![10., 1.],
+            dmatrix![1., 0.; 0., 1.],
+            dvector![f64::NEG_INFINITY, f64::NEG_INFINITY],
+            min,
+        );
+        test_case(
+            dvector![-3., 5.],
+            dmatrix![1., 0.; 0., 1.],
+            dvector![f64::INFINITY, f64::INFINITY],
+            max,
+        );
     }
 
     #[test]
     fn test_pdf() {
-        let pdf = |arg: DVector<f64>| move |x: MultivariateNormal| x.pdf(&arg);
-        test_case(vec![0., 0.], vec![1., 0., 0., 1.], 0.05854983152431917, pdf(dvec![1., 1.]));
-        test_almost(vec![0., 0.], vec![1., 0., 0., 1.], 0.013064233284684921, 1e-15, pdf(dvec![1., 2.]));
-        test_almost(vec![0., 0.], vec![1., 0., 0., 1.], 1.8618676045881531e-23, 1e-35, pdf(dvec![1., 10.]));
-        test_almost(vec![0., 0.], vec![1., 0., 0., 1.], 5.920684802611216e-45, 1e-58, pdf(dvec![10., 10.]));
-        test_almost(vec![0., 0.], vec![1., 0.9, 0.9, 1.], 1.6576716577547003e-05, 1e-18, pdf(dvec![1., -1.]));
-        test_almost(vec![0., 0.], vec![1., 0.99, 0.99, 1.], 4.1970621773477824e-44, 1e-54, pdf(dvec![1., -1.]));
-        test_almost(vec![0.5, -0.2], vec![2.0, 0.3, 0.3,  0.5], 0.0013075203140666656, 1e-15, pdf(dvec![2., 2.]));
-        test_case(vec![0., 0.], vec![f64::INFINITY, 0., 0., f64::INFINITY], 0.0, pdf(dvec![10., 10.]));
-        test_case(vec![0., 0.], vec![f64::INFINITY, 0., 0., f64::INFINITY], 0.0, pdf(dvec![100., 100.]));
+        let pdf = |arg| move |x: MultivariateNormal<_>| x.pdf(&arg);
+        test_case(
+            vector![0., 0.],
+            matrix![1., 0.; 0., 1.],
+            0.05854983152431917,
+            pdf(vector![1., 1.]),
+        );
+        test_almost(
+            vector![0., 0.],
+            matrix![1., 0.; 0., 1.],
+            0.013064233284684921,
+            1e-15,
+            pdf(vector![1., 2.]),
+        );
+        test_almost(
+            vector![0., 0.],
+            matrix![1., 0.; 0., 1.],
+            1.8618676045881531e-23,
+            1e-35,
+            pdf(vector![1., 10.]),
+        );
+        test_almost(
+            vector![0., 0.],
+            matrix![1., 0.; 0., 1.],
+            5.920684802611216e-45,
+            1e-58,
+            pdf(vector![10., 10.]),
+        );
+        test_almost(
+            vector![0., 0.],
+            matrix![1., 0.9; 0.9, 1.],
+            1.6576716577547003e-05,
+            1e-18,
+            pdf(vector![1., -1.]),
+        );
+        test_almost(
+            vector![0., 0.],
+            matrix![1., 0.99; 0.99, 1.],
+            4.1970621773477824e-44,
+            1e-54,
+            pdf(vector![1., -1.]),
+        );
+        test_almost(
+            vector![0.5, -0.2],
+            matrix![2.0, 0.3; 0.3, 0.5],
+            0.0013075203140666656,
+            1e-15,
+            pdf(vector![2., 2.]),
+        );
+        test_case(
+            vector![0., 0.],
+            matrix![f64::INFINITY, 0.; 0., f64::INFINITY],
+            0.0,
+            pdf(vector![10., 10.]),
+        );
+        test_case(
+            vector![0., 0.],
+            matrix![f64::INFINITY, 0.; 0., f64::INFINITY],
+            0.0,
+            pdf(vector![100., 100.]),
+        );
     }
 
     #[test]
     fn test_ln_pdf() {
-        let ln_pdf = |arg: DVector<_>| move |x: MultivariateNormal| x.ln_pdf(&arg);
-        test_case(vec![0., 0.], vec![1., 0., 0., 1.], (0.05854983152431917f64).ln(), ln_pdf(dvec![1., 1.]));
-        test_almost(vec![0., 0.], vec![1., 0., 0., 1.], (0.013064233284684921f64).ln(), 1e-15, ln_pdf(dvec![1., 2.]));
-        test_almost(vec![0., 0.], vec![1., 0., 0., 1.], (1.8618676045881531e-23f64).ln(), 1e-15, ln_pdf(dvec![1., 10.]));
-        test_almost(vec![0., 0.], vec![1., 0., 0., 1.], (5.920684802611216e-45f64).ln(), 1e-15, ln_pdf(dvec![10., 10.]));
-        test_almost(vec![0., 0.], vec![1., 0.9, 0.9, 1.], (1.6576716577547003e-05f64).ln(), 1e-14, ln_pdf(dvec![1., -1.]));
-        test_almost(vec![0., 0.], vec![1., 0.99, 0.99, 1.], (4.1970621773477824e-44f64).ln(), 1e-12, ln_pdf(dvec![1., -1.]));
-        test_almost(vec![0.5, -0.2], vec![2.0, 0.3, 0.3, 0.5],  (0.0013075203140666656f64).ln(), 1e-15, ln_pdf(dvec![2., 2.]));
-        test_case(vec![0., 0.], vec![f64::INFINITY, 0., 0., f64::INFINITY], f64::NEG_INFINITY, ln_pdf(dvec![10., 10.]));
-        test_case(vec![0., 0.], vec![f64::INFINITY, 0., 0., f64::INFINITY], f64::NEG_INFINITY, ln_pdf(dvec![100., 100.]));
+        let ln_pdf = |arg| move |x: MultivariateNormal<_>| x.ln_pdf(&arg);
+        test_case(
+            dvector![0., 0.],
+            dmatrix![1., 0.; 0., 1.],
+            (0.05854983152431917f64).ln(),
+            ln_pdf(dvector![1., 1.]),
+        );
+        test_almost(
+            dvector![0., 0.],
+            dmatrix![1., 0.; 0., 1.],
+            (0.013064233284684921f64).ln(),
+            1e-15,
+            ln_pdf(dvector![1., 2.]),
+        );
+        test_almost(
+            dvector![0., 0.],
+            dmatrix![1., 0.; 0., 1.],
+            (1.8618676045881531e-23f64).ln(),
+            1e-15,
+            ln_pdf(dvector![1., 10.]),
+        );
+        test_almost(
+            dvector![0., 0.],
+            dmatrix![1., 0.; 0., 1.],
+            (5.920684802611216e-45f64).ln(),
+            1e-15,
+            ln_pdf(dvector![10., 10.]),
+        );
+        test_almost(
+            dvector![0., 0.],
+            dmatrix![1., 0.9; 0.9, 1.],
+            (1.6576716577547003e-05f64).ln(),
+            1e-14,
+            ln_pdf(dvector![1., -1.]),
+        );
+        test_almost(
+            dvector![0., 0.],
+            dmatrix![1., 0.99; 0.99, 1.],
+            (4.1970621773477824e-44f64).ln(),
+            1e-12,
+            ln_pdf(dvector![1., -1.]),
+        );
+        test_almost(
+            dvector![0.5, -0.2],
+            dmatrix![2.0, 0.3; 0.3, 0.5],
+            (0.0013075203140666656f64).ln(),
+            1e-15,
+            ln_pdf(dvector![2., 2.]),
+        );
+        test_case(
+            dvector![0., 0.],
+            dmatrix![f64::INFINITY, 0.; 0., f64::INFINITY],
+            f64::NEG_INFINITY,
+            ln_pdf(dvector![10., 10.]),
+        );
+        test_case(
+            dvector![0., 0.],
+            dmatrix![f64::INFINITY, 0.; 0., f64::INFINITY],
+            f64::NEG_INFINITY,
+            ln_pdf(dvector![100., 100.]),
+        );
     }
 }