From 521c74b3e4e776b335dcb95dba0d3d2c59176e21 Mon Sep 17 00:00:00 2001 From: Henry Jacobson Date: Sat, 24 Dec 2022 00:31:15 +0100 Subject: [PATCH 01/10] feat: implement multivariate students t distribution --- src/distribution/mod.rs | 2 + src/distribution/multivariate_students_t.rs | 396 ++++++++++++++++++++ 2 files changed, 398 insertions(+) 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_students_t.rs b/src/distribution/multivariate_students_t.rs new file mode 100644 index 00000000..c8c118d0 --- /dev/null +++ b/src/distribution/multivariate_students_t.rs @@ -0,0 +1,396 @@ +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; +use nalgebra::{DMatrix, DVector}; +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, 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(test)] +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 607bc3044a68b679090b6917eaf2627494c42d1b Mon Sep 17 00:00:00 2001 From: Henry Jacobson Date: Sun, 25 Dec 2022 13:08:31 +0100 Subject: [PATCH 02/10] fix: Use ln_pdf_const instead of pdf_const --- src/distribution/multivariate_students_t.rs | 31 +++++++++++++-------- 1 file changed, 20 insertions(+), 11 deletions(-) diff --git a/src/distribution/multivariate_students_t.rs b/src/distribution/multivariate_students_t.rs index c8c118d0..fdd4a1ab 100644 --- a/src/distribution/multivariate_students_t.rs +++ b/src/distribution/multivariate_students_t.rs @@ -33,10 +33,17 @@ 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); @@ -62,10 +69,10 @@ impl MultivariateStudent { } 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 +85,7 @@ impl MultivariateStudent { scale, freedom, precision, - pdf_const, + ln_pdf_const, }) } } @@ -198,7 +205,7 @@ impl<'a> Continuous<&'a DVector, f64> for MultivariateStudent { * *(&dv.transpose() * &self.precision * &dv) .get((0, 0)) .unwrap(); - self.pdf_const * base_term.powf(-(self.freedom + self.dim as f64) / 2.) + self.ln_pdf_const.exp() * base_term.powf(-(self.freedom + self.dim as f64) / 2.) } /// Calculates the log probability density function for the multivariate @@ -210,7 +217,7 @@ impl<'a> Continuous<&'a DVector, f64> for MultivariateStudent { * *(&dv.transpose() * &self.precision * &dv) .get((0, 0)) .unwrap(); - self.pdf_const.ln() - (self.freedom + self.dim as f64) / 2. * base_term.ln() + self.ln_pdf_const - (self.freedom + self.dim as f64) / 2. * base_term.ln() } } @@ -268,6 +275,7 @@ mod tests { fn test_almost_multivariate_normal( location: Vec, scale: Vec, + freedom: f64, acc: f64, x: DVector, eval_mvs: F1, @@ -276,7 +284,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(); @@ -388,9 +396,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 517061244679e03a13c7afb81e2dcc761ed203b5 Mon Sep 17 00:00:00 2001 From: Henry Jacobson Date: Sun, 25 Dec 2022 13:39:02 +0100 Subject: [PATCH 03/10] 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 | 32 +++++++++++++++++++++ 2 files changed, 53 insertions(+), 1 deletion(-) diff --git a/src/distribution/multivariate_normal.rs b/src/distribution/multivariate_normal.rs index 9949d676..f5cefc9e 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::{Cholesky, Const, DMatrix, DVector, Dim, DimMin, Dyn, OMatrix, OVector}; @@ -117,6 +117,26 @@ where .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 std::fmt::Display for MultivariateNormal diff --git a/src/distribution/multivariate_students_t.rs b/src/distribution/multivariate_students_t.rs index fdd4a1ab..1bd785f8 100644 --- a/src/distribution/multivariate_students_t.rs +++ b/src/distribution/multivariate_students_t.rs @@ -90,6 +90,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 09db04f81f987821bd8fd01d839e390216be29ce Mon Sep 17 00:00:00 2001 From: Henry Jacobson Date: Sun, 25 Dec 2022 13:40:48 +0100 Subject: [PATCH 04/10] fix: use multivariate normal pdf when freedom = inf for multivariate student --- src/distribution/multivariate_normal.rs | 39 +++++---- src/distribution/multivariate_students_t.rs | 87 +++++++++++++-------- 2 files changed, 72 insertions(+), 54 deletions(-) diff --git a/src/distribution/multivariate_normal.rs b/src/distribution/multivariate_normal.rs index f5cefc9e..a524e440 100644 --- a/src/distribution/multivariate_normal.rs +++ b/src/distribution/multivariate_normal.rs @@ -49,6 +49,25 @@ impl MultivariateNormal { let cov = DMatrix::from_vec(mean.len(), mean.len(), cov); MultivariateNormal::new_from_nalgebra(mean, cov) } + + /// 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 { + cov_chol_decomp: mvs.scale_chol_decomp(), + mu: mvs.location(), + cov: mvs.scale(), + precision: mvs.precision(), + pdf_const, + }) + } } impl MultivariateNormal @@ -117,26 +136,6 @@ where .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 std::fmt::Display for MultivariateNormal diff --git a/src/distribution/multivariate_students_t.rs b/src/distribution/multivariate_students_t.rs index 1bd785f8..8155012f 100644 --- a/src/distribution/multivariate_students_t.rs +++ b/src/distribution/multivariate_students_t.rs @@ -1,5 +1,5 @@ use crate::distribution::Continuous; -use crate::distribution::{ChiSquared, Normal}; +use crate::distribution::{ChiSquared, MultivariateNormal, Normal}; use crate::function::gamma; use crate::statistics::{Max, MeanN, Min, Mode, VarianceN}; use crate::{Result, StatsError}; @@ -183,6 +183,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.)) @@ -227,10 +235,10 @@ impl<'a> Continuous<&'a DVector, f64> for MultivariateStudent { /// 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 @@ -243,6 +251,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 @@ -256,7 +268,6 @@ impl<'a> Continuous<&'a DVector, f64> for MultivariateStudent { #[rustfmt::skip] #[cfg(test)] mod tests { - use crate::distribution::MultivariateNormal; use core::fmt::Debug; fn try_create(location: Vec, scale: Vec, freedom: f64) -> MultivariateStudent @@ -304,26 +315,26 @@ mod tests { 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_almost_eq!(mvs_x, mvn_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_almost_eq!(mvs_x, mvn_x, acc); + // } use super::*; @@ -426,12 +437,20 @@ 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.], 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_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 0e086957caaa73de044de209f3c8c51a0cd0b68d Mon Sep 17 00:00:00 2001 From: Henry Jacobson Date: Mon, 26 Dec 2022 23:03:38 +0100 Subject: [PATCH 05/10] test: panic test for invalid pdf argument --- src/distribution/multivariate_normal.rs | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/src/distribution/multivariate_normal.rs b/src/distribution/multivariate_normal.rs index a524e440..6a5103f1 100644 --- a/src/distribution/multivariate_normal.rs +++ b/src/distribution/multivariate_normal.rs @@ -628,4 +628,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(&dvec![1.]); // x.size != mu.size + } } From 81903aea34b46d32cbcfe7033eb41588e72fb7cf Mon Sep 17 00:00:00 2001 From: Henry Jacobson Date: Sat, 2 Mar 2024 12:37:51 +0100 Subject: [PATCH 06/10] fix: tests in documentation --- src/distribution/multivariate_normal.rs | 2 +- src/distribution/multivariate_students_t.rs | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/distribution/multivariate_normal.rs b/src/distribution/multivariate_normal.rs index 6a5103f1..d3110157 100644 --- a/src/distribution/multivariate_normal.rs +++ b/src/distribution/multivariate_normal.rs @@ -633,6 +633,6 @@ mod tests { #[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 + mvn.pdf(vec![1.]); // x.size != mu.size } } diff --git a/src/distribution/multivariate_students_t.rs b/src/distribution/multivariate_students_t.rs index 8155012f..ab797bb6 100644 --- a/src/distribution/multivariate_students_t.rs +++ b/src/distribution/multivariate_students_t.rs @@ -22,8 +22,8 @@ use std::f64::consts::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 { From 9746d6983c7c7dabc587dd05942ff700d1fe0e2a Mon Sep 17 00:00:00 2001 From: Henry Jacobson Date: Sat, 2 Mar 2024 12:43:57 +0100 Subject: [PATCH 07/10] 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 ab797bb6..aab6682f 100644 --- a/src/distribution/multivariate_students_t.rs +++ b/src/distribution/multivariate_students_t.rs @@ -429,12 +429,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 612355bf8e56a4d8f7f041243a7eb7f3ad1e105e Mon Sep 17 00:00:00 2001 From: Henry Jacobson Date: Sat, 2 Mar 2024 13:02:45 +0100 Subject: [PATCH 08/10] 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 aab6682f..ec47ac24 100644 --- a/src/distribution/multivariate_students_t.rs +++ b/src/distribution/multivariate_students_t.rs @@ -42,7 +42,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(); @@ -75,7 +75,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 { @@ -129,7 +129,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, @@ -265,6 +267,7 @@ impl<'a> Continuous<&'a DVector, f64> for MultivariateStudent { } } +// TODO: Add more tests for other matrices than really straightforward symmetric positive #[rustfmt::skip] #[cfg(test)] mod tests { @@ -399,6 +402,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 b8261a0c77cbb6344b179aa045fff22ad35861e6 Mon Sep 17 00:00:00 2001 From: Henry Jacobson Date: Sun, 3 Mar 2024 14:29:06 +0100 Subject: [PATCH 09/10] tests: 3d matrices test cases for pdf. Also improves documentation for multivariate t minorly. --- src/distribution/multivariate_students_t.rs | 95 +++++++++++---------- 1 file changed, 48 insertions(+), 47 deletions(-) diff --git a/src/distribution/multivariate_students_t.rs b/src/distribution/multivariate_students_t.rs index ec47ac24..f7e8929a 100644 --- a/src/distribution/multivariate_students_t.rs +++ b/src/distribution/multivariate_students_t.rs @@ -9,9 +9,9 @@ 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 +/// 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 /// @@ -38,7 +38,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 /// @@ -91,34 +91,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 } @@ -129,8 +129,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 @@ -164,7 +164,7 @@ impl Max> for MultivariateStudent { } impl MeanN> for MultivariateStudent { - /// Returns the mean of the student distribution + /// Returns the mean of the student distribution. /// /// # Remarks /// @@ -172,11 +172,7 @@ impl MeanN> for MultivariateStudent { /// 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 } @@ -184,15 +180,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.)) @@ -203,39 +200,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 - /// - /// 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`. + /// 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(); @@ -267,7 +259,6 @@ impl<'a> Continuous<&'a DVector, f64> for MultivariateStudent { } } -// TODO: Add more tests for other matrices than really straightforward symmetric positive #[rustfmt::skip] #[cfg(test)] mod tests { @@ -364,18 +355,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] @@ -385,6 +377,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(); @@ -405,6 +398,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(); @@ -425,9 +419,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.])); } @@ -447,14 +446,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_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); // } } From f75730341762a265b1c81e83646dc538d6df9e6c Mon Sep 17 00:00:00 2001 From: Orion Yeung <11580988+orionyeung001@users.noreply.github.com> Date: Fri, 16 Aug 2024 15:01:16 -0500 Subject: [PATCH 10/10] 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 --- src/distribution/multivariate_students_t.rs | 12 +++++++++--- 1 file changed, 9 insertions(+), 3 deletions(-) diff --git a/src/distribution/multivariate_students_t.rs b/src/distribution/multivariate_students_t.rs index f7e8929a..d1348f92 100644 --- a/src/distribution/multivariate_students_t.rs +++ b/src/distribution/multivariate_students_t.rs @@ -95,28 +95,34 @@ impl MultivariateStudent { 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 { @@ -129,9 +135,9 @@ impl ::rand::distributions::Distribution> for MultivariateStudent { /// /// # Formula /// - ///```math + /// ```math /// W ⋅ L ⋅ Z + μ - ///``` + /// ``` /// /// where `W` has √(ν/Sν) distribution, Sν has Chi-squared /// distribution with ν degrees of freedom, @@ -426,7 +432,7 @@ mod tests { // 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_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.])); }