From e0dabb4208b517522e6bd1b2b597e05ae6601d78 Mon Sep 17 00:00:00 2001 From: Orion Yeung <11580988+orionyeung001@users.noreply.github.com> Date: Thu, 26 Sep 2024 13:37:38 -0500 Subject: [PATCH] feat!: replace Distribution for new traits for multivariate vector distributions --- src/distribution/multinomial.rs | 42 ++++++++++- src/distribution/multivariate_normal.rs | 79 ++++++++++----------- src/distribution/multivariate_students_t.rs | 49 +++++-------- 3 files changed, 98 insertions(+), 72 deletions(-) diff --git a/src/distribution/multinomial.rs b/src/distribution/multinomial.rs index 0794d352..c128105e 100644 --- a/src/distribution/multinomial.rs +++ b/src/distribution/multinomial.rs @@ -1,7 +1,7 @@ use crate::distribution::Discrete; use crate::function::factorial; use crate::statistics::*; -use nalgebra::{DVector, Dim, Dyn, OMatrix, OVector}; +use nalgebra::{Cholesky, DVector, Dim, DimMin, Dyn, OMatrix, OVector}; /// Implements the /// [Multinomial](https://en.wikipedia.org/wiki/Multinomial_distribution) @@ -260,6 +260,46 @@ where } } +impl CentralMoment for Multinomial +where + D: DimMin, + nalgebra::DefaultAllocator: + nalgebra::allocator::Allocator + nalgebra::allocator::Allocator, +{ + type Mu = OVector; + type Var = Cholesky; + type Kurt = (); + type Skew = (); + + fn mean(&self) -> Self::Mu { + self.p.scale(self.n as f64) + } + + fn variance(&self) -> Self::Var { + let mut cov = OMatrix::from_diagonal(&self.p.map(|x| x * (1.0 - x))); + let mut offdiag = |x: usize, y: usize| { + let elt = -self.p[x] * self.p[y]; + // cov[(x, y)] = elt; + cov[(y, x)] = elt; + }; + + for j in 0..self.p.len() { + for i in 0..j { + offdiag(i, j); + } + } + Cholesky::new_unchecked(cov.scale(self.n as f64)) + } + + fn excess_kurtosis(&self) -> Self::Kurt { + unimplemented!() + } + + fn skewness(&self) -> Self::Skew { + unimplemented!() + } +} + // impl Skewness> for Multinomial { // /// Returns the skewness of the multinomial distribution // /// diff --git a/src/distribution/multivariate_normal.rs b/src/distribution/multivariate_normal.rs index 6eb3b777..658cbb71 100644 --- a/src/distribution/multivariate_normal.rs +++ b/src/distribution/multivariate_normal.rs @@ -1,5 +1,5 @@ use crate::distribution::Continuous; -use crate::statistics::{Max, MeanN, Min, Mode, VarianceN}; +use crate::statistics::{CentralMoment, Max, Min, Mode}; use nalgebra::{Cholesky, Const, DMatrix, DVector, Dim, DimMin, Dyn, OMatrix, OVector}; use std::f64; use std::f64::consts::{E, PI}; @@ -216,14 +216,11 @@ where /// /// where `Σ` is the covariance matrix and `det` is the determinant pub fn entropy(&self) -> Option { - Some( - 0.5 * self - .variance() - .unwrap() - .scale(2. * PI * E) - .determinant() - .ln(), - ) + Some(0.5 * self.variance().scale(2. * PI * E).determinant().ln()) + } + + pub fn cov(&self) -> &OMatrix { + &self.cov } } @@ -264,6 +261,34 @@ where } } +impl CentralMoment for MultivariateNormal +where + D: DimMin, + nalgebra::DefaultAllocator: + nalgebra::allocator::Allocator + nalgebra::allocator::Allocator, +{ + type Mu = OVector; + type Var = OMatrix; + type Kurt = (); + type Skew = (); + + fn mean(&self) -> Self::Mu { + self.mu.clone() + } + + fn variance(&self) -> Self::Var { + self.cov.clone() + } + + fn excess_kurtosis(&self) -> Self::Kurt { + unimplemented!() + } + + fn skewness(&self) -> Self::Skew { + unimplemented!() + } +} + impl Min> for MultivariateNormal where D: Dim, @@ -290,34 +315,6 @@ where } } -impl MeanN> for MultivariateNormal -where - D: Dim, - nalgebra::DefaultAllocator: - nalgebra::allocator::Allocator + nalgebra::allocator::Allocator, -{ - /// Returns the mean of the normal distribution - /// - /// # Remarks - /// - /// This is the same mean used to construct the distribution - fn mean(&self) -> Option> { - Some(self.mu.clone()) - } -} - -impl VarianceN> for MultivariateNormal -where - D: Dim, - nalgebra::DefaultAllocator: - nalgebra::allocator::Allocator + nalgebra::allocator::Allocator, -{ - /// Returns the covariance matrix of the multivariate normal distribution - fn variance(&self) -> Option> { - Some(self.cov.clone()) - } -} - impl Mode> for MultivariateNormal where D: Dim, @@ -379,7 +376,7 @@ mod tests { use crate::{ distribution::{Continuous, MultivariateNormal}, - statistics::{Max, MeanN, Min, Mode, VarianceN}, + statistics::{Max, Min, Mode, CentralMoment}, }; use super::MultivariateNormalError; @@ -404,8 +401,8 @@ mod tests { + nalgebra::allocator::Allocator<(usize, usize), D>, { let mvn = try_create(mean.clone(), covariance.clone()); - assert_eq!(mean, mvn.mean().unwrap()); - assert_eq!(covariance, mvn.variance().unwrap()); + assert_eq!(mean, mvn.mean()); + assert_eq!(covariance, mvn.variance()); } fn bad_create_case(mean: OVector, covariance: OMatrix) @@ -477,7 +474,7 @@ mod tests { #[test] fn test_variance() { - let variance = |x: MultivariateNormal<_>| x.variance().unwrap(); + let variance = |x: MultivariateNormal<_>| x.variance(); test_case( vector![0., 0.], matrix![1., 0.; 0., 1.], diff --git a/src/distribution/multivariate_students_t.rs b/src/distribution/multivariate_students_t.rs index 919e3ebf..bb7e5a1e 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::function::gamma; -use crate::statistics::{Max, MeanN, Min, Mode, VarianceN}; +use crate::statistics::{CentralMoment, Max, Min, Mode}; use nalgebra::{Cholesky, Const, DMatrix, Dim, DimMin, Dyn, OMatrix, OVector}; use std::f64::consts::PI; @@ -259,52 +259,41 @@ where } } -impl MeanN> for MultivariateStudent +impl CentralMoment for MultivariateStudent where - D: Dim, + D: DimMin, nalgebra::DefaultAllocator: nalgebra::allocator::Allocator + nalgebra::allocator::Allocator, { - /// Returns the mean of the student distribution. - /// - /// # Remarks - /// - /// This is the same mean used to construct the distribution if - /// the degrees of freedom is larger than 1. - fn mean(&self) -> Option> { + type Mu = Option>; + type Var = Option>; + type Kurt = (); + type Skew = (); + + fn mean(&self) -> Self::Mu { if self.freedom > 1. { Some(self.location.clone()) } else { None } } -} -impl VarianceN> for MultivariateStudent -where - D: Dim, - nalgebra::DefaultAllocator: - nalgebra::allocator::Allocator + nalgebra::allocator::Allocator, -{ - /// Returns the covariance matrix of the multivariate student distribution. - /// - /// # Formula - /// - /// ```math - /// Σ ⋅ ν / (ν - 2) - /// ``` - /// - /// where `Σ` is the scale matrix and `ν` is the degrees of freedom. - /// Only defined if freedom is larger than 2. - fn variance(&self) -> Option> { + fn variance(&self) -> Self::Var { if self.freedom > 2. { Some(self.scale.clone() * self.freedom / (self.freedom - 2.)) } else { None } } -} + fn excess_kurtosis(&self) -> Self::Kurt { + unimplemented!() + } + + fn skewness(&self) -> Self::Skew { + unimplemented!() + } +} impl Mode> for MultivariateStudent where D: Dim, @@ -399,7 +388,7 @@ mod tests { use crate::{ distribution::{Continuous, MultivariateStudent, MultivariateNormal}, - statistics::{Max, MeanN, Min, Mode, VarianceN}, + statistics::{Max, Min, Mode, CentralMoment}, }; use super::MultivariateStudentError;