From 91036b9b2e28c05d1b4230393d13420365166b05 Mon Sep 17 00:00:00 2001 From: Henry Jacobson Date: Sat, 24 Dec 2022 00:31:15 +0100 Subject: [PATCH 1/9] feat: implement multivariate students t distribution --- src/distribution/mod.rs | 2 + src/distribution/multivariate_students_t.rs | 394 ++++++++++++++++++++ 2 files changed, 396 insertions(+) create mode 100644 src/distribution/multivariate_students_t.rs diff --git a/src/distribution/mod.rs b/src/distribution/mod.rs index 23ae7ec5..6f95bebd 100644 --- a/src/distribution/mod.rs +++ b/src/distribution/mod.rs @@ -26,6 +26,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; @@ -59,6 +60,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_students_t.rs b/src/distribution/multivariate_students_t.rs new file mode 100644 index 00000000..10e3be6e --- /dev/null +++ b/src/distribution/multivariate_students_t.rs @@ -0,0 +1,394 @@ +use crate::distribution::Continuous; +use crate::distribution::{Normal, ChiSquared}; +use crate::distribution::MultivariateNormal; +use crate::function::gamma; +use crate::statistics::{Max, MeanN, Min, Mode, VarianceN}; +use crate::{Result, StatsError}; +use nalgebra::Cholesky; +use nalgebra::{DMatrix, DVector}; +use rand::Rng; +use std::f64::consts::{E, 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, 4. * vec![1., 0., 0., 1.])); +/// assert_eq!(mvs.pdf(&DVector::from_vec(vec![1., 1.])), 0.047157020175376416); +/// ``` +#[derive(Debug, Clone, PartialEq)] +pub struct MultivariateStudent { + dim: usize, + scale_chol_decomp: DMatrix, + location: DVector, + scale: DMatrix, + freedom: f64, + precision: DMatrix, + pdf_const: f64, +} + +impl MultivariateStudent { + pub fn new(location: Vec, scale: Vec, freedom: f64) -> Result { + let dim = location.len(); + let location = DVector::from_vec(location); + let scale = DMatrix::from_vec(dim, dim, scale); + + // 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 pdf_const = gamma::gamma((freedom + (dim as f64)) / 2.) * + (gamma::gamma(freedom / 2.) * + (freedom.powi(dim as i32) * PI.powi(dim as i32) * scale_det.abs()) + .sqrt()) + .recip(); + + match Cholesky::new(scale.clone()) { + None => Err(StatsError::BadParams), + Some(cholesky_decomp) => { + let precision = cholesky_decomp.inverse(); + Ok(MultivariateStudent { + dim, + scale_chol_decomp: cholesky_decomp.unpack(), + location, + scale, + freedom, + precision, + pdf_const, + }) + } + } + } +} + +impl ::rand::distributions::Distribution> for MultivariateStudent { + /// Samples from the multivariate student distribution + /// + /// # Formula + /// + /// 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) -> DVector { + let d = Normal::new(0., 1.).unwrap(); + let s = ChiSquared::new(self.freedom).unwrap(); + let w = (self.freedom / s.sample(rng)).sqrt(); + let z = DVector::::from_distribution(self.dim, &d, rng); + (w * &self.scale_chol_decomp * z) + &self.location + } +} + +impl Min> for MultivariateStudent { + /// Returns the minimum value in the domain of the + /// multivariate student's t distribution represented by a real vector + fn min(&self) -> DVector { + DVector::from_vec(vec![f64::NEG_INFINITY; self.dim]) + } +} + +impl Max> for MultivariateStudent { + /// Returns the maximum value in the domain of the + /// multivariate student distribution represented by a real vector + fn max(&self) -> DVector { + DVector::from_vec(vec![f64::INFINITY; self.dim]) + } +} + +impl MeanN> for MultivariateStudent { + /// 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. { + let mut vec = vec![]; + for elt in self.location.clone().into_iter() { + vec.push(*elt); + } + Some(DVector::from_vec(vec)) + } else { + None + } + } +} + +impl VarianceN> for MultivariateStudent { + /// Returns the covariance matrix of the multivariate student distribution + fn variance(&self) -> Option> { + if self.freedom > 2. { + Some(self.scale.clone() * self.freedom / (self.freedom - 2.)) + } else { + None + } + } +} + +impl Mode> for MultivariateStudent { + /// Returns the mode of the multivariate student distribution + /// + /// # Formula + /// + /// ```ignore + /// μ + /// ``` + /// + /// where `μ` is the location + fn mode(&self) -> DVector { + self.location.clone() + } +} + +impl<'a> Continuous<&'a DVector, f64> for MultivariateStudent { + /// Calculates the probability density function for the multivariate + /// student distribution at `x` + /// + /// # Formula + /// + /// ```ignore + /// Gamma[(ν+p)/2] / [Gamma(ν/2) ((ν * π)^p det(Σ))^(1 / 2)] * [1 + 1/ν transpose(x - μ) inv(Σ) (x - μ)]^(-(ν+p)/2) + /// ``` + /// + /// where `ν` is the degrees of freedom, `μ` is the mean, `Gamma` + /// 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 + /// + /// TODO: Make this converge for large degrees of freedom + /// Current commented code beneath fails since `MultivariateNormal::new` accepts Vec and + /// not DVector or DMatrix. Should implement that instead of changing back to Vec, or + /// even have a constructor `MultivariateNormal::from_student`. + fn pdf(&self, x: &'a DVector) -> f64 { + // if self.freedom == f64::INFINITY { + // let mvn = MultivariateNormal::new(self.location, self.scale).unwrap(); + // return mvn.pdf(x); + // } + let dv = x - &self.location; + let base_term = 1. + 1. / self.freedom + * *(&dv.transpose() * &self.precision * &dv) + .get((0, 0)) + .unwrap(); + self.pdf_const * base_term.powf(-(self.freedom + self.dim 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 DVector) -> f64 { + let dv = x - &self.location; + let base_term = 1. + 1. / self.freedom + * *(&dv.transpose() * &self.precision * &dv) + .get((0, 0)) + .unwrap(); + self.pdf_const.ln() - (self.freedom + self.dim as f64) / 2. * base_term.ln() + } +} + +#[rustfmt::skip] +#[cfg(all(test, feature = "nightly"))] +mod tests { + use crate::distribution::MultivariateNormal; + use core::fmt::Debug; + + 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!(DVector::from_vec(location.clone()), mvs.location); + assert_eq!(DMatrix::from_vec(location.len(), location.len(), scale), mvs.scale); + } + + 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, + 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(), f64::INFINITY); + 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_almost_eq!(mvs_x, mvn_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.], 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.); + } + + #[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); + } + + #[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] + 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.], 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.])); + test_case(vec![-1., 0.], vec![f64::INFINITY, 0., 0., f64::INFINITY], 10., 0., pdf(dvec![10., 10.])); + } + + #[test] + fn test_ln_pdf() { + let 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, pdf(dvec![1., 1.])); + test_almost(vec![0., 0.], vec![1., 0., 0., 1.], 2., -4.3434030034000815, 1e-14, pdf(dvec![1., 2.])); + test_almost(vec![2., 1.], vec![5., 0., 0., 1.], 2.5, -10.542229575274265, 1e-14, pdf(dvec![1., 10.])); + test_almost(vec![-1., 0.], vec![2., 1., 1., 6.], 1.5, -9.650699521198622, 1e-14, pdf(dvec![10., 10.])); + test_case(vec![-1., 0.], vec![f64::INFINITY, 0., 0., f64::INFINITY], 10., f64::NEG_INFINITY, pdf(dvec![10., 10.])); + } + + // TODO: These tests fail because inf degrees of freedom give NaN + #[test] + fn test_pdf_freedom_inf() { + 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.], 1e-14, dvec![1., 1.], pdf_mvs, pdf_mvn) + } +} From 638f2831787635d47cceea9cf5556870a8d1739f Mon Sep 17 00:00:00 2001 From: Henry Jacobson Date: Sun, 25 Dec 2022 13:08:31 +0100 Subject: [PATCH 2/9] fix: Use ln_pdf_const instead of pdf_const --- src/distribution/multivariate_students_t.rs | 58 +++++++++++++-------- 1 file changed, 35 insertions(+), 23 deletions(-) diff --git a/src/distribution/multivariate_students_t.rs b/src/distribution/multivariate_students_t.rs index 10e3be6e..f1fa2feb 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::distribution::{Normal, ChiSquared}; use crate::distribution::MultivariateNormal; +use crate::distribution::{ChiSquared, Normal}; use crate::function::gamma; use crate::statistics::{Max, MeanN, Min, Mode, VarianceN}; use crate::{Result, StatsError}; @@ -34,13 +34,20 @@ pub struct MultivariateStudent { scale: DMatrix, freedom: f64, precision: DMatrix, - pdf_const: f64, + 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(); - let location = DVector::from_vec(location); + let location = DVector::from_vec(location); let scale = DMatrix::from_vec(dim, dim, scale); // Check that the provided scale matrix is symmetric @@ -57,15 +64,16 @@ impl MultivariateStudent { } // Check that degrees of freedom is positive if freedom <= 0. { - return Err(StatsError::ArgMustBePositive("Degrees of freedom must be positive")) + return Err(StatsError::ArgMustBePositive( + "Degrees of freedom must be positive", + )); } let scale_det = scale.determinant(); - let pdf_const = gamma::gamma((freedom + (dim as f64)) / 2.) * - (gamma::gamma(freedom / 2.) * - (freedom.powi(dim as i32) * PI.powi(dim as i32) * scale_det.abs()) - .sqrt()) - .recip(); + 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), @@ -78,7 +86,7 @@ impl MultivariateStudent { scale, freedom, precision, - pdf_const, + ln_pdf_const, }) } } @@ -193,22 +201,24 @@ impl<'a> Continuous<&'a DVector, f64> for MultivariateStudent { // return mvn.pdf(x); // } let dv = x - &self.location; - let base_term = 1. + 1. / self.freedom - * *(&dv.transpose() * &self.precision * &dv) - .get((0, 0)) - .unwrap(); - self.pdf_const * base_term.powf(-(self.freedom + self.dim as f64) / 2.) + let base_term = 1. + + 1. / self.freedom + * *(&dv.transpose() * &self.precision * &dv) + .get((0, 0)) + .unwrap(); + self.ln_pdf_const.exp() * base_term.powf(-(self.freedom + self.dim 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 DVector) -> f64 { let dv = x - &self.location; - let base_term = 1. + 1. / self.freedom - * *(&dv.transpose() * &self.precision * &dv) - .get((0, 0)) - .unwrap(); - self.pdf_const.ln() - (self.freedom + self.dim as f64) / 2. * base_term.ln() + let base_term = 1. + + 1. / self.freedom + * *(&dv.transpose() * &self.precision * &dv) + .get((0, 0)) + .unwrap(); + self.ln_pdf_const - (self.freedom + self.dim as f64) / 2. * base_term.ln() } } @@ -266,6 +276,7 @@ mod tests { fn test_almost_multivariate_normal( location: Vec, scale: Vec, + freedom: f64, acc: f64, x: DVector, eval_mvs: F1, @@ -274,7 +285,7 @@ mod tests { F1: FnOnce(MultivariateStudent, DVector) -> f64, F2: FnOnce(MultivariateNormal, DVector) -> f64, { - let mvs = try_create(location.clone(), scale.clone(), f64::INFINITY); + let mvs = try_create(location.clone(), scale.clone(), freedom); let mvn0 = MultivariateNormal::new(location, scale); assert!(mvn0.is_ok()); let mvn = mvn0.unwrap(); @@ -386,9 +397,10 @@ mod tests { // TODO: These tests fail because inf degrees of freedom give NaN #[test] - fn test_pdf_freedom_inf() { + 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.], 1e-14, 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-15, dvec![1., 1.], pdf_mvs, pdf_mvn); } } From adf85bd364a7912f1e7ef0521973c56ac51b4fa0 Mon Sep 17 00:00:00 2001 From: Henry Jacobson Date: Sun, 25 Dec 2022 13:39:02 +0100 Subject: [PATCH 3/9] feat: creation of multivariate normal distribution from same variables as multivariate students (for when freedom = inf) --- src/distribution/multivariate_normal.rs | 22 ++++++++++++- src/distribution/multivariate_students_t.rs | 34 ++++++++++++++++++++- 2 files changed, 54 insertions(+), 2 deletions(-) diff --git a/src/distribution/multivariate_normal.rs b/src/distribution/multivariate_normal.rs index ab168020..677e0225 100644 --- a/src/distribution/multivariate_normal.rs +++ b/src/distribution/multivariate_normal.rs @@ -1,5 +1,5 @@ use crate::distribution::Continuous; -use crate::distribution::Normal; +use crate::distribution::{MultivariateStudent, Normal}; use crate::statistics::{Max, MeanN, Min, Mode, VarianceN}; use crate::{Result, StatsError}; use nalgebra::{ @@ -98,6 +98,26 @@ impl MultivariateNormal { .ln(), ) } + + /// Constructs a new multivariate normal distribution from a + /// multivariate students t distribution, which have equal variables + /// when `mvs.freedom == f64::INFINITY` + pub fn from_students(mvs: MultivariateStudent) -> Result { + let mu = mvs.location(); + let scale = mvs.scale(); + let cov_det = scale.determinant(); + let pdf_const = ((2. * PI).powi(mu.nrows() as i32) * cov_det.abs()) + .recip() + .sqrt(); + Ok(MultivariateNormal { + dim: mvs.dim(), + cov_chol_decomp: mvs.scale_chol_decomp(), + mu: mvs.location(), + cov: mvs.scale(), + precision: mvs.precision(), + pdf_const, + }) + } } impl ::rand::distributions::Distribution> for MultivariateNormal { diff --git a/src/distribution/multivariate_students_t.rs b/src/distribution/multivariate_students_t.rs index f1fa2feb..8f284126 100644 --- a/src/distribution/multivariate_students_t.rs +++ b/src/distribution/multivariate_students_t.rs @@ -73,7 +73,7 @@ impl MultivariateStudent { 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(); + - 0.5 * scale_det.ln(); match Cholesky::new(scale.clone()) { None => Err(StatsError::BadParams), @@ -91,6 +91,38 @@ impl MultivariateStudent { } } } + + /// Returns the dimension of the distribution + pub fn dim(&self) -> usize { + self.dim + } + /// Returns the cholesky decomposiiton matrix of the scale matrix + /// + /// Returns A where Σ = AAᵀ + pub fn scale_chol_decomp(&self) -> DMatrix { + self.scale_chol_decomp.clone() + } + /// Returns the location of the distribution + pub fn location(&self) -> DVector { + self.location.clone() + } + /// Returns the scale matrix of the distribution + pub fn scale(&self) -> DMatrix { + self.scale.clone() + } + /// 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) -> DMatrix { + self.precision.clone() + } + /// 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 { From 6cc368b59cc5727b84f0fe1c34454be65b259ffe Mon Sep 17 00:00:00 2001 From: Henry Jacobson Date: Sun, 25 Dec 2022 13:40:48 +0100 Subject: [PATCH 4/9] fix: use multivariate normal pdf when freedom = inf for multivariate student --- src/distribution/multivariate_students_t.rs | 42 ++++++++++++++------- 1 file changed, 28 insertions(+), 14 deletions(-) diff --git a/src/distribution/multivariate_students_t.rs b/src/distribution/multivariate_students_t.rs index 8f284126..19c4fea6 100644 --- a/src/distribution/multivariate_students_t.rs +++ b/src/distribution/multivariate_students_t.rs @@ -132,7 +132,7 @@ impl ::rand::distributions::Distribution> for MultivariateStudent { /// /// W * L * Z + μ /// - /// where `W` has √(ν/Sν) distribution, Sν has Chi-squared + /// 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 @@ -167,8 +167,8 @@ impl MeanN> for MultivariateStudent { /// /// # Remarks /// - /// This is the same mean used to construct the distribution if - /// the degrees of freedom is larger than 1. + /// 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. { let mut vec = vec![]; @@ -184,6 +184,14 @@ impl MeanN> for MultivariateStudent { impl VarianceN> for MultivariateStudent { /// Returns the covariance matrix of the multivariate student distribution + /// + /// # Formula + /// ```ignore + /// Σ ⋅ ν / (ν - 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.)) @@ -219,19 +227,14 @@ impl<'a> Continuous<&'a DVector, f64> for MultivariateStudent { /// ``` /// /// where `ν` is the degrees of freedom, `μ` is the mean, `Gamma` - /// is the Gamma function, `inv(Σ)` + /// 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 - /// - /// TODO: Make this converge for large degrees of freedom - /// Current commented code beneath fails since `MultivariateNormal::new` accepts Vec and - /// not DVector or DMatrix. Should implement that instead of changing back to Vec, or - /// even have a constructor `MultivariateNormal::from_student`. fn pdf(&self, x: &'a DVector) -> f64 { - // if self.freedom == f64::INFINITY { - // let mvn = MultivariateNormal::new(self.location, self.scale).unwrap(); - // return mvn.pdf(x); - // } + if self.freedom == f64::INFINITY { + let mvn = MultivariateNormal::from_students(self.clone()).unwrap(); + return mvn.pdf(x); + } let dv = x - &self.location; let base_term = 1. + 1. / self.freedom @@ -244,6 +247,10 @@ impl<'a> Continuous<&'a DVector, f64> for MultivariateStudent { /// Calculates the log probability density function for the multivariate /// student distribution at `x`. Equivalent to pdf(x).ln(). fn ln_pdf(&self, x: &'a DVector) -> f64 { + if self.freedom == f64::INFINITY { + let mvn = MultivariateNormal::from_students(self.clone()).unwrap(); + return mvn.ln_pdf(x); + } let dv = x - &self.location; let base_term = 1. + 1. / self.freedom @@ -427,12 +434,19 @@ mod tests { test_case(vec![-1., 0.], vec![f64::INFINITY, 0., 0., f64::INFINITY], 10., f64::NEG_INFINITY, pdf(dvec![10., 10.])); } - // TODO: These tests fail because inf degrees of freedom give NaN #[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-15, dvec![1., 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.], 1e10, 1e-5, dvec![1., 1.], pdf_mvs, pdf_mvn); + test_almost_multivariate_normal(vec![0., 0.,], vec![1., 0., 0., 1.], f64::INFINITY, 1e-50, dvec![1., 1.], pdf_mvs, pdf_mvn); + } } From aeb720b447d5fd85a217e4022d805448af4d2972 Mon Sep 17 00:00:00 2001 From: Henry Jacobson Date: Mon, 26 Dec 2022 23:03:38 +0100 Subject: [PATCH 5/9] test: panic test for invalid pdf argument --- src/distribution/multivariate_normal.rs | 7 +++++++ src/distribution/multivariate_students_t.rs | 7 +++++++ 2 files changed, 14 insertions(+) diff --git a/src/distribution/multivariate_normal.rs b/src/distribution/multivariate_normal.rs index 677e0225..73dbbfe0 100644 --- a/src/distribution/multivariate_normal.rs +++ b/src/distribution/multivariate_normal.rs @@ -373,4 +373,11 @@ mod tests { 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.])); } + + #[test] + #[should_panic] + fn test_pdf_mismatched_arg_size() { + let mvn = MultivariateNormal::new(vec![0., 0.], vec![1., 0., 0., 1.,]).unwrap(); + mvn.pdf(&dvec![1.]); // x.size != mu.size + } } diff --git a/src/distribution/multivariate_students_t.rs b/src/distribution/multivariate_students_t.rs index 19c4fea6..b2f3f683 100644 --- a/src/distribution/multivariate_students_t.rs +++ b/src/distribution/multivariate_students_t.rs @@ -449,4 +449,11 @@ mod tests { test_almost_multivariate_normal(vec![0., 0.,], vec![1., 0., 0., 1.], 1e10, 1e-5, dvec![1., 1.], pdf_mvs, pdf_mvn); test_almost_multivariate_normal(vec![0., 0.,], vec![1., 0., 0., 1.], f64::INFINITY, 1e-50, dvec![1., 1.], pdf_mvs, pdf_mvn); } + + #[test] + #[should_panic] + fn test_pdf_mismatched_arg_size() { + let mvs = MultivariateStudent::new(vec![0., 0.], vec![1., 0., 0., 1.,], 3.).unwrap(); + mvs.pdf(&dvec![1.]); // x.size != mu.size + } } From 28b5c687c58118a417d795f4ca2ec4bd3080c646 Mon Sep 17 00:00:00 2001 From: Henry Jacobson Date: Sat, 2 Mar 2024 12:37:51 +0100 Subject: [PATCH 6/9] fix: tests in documentation --- src/distribution/multivariate_students_t.rs | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/distribution/multivariate_students_t.rs b/src/distribution/multivariate_students_t.rs index b2f3f683..3acdd6f4 100644 --- a/src/distribution/multivariate_students_t.rs +++ b/src/distribution/multivariate_students_t.rs @@ -23,8 +23,8 @@ use std::f64::consts::{E, PI}; /// /// 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, 4. * vec![1., 0., 0., 1.])); -/// assert_eq!(mvs.pdf(&DVector::from_vec(vec![1., 1.])), 0.047157020175376416); +/// 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 { @@ -263,7 +263,7 @@ impl<'a> Continuous<&'a DVector, f64> for MultivariateStudent { #[rustfmt::skip] #[cfg(all(test, feature = "nightly"))] -mod tests { +mod tests { use crate::distribution::MultivariateNormal; use core::fmt::Debug; From af92615778f0c23c0b32d5061347312104c40222 Mon Sep 17 00:00:00 2001 From: Henry Jacobson Date: Sat, 2 Mar 2024 12:43:57 +0100 Subject: [PATCH 7/9] fix: clearer function name in test --- src/distribution/multivariate_students_t.rs | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/distribution/multivariate_students_t.rs b/src/distribution/multivariate_students_t.rs index 3acdd6f4..d16f8537 100644 --- a/src/distribution/multivariate_students_t.rs +++ b/src/distribution/multivariate_students_t.rs @@ -426,12 +426,12 @@ mod tests { #[test] fn test_ln_pdf() { - let 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, pdf(dvec![1., 1.])); - test_almost(vec![0., 0.], vec![1., 0., 0., 1.], 2., -4.3434030034000815, 1e-14, pdf(dvec![1., 2.])); - test_almost(vec![2., 1.], vec![5., 0., 0., 1.], 2.5, -10.542229575274265, 1e-14, pdf(dvec![1., 10.])); - test_almost(vec![-1., 0.], vec![2., 1., 1., 6.], 1.5, -9.650699521198622, 1e-14, pdf(dvec![10., 10.])); - test_case(vec![-1., 0.], vec![f64::INFINITY, 0., 0., f64::INFINITY], 10., f64::NEG_INFINITY, pdf(dvec![10., 10.])); + 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] From 67ca3ac4c7d38a4f411388a1d980868bd566ad26 Mon Sep 17 00:00:00 2001 From: Henry Jacobson Date: Sat, 2 Mar 2024 13:02:45 +0100 Subject: [PATCH 8/9] fix: add documentation --- src/distribution/multivariate_students_t.rs | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/src/distribution/multivariate_students_t.rs b/src/distribution/multivariate_students_t.rs index d16f8537..16350d1c 100644 --- a/src/distribution/multivariate_students_t.rs +++ b/src/distribution/multivariate_students_t.rs @@ -43,7 +43,7 @@ impl MultivariateStudent { /// /// # Errors /// - /// Returns `StatsError::BadParams` if the scale matrix is not Symmetric-positive + /// 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(); @@ -76,7 +76,7 @@ impl MultivariateStudent { - 0.5 * scale_det.ln(); match Cholesky::new(scale.clone()) { - None => Err(StatsError::BadParams), + None => Err(StatsError::BadParams), // Scale matrix is not positive definite Some(cholesky_decomp) => { let precision = cholesky_decomp.inverse(); Ok(MultivariateStudent { @@ -130,7 +130,9 @@ impl ::rand::distributions::Distribution> for MultivariateStudent { /// /// # Formula /// + ///```ignore /// W * L * Z + μ + ///``` /// /// where `W` has √(ν/Sν) distribution, Sν has Chi-squared /// distribution with ν degrees of freedom, @@ -261,6 +263,7 @@ impl<'a> Continuous<&'a DVector, f64> for MultivariateStudent { } } +// TODO: Add more tests for other matrices than really straightforward symmetric positive #[rustfmt::skip] #[cfg(all(test, feature = "nightly"))] mod tests { @@ -396,6 +399,7 @@ mod tests { 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); } #[test] From 4f22558d175ffa6ec4523e5bfa480fe19e116a75 Mon Sep 17 00:00:00 2001 From: Henry Jacobson Date: Sun, 3 Mar 2024 14:29:06 +0100 Subject: [PATCH 9/9] tests: 3d matrices test cases for pdf. Also improves documentation for multivariate t minorly. --- src/distribution/multivariate_students_t.rs | 98 ++++++++++----------- 1 file changed, 49 insertions(+), 49 deletions(-) diff --git a/src/distribution/multivariate_students_t.rs b/src/distribution/multivariate_students_t.rs index 16350d1c..d20ff0e3 100644 --- a/src/distribution/multivariate_students_t.rs +++ b/src/distribution/multivariate_students_t.rs @@ -10,9 +10,9 @@ use rand::Rng; use std::f64::consts::{E, PI}; /// Implements the [Multivariate Student's t-distribution](https://en.wikipedia.org/wiki/Multivariate_t-distribution) -/// distribution using the "nalgebra" crate for matrix operations +/// distribution using the "nalgebra" crate for matrix operations. /// -/// Assumes all the marginal distributions have the same degree of freedom, ν +/// Assumes all the marginal distributions have the same degree of freedom, ν. /// /// # Examples /// @@ -39,7 +39,7 @@ pub struct MultivariateStudent { impl MultivariateStudent { /// Constructs a new multivariate students t distribution with a location of `location`, - /// scale matrix `scale` and `freedom` degrees of freedom + /// scale matrix `scale` and `freedom` degrees of freedom. /// /// # Errors /// @@ -92,34 +92,34 @@ impl MultivariateStudent { } } - /// Returns the dimension of the distribution + /// Returns the dimension of the distribution. pub fn dim(&self) -> usize { self.dim } - /// Returns the cholesky decomposiiton matrix of the scale matrix + /// Returns the cholesky decomposiiton matrix of the scale matrix. /// - /// Returns A where Σ = AAᵀ + /// Returns A where Σ = AAᵀ. pub fn scale_chol_decomp(&self) -> DMatrix { self.scale_chol_decomp.clone() } - /// Returns the location of the distribution + /// Returns the location of the distribution. pub fn location(&self) -> DVector { self.location.clone() } - /// Returns the scale matrix of the distribution + /// Returns the scale matrix of the distribution. pub fn scale(&self) -> DMatrix { self.scale.clone() } - /// Returns the degrees of freedom of the distribution + /// Returns the degrees of freedom of the distribution. pub fn freedom(&self) -> f64 { self.freedom } - /// Returns the inverse of the cholesky decomposition matrix + /// Returns the inverse of the cholesky decomposition matrix. pub fn precision(&self) -> DMatrix { self.precision.clone() } /// Returns the logarithmed constant part of the probability - /// distribution function + /// distribution function. pub fn ln_pdf_const(&self) -> f64 { self.ln_pdf_const } @@ -130,8 +130,8 @@ impl ::rand::distributions::Distribution> for MultivariateStudent { /// /// # Formula /// - ///```ignore - /// W * L * Z + μ + ///```math + /// W ⋅ L ⋅ Z + μ ///``` /// /// where `W` has √(ν/Sν) distribution, Sν has Chi-squared @@ -165,19 +165,15 @@ impl Max> for MultivariateStudent { } impl MeanN> for MultivariateStudent { - /// Returns the mean of the student distribution + /// 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 + /// the degrees of freedom is larger than 1. fn mean(&self) -> Option> { if self.freedom > 1. { - let mut vec = vec![]; - for elt in self.location.clone().into_iter() { - vec.push(*elt); - } - Some(DVector::from_vec(vec)) + Some(self.location.clone()) } else { None } @@ -185,15 +181,16 @@ impl MeanN> for MultivariateStudent { } impl VarianceN> for MultivariateStudent { - /// Returns the covariance matrix of the multivariate student distribution + /// Returns the covariance matrix of the multivariate student distribution. /// /// # Formula - /// ```ignore + /// + /// ```math /// Σ ⋅ ν / (ν - 2) /// ``` /// /// where `Σ` is the scale matrix and `ν` is the degrees of freedom. - /// Only defined if freedom is larger than 2 + /// 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.)) @@ -204,34 +201,34 @@ impl VarianceN> for MultivariateStudent { } impl Mode> for MultivariateStudent { - /// Returns the mode of the multivariate student distribution + /// Returns the mode of the multivariate student distribution. /// /// # Formula /// - /// ```ignore + /// ```math /// μ /// ``` /// - /// where `μ` is the location + /// where `μ` is the location. fn mode(&self) -> DVector { self.location.clone() } } impl<'a> Continuous<&'a DVector, f64> for MultivariateStudent { - /// Calculates the probability density function for the multivariate - /// student distribution at `x` + /// Calculates the probability density function for the multivariate. + /// student distribution at `x`. /// /// # Formula /// - /// ```ignore - /// Gamma[(ν+p)/2] / [Gamma(ν/2) ((ν * π)^p det(Σ))^(1 / 2)] * [1 + 1/ν transpose(x - μ) inv(Σ) (x - μ)]^(-(ν+p)/2) + /// ```math + /// [Γ(ν+p)/2] / [Γ(ν/2) ((ν * π)^p det(Σ))^(1 / 2)] * [1 + 1/ν (x - μ)ᵀ inv(Σ) (x - μ)]^(-(ν+p)/2) /// ``` /// - /// where `ν` is the degrees of freedom, `μ` is the mean, `Gamma` + /// 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 + /// of the scale matrix, and `k` is the dimension of the distribution. fn pdf(&self, x: &'a DVector) -> f64 { if self.freedom == f64::INFINITY { let mvn = MultivariateNormal::from_students(self.clone()).unwrap(); @@ -263,7 +260,6 @@ impl<'a> Continuous<&'a DVector, f64> for MultivariateStudent { } } -// TODO: Add more tests for other matrices than really straightforward symmetric positive #[rustfmt::skip] #[cfg(all(test, feature = "nightly"))] mod tests { @@ -361,18 +357,19 @@ mod tests { #[test] fn test_bad_create() { - // scale not symmetric + // scale not symmetric. bad_create_case(vec![0., 0.], vec![1., 1., 0., 1.], 1.); - // scale not positive-definite + // scale not positive-definite. bad_create_case(vec![0., 0.], vec![1., 2., 2., 1.], 1.); - // NaN in location + // NaN in location. bad_create_case(vec![0., f64::NAN], vec![1., 0., 0., 1.], 1.); - // NaN in scale Matrix + // NaN in scale Matrix. bad_create_case(vec![0., 0.], vec![1., 0., 0., f64::NAN], 1.); - // NaN in freedom + // NaN in freedom. bad_create_case(vec![0., 0.], vec![1., 0., 0., 1.], f64::NAN); - // Non-positive freedom + // 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] @@ -382,6 +379,7 @@ mod tests { 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(); @@ -402,6 +400,7 @@ mod tests { 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(); @@ -422,9 +421,14 @@ mod tests { 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.952297846610382e-16, 1e-30, pdf(dvec![0.3010, 0.1491, 0.5008])); test_case(vec![-1., 0.], vec![f64::INFINITY, 0., 0., f64::INFINITY], 10., 0., pdf(dvec![10., 10.])); } @@ -444,20 +448,16 @@ mod tests { 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-15, 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.], 1e10, 1e-5, dvec![1., 1.], pdf_mvs, pdf_mvn); - test_almost_multivariate_normal(vec![0., 0.,], vec![1., 0., 0., 1.], f64::INFINITY, 1e-50, dvec![1., 1.], pdf_mvs, pdf_mvn); - } - - #[test] - #[should_panic] - fn test_pdf_mismatched_arg_size() { - let mvs = MultivariateStudent::new(vec![0., 0.], vec![1., 0., 0., 1.,], 3.).unwrap(); - mvs.pdf(&dvec![1.]); // x.size != mu.size + 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); } }