diff --git a/src/distribution/mod.rs b/src/distribution/mod.rs index dea1f9af..b146f9d4 100644 --- a/src/distribution/mod.rs +++ b/src/distribution/mod.rs @@ -27,6 +27,7 @@ pub use self::laplace::Laplace; pub use self::log_normal::LogNormal; pub use self::multinomial::Multinomial; pub use self::multivariate_normal::MultivariateNormal; +pub use self::multivariate_students_t::MultivariateStudent; pub use self::negative_binomial::NegativeBinomial; pub use self::normal::Normal; pub use self::pareto::Pareto; @@ -60,6 +61,7 @@ mod laplace; mod log_normal; mod multinomial; mod multivariate_normal; +mod multivariate_students_t; mod negative_binomial; mod normal; mod pareto; diff --git a/src/distribution/multivariate_normal.rs b/src/distribution/multivariate_normal.rs index 9949d676..d3110157 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}; @@ -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 @@ -609,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(vec![1.]); // x.size != mu.size + } } diff --git a/src/distribution/multivariate_students_t.rs b/src/distribution/multivariate_students_t.rs new file mode 100644 index 00000000..d1348f92 --- /dev/null +++ b/src/distribution/multivariate_students_t.rs @@ -0,0 +1,467 @@ +use crate::distribution::Continuous; +use crate::distribution::{ChiSquared, MultivariateNormal, 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, vec![2., 0., 0., 2.])); +/// assert_eq!(mvs.pdf(&DVector::from_vec(vec![1., 1.])), 0.04715702017537655); +/// ``` +#[derive(Debug, Clone, PartialEq)] +pub struct MultivariateStudent { + dim: usize, + scale_chol_decomp: DMatrix, + location: DVector, + scale: DMatrix, + freedom: f64, + precision: DMatrix, + 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 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 ln_pdf_const = gamma::ln_gamma(0.5 * (freedom + dim as f64)) + - gamma::ln_gamma(0.5 * freedom) + - 0.5 * (dim as f64) * (freedom * PI).ln() + - 0.5 * scale_det.ln(); + + match Cholesky::new(scale.clone()) { + None => Err(StatsError::BadParams), // Scale matrix is not positive definite + Some(cholesky_decomp) => { + let precision = cholesky_decomp.inverse(); + Ok(MultivariateStudent { + dim, + scale_chol_decomp: cholesky_decomp.unpack(), + location, + scale, + freedom, + precision, + ln_pdf_const, + }) + } + } + } + + /// 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 { + /// Samples from the multivariate student distribution + /// + /// # Formula + /// + /// ```math + /// W ⋅ L ⋅ Z + μ + /// ``` + /// + /// where `W` has √(ν/Sν) distribution, Sν has Chi-squared + /// distribution with ν degrees of freedom, + /// `L` is the Cholesky decomposition of the scale matrix, + /// `Z` is a vector of normally distributed random variables, and + /// `μ` is the location vector + fn sample(&self, rng: &mut R) -> 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. { + Some(self.location.clone()) + } else { + None + } + } +} + +impl VarianceN> for MultivariateStudent { + /// Returns the covariance matrix of the multivariate student distribution. + /// + /// # Formula + /// + /// ```math + /// Σ ⋅ ν / (ν - 2) + /// ``` + /// + /// where `Σ` is the scale matrix and `ν` is the degrees of freedom. + /// Only defined if freedom is larger than 2. + fn variance(&self) -> Option> { + if self.freedom > 2. { + Some(self.scale.clone() * self.freedom / (self.freedom - 2.)) + } else { + None + } + } +} + +impl Mode> for MultivariateStudent { + /// Returns the mode of the multivariate student distribution. + /// + /// # Formula + /// + /// ```math + /// μ + /// ``` + /// + /// 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 + /// + /// ```math + /// [Γ(ν+p)/2] / [Γ(ν/2) ((ν * π)^p det(Σ))^(1 / 2)] * [1 + 1/ν (x - μ)ᵀ inv(Σ) (x - μ)]^(-(ν+p)/2) + /// ``` + /// + /// where `ν` is the degrees of freedom, `μ` is the mean, `Γ` + /// is the Gamma function, `inv(Σ)` + /// is the precision matrix, `det(Σ)` is the determinant + /// of the scale matrix, and `k` is the dimension of the distribution. + fn pdf(&self, x: &'a DVector) -> f64 { + 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 + * *(&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 { + 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 + * *(&dv.transpose() * &self.precision * &dv) + .get((0, 0)) + .unwrap(); + self.ln_pdf_const - (self.freedom + self.dim as f64) / 2. * base_term.ln() + } +} + +#[rustfmt::skip] +#[cfg(test)] +mod tests { + 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, + // 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::*; + + macro_rules! dvec { + ($($x:expr),*) => (DVector::from_vec(vec![$($x),*])); + } + + macro_rules! mat2 { + ($x11:expr, $x12:expr, $x21:expr, $x22:expr) => (DMatrix::from_vec(2,2,vec![$x11, $x12, $x21, $x22])); + } + + // macro_rules! mat3 { + // ($x11:expr, $x12:expr, $x13:expr, $x21:expr, $x22:expr, $x23:expr, $x31:expr, $x32:expr, $x33:expr) => (DMatrix::from_vec(3,3,vec![$x11, $x12, $x13, $x21, $x22, $x23, $x31, $x32, $x33])); + // } + + #[test] + fn test_create() { + create_case(vec![0., 0.], vec![1., 0., 0., 1.], 1.); + create_case(vec![10., 5.], vec![2., 1., 1., 2.], 3.); + create_case(vec![4., 5., 6.], vec![2., 1., 0., 1., 2., 1., 0., 1., 2.], 14.); + create_case(vec![0., f64::INFINITY], vec![1., 0., 0., 1.], f64::INFINITY); + create_case(vec![0., 0.], vec![f64::INFINITY, 0., 0., f64::INFINITY], 0.1); + } + + #[test] + fn test_bad_create() { + // scale not symmetric. + bad_create_case(vec![0., 0.], vec![1., 1., 0., 1.], 1.); + // scale not positive-definite. + bad_create_case(vec![0., 0.], vec![1., 2., 2., 1.], 1.); + // NaN in location. + bad_create_case(vec![0., f64::NAN], vec![1., 0., 0., 1.], 1.); + // NaN in scale Matrix. + bad_create_case(vec![0., 0.], vec![1., 0., 0., f64::NAN], 1.); + // NaN in freedom. + bad_create_case(vec![0., 0.], vec![1., 0., 0., 1.], f64::NAN); + // Non-positive freedom. + bad_create_case(vec![0., 0.], vec![1., 0., 0., 1.], 0.); + bad_create_case(vec![0., 0.], vec![1., 0., 0., 1.], -1.); + } + + #[test] + fn test_variance() { + let variance = |x: MultivariateStudent| x.variance().unwrap(); + test_case(vec![0., 0.], vec![1., 0., 0., 1.], 3., 3. * mat2![1., 0., 0., 1.], variance); + test_case(vec![0., 0.], vec![f64::INFINITY, 0., 0., f64::INFINITY], 3., mat2![f64::INFINITY, 0., 0., f64::INFINITY], variance); + } + + // Variance is only defined for freedom > 2. + #[test] + fn test_bad_variance() { + let variance = |x: MultivariateStudent| x.variance(); + test_case(vec![0., 0.], vec![1., 0., 0., 1.], 2., None, variance); + } + + #[test] + fn test_mode() { + let mode = |x: MultivariateStudent| x.mode(); + test_case(vec![0., 0.], vec![1., 0., 0., 1.], 1., dvec![0., 0.], mode); + test_case(vec![f64::INFINITY, f64::INFINITY], vec![1., 0., 0., 1.], 1., dvec![f64::INFINITY, f64::INFINITY], mode); + } + + #[test] + fn test_mean() { + let mean = |x: MultivariateStudent| x.mean().unwrap(); + test_case(vec![0., 0.], vec![1., 0., 0., 1.], 2., dvec![0., 0.], mean); + test_case(vec![-1., 1., 3.], vec![1., 0., 0.5, 0., 2.0, 0., 0.5, 0., 3.0], 2., dvec![-1., 1., 3.], mean); + } + + // Mean is only defined if freedom > 1. + #[test] + fn test_bad_mean() { + let mean = |x: MultivariateStudent| x.mean(); + test_case(vec![0., 0.], vec![1., 0., 0., 1.], 1., None, mean); + } + + #[test] + fn test_min_max() { + let min = |x: MultivariateStudent| x.min(); + let max = |x: MultivariateStudent| x.max(); + test_case(vec![0., 0.], vec![1., 0., 0., 1.], 1., dvec![f64::NEG_INFINITY, f64::NEG_INFINITY], min); + test_case(vec![0., 0.], vec![1., 0., 0., 1.], 1., dvec![f64::INFINITY, f64::INFINITY], max); + test_case(vec![10., 1.], vec![1., 0., 0., 1.], 1., dvec![f64::NEG_INFINITY, f64::NEG_INFINITY], min); + test_case(vec![-3., 5.], vec![1., 0., 0., 1.], 1., dvec![f64::INFINITY, f64::INFINITY], max); + } + + #[test] + fn test_pdf() { + let pdf = |arg: DVector| move |x: MultivariateStudent| x.pdf(&arg); + test_almost(vec![0., 0.], vec![1., 0., 0., 1.], 4., 0.047157020175376416, 1e-15, pdf(dvec![1., 1.])); + test_almost(vec![0., 0.], vec![1., 0., 0., 1.], 4., 0.013972450422333741737457302178882, 1e-15, pdf(dvec![1., 2.])); + test_almost(vec![0., 0.], vec![1., 0., 0., 1.], 2., 0.012992240252399619, 1e-17, pdf(dvec![1., 2.])); + test_almost(vec![2., 1.], vec![5., 0., 0., 1.], 2.5, 2.639780816598878e-5, 1e-19, pdf(dvec![1., 10.])); + test_almost(vec![-1., 0.], vec![2., 1., 1., 6.], 1.5, 6.438051574348526e-5, 1e-19, pdf(dvec![10., 10.])); + // These three are crossed checked against both python's scipy.multivariate_t.pdf and octave's mvtpdf. + test_almost(vec![-1., 1., 50.], vec![1., 0.5, 0.25, 0.5, 1., -0.1, 0.25, -0.1, 1.], 8., 6.960998836915657e-16, 1e-30, pdf(dvec![0.9718, 0.1298, 0.8134])); + test_almost(vec![-1., 1., 50.], vec![1., 0.5, 0.25, 0.5, 1., -0.1, 0.25, -0.1, 1.], 8., 7.369987979187023e-16, 1e-30, pdf(dvec![0.4922, 0.5522, 0.7185])); + test_almost(vec![-1., 1., 50.], vec![1., 0.5, 0.25, 0.5, 1., -0.1, 0.25, -0.1, 1.], 8.,6.951631724511314e-16, 1e-30, pdf(dvec![0.3020, 0.1491, 0.5008])); + test_case(vec![-1., 0.], vec![f64::INFINITY, 0., 0., f64::INFINITY], 10., 0., pdf(dvec![10., 10.])); + } + + #[test] + fn test_ln_pdf() { + let ln_pdf = |arg: DVector| move |x: MultivariateStudent| x.ln_pdf(&arg); + test_almost(vec![0., 0.], vec![1., 0., 0., 1.], 4., -3.0542723907338383, 1e-14, ln_pdf(dvec![1., 1.])); + test_almost(vec![0., 0.], vec![1., 0., 0., 1.], 2., -4.3434030034000815, 1e-14, ln_pdf(dvec![1., 2.])); + test_almost(vec![2., 1.], vec![5., 0., 0., 1.], 2.5, -10.542229575274265, 1e-14, ln_pdf(dvec![1., 10.])); + test_almost(vec![-1., 0.], vec![2., 1., 1., 6.], 1.5, -9.650699521198622, 1e-14, ln_pdf(dvec![10., 10.])); + test_case(vec![-1., 0.], vec![f64::INFINITY, 0., 0., f64::INFINITY], 10., f64::NEG_INFINITY, ln_pdf(dvec![10., 10.])); + } + + // #[test] + // fn test_pdf_freedom_large() { + // let pdf_mvs = |mv: MultivariateStudent, arg: DVector| mv.pdf(&arg); + // let pdf_mvn = |mv: MultivariateNormal, arg: DVector| mv.pdf(&arg); + // test_almost_multivariate_normal(vec![0., 0.,], vec![1., 0., 0., 1.], 1e5, 1e-6, dvec![1., 1.], pdf_mvs, pdf_mvn); + // test_almost_multivariate_normal(vec![0., 0.,], vec![1., 0., 0., 1.], 1e10, 1e-7, dvec![1., 1.], pdf_mvs, pdf_mvn); + // test_almost_multivariate_normal(vec![0., 0.,], vec![1., 0., 0., 1.], f64::INFINITY, 1e-300, dvec![1., 1.], pdf_mvs, pdf_mvn); + // test_almost_multivariate_normal(vec![5., -1.,], vec![1., 0.99, 0.99, 1.], f64::INFINITY, 1e-300, dvec![5., 1.], pdf_mvs, pdf_mvn); + // } + // #[test] + // fn test_ln_pdf_freedom_large() { + // let pdf_mvs = |mv: MultivariateStudent, arg: DVector| mv.ln_pdf(&arg); + // let pdf_mvn = |mv: MultivariateNormal, arg: DVector| mv.ln_pdf(&arg); + // test_almost_multivariate_normal(vec![0., 0.,], vec![1., 0., 0., 1.], 1e5, 1e-5, dvec![1., 1.], pdf_mvs, pdf_mvn); + // test_almost_multivariate_normal(vec![0., 0.,], vec![1., 0., 0., 1.], 1e10, 5e-6, dvec![1., 1.], pdf_mvs, pdf_mvn); + // test_almost_multivariate_normal(vec![0., 0.,], vec![1., 0., 0., 1.], f64::INFINITY, 1e-300, dvec![1., 1.], pdf_mvs, pdf_mvn); + // test_almost_multivariate_normal(vec![0., 0.,], vec![1., 0.99, 0.99, 1.], f64::INFINITY, 1e-300, dvec![1., 1.], pdf_mvs, pdf_mvn); + // } +}