From 777064ca66842a6b4b0fb540da398ec34d15768b Mon Sep 17 00:00:00 2001 From: Henry Jacobson Date: Sat, 24 Dec 2022 00:31:15 +0100 Subject: [PATCH] feat: implement multivariate students t distribution fix: Use ln_pdf_const instead of pdf_const feat: creation of multivariate normal distribution from same variables as multivariate students (for when freedom = inf) fix: use multivariate normal pdf when freedom = inf for multivariate student test: panic test for invalid pdf argument fix: tests in documentation fix: clearer function name in test fix: add documentation tests: 3d matrices test cases for pdf. Also improves documentation for multivariate t minorly. test: modify test case in multivariate_t the float chosen happens to approximate f64::LOG10_2 this leads to a linting error instead of supressing, just choosing a different value and testing against scipy also run rustfmt feat: update multivariate student to dimension generic API refactor: adds exposed Normal density for reuse in infinite DOF student distribution chore: run linting --- src/distribution/mod.rs | 2 + src/distribution/multivariate_normal.rs | 144 ++++-- src/distribution/multivariate_students_t.rs | 543 ++++++++++++++++++++ 3 files changed, 641 insertions(+), 48 deletions(-) create mode 100644 src/distribution/multivariate_students_t.rs diff --git a/src/distribution/mod.rs b/src/distribution/mod.rs index dea1f9af..b146f9d4 100644 --- a/src/distribution/mod.rs +++ b/src/distribution/mod.rs @@ -27,6 +27,7 @@ pub use self::laplace::Laplace; pub use self::log_normal::LogNormal; pub use self::multinomial::Multinomial; pub use self::multivariate_normal::MultivariateNormal; +pub use self::multivariate_students_t::MultivariateStudent; pub use self::negative_binomial::NegativeBinomial; pub use self::normal::Normal; pub use self::pareto::Pareto; @@ -60,6 +61,7 @@ mod laplace; mod log_normal; mod multinomial; mod multivariate_normal; +mod multivariate_students_t; mod negative_binomial; mod normal; mod pareto; diff --git a/src/distribution/multivariate_normal.rs b/src/distribution/multivariate_normal.rs index 9949d676..61ccfa57 100644 --- a/src/distribution/multivariate_normal.rs +++ b/src/distribution/multivariate_normal.rs @@ -1,12 +1,84 @@ use crate::distribution::Continuous; use crate::distribution::Normal; use crate::statistics::{Max, MeanN, Min, Mode, VarianceN}; -use crate::{Result, StatsError}; +use crate::StatsError; use nalgebra::{Cholesky, Const, DMatrix, DVector, Dim, DimMin, Dyn, OMatrix, OVector}; use rand::Rng; 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 +pub(super) fn density_normalization_and_exponential( + mu: &OVector, + cov: &OMatrix, + precision: &OMatrix, + x: &OVector, +) -> std::result::Result<(f64, f64), StatsError> +where + D: DimMin, + nalgebra::DefaultAllocator: nalgebra::allocator::Allocator + + nalgebra::allocator::Allocator + + nalgebra::allocator::Allocator<(usize, usize), D>, +{ + Ok(( + 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 +#[inline] +pub(super) fn density_distribution_exponential( + mu: &OVector, + precision: &OMatrix, + x: &OVector, +) -> std::result::Result +where + D: Dim, + nalgebra::DefaultAllocator: + nalgebra::allocator::Allocator + nalgebra::allocator::Allocator, +{ + if x.shape_generic().0 != precision.shape_generic().0 + || x.shape_generic().0 != mu.shape_generic().0 + || !precision.is_square() + { + return Err(StatsError::ContainersMustBeSameLength); + } + let dv = x - mu; + let exp_term: f64 = -0.5 * (precision * &dv).dot(&dv); + Ok(exp_term) + // TODO update to dimension mismatch error +} + +/// computes the argument of the normalization term in the normal distribution +/// # Errors +/// will error on dimension mismatch +#[inline] +pub(super) fn density_distribution_pdf_const( + mu: &OVector, + cov: &OMatrix, +) -> std::result::Result +where + D: DimMin, + nalgebra::DefaultAllocator: nalgebra::allocator::Allocator + + nalgebra::allocator::Allocator + + nalgebra::allocator::Allocator<(usize, usize), D>, +{ + if cov.shape_generic().0 != mu.shape_generic().0 || !cov.is_square() { + return Err(StatsError::ContainersMustBeSameLength); + } + let cov_det = cov.determinant(); + Ok(((2. * PI).powi(mu.nrows() as i32) * cov_det.abs()) + .recip() + .sqrt()) +} + /// Implements the [Multivariate Normal](https://en.wikipedia.org/wiki/Multivariate_normal_distribution) /// distribution using the "nalgebra" crate for matrix operations /// @@ -44,7 +116,7 @@ impl MultivariateNormal { /// /// Returns an error if the given covariance matrix is not /// symmetric or positive-definite - pub fn new(mean: Vec, cov: Vec) -> Result { + pub fn new(mean: Vec, cov: Vec) -> Result { let mean = DVector::from_vec(mean); let cov = DMatrix::from_vec(mean.len(), mean.len(), cov); MultivariateNormal::new_from_nalgebra(mean, cov) @@ -66,7 +138,10 @@ where /// /// Returns an error if the given covariance matrix is not /// symmetric or positive-definite - pub fn new_from_nalgebra(mean: OVector, cov: OMatrix) -> Result { + pub fn new_from_nalgebra( + mean: OVector, + cov: OMatrix, + ) -> Result { // 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 @@ -77,10 +152,6 @@ where { return Err(StatsError::BadParams); } - let cov_det = cov.determinant(); - let pdf_const = ((2. * PI).powi(mean.nrows() as i32) * cov_det.abs()) - .recip() - .sqrt(); // Store the Cholesky decomposition of the covariance matrix // for sampling match Cholesky::new(cov.clone()) { @@ -88,11 +159,11 @@ where Some(cholesky_decomp) => { let precision = cholesky_decomp.inverse(); Ok(MultivariateNormal { + pdf_const: density_distribution_pdf_const(&mean, &cov).unwrap(), cov_chol_decomp: cholesky_decomp.unpack(), mu: mean, cov, precision, - pdf_const, }) } } @@ -231,9 +302,8 @@ where impl<'a, D> Continuous<&'a OVector, f64> for MultivariateNormal where D: Dim, - nalgebra::DefaultAllocator: nalgebra::allocator::Allocator - + nalgebra::allocator::Allocator - + nalgebra::allocator::Allocator, D>, + nalgebra::DefaultAllocator: + nalgebra::allocator::Allocator + nalgebra::allocator::Allocator, { /// Calculates the probability density function for the multivariate /// normal distribution at `x` @@ -246,47 +316,18 @@ where /// /// 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 OVector) -> f64 { - let dv = x - &self.mu; - let exp_term = -0.5 - * *(&dv.transpose() * &self.precision * &dv) - .get((0, 0)) - .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 OVector) -> f64 { - let dv = x - &self.mu; - let exp_term = -0.5 - * *(&dv.transpose() * &self.precision * &dv) - .get((0, 0)) - .unwrap(); - self.pdf_const.ln() + exp_term - } -} - -impl Continuous, f64> for MultivariateNormal { - /// Calculates the probability density function for the multivariate - /// normal distribution at `x` - /// - /// # Formula - /// - /// ```text - /// (2 * π) ^ (-k / 2) * det(Σ) ^ (1 / 2) * e ^ ( -(1 / 2) * transpose(x - μ) * inv(Σ) * (x - μ)) - /// ``` - /// - /// 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: Vec) -> f64 { - self.pdf(&DVector::from(x)) + fn pdf(&self, x: &OVector) -> f64 { + self.pdf_const + * density_distribution_exponential(&self.mu, &self.precision, x) + .unwrap() + .exp() } /// 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 { - self.pdf(&DVector::from(x)) + fn ln_pdf(&self, x: &OVector) -> f64 { + self.pdf_const.ln() + + density_distribution_exponential(&self.mu, &self.precision, x).unwrap() } } @@ -609,4 +650,11 @@ mod tests { ln_pdf(dvector![100., 100.]), ); } + + #[test] + #[should_panic] + fn test_pdf_mismatched_arg_size() { + let mvn = MultivariateNormal::new(vec![0., 0.], vec![1., 0., 0., 1.,]).unwrap(); + mvn.pdf(&vec![1.].into()); // x.size != mu.size + } } diff --git a/src/distribution/multivariate_students_t.rs b/src/distribution/multivariate_students_t.rs new file mode 100644 index 00000000..4a3cd38f --- /dev/null +++ b/src/distribution/multivariate_students_t.rs @@ -0,0 +1,543 @@ +use crate::distribution::Continuous; +use crate::distribution::{ChiSquared, Normal}; +use crate::function::gamma; +use crate::statistics::{Max, MeanN, Min, Mode, VarianceN}; +use crate::{Result, StatsError}; +use nalgebra::{Cholesky, Const, DMatrix, Dim, DimMin, Dyn, OMatrix, OVector}; +use rand::Rng; +use std::f64::consts::PI; + +/// Implements the [Multivariate Student's t-distribution](https://en.wikipedia.org/wiki/Multivariate_t-distribution) +/// distribution using the "nalgebra" crate for matrix operations. +/// +/// Assumes all the marginal distributions have the same degree of freedom, ν. +/// +/// # Examples +/// +/// ``` +/// use statrs::distribution::{MultivariateStudent, Continuous}; +/// use nalgebra::{DVector, DMatrix}; +/// use statrs::statistics::{MeanN, VarianceN}; +/// +/// let mvs = MultivariateStudent::new(vec![0., 0.], vec![1., 0., 0., 1.], 4.).unwrap(); +/// 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); +/// ``` +#[derive(Debug, Clone, PartialEq)] +pub struct MultivariateStudent +where + D: Dim, + nalgebra::DefaultAllocator: + nalgebra::allocator::Allocator + nalgebra::allocator::Allocator, +{ + scale_chol_decomp: OMatrix, + location: OVector, + scale: OMatrix, + freedom: f64, + precision: OMatrix, + ln_pdf_const: f64, +} + +impl MultivariateStudent { + /// Constructs a new multivariate students t distribution with a location of `location`, + /// scale matrix `scale` and `freedom` degrees of freedom. + /// + /// # Errors + /// + /// Returns `StatsError::BadParams` if the scale matrix is not symmetric-positive + /// definite and `StatsError::ArgMustBePositive` if freedom is non-positive. + pub fn new(location: Vec, scale: Vec, freedom: f64) -> Result { + let dim = location.len(); + Self::new_from_nalgebra(location.into(), DMatrix::from_vec(dim, dim, scale), freedom) + } + + /// Returns the dimension of the distribution. + pub fn dim(&self) -> usize { + self.location.len() + } +} + +impl MultivariateStudent +where + D: DimMin, + nalgebra::DefaultAllocator: nalgebra::allocator::Allocator + + nalgebra::allocator::Allocator + + nalgebra::allocator::Allocator<(usize, usize), D>, +{ + pub fn new_from_nalgebra( + location: OVector, + scale: OMatrix, + freedom: f64, + ) -> Result { + let dim = location.len(); + + // Check that the provided scale matrix is symmetric + if scale.lower_triangle() != scale.upper_triangle().transpose() + // Check that mean and covariance do not contain NaN + || location.iter().any(|f| f.is_nan()) + || scale.iter().any(|f| f.is_nan()) + // Check that the dimensions match + || location.nrows() != scale.nrows() || scale.nrows() != scale.ncols() + // Check that the degrees of freedom is not NaN + || freedom.is_nan() + { + return Err(StatsError::BadParams); + } + // Check that degrees of freedom is positive + if freedom <= 0. { + return Err(StatsError::ArgMustBePositive( + "Degrees of freedom must be positive", + )); + } + + let scale_det = scale.determinant(); + let ln_pdf_const = gamma::ln_gamma(0.5 * (freedom + dim as f64)) + - gamma::ln_gamma(0.5 * freedom) + - 0.5 * (dim as f64) * (freedom * PI).ln() + - 0.5 * scale_det.ln(); + + match Cholesky::new(scale.clone()) { + None => Err(StatsError::BadParams), // Scale matrix is not positive definite + Some(cholesky_decomp) => { + let precision = cholesky_decomp.inverse(); + Ok(MultivariateStudent { + scale_chol_decomp: cholesky_decomp.unpack(), + location, + scale, + freedom, + precision, + ln_pdf_const, + }) + } + } + } + + /// Returns the cholesky decomposiiton matrix of the scale matrix. + /// + /// Returns A where Σ = AAᵀ. + pub fn scale_chol_decomp(&self) -> &OMatrix { + &self.scale_chol_decomp + } + + /// Returns the location of the distribution. + pub fn location(&self) -> &OVector { + &self.location + } + + /// Returns the scale matrix of the distribution. + pub fn scale(&self) -> &OMatrix { + &self.scale + } + + /// Returns the degrees of freedom of the distribution. + pub fn freedom(&self) -> f64 { + self.freedom + } + + /// Returns the inverse of the cholesky decomposition matrix. + pub fn precision(&self) -> &OMatrix { + &self.precision + } + + /// Returns the logarithmed constant part of the probability + /// distribution function. + pub fn ln_pdf_const(&self) -> f64 { + self.ln_pdf_const + } +} + +impl ::rand::distributions::Distribution> for MultivariateStudent +where + D: Dim, + nalgebra::DefaultAllocator: + nalgebra::allocator::Allocator + nalgebra::allocator::Allocator, +{ + /// Samples from the multivariate student distribution + /// + /// # Formula + /// + /// ```math + /// W ⋅ L ⋅ Z + μ + /// ``` + /// + /// where `W` has √(ν/Sν) distribution, Sν has Chi-squared + /// distribution with ν degrees of freedom, + /// `L` is the Cholesky decomposition of the scale matrix, + /// `Z` is a vector of normally distributed random variables, and + /// `μ` is the location vector + fn sample(&self, rng: &mut R) -> OVector { + let d = Normal::new(0., 1.).unwrap(); + let s = ChiSquared::new(self.freedom).unwrap(); + let w = (self.freedom / s.sample(rng)).sqrt(); + let (r, c) = self.location.shape_generic(); + let z = OVector::::from_distribution_generic(r, c, &d, rng); + (w * &self.scale_chol_decomp * z) + &self.location + } +} + +impl Min> for MultivariateStudent +where + D: Dim, + nalgebra::DefaultAllocator: + nalgebra::allocator::Allocator + nalgebra::allocator::Allocator, +{ + /// Returns the minimum value in the domain of the + /// multivariate normal distribution represented by a real vector + fn min(&self) -> OVector { + OMatrix::repeat_generic( + self.location.shape_generic().0, + Const::<1>, + f64::NEG_INFINITY, + ) + } +} + +impl Max> for MultivariateStudent +where + D: Dim, + nalgebra::DefaultAllocator: + nalgebra::allocator::Allocator + nalgebra::allocator::Allocator, +{ + /// Returns the minimum value in the domain of the + /// multivariate normal distribution represented by a real vector + fn max(&self) -> OVector { + OMatrix::repeat_generic(self.location.shape_generic().0, Const::<1>, f64::INFINITY) + } +} + +impl MeanN> for MultivariateStudent +where + D: Dim, + nalgebra::DefaultAllocator: + nalgebra::allocator::Allocator + nalgebra::allocator::Allocator, +{ + /// 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> { + if self.freedom > 1. { + Some(self.location.clone()) + } else { + None + } + } +} + +impl VarianceN> for MultivariateStudent +where + D: Dim, + nalgebra::DefaultAllocator: + nalgebra::allocator::Allocator + nalgebra::allocator::Allocator, +{ + /// 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> { + if self.freedom > 2. { + Some(self.scale.clone() * self.freedom / (self.freedom - 2.)) + } else { + None + } + } +} + +impl Mode> for MultivariateStudent +where + D: Dim, + nalgebra::DefaultAllocator: + nalgebra::allocator::Allocator + nalgebra::allocator::Allocator, +{ + /// Returns the mode of the multivariate student distribution. + /// + /// # Formula + /// + /// ```math + /// μ + /// ``` + /// + /// where `μ` is the location. + fn mode(&self) -> OVector { + self.location.clone() + } +} + +impl<'a, D> Continuous<&'a OVector, f64> for MultivariateStudent +where + D: Dim + DimMin, + nalgebra::DefaultAllocator: nalgebra::allocator::Allocator + + nalgebra::allocator::Allocator + + nalgebra::allocator::Allocator<(usize, usize), D>, +{ + /// Calculates the probability density function for the multivariate. + /// student distribution at `x`. + /// + /// # Formula + /// + /// ```math + /// [Γ(ν+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. + fn pdf(&self, x: &'a OVector) -> f64 { + if self.freedom.is_infinite() { + use super::multivariate_normal::density_normalization_and_exponential; + let (pdf_const, exp_arg) = density_normalization_and_exponential( + &self.location, + &self.scale, + &self.precision, + x, + ) + .unwrap(); + return pdf_const * exp_arg.exp(); + } + + let dv = x - &self.location; + let exp_arg: f64 = (&self.precision * &dv).dot(&dv); + let base_term = 1. + exp_arg / self.freedom; + self.ln_pdf_const.exp() * base_term.powf(-(self.freedom + self.location.len() as f64) / 2.) + } + + /// Calculates the log probability density function for the multivariate + /// student distribution at `x`. Equivalent to pdf(x).ln(). + fn ln_pdf(&self, x: &'a OVector) -> f64 { + if self.freedom.is_infinite() { + use super::multivariate_normal::density_normalization_and_exponential; + let (pdf_const, exp_arg) = density_normalization_and_exponential( + &self.location, + &self.scale, + &self.precision, + x, + ) + .unwrap(); + return pdf_const.ln() + exp_arg; + } + + let dv = x - &self.location; + let exp_arg: f64 = (&self.precision * &dv).dot(&dv); + let base_term = 1. + exp_arg / self.freedom; + self.ln_pdf_const - (self.freedom + self.location.len() as f64) / 2. * base_term.ln() + } +} + +#[rustfmt::skip] +#[cfg(test)] +mod tests { + use core::fmt::Debug; + + use approx::RelativeEq; + use nalgebra::{Dyn, DMatrix, DVector}; + + use crate::{ + distribution::{Continuous, MultivariateStudent, MultivariateNormal}, + statistics::{Max, MeanN, Min, Mode, VarianceN}, + }; + + fn try_create(location: Vec, scale: Vec, freedom: f64) -> MultivariateStudent + { + let mvs = MultivariateStudent::new(location, scale, freedom); + assert!(mvs.is_ok()); + mvs.unwrap() + } + + fn create_case(location: Vec, scale: Vec, freedom: f64) + { + let mvs = try_create(location.clone(), scale.clone(), freedom); + assert_eq!(DMatrix::from_vec(location.len(), location.len(), scale), mvs.scale); + assert_eq!(DVector::from_vec(location), mvs.location); + } + + fn bad_create_case(location: Vec, scale: Vec, freedom: f64) + { + let mvs = MultivariateStudent::new(location, scale, freedom); + assert!(mvs.is_err()); + } + + fn test_case(location: Vec, scale: Vec, freedom: f64, expected: T, eval: F) + where + T: Debug + PartialEq, + F: FnOnce(MultivariateStudent) -> T, + { + let mvs = try_create(location, scale, freedom); + let x = eval(mvs); + assert_eq!(expected, x); + } + + fn test_almost( + location: Vec, + scale: Vec, + freedom: f64, + expected: f64, + acc: f64, + eval: F, + ) where + F: FnOnce(MultivariateStudent) -> f64, + { + let mvs = try_create(location, scale, freedom); + let x = eval(mvs); + assert_almost_eq!(expected, x, acc); + } + + fn test_almost_multivariate_normal( + location: Vec, + scale: Vec, + freedom: f64, + acc: f64, + x: DVector, + eval_mvs: F1, + eval_mvn: F2, + ) where + F1: FnOnce(MultivariateStudent, DVector) -> f64, + F2: FnOnce(MultivariateNormal, DVector) -> f64, + { + let mvs = try_create(location.clone(), scale.clone(), freedom); + let mvn0 = MultivariateNormal::new(location, scale); + assert!(mvn0.is_ok()); + let mvn = mvn0.unwrap(); + let mvs_x = eval_mvs(mvs, x.clone()); + let mvn_x = eval_mvn(mvn, x.clone()); + assert!(mvs_x.relative_eq(&mvn_x, acc, acc), "mvn: {mvn_x} =/=\nmvs: {mvs_x}"); + // assert_relative_eq!(mvs_x, mvn_x, acc); + } + + + 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.], 1.); + create_case(vec![10., 5.], vec![2., 1., 1., 2.], 3.); + create_case(vec![4., 5., 6.], vec![2., 1., 0., 1., 2., 1., 0., 1., 2.], 14.); + create_case(vec![0., f64::INFINITY], vec![1., 0., 0., 1.], f64::INFINITY); + create_case(vec![0., 0.], vec![f64::INFINITY, 0., 0., f64::INFINITY], 0.1); + } + + #[test] + fn test_bad_create() { + // scale not symmetric. + bad_create_case(vec![0., 0.], vec![1., 1., 0., 1.], 1.); + // scale not positive-definite. + bad_create_case(vec![0., 0.], vec![1., 2., 2., 1.], 1.); + // NaN in location. + bad_create_case(vec![0., f64::NAN], vec![1., 0., 0., 1.], 1.); + // NaN in scale Matrix. + bad_create_case(vec![0., 0.], vec![1., 0., 0., f64::NAN], 1.); + // NaN in freedom. + bad_create_case(vec![0., 0.], vec![1., 0., 0., 1.], f64::NAN); + // Non-positive freedom. + bad_create_case(vec![0., 0.], vec![1., 0., 0., 1.], 0.); + bad_create_case(vec![0., 0.], vec![1., 0., 0., 1.], -1.); + } + + #[test] + fn test_variance() { + let variance = |x: MultivariateStudent| x.variance().unwrap(); + test_case(vec![0., 0.], vec![1., 0., 0., 1.], 3., 3. * mat2![1., 0., 0., 1.], variance); + test_case(vec![0., 0.], vec![f64::INFINITY, 0., 0., f64::INFINITY], 3., mat2![f64::INFINITY, 0., 0., f64::INFINITY], variance); + } + + // Variance is only defined for freedom > 2. + #[test] + fn test_bad_variance() { + let variance = |x: MultivariateStudent| x.variance(); + test_case(vec![0., 0.], vec![1., 0., 0., 1.], 2., None, variance); + } + + #[test] + fn test_mode() { + let mode = |x: MultivariateStudent| x.mode(); + test_case(vec![0., 0.], vec![1., 0., 0., 1.], 1., dvec![0., 0.], mode); + test_case(vec![f64::INFINITY, f64::INFINITY], vec![1., 0., 0., 1.], 1., dvec![f64::INFINITY, f64::INFINITY], mode); + } + + #[test] + fn test_mean() { + let mean = |x: MultivariateStudent| x.mean().unwrap(); + test_case(vec![0., 0.], vec![1., 0., 0., 1.], 2., dvec![0., 0.], mean); + test_case(vec![-1., 1., 3.], vec![1., 0., 0.5, 0., 2.0, 0., 0.5, 0., 3.0], 2., dvec![-1., 1., 3.], mean); + } + + // Mean is only defined if freedom > 1. + #[test] + fn test_bad_mean() { + let mean = |x: MultivariateStudent| x.mean(); + test_case(vec![0., 0.], vec![1., 0., 0., 1.], 1., None, mean); + } + + #[test] + fn test_min_max() { + let min = |x: MultivariateStudent| x.min(); + let max = |x: MultivariateStudent| x.max(); + test_case(vec![0., 0.], vec![1., 0., 0., 1.], 1., dvec![f64::NEG_INFINITY, f64::NEG_INFINITY], min); + test_case(vec![0., 0.], vec![1., 0., 0., 1.], 1., dvec![f64::INFINITY, f64::INFINITY], max); + test_case(vec![10., 1.], vec![1., 0., 0., 1.], 1., dvec![f64::NEG_INFINITY, f64::NEG_INFINITY], min); + test_case(vec![-3., 5.], vec![1., 0., 0., 1.], 1., dvec![f64::INFINITY, f64::INFINITY], max); + } + + #[test] + fn test_pdf() { + let pdf = |arg: DVector| move |x: MultivariateStudent| x.pdf(&arg); + test_almost(vec![0., 0.], vec![1., 0., 0., 1.], 4., 0.047157020175376416, 1e-15, pdf(dvec![1., 1.])); + test_almost(vec![0., 0.], vec![1., 0., 0., 1.], 4., 0.013972450422333741737457302178882, 1e-15, pdf(dvec![1., 2.])); + test_almost(vec![0., 0.], vec![1., 0., 0., 1.], 2., 0.012992240252399619, 1e-17, pdf(dvec![1., 2.])); + test_almost(vec![2., 1.], vec![5., 0., 0., 1.], 2.5, 2.639780816598878e-5, 1e-19, pdf(dvec![1., 10.])); + test_almost(vec![-1., 0.], vec![2., 1., 1., 6.], 1.5, 6.438051574348526e-5, 1e-19, pdf(dvec![10., 10.])); + // These three are crossed checked against both python's scipy.multivariate_t.pdf and octave's mvtpdf. + test_almost(vec![-1., 1., 50.], vec![1., 0.5, 0.25, 0.5, 1., -0.1, 0.25, -0.1, 1.], 8., 6.960998836915657e-16, 1e-30, pdf(dvec![0.9718, 0.1298, 0.8134])); + test_almost(vec![-1., 1., 50.], vec![1., 0.5, 0.25, 0.5, 1., -0.1, 0.25, -0.1, 1.], 8., 7.369987979187023e-16, 1e-30, pdf(dvec![0.4922, 0.5522, 0.7185])); + test_almost(vec![-1., 1., 50.], vec![1., 0.5, 0.25, 0.5, 1., -0.1, 0.25, -0.1, 1.], 8.,6.951631724511314e-16, 1e-30, pdf(dvec![0.3020, 0.1491, 0.5008])); + test_case(vec![-1., 0.], vec![f64::INFINITY, 0., 0., f64::INFINITY], 10., 0., pdf(dvec![10., 10.])); + } + + #[test] + fn test_ln_pdf() { + let ln_pdf = |arg: DVector| move |x: MultivariateStudent| x.ln_pdf(&arg); + test_almost(vec![0., 0.], vec![1., 0., 0., 1.], 4., -3.0542723907338383, 1e-14, ln_pdf(dvec![1., 1.])); + test_almost(vec![0., 0.], vec![1., 0., 0., 1.], 2., -4.3434030034000815, 1e-14, ln_pdf(dvec![1., 2.])); + test_almost(vec![2., 1.], vec![5., 0., 0., 1.], 2.5, -10.542229575274265, 1e-14, ln_pdf(dvec![1., 10.])); + test_almost(vec![-1., 0.], vec![2., 1., 1., 6.], 1.5, -9.650699521198622, 1e-14, ln_pdf(dvec![10., 10.])); + // test_case(vec![-1., 0.], vec![f64::INFINITY, 0., 0., f64::INFINITY], 10., f64::NEG_INFINITY, ln_pdf(dvec![10., 10.])); + } + + #[test] + fn test_pdf_freedom_large() { + let pdf_mvs = |mv: MultivariateStudent, arg: DVector| mv.pdf(&arg); + let pdf_mvn = |mv: MultivariateNormal, arg: DVector| mv.pdf(&arg); + test_almost_multivariate_normal(vec![0., 0.,], vec![1., 0., 0., 1.], 1e5, 1e-6, dvec![1., 1.], pdf_mvs, pdf_mvn); + test_almost_multivariate_normal(vec![0., 0.,], vec![1., 0., 0., 1.], 1e10, 1e-7, dvec![1., 1.], pdf_mvs, pdf_mvn); + test_almost_multivariate_normal(vec![0., 0.,], vec![1., 0., 0., 1.], f64::INFINITY, 1e-300, dvec![1., 1.], pdf_mvs, pdf_mvn); + test_almost_multivariate_normal(vec![5., -1.,], vec![1., 0.99, 0.99, 1.], f64::INFINITY, 1e-300, dvec![5., 1.], pdf_mvs, pdf_mvn); + } + #[test] + fn test_ln_pdf_freedom_large() { + let pdf_mvs = |mv: MultivariateStudent, arg: DVector| mv.ln_pdf(&arg); + let pdf_mvn = |mv: MultivariateNormal, arg: DVector| mv.ln_pdf(&arg); + test_almost_multivariate_normal(vec![0., 0.,], vec![1., 0., 0., 1.], 1e5, 1e-5, dvec![1., 1.], pdf_mvs, pdf_mvn); + test_almost_multivariate_normal(vec![0., 0.,], vec![1., 0., 0., 1.], 1e10, 5e-6, dvec![1., 1.], pdf_mvs, pdf_mvn); + test_almost_multivariate_normal(vec![0., 0.,], vec![1., 0., 0., 1.], f64::INFINITY, 1e-300, dvec![1., 1.], pdf_mvs, pdf_mvn); + test_almost_multivariate_normal(vec![0., 0.,], vec![1., 0.99, 0.99, 1.], f64::INFINITY, 1e-300, dvec![1., 1.], pdf_mvs, pdf_mvn); + } +}