diff --git a/src/consts.rs b/src/consts.rs index 32e83a64..a666c9f4 100644 --- a/src/consts.rs +++ b/src/consts.rs @@ -24,5 +24,8 @@ pub const TWO_SQRT_E_OVER_PI: f64 = 1.860382734205265717336249247266663112059421 pub const EULER_MASCHERONI: f64 = 0.5772156649015328606065120900824024310421593359399235988057672348849; +/// Constant value for `zeta(3)` +pub const ZETA_3: f64 = 1.2020569031595942853997381615114499907649862923404988817922715553; + /// Targeted accuracy instantiated over `f64` pub const ACC: f64 = 10e-11; diff --git a/src/distribution/bernoulli.rs b/src/distribution/bernoulli.rs index 21e3e27f..c4b413c6 100644 --- a/src/distribution/bernoulli.rs +++ b/src/distribution/bernoulli.rs @@ -10,13 +10,14 @@ use crate::statistics::*; /// # Examples /// /// ``` -/// use statrs::distribution::{Bernoulli, Discrete}; -/// use statrs::statistics::Distribution; +/// use statrs::distribution::{Bernoulli, BinomialError, Discrete}; +/// use statrs::statistics::*; /// -/// let n = Bernoulli::new(0.5).unwrap(); -/// assert_eq!(n.mean().unwrap(), 0.5); +/// let n = Bernoulli::new(0.5)?; +/// assert_eq!(n.mean(), 0.5); /// assert_eq!(n.pmf(0), 0.5); /// assert_eq!(n.pmf(1), 0.5); +/// # Ok::<(), BinomialError>(()) /// ``` #[derive(Copy, Clone, PartialEq, Debug)] pub struct Bernoulli { @@ -53,10 +54,11 @@ impl Bernoulli { /// # Examples /// /// ``` - /// use statrs::distribution::Bernoulli; + /// use statrs::distribution::{Bernoulli, BinomialError}; /// - /// let n = Bernoulli::new(0.5).unwrap(); + /// let n = Bernoulli::new(0.5)?; /// assert_eq!(n.p(), 0.5); + /// # Ok::<(), BinomialError>(()) /// ``` pub fn p(&self) -> f64 { self.b.p() @@ -68,10 +70,11 @@ impl Bernoulli { /// # Examples /// /// ``` - /// use statrs::distribution::Bernoulli; + /// use statrs::distribution::{Bernoulli, BinomialError}; /// - /// let n = Bernoulli::new(0.5).unwrap(); + /// let n = Bernoulli::new(0.5)?; /// assert_eq!(n.n(), 1); + /// # Ok::<(), BinomialError>(()) /// ``` pub fn n(&self) -> u64 { 1 @@ -160,58 +163,6 @@ impl Max for Bernoulli { } } -impl Distribution for Bernoulli { - /// Returns the mean of the bernoulli - /// distribution - /// - /// # Formula - /// - /// ```text - /// p - /// ``` - fn mean(&self) -> Option { - self.b.mean() - } - - /// Returns the variance of the bernoulli - /// distribution - /// - /// # Formula - /// - /// ```text - /// p * (1 - p) - /// ``` - fn variance(&self) -> Option { - self.b.variance() - } - - /// Returns the entropy of the bernoulli - /// distribution - /// - /// # Formula - /// - /// ```text - /// q = (1 - p) - /// -q * ln(q) - p * ln(p) - /// ``` - fn entropy(&self) -> Option { - self.b.entropy() - } - - /// Returns the skewness of the bernoulli - /// distribution - /// - /// # Formula - /// - /// ```text - /// q = (1 - p) - /// (1 - 2p) / sqrt(p * q) - /// ``` - fn skewness(&self) -> Option { - self.b.skewness() - } -} - impl Median for Bernoulli { /// Returns the median of the bernoulli /// distribution @@ -270,6 +221,65 @@ impl Discrete for Bernoulli { } } +/// returns the mean of a bernoulli variable +/// +/// this is also it's probability parameter + +impl Mean for Bernoulli { + type Mu = f64; + fn mean(&self) -> Self::Mu { + self.p() + } +} + +/// returns the variance of a bernoulli variable +/// +/// # Formula +/// ```text +/// p (1 - p) +/// ``` + +impl Variance for Bernoulli { + type Var = f64; + fn variance(&self) -> Self::Var { + let p = self.p(); + p.mul_add(-p, p) + } +} + +/// Returns the skewness of a bernoulli variable +/// +/// # Formula +/// +/// ```text +/// (1 - 2p) / sqrt(p * (1 - p))) +/// ``` + +impl Skewness for Bernoulli { + type Skew = f64; + fn skewness(&self) -> Self::Skew { + let d = 0.5 - self.p(); + 2.0 * d / (0.25 - d * d).sqrt() + } +} + +/// Returns the excess kurtosis of a bernoulli variable +/// +/// # Formula +/// +/// ```text +/// pq^-1 - 6; pq = p (1-p) +/// ``` + +impl ExcessKurtosis for Bernoulli { + type Kurt = f64; + fn excess_kurtosis(&self) -> Self::Kurt { + let p = self.p(); + let pq = p.mul_add(-p, p); + pq.recip() - 6.0 + } +} + #[rustfmt::skip] #[cfg(test)] mod testing { diff --git a/src/distribution/beta.rs b/src/distribution/beta.rs index 81856468..dd0d1cee 100644 --- a/src/distribution/beta.rs +++ b/src/distribution/beta.rs @@ -8,13 +8,14 @@ use crate::statistics::*; /// # Examples /// /// ``` -/// use statrs::distribution::{Beta, Continuous}; +/// use statrs::distribution::{Beta, Continuous, BetaError}; /// use statrs::statistics::*; /// use statrs::prec; /// -/// let n = Beta::new(2.0, 2.0).unwrap(); -/// assert_eq!(n.mean().unwrap(), 0.5); +/// let n = Beta::new(2.0, 2.0)?; +/// assert_eq!(n.mean(), 0.5); /// assert!(prec::almost_eq(n.pdf(0.5), 1.5, 1e-14)); +/// # Ok::<(), BetaError>(()) /// ``` #[derive(Copy, Clone, PartialEq, Debug)] pub struct Beta { @@ -82,10 +83,11 @@ impl Beta { /// # Examples /// /// ``` - /// use statrs::distribution::Beta; + /// use statrs::distribution::{Beta, BetaError}; /// - /// let n = Beta::new(1.0, 2.0).unwrap(); + /// let n = Beta::new(1.0, 2.0)?; /// assert_eq!(n.shape_a(), 1.0); + /// # Ok::<(), BetaError>(()) /// ``` pub fn shape_a(&self) -> f64 { self.shape_a @@ -96,10 +98,11 @@ impl Beta { /// # Examples /// /// ``` - /// use statrs::distribution::Beta; + /// use statrs::distribution::{Beta, BetaError}; /// - /// let n = Beta::new(1.0, 2.0).unwrap(); + /// let n = Beta::new(1.0, 2.0)?; /// assert_eq!(n.shape_b(), 2.0); + /// # Ok::<(), BetaError>(()) /// ``` pub fn shape_b(&self) -> f64 { self.shape_b @@ -221,38 +224,70 @@ impl Max for Beta { } } -impl Distribution for Beta { - /// Returns the mean of the beta distribution. - /// - /// # Formula - /// - /// ```text - /// α / (α + β) - /// ``` - /// - /// where `α` is shapeA and `β` is shapeB. - fn mean(&self) -> Option { - Some(self.shape_a / (self.shape_a + self.shape_b)) +/// Returns the mean of the beta distribution. +/// +/// # Formula +/// +/// ```text +/// α / (α + β) +/// ``` +/// +/// where `α` is shapeA and `β` is shapeB. +impl Mean for Beta { + type Mu = f64; + fn mean(&self) -> Self::Mu { + self.shape_a / (self.shape_a + self.shape_b) } +} - /// Returns the variance of the beta distribution. - /// - /// # Formula - /// - /// ```text - /// (α * β) / ((α + β)^2 * (α + β + 1)) - /// ``` - /// - /// where `α` is shapeA and `β` is shapeB. - fn variance(&self) -> Option { - Some( - self.shape_a * self.shape_b - / ((self.shape_a + self.shape_b) - * (self.shape_a + self.shape_b) - * (self.shape_a + self.shape_b + 1.0)), - ) +/// Returns the variance of the beta distribution. +/// +/// # Formula +/// +/// ```text +/// (α * β) / ((α + β)^2 * (α + β + 1)) +/// ``` +/// +/// where `α` is shapeA and `β` is shapeB. +impl Variance for Beta { + type Var = f64; + fn variance(&self) -> Self::Var { + self.shape_a * self.shape_b + / ((self.shape_a + self.shape_b) + * (self.shape_a + self.shape_b) + * (self.shape_a + self.shape_b + 1.0)) } +} +/// Returns the skewness of the Beta distribution. +/// +/// # Formula +/// +/// ```text +/// 2(β - α) * sqrt(α + β + 1) / ((α + β + 2) * sqrt(αβ)) +/// ``` +/// +/// where `α` is shapeA and `β` is shapeB. +impl Skewness for Beta { + type Skew = f64; + fn skewness(&self) -> Self::Skew { + 2.0 * (self.shape_b - self.shape_a) * (self.shape_a + self.shape_b + 1.0).sqrt() + / ((self.shape_a + self.shape_b + 2.0) * (self.shape_a * self.shape_b).sqrt()) + } +} + +/// Returns the excess kurtosis of the Beta distribution +impl ExcessKurtosis for Beta { + type Kurt = f64; + fn excess_kurtosis(&self) -> Self::Kurt { + let a = self.shape_a; + let b = self.shape_b; + let numer = 6. * ((a - b).powi(2) * (a + b + 1.) - a * b * (a + b + 2.)); + let denom = a * b * (a + b + 2.) * (a + b + 3.); + numer / denom + } +} +impl Entropy for Beta { /// Returns the entropy of the beta distribution. /// /// # Formula @@ -262,32 +297,13 @@ impl Distribution for Beta { /// ``` /// /// where `α` is shapeA, `β` is shapeB and `ψ` is the digamma function. - fn entropy(&self) -> Option { - Some( - beta::ln_beta(self.shape_a, self.shape_b) - - (self.shape_a - 1.0) * gamma::digamma(self.shape_a) - - (self.shape_b - 1.0) * gamma::digamma(self.shape_b) - + (self.shape_a + self.shape_b - 2.0) * gamma::digamma(self.shape_a + self.shape_b), - ) - } - - /// Returns the skewness of the Beta distribution. - /// - /// # Formula - /// - /// ```text - /// 2(β - α) * sqrt(α + β + 1) / ((α + β + 2) * sqrt(αβ)) - /// ``` - /// - /// where `α` is shapeA and `β` is shapeB. - fn skewness(&self) -> Option { - Some( - 2.0 * (self.shape_b - self.shape_a) * (self.shape_a + self.shape_b + 1.0).sqrt() - / ((self.shape_a + self.shape_b + 2.0) * (self.shape_a * self.shape_b).sqrt()), - ) + fn entropy(&self) -> f64 { + beta::ln_beta(self.shape_a, self.shape_b) + - (self.shape_a - 1.0) * gamma::digamma(self.shape_a) + - (self.shape_b - 1.0) * gamma::digamma(self.shape_b) + + (self.shape_a + self.shape_b - 2.0) * gamma::digamma(self.shape_a + self.shape_b) } } - impl Mode> for Beta { /// Returns the mode of the Beta distribution. Returns `None` if `α <= 1` /// or `β <= 1`. @@ -423,7 +439,7 @@ mod tests { #[test] fn test_mean() { - let f = |x: Beta| x.mean().unwrap(); + let f = |x: Beta| x.mean(); let test = [ ((1.0, 1.0), 0.5), ((9.0, 1.0), 0.9), @@ -436,7 +452,7 @@ mod tests { #[test] fn test_variance() { - let f = |x: Beta| x.variance().unwrap(); + let f = |x: Beta| x.variance(); let test = [ ((1.0, 1.0), 1.0 / 12.0), ((9.0, 1.0), 9.0 / 1100.0), @@ -449,7 +465,7 @@ mod tests { #[test] fn test_entropy() { - let f = |x: Beta| x.entropy().unwrap(); + let f = |x: Beta| x.entropy(); let test = [ ((9.0, 1.0), -1.3083356884473304939016015), ((5.0, 100.0), -2.52016231876027436794592), @@ -462,7 +478,7 @@ mod tests { #[test] fn test_skewness() { - let skewness = |x: Beta| x.skewness().unwrap(); + let skewness = |x: Beta| x.skewness(); test_relative(1.0, 1.0, 0.0, skewness); test_relative(9.0, 1.0, -1.4740554623801777107177478829, skewness); test_relative(5.0, 100.0, 0.817594109275534303545831591, skewness); diff --git a/src/distribution/binomial.rs b/src/distribution/binomial.rs index 1c5e2be2..4b8a3d66 100644 --- a/src/distribution/binomial.rs +++ b/src/distribution/binomial.rs @@ -10,13 +10,14 @@ use std::f64; /// # Examples /// /// ``` -/// use statrs::distribution::{Binomial, Discrete}; -/// use statrs::statistics::Distribution; +/// use statrs::distribution::{Binomial, BinomialError, Discrete}; +/// use statrs::statistics::*; /// -/// let n = Binomial::new(0.5, 5).unwrap(); -/// assert_eq!(n.mean().unwrap(), 2.5); +/// let n = Binomial::new(0.5, 5)?; +/// assert_eq!(n.mean(), 2.5); /// assert_eq!(n.pmf(0), 0.03125); /// assert_eq!(n.pmf(3), 0.3125); +/// # Ok::<(), BinomialError>(()) /// ``` #[derive(Copy, Clone, PartialEq, Debug)] pub struct Binomial { @@ -78,10 +79,11 @@ impl Binomial { /// # Examples /// /// ``` - /// use statrs::distribution::Binomial; + /// use statrs::distribution::{Binomial, BinomialError}; /// - /// let n = Binomial::new(0.5, 5).unwrap(); + /// let n = Binomial::new(0.5, 5)?; /// assert_eq!(n.p(), 0.5); + /// # Ok::<(), BinomialError>(()) /// ``` pub fn p(&self) -> f64 { self.p @@ -93,10 +95,11 @@ impl Binomial { /// # Examples /// /// ``` - /// use statrs::distribution::Binomial; + /// use statrs::distribution::{Binomial, BinomialError}; /// - /// let n = Binomial::new(0.5, 5).unwrap(); + /// let n = Binomial::new(0.5, 5)?; /// assert_eq!(n.n(), 5); + /// # Ok::<(), BinomialError>(()) /// ``` pub fn n(&self) -> u64 { self.n @@ -202,57 +205,90 @@ impl Max for Binomial { } } -impl Distribution for Binomial { - /// Returns the mean of the binomial distribution - /// - /// # Formula - /// - /// ```text - /// p * n - /// ``` - fn mean(&self) -> Option { - Some(self.p * self.n as f64) +/// Returns the mean of the binomial distribution +/// +/// # Formula +/// +/// ```text +/// p * n +/// ``` + +impl Mean for Binomial { + type Mu = f64; + fn mean(&self) -> Self::Mu { + self.n() as f64 * self.p() } +} - /// Returns the variance of the binomial distribution - /// - /// # Formula - /// - /// ```text - /// n * p * (1 - p) - /// ``` - fn variance(&self) -> Option { - Some(self.p * (1.0 - self.p) * self.n as f64) +/// Returns the variance of the binomial distribution +/// +/// # Formula +/// +/// ```text +/// n * p * (1 - p) +/// ``` + +impl Variance for Binomial { + type Var = f64; + fn variance(&self) -> Self::Var { + let n = self.n() as f64; + let p = self.p(); + let pq = p.mul_add(-p, p); + n * pq } +} +/// Returns the skewness of the binomial distribution +/// +/// # Formula +/// +/// ```text +/// (1 - 2p) / sqrt(n * p * (1 - p))) +/// ``` + +impl Skewness for Binomial { + type Skew = f64; + fn skewness(&self) -> Self::Skew { + let n = self.n as f64; + (1.0 - 2.0 * self.p) / (n * self.p * (1.0 - self.p)).sqrt() + } +} + +/// Returns the excess kurtosis of the binomial distribution +/// +/// # Formula +/// +/// ```text +/// (1 - 6 pq) / (n pq); pq = p - p^2 +/// ``` + +impl ExcessKurtosis for Binomial { + type Kurt = f64; + fn excess_kurtosis(&self) -> Self::Kurt { + let n = self.n() as f64; + let pq = self.p() - (1.0 - self.p()); + (1.0 - 6.0 * pq) / (n * pq) + } +} + +impl Entropy for Binomial { /// Returns the entropy of the binomial distribution /// /// # Formula + /// via Stirling approximation to O(1/n) /// /// ```text /// (1 / 2) * ln (2 * π * e * n * p * (1 - p)) /// ``` - fn entropy(&self) -> Option { - let entr = if self.p == 0.0 || ulps_eq!(self.p, 1.0) { + fn entropy(&self) -> f64 { + if self.p == 0.0 || ulps_eq!(self.p, 1.0) { 0.0 } else { (0..self.n + 1).fold(0.0, |acc, x| { let p = self.pmf(x); acc - p * p.ln() }) - }; - Some(entr) - } - - /// Returns the skewness of the binomial distribution - /// - /// # Formula - /// - /// ```text - /// (1 - 2p) / sqrt(n * p * (1 - p))) - /// ``` - fn skewness(&self) -> Option { - Some((1.0 - 2.0 * self.p) / (self.n as f64 * self.p * (1.0 - self.p)).sqrt()) + } } } @@ -377,7 +413,7 @@ mod tests { #[test] fn test_mean() { - let mean = |x: Binomial| x.mean().unwrap(); + let mean = |x: Binomial| x.mean(); test_exact(0.0, 4, 0.0, mean); test_absolute(0.3, 3, 0.9, 1e-15, mean); test_exact(1.0, 2, 2.0, mean); @@ -385,7 +421,7 @@ mod tests { #[test] fn test_variance() { - let variance = |x: Binomial| x.variance().unwrap(); + let variance = |x: Binomial| x.variance(); test_exact(0.0, 4, 0.0, variance); test_exact(0.3, 3, 0.63, variance); test_exact(1.0, 2, 0.0, variance); @@ -393,7 +429,7 @@ mod tests { #[test] fn test_entropy() { - let entropy = |x: Binomial| x.entropy().unwrap(); + let entropy = |x: Binomial| x.entropy(); test_exact(0.0, 4, 0.0, entropy); test_absolute(0.3, 3, 1.1404671643037712668976423399228972051669206536461, 1e-15, entropy); test_exact(1.0, 2, 0.0, entropy); @@ -401,7 +437,7 @@ mod tests { #[test] fn test_skewness() { - let skewness = |x: Binomial| x.skewness().unwrap(); + let skewness = |x: Binomial| x.skewness(); test_exact(0.0, 4, f64::INFINITY, skewness); test_exact(0.3, 3, 0.503952630678969636286, skewness); test_exact(1.0, 2, f64::NEG_INFINITY, skewness); diff --git a/src/distribution/categorical.rs b/src/distribution/categorical.rs index 008e89f9..e26cf626 100644 --- a/src/distribution/categorical.rs +++ b/src/distribution/categorical.rs @@ -10,13 +10,14 @@ use std::f64; /// # Examples /// /// ``` -/// use statrs::distribution::{Categorical, Discrete}; -/// use statrs::statistics::Distribution; +/// use statrs::distribution::{Categorical, CategoricalError, Discrete}; +/// use statrs::statistics::*; /// use statrs::prec; /// -/// let n = Categorical::new(&[0.0, 1.0, 2.0]).unwrap(); -/// assert!(prec::almost_eq(n.mean().unwrap(), 5.0 / 3.0, 1e-15)); +/// let n = Categorical::new(&[0.0, 1.0, 2.0])?; +/// assert!(prec::almost_eq(n.mean(), 5.0 / 3.0, 1e-15)); /// assert_eq!(n.pmf(1), 1.0 / 3.0); +/// # Ok::<(), CategoricalError>(()) /// ``` #[derive(Clone, PartialEq, Debug)] pub struct Categorical { @@ -238,51 +239,53 @@ impl Max for Categorical { } } -impl Distribution for Categorical { - /// Returns the mean of the categorical distribution - /// - /// # Formula - /// - /// ```text - /// Σ(j * p_j) - /// ``` - /// - /// where `p_j` is the `j`th probability mass, - /// `Σ` is the sum from `0` to `k - 1`, - /// and `k` is the number of categories - fn mean(&self) -> Option { - Some( - self.norm_pmf - .iter() - .enumerate() - .fold(0.0, |acc, (idx, &val)| acc + idx as f64 * val), - ) +/// Returns the mean of the categorical distribution +/// +/// # Formula +/// +/// ```text +/// Σ(j * p_j) +/// ``` +/// +/// where `p_j` is the `j`th probability mass, +/// `Σ` is the sum from `0` to `k - 1`, +/// and `k` is the number of categories +impl Mean for Categorical { + type Mu = f64; + fn mean(&self) -> Self::Mu { + self.norm_pmf + .iter() + .enumerate() + .fold(0.0, |acc, (idx, &val)| acc + idx as f64 * val) } +} - /// Returns the variance of the categorical distribution - /// - /// # Formula - /// - /// ```text - /// Σ(p_j * (j - μ)^2) - /// ``` - /// - /// where `p_j` is the `j`th probability mass, `μ` is the mean, - /// `Σ` is the sum from `0` to `k - 1`, - /// and `k` is the number of categories - fn variance(&self) -> Option { - let mu = self.mean()?; - let var = self - .norm_pmf +/// Returns the variance of the categorical distribution, 0-indexed +/// +/// Preserved for backward compatibility +/// +/// # Formula +/// +/// ```text +/// Σ(p_j * (j - μ)^2) +/// ``` +/// +/// where `p_j` is the `j`th probability mass, `μ` is the mean, +/// `Σ` is the sum from `0` to `k - 1`, +/// and `k` is the number of categories +impl Variance for Categorical { + type Var = f64; + fn variance(&self) -> Self::Var { + let mu = self.mean(); + self.norm_pmf .iter() .enumerate() - .fold(0.0, |acc, (idx, &val)| { - let r = idx as f64 - mu; - acc + r * r * val - }); - Some(var) + .map(|(k, &p)| (k as f64 - mu).powi(2) * p) + .sum() } +} +impl Entropy for Categorical { /// Returns the entropy of the categorical distribution /// /// # Formula @@ -294,14 +297,13 @@ impl Distribution for Categorical { /// where `p_j` is the `j`th probability mass, /// `Σ` is the sum from `0` to `k - 1`, /// and `k` is the number of categories - fn entropy(&self) -> Option { - let entr = -self + fn entropy(&self) -> f64 { + -self .norm_pmf .iter() .filter(|&&p| p > 0.0) .map(|p| p * p.ln()) - .sum::(); - Some(entr) + .sum::() } } impl Median for Categorical { @@ -434,7 +436,7 @@ mod tests { #[test] fn test_mean() { - let mean = |x: Categorical| x.mean().unwrap(); + let mean = |x: Categorical| x.mean(); test_exact(&[0.0, 0.25, 0.5, 0.25], 2.0, mean); test_exact(&[0.0, 1.0, 2.0, 1.0], 2.0, mean); test_exact(&[0.0, 0.5, 0.5], 1.5, mean); @@ -444,7 +446,7 @@ mod tests { #[test] fn test_variance() { - let variance = |x: Categorical| x.variance().unwrap(); + let variance = |x: Categorical| x.variance(); test_exact(&[0.0, 0.25, 0.5, 0.25], 0.5, variance); test_exact(&[0.0, 1.0, 2.0, 1.0], 0.5, variance); test_exact(&[0.0, 0.5, 0.5], 0.25, variance); @@ -454,7 +456,7 @@ mod tests { #[test] fn test_entropy() { - let entropy = |x: Categorical| x.entropy().unwrap(); + let entropy = |x: Categorical| x.entropy(); test_exact(&[0.0, 1.0], 0.0, entropy); test_absolute(&[0.0, 1.0, 1.0], 2f64.ln(), 1e-15, entropy); test_absolute(&[1.0, 1.0, 1.0], 3f64.ln(), 1e-15, entropy); diff --git a/src/distribution/cauchy.rs b/src/distribution/cauchy.rs index 7353f02d..32e5b37d 100644 --- a/src/distribution/cauchy.rs +++ b/src/distribution/cauchy.rs @@ -8,12 +8,13 @@ use std::f64; /// # Examples /// /// ``` -/// use statrs::distribution::{Cauchy, Continuous}; -/// use statrs::statistics::Mode; +/// use statrs::distribution::{Cauchy, Continuous, CauchyError}; +/// use statrs::statistics::*; /// -/// let n = Cauchy::new(0.0, 1.0).unwrap(); +/// let n = Cauchy::new(0.0, 1.0)?; /// assert_eq!(n.mode().unwrap(), 0.0); /// assert_eq!(n.pdf(1.0), 0.1591549430918953357689); +/// # Ok::<(), CauchyError>(()) /// ``` #[derive(Copy, Clone, PartialEq, Debug)] pub struct Cauchy { @@ -80,10 +81,11 @@ impl Cauchy { /// # Examples /// /// ``` - /// use statrs::distribution::Cauchy; + /// use statrs::distribution::{Cauchy, CauchyError}; /// - /// let n = Cauchy::new(0.0, 1.0).unwrap(); + /// let n = Cauchy::new(0.0, 1.0)?; /// assert_eq!(n.location(), 0.0); + /// # Ok::<(), CauchyError>(()) /// ``` pub fn location(&self) -> f64 { self.location @@ -94,10 +96,11 @@ impl Cauchy { /// # Examples /// /// ``` - /// use statrs::distribution::Cauchy; + /// use statrs::distribution::{Cauchy, CauchyError}; /// - /// let n = Cauchy::new(0.0, 1.0).unwrap(); + /// let n = Cauchy::new(0.0, 1.0)?; /// assert_eq!(n.scale(), 1.0); + /// # Ok::<(), CauchyError>(()) /// ``` pub fn scale(&self) -> f64 { self.scale @@ -196,7 +199,13 @@ impl Max for Cauchy { } } -impl Distribution for Cauchy { +/// This is implemented so that there are docs stating that Cauchy does not have the convergence for a mean. +impl Mean for Cauchy { + type Mu = (); + fn mean(&self) -> Self::Mu {} +} + +impl Entropy for Cauchy { /// Returns the entropy of the cauchy distribution /// /// # Formula @@ -206,8 +215,8 @@ impl Distribution for Cauchy { /// ``` /// /// where `γ` is the scale - fn entropy(&self) -> Option { - Some((4.0 * f64::consts::PI * self.scale).ln()) + fn entropy(&self) -> f64 { + (4.0 * f64::consts::PI * self.scale).ln() } } @@ -311,7 +320,7 @@ mod tests { #[test] fn test_entropy() { - let entropy = |x: Cauchy| x.entropy().unwrap(); + let entropy = |x: Cauchy| x.entropy(); test_exact(0.0, 2.0, 3.224171427529236102395, entropy); test_exact(0.1, 4.0, 3.917318608089181411812, entropy); test_exact(1.0, 10.0, 4.833609339963336476996, entropy); diff --git a/src/distribution/chi.rs b/src/distribution/chi.rs index 8eb50603..75729b3e 100644 --- a/src/distribution/chi.rs +++ b/src/distribution/chi.rs @@ -9,13 +9,14 @@ use std::f64; /// # Examples /// /// ``` -/// use statrs::distribution::{Chi, Continuous}; -/// use statrs::statistics::Distribution; +/// use statrs::distribution::{Chi, Continuous, ChiError}; +/// use statrs::statistics::*; /// use statrs::prec; /// -/// let n = Chi::new(2.0).unwrap(); +/// let n = Chi::new(2.0)?; /// assert!(prec::almost_eq(n.mean().unwrap(), 1.25331413731550025121, 1e-14)); /// assert!(prec::almost_eq(n.pdf(1.0), 0.60653065971263342360, 1e-15)); +/// # Ok::<(), ChiError>(()) /// ``` #[derive(Copy, Clone, PartialEq, Debug)] pub struct Chi { @@ -77,10 +78,12 @@ impl Chi { /// # Examples /// /// ``` - /// use statrs::distribution::Chi; + /// use statrs::distribution::{Chi, ChiError}; + /// /// - /// let n = Chi::new(2.0).unwrap(); + /// let n = Chi::new(2.0)?; /// assert_eq!(n.freedom(), 2.0); + /// # Ok::<(), ChiError>(()) /// ``` pub fn freedom(&self) -> f64 { self.freedom @@ -177,21 +180,23 @@ impl Max for Chi { } } -impl Distribution for Chi { - /// Returns the mean of the chi distribution - /// - /// # Remarks - /// - /// Returns `NaN` if `freedom` is `INF` - /// - /// # Formula - /// - /// ```text - /// sqrt2 * Γ((k + 1) / 2) / Γ(k / 2) - /// ``` - /// - /// where `k` is degrees of freedom and `Γ` is the gamma function - fn mean(&self) -> Option { +/// Returns the mean of the chi distribution +/// +/// # Remarks +/// +/// Returns `None` if `freedom` is `INF` +/// +/// # Formula +/// +/// ```text +/// sqrt2 * Γ((k + 1) / 2) / Γ(k / 2) +/// ``` +/// +/// where `k` is degrees of freedom and `Γ` is the gamma function + +impl Mean for Chi { + type Mu = Option; + fn mean(&self) -> Self::Mu { if self.freedom.is_infinite() { None } else if self.freedom > 300.0 { @@ -212,26 +217,78 @@ impl Distribution for Chi { Some(mean) } } +} - /// Returns the variance of the chi distribution - /// - /// # Remarks - /// - /// Returns `NaN` if `freedom` is `INF` - /// - /// # Formula - /// - /// ```text - /// k - μ^2 - /// ``` - /// - /// where `k` is degrees of freedom and `μ` is the mean - /// of the distribution - fn variance(&self) -> Option { +/// Returns the variance of the chi distribution +/// +/// # Remarks +/// +/// Returns `None` if `freedom` is `INF` +/// +/// # Formula +/// +/// ```text +/// k - μ^2 +/// ``` +/// +/// where `k` is degrees of freedom and `μ` is the mean +/// of the distribution + +impl Variance for Chi { + type Var = Option; + fn variance(&self) -> Self::Var { let mean = self.mean()?; Some(self.freedom - mean * mean) } +} + +/// Returns the skewness of the chi distribution +/// +/// # Remarks +/// +/// Returns `None` if `freedom` is `INF` +/// +/// # Formula +/// +/// ```text +/// (μ / σ^3) * (1 - 2σ^2) +/// ``` +/// where `μ` is the mean and `σ` the standard deviation +/// of the distribution + +impl Skewness for Chi { + type Skew = Option; + fn skewness(&self) -> Self::Skew { + let sigma = self.variance()?.sqrt(); + let skew = self.mean()? * (1.0 - 2.0 * sigma * sigma) / (sigma * sigma * sigma); + Some(skew) + } +} + +/// Returns the excess kurtosis of the chi distribution +/// +/// # Remarks +/// +/// Returns `None` if `freedom` is `INF` +/// +/// # Formula +/// +/// ```text +/// 2 / σ^2 (1 - μσγ - σ^2) +/// ``` +/// where `μ` is the mean, `σ` the standard deviation, and `γ` the skewness +/// of the distribution + +impl ExcessKurtosis for Chi { + type Kurt = Option; + fn excess_kurtosis(&self) -> Self::Kurt { + let sig_sqr = self.variance()?; + let mu_over_var = self.mean()? / sig_sqr; + Some(2. / sig_sqr - 2. * mu_over_var.powi(2) * (1. - 2. * sig_sqr) - 2.) + } +} +impl Entropy> for Chi { /// Returns the entropy of the chi distribution /// /// # Remarks @@ -257,25 +314,6 @@ impl Distribution for Chi { / 2.0; Some(entr) } - - /// Returns the skewness of the chi distribution - /// - /// # Remarks - /// - /// Returns `NaN` if `freedom` is `INF` - /// - /// # Formula - /// - /// ```text - /// (μ / σ^3) * (1 - 2σ^2) - /// ``` - /// where `μ` is the mean and `σ` the standard deviation - /// of the distribution - fn skewness(&self) -> Option { - let sigma = self.std_dev()?; - let skew = self.mean()? * (1.0 - 2.0 * sigma * sigma) / (sigma * sigma * sigma); - Some(skew) - } } impl Mode> for Chi { diff --git a/src/distribution/chi_squared.rs b/src/distribution/chi_squared.rs index bd31bcf8..d0cb939d 100644 --- a/src/distribution/chi_squared.rs +++ b/src/distribution/chi_squared.rs @@ -11,13 +11,14 @@ use std::f64; /// # Examples /// /// ``` -/// use statrs::distribution::{ChiSquared, Continuous}; -/// use statrs::statistics::Distribution; +/// use statrs::distribution::{ChiSquared, Continuous, GammaError}; +/// use statrs::statistics::*; /// use statrs::prec; /// -/// let n = ChiSquared::new(3.0).unwrap(); -/// assert_eq!(n.mean().unwrap(), 3.0); +/// let n = ChiSquared::new(3.0)?; +/// assert_eq!(n.mean(), 3.0); /// assert!(prec::almost_eq(n.pdf(4.0), 0.107981933026376103901, 1e-15)); +/// # Ok::<(), GammaError>(()) /// ``` #[derive(Copy, Clone, PartialEq, Debug)] pub struct ChiSquared { @@ -56,10 +57,11 @@ impl ChiSquared { /// # Examples /// /// ``` - /// use statrs::distribution::ChiSquared; + /// use statrs::distribution::{ChiSquared, GammaError}; /// - /// let n = ChiSquared::new(3.0).unwrap(); + /// let n = ChiSquared::new(3.0)?; /// assert_eq!(n.freedom(), 3.0); + /// # Ok::<(), GammaError>(()) /// ``` pub fn freedom(&self) -> f64 { self.freedom @@ -70,10 +72,11 @@ impl ChiSquared { /// # Examples /// /// ``` - /// use statrs::distribution::ChiSquared; + /// use statrs::distribution::{ChiSquared, GammaError}; /// - /// let n = ChiSquared::new(3.0).unwrap(); + /// let n = ChiSquared::new(3.0)?; /// assert_eq!(n.shape(), 3.0 / 2.0); + /// # Ok::<(), GammaError>(()) /// ``` pub fn shape(&self) -> f64 { self.g.shape() @@ -84,10 +87,11 @@ impl ChiSquared { /// # Examples /// /// ``` - /// use statrs::distribution::ChiSquared; + /// use statrs::distribution::{ChiSquared, GammaError}; /// - /// let n = ChiSquared::new(3.0).unwrap(); + /// let n = ChiSquared::new(3.0)?; /// assert_eq!(n.rate(), 0.5); + /// # Ok::<(), GammaError>(()) /// ``` pub fn rate(&self) -> f64 { self.g.rate() @@ -185,33 +189,62 @@ impl Max for ChiSquared { } } -impl Distribution for ChiSquared { - /// Returns the mean of the chi-squared distribution - /// - /// # Formula - /// - /// ```text - /// k - /// ``` - /// - /// where `k` is the degrees of freedom - fn mean(&self) -> Option { +/// Returns the mean of the chi-squared distribution +/// +/// # Formula +/// +/// ```text +/// k +/// ``` +/// +/// where `k` is the degrees of freedom +impl Mean for ChiSquared { + type Mu = ::Mu; + fn mean(&self) -> Self::Mu { self.g.mean() } +} - /// Returns the variance of the chi-squared distribution - /// - /// # Formula - /// - /// ```text - /// 2k - /// ``` - /// - /// where `k` is the degrees of freedom - fn variance(&self) -> Option { +/// Returns the variance of the chi-squared distribution +/// +/// # Formula +/// +/// ```text +/// 2k +/// ``` +/// +/// where `k` is the degrees of freedom +impl Variance for ChiSquared { + type Var = ::Var; + fn variance(&self) -> Self::Var { self.g.variance() } +} + +/// Returns the skewness of the chi-squared distribution +/// +/// # Formula +/// +/// ```text +/// sqrt(8 / k) +/// ``` +/// +/// where `k` is the degrees of freedom +impl Skewness for ChiSquared { + type Skew = ::Skew; + fn skewness(&self) -> Self::Skew { + self.g.skewness() + } +} + +impl ExcessKurtosis for ChiSquared { + type Kurt = ::Kurt; + fn excess_kurtosis(&self) -> Self::Kurt { + self.g.excess_kurtosis() + } +} +impl Entropy for ChiSquared { /// Returns the entropy of the chi-squared distribution /// /// # Formula @@ -222,22 +255,9 @@ impl Distribution for ChiSquared { /// /// where `k` is the degrees of freedom, `Γ` is the gamma function, /// and `ψ` is the digamma function - fn entropy(&self) -> Option { + fn entropy(&self) -> f64 { self.g.entropy() } - - /// Returns the skewness of the chi-squared distribution - /// - /// # Formula - /// - /// ```text - /// sqrt(8 / k) - /// ``` - /// - /// where `k` is the degrees of freedom - fn skewness(&self) -> Option { - self.g.skewness() - } } impl Median for ChiSquared { diff --git a/src/distribution/dirac.rs b/src/distribution/dirac.rs index 8bd782f6..056213b4 100644 --- a/src/distribution/dirac.rs +++ b/src/distribution/dirac.rs @@ -7,11 +7,12 @@ use crate::statistics::*; /// # Examples /// /// ``` -/// use statrs::distribution::{Dirac, Continuous}; -/// use statrs::statistics::Distribution; +/// use statrs::distribution::{Dirac, Continuous, DiracError}; +/// use statrs::statistics::*; /// -/// let n = Dirac::new(3.0).unwrap(); -/// assert_eq!(n.mean().unwrap(), 3.0); +/// let n = Dirac::new(3.0)?; +/// assert_eq!(n.mean(), 3.0); +/// # Ok::<(), DiracError>(()) /// ``` #[derive(Debug, Copy, Clone, PartialEq)] pub struct Dirac(f64); @@ -130,52 +131,65 @@ impl Max for Dirac { } } -impl Distribution for Dirac { - /// Returns the mean of the dirac distribution - /// - /// # Remarks - /// - /// Since the only value that can be produced by this distribution is `v` with probability - /// 1, it is just `v`. - fn mean(&self) -> Option { - Some(self.0) - } +/// Returns the mean of the dirac distribution +/// +/// # Remarks +/// +/// Since the only value that can be produced by this distribution is `v` with probability +/// 1, it is just `v`. - /// Returns the variance of the dirac distribution - /// - /// # Formula - /// - /// ```text - /// 0 - /// ``` - /// - /// Since only one value can be produced there is no variance. - fn variance(&self) -> Option { - Some(0.0) +impl Mean for Dirac { + type Mu = f64; + fn mean(&self) -> Self::Mu { + self.0 } +} - /// Returns the entropy of the dirac distribution - /// - /// # Formula - /// - /// ```text - /// 0 - /// ``` - /// - /// Since this distribution has full certainty, it encodes no information - fn entropy(&self) -> Option { - Some(0.0) +/// Returns the variance of the dirac distribution +/// +/// # Formula +/// +/// ```text +/// 0 +/// ``` +/// +/// Since only one value can be produced there is no variance. + +impl Variance for Dirac { + type Var = f64; + fn variance(&self) -> Self::Var { + 0.0 } +} + +/// Returns the skewness of the dirac distribution +/// +/// # Formula +/// +/// ```text +/// 0 +/// ``` + +impl Skewness for Dirac { + type Skew = (); + fn skewness(&self) -> Self::Skew {} +} + +impl ExcessKurtosis for Dirac { + type Kurt = (); + fn excess_kurtosis(&self) -> Self::Kurt {} +} - /// Returns the skewness of the dirac distribution +impl Entropy for Dirac { + /// Returns the entropy of the dirac distribution /// /// # Formula /// /// ```text /// 0 /// ``` - fn skewness(&self) -> Option { - Some(0.0) + fn entropy(&self) -> f64 { + 0.0 } } @@ -233,7 +247,7 @@ mod tests { #[test] fn test_variance() { - let variance = |x: Dirac| x.variance().unwrap(); + let variance = |x: Dirac| x.variance(); test_exact(0.0, 0.0, variance); test_exact(-5.0, 0.0, variance); test_exact(f64::INFINITY, 0.0, variance); @@ -241,20 +255,11 @@ mod tests { #[test] fn test_entropy() { - let entropy = |x: Dirac| x.entropy().unwrap(); + let entropy = |x: Dirac| x.entropy(); test_exact(0.0, 0.0, entropy); test_exact(f64::INFINITY, 0.0, entropy); } - #[test] - fn test_skewness() { - let skewness = |x: Dirac| x.skewness().unwrap(); - test_exact(0.0, 0.0, skewness); - test_exact(4.0, 0.0, skewness); - test_exact(0.3, 0.0, skewness); - test_exact(f64::INFINITY, 0.0, skewness); - } - #[test] fn test_mode() { let mode = |x: Dirac| x.mode().unwrap(); diff --git a/src/distribution/dirichlet.rs b/src/distribution/dirichlet.rs index 5e08d0a4..69a424ef 100644 --- a/src/distribution/dirichlet.rs +++ b/src/distribution/dirichlet.rs @@ -12,14 +12,14 @@ use std::f64; /// # Examples /// /// ``` -/// use statrs::distribution::{Dirichlet, Continuous}; -/// use statrs::statistics::Distribution; +/// use statrs::distribution::{Dirichlet, Continuous, DirichletError}; +/// use statrs::statistics::*; /// use nalgebra::DVector; -/// use statrs::statistics::MeanN; /// -/// let n = Dirichlet::new(vec![1.0, 2.0, 3.0]).unwrap(); -/// assert_eq!(n.mean().unwrap(), DVector::from_vec(vec![1.0 / 6.0, 1.0 / 3.0, 0.5])); +/// let n = Dirichlet::new(vec![1.0, 2.0, 3.0])?; +/// assert_eq!(n.mean(), DVector::from_vec(vec![1.0 / 6.0, 1.0 / 3.0, 0.5])); /// assert_eq!(n.pdf(&DVector::from_vec(vec![0.33333, 0.33333, 0.33333])), 2.222155556222205); +/// # Ok::<(), DirichletError>(()) /// ``` #[derive(Clone, PartialEq, Debug)] pub struct Dirichlet @@ -138,11 +138,12 @@ where /// # Examples /// /// ``` - /// use statrs::distribution::Dirichlet; + /// use statrs::distribution::{Dirichlet, DirichletError}; /// use nalgebra::DVector; /// - /// let n = Dirichlet::new(vec![1.0, 2.0, 3.0]).unwrap(); + /// let n = Dirichlet::new(vec![1.0, 2.0, 3.0])?; /// assert_eq!(n.alpha(), &DVector::from_vec(vec![1.0, 2.0, 3.0])); + /// # Ok::<(), DirichletError>(()) /// ``` pub fn alpha(&self) -> &nalgebra::OVector { &self.alpha @@ -212,11 +213,12 @@ where } } -impl MeanN> for Dirichlet +impl Mean for Dirichlet where D: Dim, nalgebra::DefaultAllocator: nalgebra::allocator::Allocator, { + type Mu = OVector; /// Returns the means of the dirichlet distribution /// /// # Formula @@ -227,18 +229,19 @@ where /// /// for the `i`th element where `α_i` is the `i`th concentration parameter /// and `α_0` is the sum of all concentration parameters - fn mean(&self) -> Option> { + fn mean(&self) -> Self::Mu { let sum = self.alpha_sum(); - Some(self.alpha.map(|x| x / sum)) + self.alpha.map(|x| x / sum) } } -impl VarianceN> for Dirichlet +impl Variance for Dirichlet where D: Dim, nalgebra::DefaultAllocator: nalgebra::allocator::Allocator + nalgebra::allocator::Allocator, { + type Var = OMatrix; /// Returns the variances of the dirichlet distribution /// /// # Formula @@ -249,7 +252,7 @@ where /// /// for the `i`th element where `α_i` is the `i`th concentration parameter /// and `α_0` is the sum of all concentration parameters - fn variance(&self) -> Option> { + fn variance(&self) -> Self::Var { let sum = self.alpha_sum(); let normalizing = sum * sum * (sum + 1.0); let mut cov = OMatrix::from_diagonal(&self.alpha.map(|x| x * (sum - x) / normalizing)); @@ -263,7 +266,7 @@ where offdiag(i, j); } } - Some(cov) + cov } } @@ -426,7 +429,7 @@ mod tests { #[test] fn test_mean() { - let mean = |dd: Dirichlet<_>| dd.mean().unwrap(); + let mean = |dd: Dirichlet<_>| dd.mean(); test_almost(vec![0.5; 5].into(), vec![1.0 / 5.0; 5].into(), 1e-15, mean); @@ -447,7 +450,7 @@ mod tests { #[test] fn test_variance() { - let variance = |dd: Dirichlet<_>| dd.variance().unwrap(); + let variance = |dd: Dirichlet<_>| dd.variance(); test_almost( dvector![1.0, 2.0], diff --git a/src/distribution/discrete_uniform.rs b/src/distribution/discrete_uniform.rs index a9760808..382c5104 100644 --- a/src/distribution/discrete_uniform.rs +++ b/src/distribution/discrete_uniform.rs @@ -8,12 +8,13 @@ use crate::statistics::*; /// # Examples /// /// ``` -/// use statrs::distribution::{DiscreteUniform, Discrete}; -/// use statrs::statistics::Distribution; +/// use statrs::distribution::{DiscreteUniform, Discrete, DiscreteUniformError}; +/// use statrs::statistics::*; /// -/// let n = DiscreteUniform::new(0, 5).unwrap(); -/// assert_eq!(n.mean().unwrap(), 2.5); +/// let n = DiscreteUniform::new(0, 5)?; +/// assert_eq!(n.mean(), 2.5); /// assert_eq!(n.pmf(3), 1.0 / 6.0); +/// # Ok::<(), DiscreteUniformError>(()) /// ``` #[derive(Copy, Clone, PartialEq, Debug)] pub struct DiscreteUniform { @@ -66,6 +67,11 @@ impl DiscreteUniform { Ok(DiscreteUniform { min, max }) } } + + pub fn std_dev(&self) -> f64 { + let dist = self.max - self.min + 1; + (dist * dist) as f64 / 12. + } } impl std::fmt::Display for DiscreteUniform { @@ -159,30 +165,62 @@ impl Max for DiscreteUniform { } } -impl Distribution for DiscreteUniform { - /// Returns the mean of the discrete uniform distribution - /// - /// # Formula - /// - /// ```text - /// (min + max) / 2 - /// ``` - fn mean(&self) -> Option { - Some((self.min + self.max) as f64 / 2.0) +/// Returns the mean of the discrete uniform distribution +/// +/// # Formula +/// +/// ```text +/// (min + max) / 2 +/// ``` + +impl Mean for DiscreteUniform { + type Mu = f64; + fn mean(&self) -> Self::Mu { + (self.min + self.max) as f64 / 2.0 } +} - /// Returns the variance of the discrete uniform distribution - /// - /// # Formula - /// - /// ```text - /// ((max - min + 1)^2 - 1) / 12 - /// ``` - fn variance(&self) -> Option { +/// Returns the variance of the discrete uniform distribution +/// +/// # Formula +/// +/// ```text +/// ((max - min + 1)^2 - 1) / 12 +/// ``` + +impl Variance for DiscreteUniform { + type Var = f64; + fn variance(&self) -> Self::Var { let diff = (self.max - self.min) as f64; - Some(((diff + 1.0) * (diff + 1.0) - 1.0) / 12.0) + ((diff + 1.0) * (diff + 1.0) - 1.0) / 12.0 + } +} + +/// Returns the skewness of the discrete uniform distribution +/// +/// # Formula +/// +/// ```text +/// 0 +/// ``` + +impl Skewness for DiscreteUniform { + type Skew = f64; + fn skewness(&self) -> Self::Skew { + 0.0 + } +} + +/// docs + +impl ExcessKurtosis for DiscreteUniform { + type Kurt = f64; + fn excess_kurtosis(&self) -> Self::Kurt { + -1.2 } +} +impl Entropy for DiscreteUniform { /// Returns the entropy of the discrete uniform distribution /// /// # Formula @@ -190,20 +228,9 @@ impl Distribution for DiscreteUniform { /// ```text /// ln(max - min + 1) /// ``` - fn entropy(&self) -> Option { + fn entropy(&self) -> f64 { let diff = (self.max - self.min) as f64; - Some((diff + 1.0).ln()) - } - - /// Returns the skewness of the discrete uniform distribution - /// - /// # Formula - /// - /// ```text - /// 0 - /// ``` - fn skewness(&self) -> Option { - Some(0.0) + (diff + 1.0).ln() } } @@ -304,7 +331,7 @@ mod tests { #[test] fn test_mean() { - let mean = |x: DiscreteUniform| x.mean().unwrap(); + let mean = |x: DiscreteUniform| x.mean(); test_exact(-10, 10, 0.0, mean); test_exact(0, 4, 2.0, mean); test_exact(10, 20, 15.0, mean); @@ -313,7 +340,7 @@ mod tests { #[test] fn test_variance() { - let variance = |x: DiscreteUniform| x.variance().unwrap(); + let variance = |x: DiscreteUniform| x.variance(); test_exact(-10, 10, 36.66666666666666666667, variance); test_exact(0, 4, 2.0, variance); test_exact(10, 20, 10.0, variance); @@ -322,7 +349,7 @@ mod tests { #[test] fn test_entropy() { - let entropy = |x: DiscreteUniform| x.entropy().unwrap(); + let entropy = |x: DiscreteUniform| x.entropy(); test_exact(-10, 10, 3.0445224377234229965005979803657054342845752874046093, entropy); test_exact(0, 4, 1.6094379124341003746007593332261876395256013542685181, entropy); test_exact(10, 20, 2.3978952727983705440619435779651292998217068539374197, entropy); @@ -331,7 +358,7 @@ mod tests { #[test] fn test_skewness() { - let skewness = |x: DiscreteUniform| x.skewness().unwrap(); + let skewness = |x: DiscreteUniform| x.skewness(); test_exact(-10, 10, 0.0, skewness); test_exact(0, 4, 0.0, skewness); test_exact(10, 20, 0.0, skewness); diff --git a/src/distribution/empirical.rs b/src/distribution/empirical.rs index b2ab498d..9545e637 100644 --- a/src/distribution/empirical.rs +++ b/src/distribution/empirical.rs @@ -27,7 +27,7 @@ impl Ord for NonNan { /// /// ``` /// use statrs::distribution::{Continuous, Empirical}; -/// use statrs::statistics::Distribution; +/// use statrs::statistics::*; /// /// let samples = vec![0.0, 5.0, 10.0]; /// @@ -199,12 +199,17 @@ impl Min for Empirical { } } -impl Distribution for Empirical { - fn mean(&self) -> Option { +impl Mean for Empirical { + type Mu = Option; + fn mean(&self) -> Self::Mu { self.mean_and_var.map(|(mean, _)| mean) } +} + +impl Variance for Empirical { + type Var = Option; - fn variance(&self) -> Option { + fn variance(&self) -> Self::Var { self.mean_and_var.map(|(_, var)| var / (self.sum - 1.)) } } diff --git a/src/distribution/erlang.rs b/src/distribution/erlang.rs index 89e5e4ec..97d1b264 100644 --- a/src/distribution/erlang.rs +++ b/src/distribution/erlang.rs @@ -10,13 +10,14 @@ use crate::statistics::*; /// # Examples /// /// ``` -/// use statrs::distribution::{Erlang, Continuous}; -/// use statrs::statistics::Distribution; +/// use statrs::distribution::{Erlang, Continuous, GammaError}; +/// use statrs::statistics::*; /// use statrs::prec; /// -/// let n = Erlang::new(3, 1.0).unwrap(); -/// assert_eq!(n.mean().unwrap(), 3.0); +/// let n = Erlang::new(3, 1.0)?; +/// assert_eq!(n.mean(), 3.0); /// assert!(prec::almost_eq(n.pdf(2.0), 0.270670566473225383788, 1e-15)); +/// # Ok::<(), GammaError>(()) /// ``` #[derive(Copy, Clone, PartialEq, Debug)] pub struct Erlang { @@ -52,10 +53,11 @@ impl Erlang { /// # Examples /// /// ``` - /// use statrs::distribution::Erlang; + /// use statrs::distribution::{Erlang, GammaError}; /// - /// let n = Erlang::new(3, 1.0).unwrap(); + /// let n = Erlang::new(3, 1.0)?; /// assert_eq!(n.shape(), 3); + /// # Ok::<(), GammaError>(()) /// ``` pub fn shape(&self) -> u64 { self.g.shape() as u64 @@ -66,10 +68,11 @@ impl Erlang { /// # Examples /// /// ``` - /// use statrs::distribution::Erlang; + /// use statrs::distribution::{Erlang, GammaError}; /// /// let n = Erlang::new(3, 1.0).unwrap(); /// assert_eq!(n.rate(), 1.0); + /// # Ok::<(), GammaError>(()) /// ``` pub fn rate(&self) -> f64 { self.g.rate() @@ -169,38 +172,68 @@ impl Max for Erlang { } } -impl Distribution for Erlang { - /// Returns the mean of the erlang distribution - /// - /// # Remarks - /// - /// Returns `shape` if `rate == f64::INFINITY`. This behavior - /// is borrowed from the Math.NET implementation - /// - /// # Formula - /// - /// ```text - /// k / λ - /// ``` - /// - /// where `k` is the shape and `λ` is the rate - fn mean(&self) -> Option { +/// Returns the mean of the erlang distribution +/// +/// # Remarks +/// +/// Returns `shape` if `rate == f64::INFINITY`. This behavior +/// is borrowed from the Math.NET implementation +/// +/// # Formula +/// +/// ```text +/// k / λ +/// ``` +/// +/// where `k` is the shape and `λ` is the rate +impl Mean for Erlang { + type Mu = ::Mu; + fn mean(&self) -> Self::Mu { self.g.mean() } +} - /// Returns the variance of the erlang distribution - /// - /// # Formula - /// - /// ```text - /// k / λ^2 - /// ``` - /// - /// where `α` is the shape and `λ` is the rate - fn variance(&self) -> Option { +/// Returns the variance of the erlang distribution +/// +/// # Formula +/// +/// ```text +/// k / λ^2 +/// ``` +/// +/// where `α` is the shape and `λ` is the rate +impl Variance for Erlang { + type Var = ::Var; + fn variance(&self) -> Self::Var { self.g.variance() } +} + +/// Returns the skewness of the erlang distribution +/// +/// # Formula +/// +/// ```text +/// 2 / sqrt(k) +/// ``` +/// +/// where `k` is the shape +impl Skewness for Erlang { + type Skew = ::Skew; + fn skewness(&self) -> Self::Skew { + self.g.skewness() + } +} + +/// docs +impl ExcessKurtosis for Erlang { + type Kurt = ::Kurt; + fn excess_kurtosis(&self) -> Self::Kurt { + self.g.excess_kurtosis() + } +} +impl Entropy for Erlang { /// Returns the entropy of the erlang distribution /// /// # Formula @@ -211,22 +244,9 @@ impl Distribution for Erlang { /// /// where `k` is the shape, `λ` is the rate, `Γ` is the gamma function, /// and `ψ` is the digamma function - fn entropy(&self) -> Option { + fn entropy(&self) -> f64 { self.g.entropy() } - - /// Returns the skewness of the erlang distribution - /// - /// # Formula - /// - /// ```text - /// 2 / sqrt(k) - /// ``` - /// - /// where `k` is the shape - fn skewness(&self) -> Option { - self.g.skewness() - } } impl Mode> for Erlang { diff --git a/src/distribution/exponential.rs b/src/distribution/exponential.rs index 556812a5..18364a93 100644 --- a/src/distribution/exponential.rs +++ b/src/distribution/exponential.rs @@ -11,13 +11,15 @@ use std::f64; /// # Examples /// /// ``` -/// use statrs::distribution::{Exp, Continuous}; -/// use statrs::statistics::Distribution; +/// use statrs::distribution::{Exp, Continuous, ExpError}; +/// use statrs::statistics::*; /// -/// let n = Exp::new(1.0).unwrap(); -/// assert_eq!(n.mean().unwrap(), 1.0); +/// let n = Exp::new(1.0)?; +/// assert_eq!(n.mean(), 1.0); /// assert_eq!(n.pdf(1.0), 0.3678794411714423215955); +/// # Ok::<(), ExpError>(()) /// ``` +#[doc(alias = "Exponential")] #[derive(Copy, Clone, PartialEq, Debug)] pub struct Exp { rate: f64, @@ -74,10 +76,11 @@ impl Exp { /// # Examples /// /// ``` - /// use statrs::distribution::Exp; + /// use statrs::distribution::{Exp, ExpError}; /// - /// let n = Exp::new(1.0).unwrap(); + /// let n = Exp::new(1.0)?; /// assert_eq!(n.rate(), 1.0); + /// # Ok::<(), ExpError>(()) /// ``` pub fn rate(&self) -> f64 { self.rate @@ -179,33 +182,71 @@ impl Max for Exp { } } -impl Distribution for Exp { - /// Returns the mean of the exponential distribution - /// - /// # Formula - /// - /// ```text - /// 1 / λ - /// ``` - /// - /// where `λ` is the rate - fn mean(&self) -> Option { - Some(1.0 / self.rate) +/// Returns the mean of the exponential distribution +/// +/// # Formula +/// +/// ```text +/// 1 / λ +/// ``` +/// +/// where `λ` is the rate + +impl Mean for Exp { + type Mu = f64; + fn mean(&self) -> Self::Mu { + self.rate.recip() } +} - /// Returns the variance of the exponential distribution - /// - /// # Formula - /// - /// ```text - /// 1 / λ^2 - /// ``` - /// - /// where `λ` is the rate - fn variance(&self) -> Option { - Some(1.0 / (self.rate * self.rate)) +/// Returns the variance of the exponential distribution +/// +/// # Formula +/// +/// ```text +/// 1 / λ^2 +/// ``` +/// +/// where `λ` is the rate + +impl Variance for Exp { + type Var = f64; + fn variance(&self) -> Self::Var { + self.rate.powi(-2) + } +} + +/// Returns the skewness of the exponential distribution +/// +/// # Formula +/// +/// ```text +/// 2 +/// ``` + +impl Skewness for Exp { + type Skew = f64; + fn skewness(&self) -> Self::Skew { + 2.0 + } +} + +/// Returns the excess kurtosis of the exponential distribution +/// +/// # Formula +/// +/// ```text +/// 6 +/// ``` + +impl ExcessKurtosis for Exp { + type Kurt = f64; + fn excess_kurtosis(&self) -> Self::Kurt { + 6.0 } +} +impl Entropy for Exp { /// Returns the entropy of the exponential distribution /// /// # Formula @@ -215,19 +256,8 @@ impl Distribution for Exp { /// ``` /// /// where `λ` is the rate - fn entropy(&self) -> Option { - Some(1.0 - self.rate.ln()) - } - - /// Returns the skewness of the exponential distribution - /// - /// # Formula - /// - /// ```text - /// 2 - /// ``` - fn skewness(&self) -> Option { - Some(2.0) + fn entropy(&self) -> f64 { + 1.0 - self.rate.ln() } } @@ -323,7 +353,7 @@ mod tests { #[test] fn test_mean() { - let mean = |x: Exp| x.mean().unwrap(); + let mean = |x: Exp| x.mean(); test_exact(0.1, 10.0, mean); test_exact(1.0, 1.0, mean); test_exact(10.0, 0.1, mean); @@ -331,7 +361,7 @@ mod tests { #[test] fn test_variance() { - let variance = |x: Exp| x.variance().unwrap(); + let variance = |x: Exp| x.variance(); test_absolute(0.1, 100.0, 1e-13, variance); test_exact(1.0, 1.0, variance); test_exact(10.0, 0.01, variance); @@ -339,7 +369,7 @@ mod tests { #[test] fn test_entropy() { - let entropy = |x: Exp| x.entropy().unwrap(); + let entropy = |x: Exp| x.entropy(); test_absolute(0.1, 3.302585092994045684018, 1e-15, entropy); test_exact(1.0, 1.0, entropy); test_absolute(10.0, -1.302585092994045684018, 1e-15, entropy); @@ -347,7 +377,7 @@ mod tests { #[test] fn test_skewness() { - let skewness = |x: Exp| x.skewness().unwrap(); + let skewness = |x: Exp| x.skewness(); test_exact(0.1, 2.0, skewness); test_exact(1.0, 2.0, skewness); test_exact(10.0, 2.0, skewness); diff --git a/src/distribution/fisher_snedecor.rs b/src/distribution/fisher_snedecor.rs index 4e204225..12d1c3ca 100644 --- a/src/distribution/fisher_snedecor.rs +++ b/src/distribution/fisher_snedecor.rs @@ -1,5 +1,5 @@ use crate::distribution::{Continuous, ContinuousCDF}; -use crate::function::beta; +use crate::function::{beta, gamma}; use crate::statistics::*; use std::f64; @@ -10,13 +10,14 @@ use std::f64; /// # Examples /// /// ``` -/// use statrs::distribution::{FisherSnedecor, Continuous}; -/// use statrs::statistics::Distribution; +/// use statrs::distribution::{FisherSnedecor, Continuous, FisherSnedecorError}; +/// use statrs::statistics::*; /// use statrs::prec; /// -/// let n = FisherSnedecor::new(3.0, 3.0).unwrap(); +/// let n = FisherSnedecor::new(3.0, 3.0)?; /// assert_eq!(n.mean().unwrap(), 3.0); /// assert!(prec::almost_eq(n.pdf(1.0), 0.318309886183790671538, 1e-15)); +/// # Ok::<(), FisherSnedecorError>(()) /// ``` #[derive(Debug, Copy, Clone, PartialEq)] pub struct FisherSnedecor { @@ -92,10 +93,11 @@ impl FisherSnedecor { /// # Examples /// /// ``` - /// use statrs::distribution::FisherSnedecor; + /// use statrs::distribution::{FisherSnedecor, FisherSnedecorError}; /// - /// let n = FisherSnedecor::new(2.0, 3.0).unwrap(); + /// let n = FisherSnedecor::new(2.0, 3.0)?; /// assert_eq!(n.freedom_1(), 2.0); + /// # Ok::<(), FisherSnedecorError>(()) /// ``` pub fn freedom_1(&self) -> f64 { self.freedom_1 @@ -107,10 +109,11 @@ impl FisherSnedecor { /// # Examples /// /// ``` - /// use statrs::distribution::FisherSnedecor; + /// use statrs::distribution::{FisherSnedecor, FisherSnedecorError}; /// /// let n = FisherSnedecor::new(2.0, 3.0).unwrap(); /// assert_eq!(n.freedom_2(), 3.0); + /// # Ok::<(), FisherSnedecorError>(()) /// ``` pub fn freedom_2(&self) -> f64 { self.freedom_2 @@ -239,51 +242,50 @@ impl Max for FisherSnedecor { } } -impl Distribution for FisherSnedecor { - /// Returns the mean of the fisher-snedecor distribution - /// - /// # Panics - /// - /// If `freedom_2 <= 2.0` - /// - /// # Remarks - /// - /// Returns `NaN` if `freedom_2` is `INF` - /// - /// # Formula - /// - /// ```text - /// d2 / (d2 - 2) - /// ``` - /// - /// where `d2` is the second degree of freedom - fn mean(&self) -> Option { +/// Returns the mean of the fisher-snedecor distribution as `Some` if defined +/// +/// # Remarks +/// +/// Returns `NaN` if `freedom_2` is `INF` +/// +/// # Formula +/// +/// ```text +/// d2 / (d2 - 2) +/// ``` +/// +/// where `d2` is the second degree of freedom + +impl Mean for FisherSnedecor { + type Mu = Option; + fn mean(&self) -> Self::Mu { if self.freedom_2 <= 2.0 { None } else { Some(self.freedom_2 / (self.freedom_2 - 2.0)) } } +} - /// Returns the variance of the fisher-snedecor distribution - /// - /// # Panics - /// - /// If `freedom_2 <= 4.0` - /// - /// # Remarks - /// - /// Returns `NaN` if `freedom_1` or `freedom_2` is `INF` - /// - /// # Formula - /// - /// ```text - /// (2 * d2^2 * (d1 + d2 - 2)) / (d1 * (d2 - 2)^2 * (d2 - 4)) - /// ``` - /// - /// where `d1` is the first degree of freedom and `d2` is - /// the second degree of freedom - fn variance(&self) -> Option { +/// Returns the variance of the fisher-snedecor distribution +/// +/// # Remarks +/// +/// Returns `NaN` if `freedom_1` or `freedom_2` is `INF` +/// Return `None` if `freedom_2` <= 4 +/// +/// # Formula +/// +/// ```text +/// (2 * d2^2 * (d1 + d2 - 2)) / (d1 * (d2 - 2)^2 * (d2 - 4)) +/// ``` +/// +/// where `d1` is the first degree of freedom and `d2` is +/// the second degree of freedom + +impl Variance for FisherSnedecor { + type Var = Option; + fn variance(&self) -> Self::Var { if self.freedom_2 <= 4.0 { None } else { @@ -296,27 +298,27 @@ impl Distribution for FisherSnedecor { Some(val) } } +} - /// Returns the skewness of the fisher-snedecor distribution - /// - /// # Panics - /// - /// If `freedom_2 <= 6.0` - /// - /// # Remarks - /// - /// Returns `NaN` if `freedom_1` or `freedom_2` is `INF` - /// - /// # Formula - /// - /// ```text - /// ((2d1 + d2 - 2) * sqrt(8 * (d2 - 4))) / ((d2 - 6) * sqrt(d1 * (d1 + d2 - /// - 2))) - /// ``` - /// - /// where `d1` is the first degree of freedom and `d2` is - /// the second degree of freedom - fn skewness(&self) -> Option { +/// Returns the skewness of the fisher-snedecor distribution +/// +/// # Remarks +/// +/// Returns `NaN` if `freedom_1` or `freedom_2` is `INF` +/// +/// # Formula +/// +/// ```text +/// ((2d1 + d2 - 2) * sqrt(8 * (d2 - 4))) / ((d2 - 6) * sqrt(d1 * (d1 + d2 +/// - 2))) +/// ``` +/// +/// where `d1` is the first degree of freedom and `d2` is +/// the second degree of freedom + +impl Skewness for FisherSnedecor { + type Skew = Option; + fn skewness(&self) -> Self::Skew { if self.freedom_2 <= 6.0 { None } else { @@ -329,6 +331,33 @@ impl Distribution for FisherSnedecor { } } +impl ExcessKurtosis for FisherSnedecor { + type Kurt = Option; + fn excess_kurtosis(&self) -> Self::Kurt { + if self.freedom_2 <= 8.0 { + None + } else { + let d1 = self.freedom_1; + let d2 = self.freedom_2; + let numer = (5. * d1 - 22.) * (d1 + d2 - 2.) + (d2 - 4.) * (d2 - 2.) * (d2 - 2.); + let denom = d1 * (d2 - 6.) * (d2 - 8.) * (d1 + d2 - 2.); + Some(12. * numer / denom) + } + } +} + +impl Entropy for FisherSnedecor { + fn entropy(&self) -> f64 { + let d1 = self.freedom_1; + let d2 = self.freedom_2; + gamma::ln_gamma(d1 / 2.) + gamma::ln_gamma(d2 / 2.) - gamma::ln_gamma((d1 + d2) / 2.) + + (1. - d1 / 2.) * gamma::digamma(1. + d1 / 2.) + - (1. + d2 / 2.) * gamma::digamma(1. + d2 / 2.) + + (d1 + d2) / 2. * gamma::digamma((d1 + d2) / 2.) + + (d2 / d1).ln() + } +} + impl Mode> for FisherSnedecor { /// Returns the mode for the fisher-snedecor distribution /// diff --git a/src/distribution/gamma.rs b/src/distribution/gamma.rs index a403355d..7542b1ec 100644 --- a/src/distribution/gamma.rs +++ b/src/distribution/gamma.rs @@ -9,13 +9,14 @@ use crate::statistics::*; /// # Examples /// /// ``` -/// use statrs::distribution::{Gamma, Continuous}; -/// use statrs::statistics::Distribution; +/// use statrs::distribution::{Gamma, Continuous, GammaError}; +/// use statrs::statistics::*; /// use statrs::prec; /// -/// let n = Gamma::new(3.0, 1.0).unwrap(); -/// assert_eq!(n.mean().unwrap(), 3.0); +/// let n = Gamma::new(3.0, 1.0)?; +/// assert_eq!(n.mean(), 3.0); /// assert!(prec::almost_eq(n.pdf(2.0), 0.270670566473225383788, 1e-15)); +/// # Ok::<(), GammaError>(()) /// ``` #[derive(Copy, Clone, PartialEq, Debug)] pub struct Gamma { @@ -91,10 +92,11 @@ impl Gamma { /// # Examples /// /// ``` - /// use statrs::distribution::Gamma; + /// use statrs::distribution::{Gamma, GammaError}; /// - /// let n = Gamma::new(3.0, 1.0).unwrap(); + /// let n = Gamma::new(3.0, 1.0)?; /// assert_eq!(n.shape(), 3.0); + /// # Ok::<(), GammaError>(()) /// ``` pub fn shape(&self) -> f64 { self.shape @@ -105,10 +107,11 @@ impl Gamma { /// # Examples /// /// ``` - /// use statrs::distribution::Gamma; + /// use statrs::distribution::{Gamma, GammaError}; /// /// let n = Gamma::new(3.0, 1.0).unwrap(); /// assert_eq!(n.rate(), 1.0); + /// # Ok::<(), GammaError>(()) /// ``` pub fn rate(&self) -> f64 { self.rate @@ -256,33 +259,63 @@ impl Max for Gamma { } } -impl Distribution for Gamma { - /// Returns the mean of the gamma distribution - /// - /// # Formula - /// - /// ```text - /// α / β - /// ``` - /// - /// where `α` is the shape and `β` is the rate - fn mean(&self) -> Option { - Some(self.shape / self.rate) +/// Returns the mean of the gamma distribution +/// +/// # Formula +/// +/// ```text +/// α / β +/// ``` +/// +/// where `α` is the shape and `β` is the rate +impl Mean for Gamma { + type Mu = f64; + fn mean(&self) -> Self::Mu { + self.shape / self.rate } +} - /// Returns the variance of the gamma distribution - /// - /// # Formula - /// - /// ```text - /// α / β^2 - /// ``` - /// - /// where `α` is the shape and `β` is the rate - fn variance(&self) -> Option { - Some(self.shape / (self.rate * self.rate)) +/// Returns the variance of the gamma distribution +/// +/// # Formula +/// +/// ```text +/// α / β^2 +/// ``` +/// +/// where `α` is the shape and `β` is the rate +impl Variance for Gamma { + type Var = f64; + fn variance(&self) -> Self::Var { + self.shape / (self.rate * self.rate) + } +} + +/// Returns the skewness of the gamma distribution +/// +/// # Formula +/// +/// ```text +/// 2 / sqrt(α) +/// ``` +/// +/// where `α` is the shape +impl Skewness for Gamma { + type Skew = f64; + fn skewness(&self) -> Self::Skew { + 2.0 / self.shape.sqrt() + } +} + +/// docs +impl ExcessKurtosis for Gamma { + type Kurt = (); + fn excess_kurtosis(&self) -> Self::Kurt { + unimplemented!() } +} +impl Entropy for Gamma { /// Returns the entropy of the gamma distribution /// /// # Formula @@ -293,24 +326,10 @@ impl Distribution for Gamma { /// /// where `α` is the shape, `β` is the rate, `Γ` is the gamma function, /// and `ψ` is the digamma function - fn entropy(&self) -> Option { - let entr = self.shape - self.rate.ln() + fn entropy(&self) -> f64 { + self.shape - self.rate.ln() + gamma::ln_gamma(self.shape) - + (1.0 - self.shape) * gamma::digamma(self.shape); - Some(entr) - } - - /// Returns the skewness of the gamma distribution - /// - /// # Formula - /// - /// ```text - /// 2 / sqrt(α) - /// ``` - /// - /// where `α` is the shape - fn skewness(&self) -> Option { - Some(2.0 / self.shape.sqrt()) + + (1.0 - self.shape) * gamma::digamma(self.shape) } } @@ -478,7 +497,7 @@ mod tests { #[test] fn test_mean() { - let f = |x: Gamma| x.mean().unwrap(); + let f = |x: Gamma| x.mean(); let test = [ ((1.0, 0.1), 10.0), ((1.0, 1.0), 1.0), @@ -493,7 +512,7 @@ mod tests { #[test] fn test_variance() { - let f = |x: Gamma| x.variance().unwrap(); + let f = |x: Gamma| x.variance(); let test = [ ((1.0, 0.1), 100.0), ((1.0, 1.0), 1.0), @@ -508,7 +527,7 @@ mod tests { #[test] fn test_entropy() { - let f = |x: Gamma| x.entropy().unwrap(); + let f = |x: Gamma| x.entropy(); let test = [ ((1.0, 0.1), 3.302585092994045628506840223), ((1.0, 1.0), 1.0), @@ -523,7 +542,7 @@ mod tests { #[test] fn test_skewness() { - let f = |x: Gamma| x.skewness().unwrap(); + let f = |x: Gamma| x.skewness(); let test = [ ((1.0, 0.1), 2.0), ((1.0, 1.0), 2.0), diff --git a/src/distribution/geometric.rs b/src/distribution/geometric.rs index 4ec1b30b..e1477101 100644 --- a/src/distribution/geometric.rs +++ b/src/distribution/geometric.rs @@ -9,13 +9,14 @@ use std::f64; /// # Examples /// /// ``` -/// use statrs::distribution::{Geometric, Discrete}; -/// use statrs::statistics::Distribution; +/// use statrs::distribution::{Geometric, GeometricError, Discrete}; +/// use statrs::statistics::*; /// -/// let n = Geometric::new(0.3).unwrap(); -/// assert_eq!(n.mean().unwrap(), 1.0 / 0.3); +/// let n = Geometric::new(0.3)?; +/// assert_eq!(n.mean(), 1.0 / 0.3); /// assert_eq!(n.pmf(1), 0.3); /// assert_eq!(n.pmf(2), 0.21); +/// # Ok::<(), GeometricError>(()) /// ``` #[derive(Copy, Clone, PartialEq, Debug)] pub struct Geometric { @@ -75,9 +76,11 @@ impl Geometric { /// /// ``` /// use statrs::distribution::Geometric; + /// # use statrs::distribution::GeometricError; /// - /// let n = Geometric::new(0.5).unwrap(); + /// let n = Geometric::new(0.5)?; /// assert_eq!(n.p(), 0.5); + /// # Ok::<(), GeometricError>(()) /// ``` pub fn p(&self) -> f64 { self.p @@ -183,53 +186,74 @@ impl Max for Geometric { } } -impl Distribution for Geometric { - /// Returns the mean of the geometric distribution - /// - /// # Formula - /// - /// ```text - /// 1 / p - /// ``` - fn mean(&self) -> Option { - Some(1.0 / self.p) +/// Returns the mean of the geometric distribution +/// +/// # Formula +/// +/// ```text +/// 1 / p +/// ``` + +impl Mean for Geometric { + type Mu = f64; + fn mean(&self) -> Self::Mu { + self.p.recip() } +} - /// Returns the standard deviation of the geometric distribution - /// - /// # Formula - /// - /// ```text - /// (1 - p) / p^2 - /// ``` - fn variance(&self) -> Option { - Some((1.0 - self.p) / (self.p * self.p)) +/// Returns the standard deviation of the geometric distribution +/// +/// # Formula +/// +/// ```text +/// (1 - p) / p^2 +/// ``` + +impl Variance for Geometric { + type Var = f64; + fn variance(&self) -> Self::Var { + (1.0 - self.p) / (self.p * self.p) } +} - /// Returns the entropy of the geometric distribution - /// - /// # Formula - /// - /// ```text - /// (-(1 - p) * log_2(1 - p) - p * log_2(p)) / p - /// ``` - fn entropy(&self) -> Option { - let inv = 1.0 / self.p; - Some(-inv * (1. - self.p).log(2.0) + (inv - 1.).log(2.0)) +/// Returns the skewness of the geometric distribution +/// +/// # Formula +/// +/// ```text +/// (2 - p) / sqrt(1 - p) +/// ``` + +impl Skewness for Geometric { + type Skew = f64; + fn skewness(&self) -> Self::Skew { + if ulps_eq!(self.p, 1.0) { + return f64::INFINITY; + }; + (2.0 - self.p) / (1.0 - self.p).sqrt() + } +} + +/// Returns the excess kurtosis of the geometric distribution + +impl ExcessKurtosis for Geometric { + type Kurt = f64; + fn excess_kurtosis(&self) -> Self::Kurt { + 6. + self.p.powi(2) / (1. - self.p) } +} - /// Returns the skewness of the geometric distribution +impl Entropy for Geometric { + /// Returns the entropy of the geometric distribution /// /// # Formula /// /// ```text - /// (2 - p) / sqrt(1 - p) + /// (-(1 - p) * log_2(1 - p) - p * log_2(p)) / p /// ``` - fn skewness(&self) -> Option { - if ulps_eq!(self.p, 1.0) { - return Some(f64::INFINITY); - }; - Some((2.0 - self.p) / (1.0 - self.p).sqrt()) + fn entropy(&self) -> f64 { + let inv = self.p.recip(); + -inv * (1. - self.p).log(2.0) + (inv - 1.).log(2.0) } } @@ -324,28 +348,28 @@ mod tests { #[test] fn test_mean() { - let mean = |x: Geometric| x.mean().unwrap(); + let mean = |x: Geometric| x.mean(); test_exact(0.3, 1.0 / 0.3, mean); test_exact(1.0, 1.0, mean); } #[test] fn test_variance() { - let variance = |x: Geometric| x.variance().unwrap(); + let variance = |x: Geometric| x.variance(); test_exact(0.3, 0.7 / (0.3 * 0.3), variance); test_exact(1.0, 0.0, variance); } #[test] fn test_entropy() { - let entropy = |x: Geometric| x.entropy().unwrap(); + let entropy = |x: Geometric| x.entropy(); test_absolute(0.3, 2.937636330768973333333, 1e-14, entropy); test_is_nan(1.0, entropy); } #[test] fn test_skewness() { - let skewness = |x: Geometric| x.skewness().unwrap(); + let skewness = |x: Geometric| x.skewness(); test_absolute(0.3, 2.031888635868469187947, 1e-15, skewness); test_exact(1.0, f64::INFINITY, skewness); } diff --git a/src/distribution/gumbel.rs b/src/distribution/gumbel.rs index ddd17e4e..42a8dc45 100644 --- a/src/distribution/gumbel.rs +++ b/src/distribution/gumbel.rs @@ -1,5 +1,5 @@ use super::{Continuous, ContinuousCDF}; -use crate::consts::EULER_MASCHERONI; +use crate::consts::{self, EULER_MASCHERONI}; use crate::statistics::*; use std::f64::{self, consts::PI}; @@ -9,13 +9,14 @@ use std::f64::{self, consts::PI}; /// # Examples /// /// ``` -/// use statrs::distribution::{Gumbel, Continuous}; -/// use statrs::{consts::EULER_MASCHERONI, statistics::Distribution}; +/// use statrs::distribution::{Gumbel, Continuous, GumbelError}; +/// use statrs::{consts::EULER_MASCHERONI, statistics::*}; /// -/// let n = Gumbel::new(0.0, 1.0).unwrap(); +/// let n = Gumbel::new(0.0, 1.0)?; /// assert_eq!(n.location(), 0.0); -/// assert_eq!(n.skewness().unwrap(), 1.13955); +/// assert_eq!(n.mean(), 0.5772156649015329); /// assert_eq!(n.pdf(0.0), 0.36787944117144233); +/// # Ok::<(), GumbelError>(()) /// ``` #[derive(Copy, Clone, PartialEq, Debug)] pub struct Gumbel { @@ -82,10 +83,11 @@ impl Gumbel { /// # Examples /// /// ``` - /// use statrs::distribution::Gumbel; + /// use statrs::distribution::{Gumbel, GumbelError}; /// - /// let n = Gumbel::new(0.0, 1.0).unwrap(); + /// let n = Gumbel::new(0.0, 1.0)?; /// assert_eq!(n.location(), 0.0); + /// # Ok::<(), GumbelError>(()) /// ``` pub fn location(&self) -> f64 { self.location @@ -96,14 +98,19 @@ impl Gumbel { /// # Examples /// /// ``` - /// use statrs::distribution::Gumbel; + /// use statrs::distribution::{Gumbel, GumbelError}; /// - /// let n = Gumbel::new(0.0, 1.0).unwrap(); + /// let n = Gumbel::new(0.0, 1.0)?; /// assert_eq!(n.scale(), 1.0); + /// # Ok::<(), GumbelError>(()) /// ``` pub fn scale(&self) -> f64 { self.scale } + + pub fn std_dev(&self) -> f64 { + (PI * self.scale) / 6_f64.sqrt() + } } impl std::fmt::Display for Gumbel { @@ -199,74 +206,88 @@ impl Max for Gumbel { } } -impl Distribution for Gumbel { - /// Returns the entropy of the Gumbel distribution - /// - /// # Formula - /// - /// ```text - /// ln(β) + γ + 1 - /// ``` - /// - /// where `β` is the scale - /// and `γ` is the Euler-Mascheroni constant (approx 0.57721) - fn entropy(&self) -> Option { - Some(1.0 + EULER_MASCHERONI + (self.scale).ln()) +/// Returns the mean of the Gumbel distribution +/// +/// # Formula +/// +/// ```text +/// μ + γβ +/// ``` +/// +/// where `μ` is the location, `β` is the scale +/// and `γ` is the Euler-Mascheroni constant (approx 0.57721) + +impl Mean for Gumbel { + type Mu = f64; + fn mean(&self) -> Self::Mu { + self.location + (EULER_MASCHERONI * self.scale) } +} - /// Returns the mean of the Gumbel distribution - /// - /// # Formula - /// - /// ```text - /// μ + γβ - /// ``` - /// - /// where `μ` is the location, `β` is the scale - /// and `γ` is the Euler-Mascheroni constant (approx 0.57721) - fn mean(&self) -> Option { - Some(self.location + (EULER_MASCHERONI * self.scale)) +/// Returns the variance of the Gumbel distribution +/// +/// # Formula +/// +/// ```text +/// (π^2 / 6) * β^2 +/// ``` +/// +/// where `β` is the scale and `π` is the constant PI (approx 3.14159) + +impl Variance for Gumbel { + type Var = f64; + fn variance(&self) -> Self::Var { + ((PI * PI) / 6.0) * self.scale * self.scale } +} - /// Returns the skewness of the Gumbel distribution - /// - /// # Formula - /// - /// ```text - /// 12 * sqrt(6) * ζ(3) / π^3 ≈ 1.13955 - /// ``` - /// ζ(3) is the Riemann zeta function evaluated at 3 (approx 1.20206) - /// and π is the constant PI (approx 3.14159) - /// - /// This approximately evaluates to 1.13955 - fn skewness(&self) -> Option { - Some(1.13955) +/// Returns the skewness of the Gumbel distribution +/// +/// # Formula +/// +/// ```text +/// 12 * sqrt(6) * ζ(3) / π^3 ≈ 1.13955 +/// ``` +/// ζ(3) is the Riemann zeta function evaluated at 3 (approx 1.20206) +/// and π is the constant PI (approx 3.14159) +/// +/// This approximately evaluates to 1.13955 + +impl Skewness for Gumbel { + type Skew = f64; + fn skewness(&self) -> Self::Skew { + 12. * 6_f64.sqrt() * consts::ZETA_3 * PI.powi(-3) } +} - /// Returns the variance of the Gumbel distribution - /// - /// # Formula - /// - /// ```text - /// (π^2 / 6) * β^2 - /// ``` - /// - /// where `β` is the scale and `π` is the constant PI (approx 3.14159) - fn variance(&self) -> Option { - Some(((PI * PI) / 6.0) * self.scale * self.scale) +/// Return the excess kurtosis of the Gumbel Distribution +/// +/// # Formula +/// +/// ```text +/// 12 / 5 +/// ``` + +impl ExcessKurtosis for Gumbel { + type Kurt = f64; + fn excess_kurtosis(&self) -> Self::Kurt { + 2.4 } +} - /// Returns the standard deviation of the Gumbel distribution +impl Entropy for Gumbel { + /// Returns the entropy of the Gumbel distribution /// /// # Formula /// /// ```text - /// β * π / sqrt(6) + /// ln(β) + γ + 1 /// ``` /// - /// where `β` is the scale and `π` is the constant PI (approx 3.14159) - fn std_dev(&self) -> Option { - Some(self.scale * PI / 6.0_f64.sqrt()) + /// where `β` is the scale + /// and `γ` is the Euler-Mascheroni constant (approx 0.57721) + fn entropy(&self) -> f64 { + 1.0 + EULER_MASCHERONI + (self.scale).ln() } } @@ -379,7 +400,7 @@ mod tests { #[test] fn test_entropy() { - let entropy = |x: Gumbel| x.entropy().unwrap(); + let entropy = |x: Gumbel| x.entropy(); test_exact(0.0, 2.0, 2.270362845461478, entropy); test_exact(0.1, 4.0, 2.9635100260214235, entropy); test_exact(1.0, 10.0, 3.8798007578955787, entropy); @@ -388,7 +409,7 @@ mod tests { #[test] fn test_mean() { - let mean = |x: Gumbel| x.mean().unwrap(); + let mean = |x: Gumbel| x.mean(); test_exact(0.0, 2.0, 1.1544313298030658, mean); test_exact(0.1, 4.0, 2.4088626596061316, mean); test_exact(1.0, 10.0, 6.772156649015328, mean); @@ -398,17 +419,19 @@ mod tests { #[test] fn test_skewness() { - let skewness = |x: Gumbel| x.skewness().unwrap(); - test_exact(0.0, 2.0, 1.13955, skewness); - test_exact(0.1, 4.0, 1.13955, skewness); - test_exact(1.0, 10.0, 1.13955, skewness); - test_exact(10.0, 11.0, 1.13955, skewness); - test_exact(10.0, f64::INFINITY, 1.13955, skewness); + let skewness = |x: Gumbel| x.skewness(); + let value = skewness(Gumbel{location: 0.0, scale: 1.0}); + + test_exact(0.0, 2.0, value, skewness); + test_exact(0.1, 4.0, value, skewness); + test_exact(1.0, 10.0, value, skewness); + test_exact(10.0, 11.0, value, skewness); + test_exact(10.0, f64::INFINITY, value, skewness); } #[test] fn test_variance() { - let variance = |x: Gumbel| x.variance().unwrap(); + let variance = |x: Gumbel| x.variance(); test_exact(0.0, 2.0, 6.579736267392906, variance); test_exact(0.1, 4.0, 26.318945069571624, variance); test_exact(1.0, 10.0, 164.49340668482265, variance); @@ -417,7 +440,7 @@ mod tests { #[test] fn test_std_dev() { - let std_dev = |x: Gumbel| x.std_dev().unwrap(); + let std_dev = |x: Gumbel| x.std_dev(); test_exact(0.0, 2.0, 2.565099660323728, std_dev); test_exact(0.1, 4.0, 5.130199320647456, std_dev); test_exact(1.0, 10.0, 12.82549830161864, std_dev); diff --git a/src/distribution/hypergeometric.rs b/src/distribution/hypergeometric.rs index 4c03bebe..409ee548 100644 --- a/src/distribution/hypergeometric.rs +++ b/src/distribution/hypergeometric.rs @@ -12,13 +12,13 @@ use std::f64; /// /// ``` /// use statrs::distribution::{Hypergeometric, Discrete}; -/// use statrs::statistics::Distribution; -/// use statrs::prec; +/// use statrs::statistics::*; +/// use approx::assert_relative_eq; /// /// let n = Hypergeometric::new(500, 50, 100).unwrap(); /// assert_eq!(n.mean().unwrap(), 10.); -/// assert!(prec::almost_eq(n.pmf(10), 0.14736784, 1e-8)); -/// assert!(prec::almost_eq(n.pmf(25), 3.537e-7, 1e-10)); +/// assert_relative_eq!(n.pmf(10), 0.14736784, epsilon=1e-8); +/// assert_relative_eq!(n.pmf(25), 3.537e-7, epsilon=1e-10); /// ``` #[derive(Copy, Clone, PartialEq, Eq, Debug)] pub struct Hypergeometric { @@ -292,42 +292,48 @@ impl Max for Hypergeometric { } } -impl Distribution for Hypergeometric { - /// Returns the mean of the hypergeometric distribution - /// - /// # None - /// - /// If `N` is `0` - /// - /// # Formula - /// - /// ```text - /// K * n / N - /// ``` - /// - /// where `N` is population, `K` is successes, and `n` is draws - fn mean(&self) -> Option { +/// Returns the mean of the hypergeometric distribution +/// +/// # None +/// +/// If `N` is `0` +/// +/// # Formula +/// +/// ```text +/// K * n / N +/// ``` +/// +/// where `N` is population, `K` is successes, and `n` is draws + +impl Mean for Hypergeometric { + type Mu = Option; + fn mean(&self) -> Self::Mu { if self.population == 0 { None } else { Some(self.successes as f64 * self.draws as f64 / self.population as f64) } } +} - /// Returns the variance of the hypergeometric distribution - /// - /// # None - /// - /// If `N <= 1` - /// - /// # Formula - /// - /// ```text - /// n * (K / N) * ((N - K) / N) * ((N - n) / (N - 1)) - /// ``` - /// - /// where `N` is population, `K` is successes, and `n` is draws - fn variance(&self) -> Option { +/// Returns the variance of the hypergeometric distribution +/// +/// # None +/// +/// If `N <= 1` +/// +/// # Formula +/// +/// ```text +/// n * (K / N) * ((N - K) / N) * ((N - n) / (N - 1)) +/// ``` +/// +/// where `N` is population, `K` is successes, and `n` is draws + +impl Variance for Hypergeometric { + type Var = Option; + fn variance(&self) -> Self::Var { if self.population <= 1 { None } else { @@ -337,22 +343,26 @@ impl Distribution for Hypergeometric { Some(val) } } +} - /// Returns the skewness of the hypergeometric distribution - /// - /// # None - /// - /// If `N <= 2` - /// - /// # Formula - /// - /// ```text - /// ((N - 2K) * (N - 1)^(1 / 2) * (N - 2n)) / ([n * K * (N - K) * (N - - /// n)]^(1 / 2) * (N - 2)) - /// ``` - /// - /// where `N` is population, `K` is successes, and `n` is draws - fn skewness(&self) -> Option { +/// Returns the skewness of the hypergeometric distribution +/// +/// # None +/// +/// If `N <= 2` +/// +/// # Formula +/// +/// ```text +/// ((N - 2K) * (N - 1)^(1 / 2) * (N - 2n)) / ([n * K * (N - K) * (N - +/// n)]^(1 / 2) * (N - 2)) +/// ``` +/// +/// where `N` is population, `K` is successes, and `n` is draws + +impl Skewness for Hypergeometric { + type Skew = Option; + fn skewness(&self) -> Self::Skew { if self.population <= 2 { None } else { @@ -367,6 +377,41 @@ impl Distribution for Hypergeometric { } } +/// Returns the excess kurtosis of the hypergeometric distribution +/// +/// # Formula +/// if N >= 3 +/// ```text +/// $$ +/// \frac{1}{n K (N - K)(N-n)(N-2)(N-3)}\left[N^2(N-1)\left(N(N+1) - 6K(N-K) - 6n(N-n)\right) + 6 n K (N-K)(N-n)(5N-6)\right] +/// $$ +/// ``` +/// `None` otherwise + +impl ExcessKurtosis for Hypergeometric { + type Kurt = Option; + fn excess_kurtosis(&self) -> Self::Kurt { + if self.population < 3 { + return None; + } + let pop = self.population as f64; + let succ = self.successes as f64; + let n = self.draws as f64; + + let pop_sqr = (self.population * self.population) as f64; + let pop_m2 = (self.population - 2) as f64; + let pop_m3 = (self.population - 3) as f64; + let pop_mn = (self.population - self.draws) as f64; + let pop_msucc = (self.population - self.successes) as f64; + + let kurt_a_inv = n * succ * pop_msucc * pop_mn * pop_m2 * pop_m3; + let kurt_b = + pop * (pop_sqr - pop) * (pop_sqr + pop - 6. * succ * pop_msucc - 6. * n * pop_mn) + + 6. * n * succ * pop_msucc * pop_mn * (5. * pop - 6.); + Some(kurt_b / kurt_a_inv) + } +} + impl Mode> for Hypergeometric { /// Returns the mode of the hypergeometric distribution /// diff --git a/src/distribution/internal.rs b/src/distribution/internal.rs index 9075669e..b4fc13e9 100644 --- a/src/distribution/internal.rs +++ b/src/distribution/internal.rs @@ -335,7 +335,7 @@ pub mod test { #[test] #[should_panic] fn test_is_nan_failure() { - test_is_nan(0.8, 1.2, |dist| dist.mean().unwrap()); + test_is_nan(0.8, 1.2, |dist| dist.mean()); } #[test] @@ -346,7 +346,7 @@ pub mod test { #[test] #[should_panic] fn test_is_none_failure() { - test_none(0.8, 1.2, |dist| dist.mean()); + test_none(0.8, 1.2, |dist| Some(dist.mean())); } } diff --git a/src/distribution/inverse_gamma.rs b/src/distribution/inverse_gamma.rs index c34f99fa..0cb40da7 100644 --- a/src/distribution/inverse_gamma.rs +++ b/src/distribution/inverse_gamma.rs @@ -11,11 +11,11 @@ use std::f64; /// /// ``` /// use statrs::distribution::{InverseGamma, Continuous}; -/// use statrs::statistics::Distribution; -/// use statrs::prec; +/// use statrs::statistics::*; +/// use approx::assert_relative_eq; /// /// let n = InverseGamma::new(1.1, 0.1).unwrap(); -/// assert!(prec::almost_eq(n.mean().unwrap(), 1.0, 1e-14)); +/// assert_relative_eq!(n.mean().unwrap(), 1.0, epsilon=1e-14); /// assert_eq!(n.pdf(1.0), 0.07554920138253064); /// ``` #[derive(Copy, Clone, PartialEq, Debug)] @@ -202,42 +202,48 @@ impl Max for InverseGamma { } } -impl Distribution for InverseGamma { - /// Returns the mean of the inverse distribution - /// - /// # None - /// - /// If `shape <= 1.0` - /// - /// # Formula - /// - /// ```text - /// β / (α - 1) - /// ``` - /// - /// where `α` is the shape and `β` is the rate - fn mean(&self) -> Option { +/// Returns the mean of the inverse distribution +/// +/// # None +/// +/// If `shape <= 1.0` +/// +/// # Formula +/// +/// ```text +/// β / (α - 1) +/// ``` +/// +/// where `α` is the shape and `β` is the rate + +impl Mean for InverseGamma { + type Mu = Option; + fn mean(&self) -> Self::Mu { if self.shape <= 1.0 { None } else { Some(self.rate / (self.shape - 1.0)) } } +} - /// Returns the variance of the inverse gamma distribution - /// - /// # None - /// - /// If `shape <= 2.0` - /// - /// # Formula - /// - /// ```text - /// β^2 / ((α - 1)^2 * (α - 2)) - /// ``` - /// - /// where `α` is the shape and `β` is the rate - fn variance(&self) -> Option { +/// Returns the variance of the inverse gamma distribution +/// +/// # None +/// +/// If `shape <= 2.0` +/// +/// # Formula +/// +/// ```text +/// β^2 / ((α - 1)^2 * (α - 2)) +/// ``` +/// +/// where `α` is the shape and `β` is the rate + +impl Variance for InverseGamma { + type Var = Option; + fn variance(&self) -> Self::Var { if self.shape <= 2.0 { None } else { @@ -246,7 +252,59 @@ impl Distribution for InverseGamma { Some(val) } } +} +/// Returns the skewness of the inverse gamma distribution +/// +/// # None +/// +/// If `shape <= 3` +/// +/// # Formula +/// +/// ```text +/// 4 * sqrt(α - 2) / (α - 3) +/// ``` +/// +/// where `α` is the shape + +impl Skewness for InverseGamma { + type Skew = Option; + fn skewness(&self) -> Self::Skew { + if self.shape <= 3.0 { + None + } else { + Some(4.0 * (self.shape - 2.0).sqrt() / (self.shape - 3.0)) + } + } +} + +/// Returns the excess kurtosis of the inverse gamma distribution +/// +/// # None +/// +/// If `shape <= 3` +/// +/// # Formula +/// +/// ```text +/// 6 (5α - 11) / (α - 3) / (α - 4) +/// ``` +/// +/// where `α` is the shape + +impl ExcessKurtosis for InverseGamma { + type Kurt = Option; + fn excess_kurtosis(&self) -> Self::Kurt { + if self.shape <= 4. { + None + } else { + Some(6. * (5. * self.shape - 11.) / (self.shape - 3.) / (self.shape - 4.)) + } + } +} + +impl Entropy for InverseGamma { /// Returns the entropy of the inverse gamma distribution /// /// # Formula @@ -257,31 +315,9 @@ impl Distribution for InverseGamma { /// /// where `α` is the shape, `β` is the rate, `Γ` is the gamma function, /// and `ψ` is the digamma function - fn entropy(&self) -> Option { - let entr = self.shape + self.rate.ln() + gamma::ln_gamma(self.shape) - - (1.0 + self.shape) * gamma::digamma(self.shape); - Some(entr) - } - - /// Returns the skewness of the inverse gamma distribution - /// - /// # None - /// - /// If `shape <= 3` - /// - /// # Formula - /// - /// ```text - /// 4 * sqrt(α - 2) / (α - 3) - /// ``` - /// - /// where `α` is the shape - fn skewness(&self) -> Option { - if self.shape <= 3.0 { - None - } else { - Some(4.0 * (self.shape - 2.0).sqrt() / (self.shape - 3.0)) - } + fn entropy(&self) -> f64 { + self.shape + self.rate.ln() + gamma::ln_gamma(self.shape) + - (1.0 + self.shape) * gamma::digamma(self.shape) } } @@ -395,7 +431,7 @@ mod tests { #[test] fn test_entropy() { - let entropy = |x: InverseGamma| x.entropy().unwrap(); + let entropy = |x: InverseGamma| x.entropy(); test_absolute(0.1, 0.1, 11.51625799319234475054, 1e-14, entropy); test_absolute(1.0, 1.0, 2.154431329803065721213, 1e-14, entropy); } diff --git a/src/distribution/laplace.rs b/src/distribution/laplace.rs index f0d5b8ea..13f4e7b6 100644 --- a/src/distribution/laplace.rs +++ b/src/distribution/laplace.rs @@ -1,5 +1,5 @@ use crate::distribution::{Continuous, ContinuousCDF}; -use crate::statistics::{Distribution, Max, Median, Min, Mode}; +use crate::statistics::*; use std::f64; /// Implements the [Laplace](https://en.wikipedia.org/wiki/Laplace_distribution) @@ -213,33 +213,55 @@ impl Max for Laplace { } } -impl Distribution for Laplace { - /// Returns the mode of the laplace distribution - /// - /// # Formula - /// - /// ```text - /// μ - /// ``` - /// - /// where `μ` is the location - fn mean(&self) -> Option { - Some(self.location) +/// Returns the mode of the laplace distribution +/// +/// # Formula +/// +/// ```text +/// μ +/// ``` +/// +/// where `μ` is the location + +impl Mean for Laplace { + type Mu = f64; + fn mean(&self) -> Self::Mu { + self.location } +} - /// Returns the variance of the laplace distribution - /// - /// # Formula - /// - /// ```text - /// 2*b^2 - /// ``` - /// - /// where `b` is the scale - fn variance(&self) -> Option { - Some(2. * self.scale * self.scale) +/// Returns the variance of the laplace distribution +/// +/// # Formula +/// +/// ```text +/// 2*b^2 +/// ``` +/// +/// where `b` is the scale + +impl Variance for Laplace { + type Var = f64; + fn variance(&self) -> Self::Var { + 2. * self.scale * self.scale } +} + +impl Skewness for Laplace { + type Skew = f64; + fn skewness(&self) -> Self::Skew { + 0. + } +} +impl ExcessKurtosis for Laplace { + type Kurt = f64; + fn excess_kurtosis(&self) -> Self::Kurt { + 3.0 + } +} + +impl Entropy for Laplace { /// Returns the entropy of the laplace distribution /// /// # Formula @@ -249,19 +271,8 @@ impl Distribution for Laplace { /// ``` /// /// where `b` is the scale - fn entropy(&self) -> Option { - Some((2. * self.scale).ln() + 1.) - } - - /// Returns the skewness of the laplace distribution - /// - /// # Formula - /// - /// ```text - /// 0 - /// ``` - fn skewness(&self) -> Option { - Some(0.) + fn entropy(&self) -> f64 { + (2. * self.scale).ln() + 1. } } @@ -365,7 +376,7 @@ mod tests { #[test] fn test_mean() { - let mean = |x: Laplace| x.mean().unwrap(); + let mean = |x: Laplace| x.mean(); test_exact(f64::NEG_INFINITY, 0.1, f64::NEG_INFINITY, mean); test_exact(-5.0 - 1.0, 1.0, -6.0, mean); test_exact(0.0, 5.0, 0.0, mean); @@ -375,7 +386,7 @@ mod tests { #[test] fn test_variance() { - let variance = |x: Laplace| x.variance().unwrap(); + let variance = |x: Laplace| x.variance(); test_absolute(f64::NEG_INFINITY, 0.1, 0.02, 1E-12, variance); test_absolute(-5.0 - 1.0, 1.0, 2.0, 1E-12, variance); test_absolute(0.0, 5.0, 50.0, 1E-12, variance); @@ -385,7 +396,7 @@ mod tests { } #[test] fn test_entropy() { - let entropy = |x: Laplace| x.entropy().unwrap(); + let entropy = |x: Laplace| x.entropy(); test_absolute( f64::NEG_INFINITY, 0.1, @@ -401,7 +412,7 @@ mod tests { #[test] fn test_skewness() { - let skewness = |x: Laplace| x.skewness().unwrap(); + let skewness = |x: Laplace| x.skewness(); test_exact(f64::NEG_INFINITY, 0.1, 0.0, skewness); test_exact(-6.0, 1.0, 0.0, skewness); test_exact(1.0, 7.0, 0.0, skewness); @@ -554,8 +565,8 @@ mod tests { #[cfg(feature = "rand")] #[test] fn test_sample() { - use ::rand::distributions::Distribution; - use ::rand::thread_rng; + use rand::distributions::Distribution; + use rand::thread_rng; let l = create_ok(0.1, 0.5); l.sample(&mut thread_rng()); @@ -564,9 +575,9 @@ mod tests { #[cfg(feature = "rand")] #[test] fn test_sample_distribution() { - use ::rand::distributions::Distribution; - use ::rand::rngs::StdRng; - use ::rand::SeedableRng; + use rand::distributions::Distribution; + use rand::rngs::StdRng; + use rand::SeedableRng; // sanity check sampling let location = 0.0; diff --git a/src/distribution/log_normal.rs b/src/distribution/log_normal.rs index 4ddbb43e..2bcb5960 100644 --- a/src/distribution/log_normal.rs +++ b/src/distribution/log_normal.rs @@ -12,12 +12,12 @@ use std::f64; /// /// ``` /// use statrs::distribution::{LogNormal, Continuous}; -/// use statrs::statistics::Distribution; -/// use statrs::prec; +/// use statrs::statistics::*; +/// use approx::assert_relative_eq; /// /// let n = LogNormal::new(0.0, 1.0).unwrap(); -/// assert_eq!(n.mean().unwrap(), (0.5f64).exp()); -/// assert!(prec::almost_eq(n.pdf(1.0), 0.3989422804014326779399, 1e-16)); +/// assert_eq!(n.mean(), (0.5f64).exp()); +/// assert_relative_eq!(n.pdf(1.0), 0.3989422804014326779399, epsilon=1e-16); /// ``` #[derive(Copy, Clone, PartialEq, Debug)] pub struct LogNormal { @@ -204,59 +204,89 @@ impl Max for LogNormal { } } -impl Distribution for LogNormal { - /// Returns the mean of the log-normal distribution - /// - /// # Formula - /// - /// ```text - /// e^(μ + σ^2 / 2) - /// ``` - /// - /// where `μ` is the location and `σ` is the scale - fn mean(&self) -> Option { - Some((self.location + self.scale * self.scale / 2.0).exp()) +/// Returns the mean of the log-normal distribution +/// +/// # Formula +/// +/// ```text +/// e^(μ + σ^2 / 2) +/// ``` +/// +/// where `μ` is the location and `σ` is the scale + +impl Mean for LogNormal { + type Mu = f64; + fn mean(&self) -> Self::Mu { + (self.location + self.scale * self.scale / 2.0).exp() } +} - /// Returns the variance of the log-normal distribution - /// - /// # Formula - /// - /// ```text - /// (e^(σ^2) - 1) * e^(2μ + σ^2) - /// ``` - /// - /// where `μ` is the location and `σ` is the scale - fn variance(&self) -> Option { +/// Returns the variance of the log-normal distribution +/// +/// # Formula +/// +/// ```text +/// (e^(σ^2) - 1) * e^(2μ + σ^2) +/// ``` +/// +/// where `μ` is the location and `σ` is the scale + +impl Variance for LogNormal { + type Var = f64; + fn variance(&self) -> Self::Var { let sigma2 = self.scale * self.scale; - Some((sigma2.exp() - 1.0) * (self.location + self.location + sigma2).exp()) + (sigma2.exp() - 1.0) * (self.location + self.location + sigma2).exp() } +} - /// Returns the entropy of the log-normal distribution - /// - /// # Formula - /// - /// ```text - /// ln(σe^(μ + 1 / 2) * sqrt(2π)) - /// ``` - /// - /// where `μ` is the location and `σ` is the scale - fn entropy(&self) -> Option { - Some(0.5 + self.scale.ln() + self.location + consts::LN_SQRT_2PI) +/// Returns the skewness of the log-normal distribution +/// +/// # Formula +/// +/// ```text +/// (e^(σ^2) + 2) * sqrt(e^(σ^2) - 1) +/// ``` +/// +/// where `μ` is the location and `σ` is the scale + +impl Skewness for LogNormal { + type Skew = f64; + fn skewness(&self) -> Self::Skew { + let expsigma2 = (self.scale * self.scale).exp(); + (expsigma2 + 2.0) * (expsigma2 - 1.0).sqrt() } +} - /// Returns the skewness of the log-normal distribution +/// Returns the excess kurtosis of the log-normal distribution +/// +/// # Formula +/// +/// ```text +/// +/// ``` +/// +/// where `μ` is the location and `σ` is the scale + +impl ExcessKurtosis for LogNormal { + type Kurt = f64; + fn excess_kurtosis(&self) -> Self::Kurt { + let expsigma2 = (self.scale * self.scale).exp(); + expsigma2.powi(4) + 2. * expsigma2.powi(3) + 3. * expsigma2.powi(2) - 6. + } +} + +impl Entropy for LogNormal { + /// Returns the entropy of the log-normal distribution /// /// # Formula /// /// ```text - /// (e^(σ^2) + 2) * sqrt(e^(σ^2) - 1) + /// ln(σe^(μ + 1 / 2) * sqrt(2π)) /// ``` /// /// where `μ` is the location and `σ` is the scale - fn skewness(&self) -> Option { - let expsigma2 = (self.scale * self.scale).exp(); - Some((expsigma2 + 2.0) * (expsigma2 - 1.0).sqrt()) + fn entropy(&self) -> f64 { + 0.5 + self.scale.ln() + self.location + consts::LN_SQRT_2PI } } @@ -359,7 +389,7 @@ mod tests { #[test] fn test_mean() { - let mean = |x: LogNormal| x.mean().unwrap(); + let mean = |x: LogNormal| x.mean(); test_exact(-1.0, 0.1, 0.369723444544058982601, mean); test_exact(-1.0, 1.5, 1.133148453066826316829, mean); test_exact(-1.0, 2.5, 8.372897488127264663205, mean); @@ -388,7 +418,7 @@ mod tests { #[test] fn test_variance() { - let variance = |x: LogNormal| x.variance().unwrap(); + let variance = |x: LogNormal| x.variance(); test_absolute(-1.0, 0.1, 0.001373811865368952608715, 1e-16, variance); test_exact(-1.0, 1.5, 10.898468544015731954, variance); test_exact(-1.0, 2.5, 36245.39726189994988081, variance); @@ -417,7 +447,7 @@ mod tests { #[test] fn test_entropy() { - let entropy = |x: LogNormal| x.entropy().unwrap(); + let entropy = |x: LogNormal| x.entropy(); test_exact(-1.0, 0.1, -1.8836465597893728867265104870209210873020761202386, entropy); test_exact(-1.0, 1.5, 0.82440364131283712375834285186996677643338789710028, entropy); test_exact(-1.0, 2.5, 1.335229265078827806963856948173628711311498693546, entropy); @@ -446,7 +476,7 @@ mod tests { #[test] fn test_skewness() { - let skewness = |x: LogNormal| x.skewness().unwrap(); + let skewness = |x: LogNormal| x.skewness(); test_absolute(-1.0, 0.1, 0.30175909933883402945387113824982918009810212213629, 1e-14, skewness); test_exact(-1.0, 1.5, 33.46804679732172529147579024311650645764144530123, skewness); test_absolute(-1.0, 2.5, 11824.007933610287521341659465200553739278936344799, 1e-11, skewness); diff --git a/src/distribution/mod.rs b/src/distribution/mod.rs index 5aa5ff5d..fc573229 100644 --- a/src/distribution/mod.rs +++ b/src/distribution/mod.rs @@ -2,8 +2,8 @@ //! and provides //! concrete implementations for a variety of distributions. use super::statistics::{Max, Min}; -use ::num_traits::{Float, Num}; use num_traits::NumAssignOps; +use num_traits::{Float, Num}; pub use self::bernoulli::Bernoulli; pub use self::beta::{Beta, BetaError}; @@ -265,8 +265,7 @@ pub trait Continuous { /// # Remarks /// /// All methods provided by the `Discrete` trait are unchecked, meaning -/// they can panic if in an invalid state or encountering invalid input -/// depending on the implementing distribution. +/// they can panic if in an invalid state or encountering invalid input depending on the implementing distribution. pub trait Discrete { /// Returns the probability mass function calculated at `x` for a given /// distribution. diff --git a/src/distribution/multinomial.rs b/src/distribution/multinomial.rs index 0794d352..40532919 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, Dim, DimMin, Dyn, OMatrix, OVector}; /// Implements the /// [Multinomial](https://en.wikipedia.org/wiki/Multinomial_distribution) @@ -13,11 +13,11 @@ use nalgebra::{DVector, Dim, Dyn, OMatrix, OVector}; /// /// ``` /// use statrs::distribution::Multinomial; -/// use statrs::statistics::MeanN; +/// use statrs::statistics::*; /// use nalgebra::vector; /// /// let n = Multinomial::new_from_nalgebra(vector![0.3, 0.7], 5).unwrap(); -/// assert_eq!(n.mean().unwrap(), (vector![1.5, 3.5])); +/// assert_eq!(n.mean(), (vector![1.5, 3.5])); /// ``` #[derive(Debug, Clone, PartialEq)] pub struct Multinomial @@ -204,45 +204,26 @@ where res } -impl MeanN> for Multinomial +impl Mean for Multinomial where - D: Dim, - nalgebra::DefaultAllocator: nalgebra::allocator::Allocator, + D: DimMin, + nalgebra::DefaultAllocator: + nalgebra::allocator::Allocator + nalgebra::allocator::Allocator, { - /// Returns the mean of the multinomial distribution - /// - /// # Formula - /// - /// ```text - /// n * p_i for i in 1...k - /// ``` - /// - /// where `n` is the number of trials, `p_i` is the `i`th probability, - /// and `k` is the total number of probabilities - fn mean(&self) -> Option> { - Some(DVector::from_vec( - self.p.iter().map(|x| x * self.n as f64).collect(), - )) + type Mu = OVector; + fn mean(&self) -> Self::Mu { + self.p.scale(self.n as f64) } } -impl VarianceN> for Multinomial +impl Variance for Multinomial where - D: Dim, + D: DimMin, nalgebra::DefaultAllocator: nalgebra::allocator::Allocator + nalgebra::allocator::Allocator, { - /// Returns the variance of the multinomial distribution - /// - /// # Formula - /// - /// ```text - /// n * p_i * (1 - p_i) for i in 1...k - /// ``` - /// - /// where `n` is the number of trials, `p_i` is the `i`th probability, - /// and `k` is the total number of probabilities - fn variance(&self) -> Option> { + type Var = Cholesky; + 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]; @@ -250,13 +231,12 @@ where cov[(y, x)] = elt; }; - for i in 0..self.p.len() { - for j in 0..i { + for j in 0..self.p.len() { + for i in 0..j { offdiag(i, j); } } - cov.fill_lower_triangle_with_upper_triangle(); - Some(cov.scale(self.n as f64)) + Cholesky::new_unchecked(cov.scale(self.n as f64)) } } @@ -364,7 +344,7 @@ where mod tests { use crate::{ distribution::{Discrete, Multinomial, MultinomialError}, - statistics::{MeanN, VarianceN}, + statistics::{Mean, Variance}, }; use nalgebra::{dmatrix, dvector, vector, DimMin, Dyn, OVector}; use std::fmt::{Debug, Display}; @@ -434,7 +414,7 @@ mod tests { #[test] fn test_mean() { - let mean = |x: Multinomial<_>| x.mean().unwrap(); + let mean = |x: Multinomial<_>| x.mean(); test_almost(dvector![0.3, 0.7], 5, dvector![1.5, 3.5], 1e-12, mean); test_almost( dvector![0.1, 0.3, 0.6], @@ -461,7 +441,10 @@ mod tests { #[test] fn test_variance() { - let variance = |x: Multinomial<_>| x.variance().unwrap(); + let variance = |x: Multinomial| { + let l = x.variance().l(); + l.clone()*l.transpose() + }; test_almost( dvector![0.3, 0.7], 5, diff --git a/src/distribution/multivariate_normal.rs b/src/distribution/multivariate_normal.rs index 6eb3b777..1c52d375 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::*; use nalgebra::{Cholesky, Const, DMatrix, DVector, Dim, DimMin, Dyn, OMatrix, OVector}; use std::f64; use std::f64::consts::{E, PI}; @@ -78,11 +78,11 @@ where /// ``` /// use statrs::distribution::{MultivariateNormal, Continuous}; /// use nalgebra::{matrix, vector}; -/// use statrs::statistics::{MeanN, VarianceN}; +/// use statrs::statistics::*; /// /// let mvn = MultivariateNormal::new_from_nalgebra(vector![0., 0.], matrix![1., 0.; 0., 1.]).unwrap(); -/// assert_eq!(mvn.mean().unwrap(), vector![0., 0.]); -/// assert_eq!(mvn.variance().unwrap(), matrix![1., 0.; 0., 1.]); +/// assert_eq!(mvn.mean(), vector![0., 0.]); +/// assert_eq!(mvn.variance(), matrix![1., 0.; 0., 1.]); /// assert_eq!(mvn.pdf(&vector![1., 1.]), 0.05854983152431917); /// ``` #[derive(Clone, PartialEq, Debug)] @@ -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,57 +261,53 @@ where } } -impl Min> for MultivariateNormal +impl Mean for MultivariateNormal where - D: Dim, + D: DimMin, nalgebra::DefaultAllocator: nalgebra::allocator::Allocator + nalgebra::allocator::Allocator, { - /// Returns the minimum value in the domain of the - /// multivariate normal distribution represented by a real vector - fn min(&self) -> OVector { - OMatrix::repeat_generic(self.mu.shape_generic().0, Const::<1>, f64::NEG_INFINITY) + type Mu = OVector; + fn mean(&self) -> Self::Mu { + self.mu.clone() } } -impl Max> for MultivariateNormal +impl Variance for MultivariateNormal where - D: Dim, + D: DimMin, nalgebra::DefaultAllocator: nalgebra::allocator::Allocator + nalgebra::allocator::Allocator, { - /// Returns the maximum value in the domain of the - /// multivariate normal distribution represented by a real vector - fn max(&self) -> OVector { - OMatrix::repeat_generic(self.mu.shape_generic().0, Const::<1>, f64::INFINITY) + type Var = OMatrix; + fn variance(&self) -> Self::Var { + self.cov.clone() } } -impl MeanN> for MultivariateNormal +impl Min> 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()) + /// Returns the minimum value in the domain of the + /// multivariate normal distribution represented by a real vector + fn min(&self) -> OVector { + OMatrix::repeat_generic(self.mu.shape_generic().0, Const::<1>, f64::NEG_INFINITY) } } -impl VarianceN> for MultivariateNormal +impl Max> 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()) + /// Returns the maximum value in the domain of the + /// multivariate normal distribution represented by a real vector + fn max(&self) -> OVector { + OMatrix::repeat_generic(self.mu.shape_generic().0, Const::<1>, f64::INFINITY) } } @@ -379,7 +372,7 @@ mod tests { use crate::{ distribution::{Continuous, MultivariateNormal}, - statistics::{Max, MeanN, Min, Mode, VarianceN}, + statistics::*, }; use super::MultivariateNormalError; @@ -404,8 +397,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 +470,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..3f169428 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::*; use nalgebra::{Cholesky, Const, DMatrix, Dim, DimMin, Dyn, OMatrix, OVector}; use std::f64::consts::PI; @@ -14,11 +14,11 @@ use std::f64::consts::PI; /// ``` /// use statrs::distribution::{MultivariateStudent, Continuous}; /// use nalgebra::{DVector, DMatrix}; -/// use statrs::statistics::{MeanN, VarianceN}; +/// use statrs::statistics::*; /// /// 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.mean(), Some(DVector::from_vec(vec![0., 0.]))); +/// assert_eq!(mvs.variance(), Some(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)] @@ -259,19 +259,14 @@ where } } -impl MeanN> for MultivariateStudent +impl Mean 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>; + fn mean(&self) -> Self::Mu { if self.freedom > 1. { Some(self.location.clone()) } else { @@ -280,23 +275,14 @@ where } } -impl VarianceN> for MultivariateStudent +impl Variance for MultivariateStudent where - D: Dim, + D: DimMin, 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> { + type Var = Option>; + fn variance(&self) -> Self::Var { if self.freedom > 2. { Some(self.scale.clone() * self.freedom / (self.freedom - 2.)) } else { @@ -399,7 +385,7 @@ mod tests { use crate::{ distribution::{Continuous, MultivariateStudent, MultivariateNormal}, - statistics::{Max, MeanN, Min, Mode, VarianceN}, + statistics::*, }; use super::MultivariateStudentError; diff --git a/src/distribution/negative_binomial.rs b/src/distribution/negative_binomial.rs index 27b90aa4..f42252dd 100644 --- a/src/distribution/negative_binomial.rs +++ b/src/distribution/negative_binomial.rs @@ -25,13 +25,13 @@ use std::f64; /// /// ``` /// use statrs::distribution::{NegativeBinomial, Discrete}; -/// use statrs::statistics::DiscreteDistribution; -/// use statrs::prec::almost_eq; +/// use statrs::statistics::*; +/// use approx::assert_relative_eq; /// /// let r = NegativeBinomial::new(4.0, 0.5).unwrap(); -/// assert_eq!(r.mean().unwrap(), 4.0); -/// assert!(almost_eq(r.pmf(0), 0.0625, 1e-8)); -/// assert!(almost_eq(r.pmf(3), 0.15625, 1e-8)); +/// assert_eq!(r.mean(), 4.0); +/// assert_relative_eq!(r.pmf(0), 0.0625, epsilon=1e-8); +/// assert_relative_eq!(r.pmf(3), 0.15625, epsilon=1e-8); /// ``` #[derive(Copy, Clone, PartialEq, Debug)] pub struct NegativeBinomial { @@ -212,38 +212,65 @@ impl Max for NegativeBinomial { } } -impl DiscreteDistribution for NegativeBinomial { - /// Returns the mean of the negative binomial distribution. - /// - /// # Formula - /// - /// ```text - /// r * (1-p) / p - /// ``` - fn mean(&self) -> Option { - Some(self.r * (1.0 - self.p) / self.p) +/// Returns the mean of the negative binomial distribution. +/// +/// # Formula +/// +/// ```text +/// r * (1-p) / p +/// ``` + +impl Mean for NegativeBinomial { + type Mu = f64; + fn mean(&self) -> Self::Mu { + self.r * (1.0 - self.p) / self.p } +} - /// Returns the variance of the negative binomial distribution. - /// - /// # Formula - /// - /// ```text - /// r * (1-p) / p^2 - /// ``` - fn variance(&self) -> Option { - Some(self.r * (1.0 - self.p) / (self.p * self.p)) +/// Returns the variance of the negative binomial distribution. +/// +/// # Formula +/// +/// ```text +/// r * (1-p) / p^2 +/// ``` + +impl Variance for NegativeBinomial { + type Var = f64; + fn variance(&self) -> Self::Var { + // use logarithm as p on (0,1) as a way to multiply without relying on subtraction in 1-p + let ln_var = self.r.ln() + (-self.p).ln_1p() - 2. * self.p.ln(); + ln_var.exp() } +} - /// Returns the skewness of the negative binomial distribution. - /// - /// # Formula - /// - /// ```text - /// (2-p) / sqrt(r * (1-p)) - /// ``` - fn skewness(&self) -> Option { - Some((2.0 - self.p) / f64::sqrt(self.r * (1.0 - self.p))) +/// Returns the skewness of the negative binomial distribution. +/// +/// # Formula +/// +/// ```text +/// (2-p) / sqrt(r * (1-p)) +/// ``` + +impl Skewness for NegativeBinomial { + type Skew = f64; + fn skewness(&self) -> Self::Skew { + (2.0 - self.p) / f64::sqrt(self.r * (1.0 - self.p)) + } +} + +/// Returns the skewness of the negative binomial distribution. +/// +/// # Formula +/// +/// ```text +/// 6/r + p^2 / (r * (1-p)) +/// ``` + +impl ExcessKurtosis for NegativeBinomial { + type Kurt = f64; + fn excess_kurtosis(&self) -> Self::Kurt { + 6. / self.r + self.variance().recip() } } @@ -343,7 +370,7 @@ mod tests { #[test] fn test_mean() { - let mean = |x: NegativeBinomial| x.mean().unwrap(); + let mean = |x: NegativeBinomial| x.mean(); test_exact(4.0, 0.0, f64::INFINITY, mean); test_absolute(3.0, 0.3, 7.0, 1e-15 , mean); test_exact(2.0, 1.0, 0.0, mean); @@ -351,7 +378,7 @@ mod tests { #[test] fn test_variance() { - let variance = |x: NegativeBinomial| x.variance().unwrap(); + let variance = |x: NegativeBinomial| x.variance(); test_exact(4.0, 0.0, f64::INFINITY, variance); test_absolute(3.0, 0.3, 23.333333333333, 1e-12, variance); test_exact(2.0, 1.0, 0.0, variance); @@ -359,7 +386,7 @@ mod tests { #[test] fn test_skewness() { - let skewness = |x: NegativeBinomial| x.skewness().unwrap(); + let skewness = |x: NegativeBinomial| x.skewness(); test_exact(0.0, 0.0, f64::INFINITY, skewness); test_absolute(0.1, 0.3, 6.425396041, 1e-09, skewness); test_exact(1.0, 1.0, f64::INFINITY, skewness); diff --git a/src/distribution/normal.rs b/src/distribution/normal.rs index 653b1cc8..0407ec7f 100644 --- a/src/distribution/normal.rs +++ b/src/distribution/normal.rs @@ -11,11 +11,12 @@ use std::f64; /// /// ``` /// use statrs::distribution::{Normal, Continuous}; -/// use statrs::statistics::Distribution; +/// use statrs::statistics::*; +/// use approx::assert_relative_eq; /// /// let n = Normal::new(0.0, 1.0).unwrap(); -/// assert_eq!(n.mean().unwrap(), 0.0); -/// assert_eq!(n.pdf(1.0), 0.2419707245191433497978); +/// assert_eq!(n.mean(), 0.0); +/// assert_relative_eq!(n.pdf(1.0), 0.2419707245191433497978); /// ``` #[derive(Copy, Clone, PartialEq, Debug)] pub struct Normal { @@ -97,6 +98,10 @@ impl Normal { std_dev: 1.0, } } + + pub fn std_dev(&self) -> f64 { + self.std_dev + } } impl std::fmt::Display for Normal { @@ -206,36 +211,51 @@ impl Max for Normal { } } -impl Distribution for Normal { - /// Returns the mean of the normal distribution - /// - /// # Remarks - /// - /// This is the same mean used to construct the distribution - fn mean(&self) -> Option { - Some(self.mean) +/// Returns the mean of the normal distribution +/// +/// # Remarks +/// +/// This is the same mean used to construct the distribution + +impl Mean for Normal { + type Mu = f64; + fn mean(&self) -> Self::Mu { + self.mean } +} - /// Returns the variance of the normal distribution - /// - /// # Formula - /// - /// ```text - /// σ^2 - /// ``` - /// - /// where `σ` is the standard deviation - fn variance(&self) -> Option { - Some(self.std_dev * self.std_dev) +/// Returns the variance of the normal distribution +/// +/// # Formula +/// +/// ```text +/// σ^2 +/// ``` +/// +/// where `σ` is the standard deviation + +impl Variance for Normal { + type Var = f64; + fn variance(&self) -> Self::Var { + self.std_dev * self.std_dev } +} - /// Returns the standard deviation of the normal distribution - /// # Remarks - /// This is the same standard deviation used to construct the distribution - fn std_dev(&self) -> Option { - Some(self.std_dev) +impl Skewness for Normal { + type Skew = f64; + fn skewness(&self) -> Self::Skew { + 0.0 } +} +impl ExcessKurtosis for Normal { + type Kurt = f64; + fn excess_kurtosis(&self) -> Self::Kurt { + 0.0 + } +} + +impl Entropy for Normal { /// Returns the entropy of the normal distribution /// /// # Formula @@ -245,19 +265,8 @@ impl Distribution for Normal { /// ``` /// /// where `σ` is the standard deviation - fn entropy(&self) -> Option { - Some(self.std_dev.ln() + consts::LN_SQRT_2PIE) - } - - /// Returns the skewness of the normal distribution - /// - /// # Formula - /// - /// ```text - /// 0 - /// ``` - fn skewness(&self) -> Option { - Some(0.0) + fn entropy(&self) -> f64 { + self.std_dev.ln() + consts::LN_SQRT_2PIE } } @@ -394,7 +403,7 @@ mod tests { #[test] fn test_variance() { - let variance = |x: Normal| x.variance().unwrap(); + let variance = |x: Normal| x.variance(); test_exact(0.0, 0.1, 0.1 * 0.1, variance); test_exact(0.0, 1.0, 1.0, variance); test_exact(0.0, 10.0, 100.0, variance); @@ -403,7 +412,7 @@ mod tests { #[test] fn test_entropy() { - let entropy = |x: Normal| x.entropy().unwrap(); + let entropy = |x: Normal| x.entropy(); test_absolute(0.0, 0.1, -0.8836465597893729422377, 1e-15, entropy); test_exact(0.0, 1.0, 1.41893853320467274178, entropy); test_exact(0.0, 10.0, 3.721523626198718425798, entropy); @@ -412,7 +421,7 @@ mod tests { #[test] fn test_skewness() { - let skewness = |x: Normal| x.skewness().unwrap(); + let skewness = |x: Normal| x.skewness(); test_exact(0.0, 0.1, 0.0, skewness); test_exact(4.0, 1.0, 0.0, skewness); test_exact(0.3, 10.0, 0.0, skewness); @@ -558,8 +567,8 @@ mod tests { fn test_default() { let n = Normal::default(); - let n_mean = n.mean().unwrap(); - let n_std = n.std_dev().unwrap(); + let n_mean = n.mean(); + let n_std = n.std_dev(); // Check that the mean of the distribution is close to 0 assert_almost_eq!(n_mean, 0.0, 1e-15); diff --git a/src/distribution/pareto.rs b/src/distribution/pareto.rs index d3dd9681..e056d8f9 100644 --- a/src/distribution/pareto.rs +++ b/src/distribution/pareto.rs @@ -9,12 +9,12 @@ use std::f64; /// /// ``` /// use statrs::distribution::{Pareto, Continuous}; -/// use statrs::statistics::Distribution; -/// use statrs::prec; +/// use statrs::statistics::*; +/// use approx::assert_relative_eq; /// /// let p = Pareto::new(1.0, 2.0).unwrap(); -/// assert_eq!(p.mean().unwrap(), 2.0); -/// assert!(prec::almost_eq(p.pdf(2.0), 0.25, 1e-15)); +/// assert_eq!(p.mean(), Some(2.0)); +/// assert_relative_eq!(p.pdf(2.0), 0.25); /// ``` #[derive(Copy, Clone, PartialEq, Debug)] pub struct Pareto { @@ -218,42 +218,44 @@ impl Max for Pareto { } } -impl Distribution for Pareto { - /// Returns the mean of the Pareto distribution - /// - /// # Formula - /// - /// ```text - /// if α <= 1 { - /// f64::INFINITY - /// } else { - /// (α * x_m)/(α - 1) - /// } - /// ``` - /// - /// where `x_m` is the scale and `α` is the shape - fn mean(&self) -> Option { +/// Returns the mean of the Pareto distribution +/// +/// # Formula +/// +/// ```text +/// (α * x_m)/(α - 1), where α > 1 +/// ``` +/// +/// where `x_m` is the scale and `α` is the shape + +impl Mean for Pareto { + type Mu = Option; + fn mean(&self) -> Self::Mu { if self.shape <= 1.0 { None } else { Some((self.shape * self.scale) / (self.shape - 1.0)) } } +} - /// Returns the variance of the Pareto distribution - /// - /// # Formula - /// - /// ```text - /// if α <= 2 { - /// f64::INFINITY - /// } else { - /// (x_m/(α - 1))^2 * (α/(α - 2)) - /// } - /// ``` - /// - /// where `x_m` is the scale and `α` is the shape - fn variance(&self) -> Option { +/// Returns the variance of the Pareto distribution +/// +/// # Formula +/// +/// ```text +/// if α <= 2 { +/// f64::INFINITY +/// } else { +/// (x_m/(α - 1))^2 * (α/(α - 2)) +/// } +/// ``` +/// +/// where `x_m` is the scale and `α` is the shape + +impl Variance for Pareto { + type Var = Option; + fn variance(&self) -> Self::Var { if self.shape <= 2.0 { None } else { @@ -261,36 +263,27 @@ impl Distribution for Pareto { Some(a * a * self.shape / (self.shape - 2.0)) } } +} - /// Returns the entropy for the Pareto distribution - /// - /// # Formula - /// - /// ```text - /// ln(α/x_m) - 1/α - 1 - /// ``` - /// - /// where `x_m` is the scale and `α` is the shape - fn entropy(&self) -> Option { - Some(self.shape.ln() - self.scale.ln() - (1.0 / self.shape) - 1.0) - } +/// Returns the skewness of the Pareto distribution +/// +/// # Panics +/// +/// If `α <= 3.0` +/// +/// where `α` is the shape +/// +/// # Formula +/// +/// ```text +/// (2*(α + 1)/(α - 3))*sqrt((α - 2)/α) +/// ``` +/// +/// where `α` is the shape - /// Returns the skewness of the Pareto distribution - /// - /// # Panics - /// - /// If `α <= 3.0` - /// - /// where `α` is the shape - /// - /// # Formula - /// - /// ```text - /// (2*(α + 1)/(α - 3))*sqrt((α - 2)/α) - /// ``` - /// - /// where `α` is the shape - fn skewness(&self) -> Option { +impl Skewness for Pareto { + type Skew = Option; + fn skewness(&self) -> Self::Skew { if self.shape <= 3.0 { None } else { @@ -302,6 +295,39 @@ impl Distribution for Pareto { } } +/// Returns the excess kurtosis of the Pareto distribution + +impl ExcessKurtosis for Pareto { + type Kurt = Option; + fn excess_kurtosis(&self) -> Self::Kurt { + if self.shape <= 4.0 { + None + } else { + Some( + 6.0 * (self.shape.powi(3) + self.shape.powi(2) - 6.0 * self.shape - 2.0) + / self.shape + / (self.shape - 3.0) + / (self.shape - 4.0), + ) + } + } +} + +impl Entropy for Pareto { + /// Returns the entropy for the Pareto distribution + /// + /// # Formula + /// + /// ```text + /// ln(α/x_m) - 1/α - 1 + /// ``` + /// + /// where `x_m` is the scale and `α` is the shape + fn entropy(&self) -> f64 { + self.shape.ln() - self.scale.ln() - (1.0 / self.shape) - 1.0 + } +} + impl Median for Pareto { /// Returns the median of the Pareto distribution /// @@ -422,7 +448,7 @@ mod tests { #[test] fn test_entropy() { - let entropy = |x: Pareto| x.entropy().unwrap(); + let entropy = |x: Pareto| x.entropy(); test_exact(0.1, 0.1, -11.0, entropy); test_exact(1.0, 1.0, -2.0, entropy); test_exact(10.0, 10.0, -1.1, entropy); diff --git a/src/distribution/poisson.rs b/src/distribution/poisson.rs index 1c0598c4..36f2c877 100644 --- a/src/distribution/poisson.rs +++ b/src/distribution/poisson.rs @@ -10,12 +10,12 @@ use std::f64; /// /// ``` /// use statrs::distribution::{Poisson, Discrete}; -/// use statrs::statistics::Distribution; -/// use statrs::prec; +/// use statrs::statistics::*; +/// use approx::assert_relative_eq; /// /// let n = Poisson::new(1.0).unwrap(); -/// assert_eq!(n.mean().unwrap(), 1.0); -/// assert!(prec::almost_eq(n.pmf(1), 0.367879441171442, 1e-15)); +/// assert_eq!(n.mean(), 1.0); +/// assert_relative_eq!(n.pmf(1), 0.367879441171442, epsilon=1e-14); /// ``` #[derive(Copy, Clone, PartialEq, Debug)] pub struct Poisson { @@ -173,62 +173,82 @@ impl Max for Poisson { } } -impl Distribution for Poisson { - /// Returns the mean of the poisson distribution - /// - /// # Formula - /// - /// ```text - /// λ - /// ``` - /// - /// where `λ` is the rate - fn mean(&self) -> Option { - Some(self.lambda) +/// Returns the mean of the poisson distribution +/// +/// # Formula +/// +/// ```text +/// λ +/// ``` +/// +/// where `λ` is the rate + +impl Mean for Poisson { + type Mu = f64; + fn mean(&self) -> Self::Mu { + self.lambda } +} - /// Returns the variance of the poisson distribution - /// - /// # Formula - /// - /// ```text - /// λ - /// ``` - /// - /// where `λ` is the rate - fn variance(&self) -> Option { - Some(self.lambda) +/// Returns the variance of the poisson distribution +/// +/// # Formula +/// +/// ```text +/// λ +/// ``` +/// +/// where `λ` is the rate + +impl Variance for Poisson { + type Var = f64; + fn variance(&self) -> Self::Var { + self.lambda } +} - /// Returns the entropy of the poisson distribution - /// - /// # Formula - /// - /// ```text - /// (1 / 2) * ln(2πeλ) - 1 / (12λ) - 1 / (24λ^2) - 19 / (360λ^3) - /// ``` - /// - /// where `λ` is the rate - fn entropy(&self) -> Option { - Some( - 0.5 * (2.0 * f64::consts::PI * f64::consts::E * self.lambda).ln() - - 1.0 / (12.0 * self.lambda) - - 1.0 / (24.0 * self.lambda * self.lambda) - - 19.0 / (360.0 * self.lambda * self.lambda * self.lambda), - ) +/// Returns the skewness of the poisson distribution +/// +/// # Formula +/// +/// ```text +/// λ^(-1/2) +/// ``` +/// +/// where `λ` is the rate + +impl Skewness for Poisson { + type Skew = f64; + fn skewness(&self) -> Self::Skew { + self.lambda.sqrt().recip() + } +} + +/// Returns the excess kurtosis for the Poisson distribution + +impl ExcessKurtosis for Poisson { + type Kurt = f64; + fn excess_kurtosis(&self) -> Self::Kurt { + self.lambda.recip() } +} - /// Returns the skewness of the poisson distribution +impl Entropy for Poisson { + /// Returns the entropy of the poisson distribution /// /// # Formula /// /// ```text - /// λ^(-1/2) + /// (1 / 2) * ln(2πeλ) - 1 / (12λ) - 1 / (24λ^2) - 19 / (360λ^3) /// ``` /// /// where `λ` is the rate - fn skewness(&self) -> Option { - Some(1.0 / self.lambda.sqrt()) + fn entropy(&self) -> f64 { + // TODO: this is only valid for large lambda + 0.5 * (2.0 * f64::consts::PI * f64::consts::E * self.lambda).ln() + - 1.0 / (12.0 * self.lambda) + - 1.0 / (24.0 * self.lambda * self.lambda) + - 19.0 / (360.0 * self.lambda * self.lambda * self.lambda) } } @@ -361,7 +381,7 @@ mod tests { #[test] fn test_mean() { - let mean = |x: Poisson| x.mean().unwrap(); + let mean = |x: Poisson| x.mean(); test_exact(1.5, 1.5, mean); test_exact(5.4, 5.4, mean); test_exact(10.8, 10.8, mean); @@ -369,7 +389,7 @@ mod tests { #[test] fn test_variance() { - let variance = |x: Poisson| x.variance().unwrap(); + let variance = |x: Poisson| x.variance(); test_exact(1.5, 1.5, variance); test_exact(5.4, 5.4, variance); test_exact(10.8, 10.8, variance); @@ -377,7 +397,7 @@ mod tests { #[test] fn test_entropy() { - let entropy = |x: Poisson| x.entropy().unwrap(); + let entropy = |x: Poisson| x.entropy(); test_absolute(1.5, 1.531959153102376331946, 1e-15, entropy); test_absolute(5.4, 2.244941839577643504608, 1e-15, entropy); test_exact(10.8, 2.600596429676975222694, entropy); @@ -385,7 +405,7 @@ mod tests { #[test] fn test_skewness() { - let skewness = |x: Poisson| x.skewness().unwrap(); + let skewness = |x: Poisson| x.skewness(); test_absolute(1.5, 0.8164965809277260327324, 1e-15, skewness); test_absolute(5.4, 0.4303314829119352094644, 1e-16, skewness); test_absolute(10.8, 0.3042903097250922852539, 1e-16, skewness); diff --git a/src/distribution/students_t.rs b/src/distribution/students_t.rs index 485092b5..f11e6e8a 100644 --- a/src/distribution/students_t.rs +++ b/src/distribution/students_t.rs @@ -10,12 +10,12 @@ use std::f64; /// /// ``` /// use statrs::distribution::{StudentsT, Continuous}; -/// use statrs::statistics::Distribution; -/// use statrs::prec; +/// use statrs::statistics::*; +/// use approx::assert_relative_eq; /// /// let n = StudentsT::new(0.0, 1.0, 2.0).unwrap(); -/// assert_eq!(n.mean().unwrap(), 0.0); -/// assert!(prec::almost_eq(n.pdf(0.0), 0.353553390593274, 1e-15)); +/// assert_eq!(n.mean(), Some(0.0)); +/// assert_relative_eq!(n.pdf(0.0), 0.353553390593274, epsilon=1e-14); /// ``` #[derive(Copy, Clone, PartialEq, Debug)] pub struct StudentsT { @@ -271,48 +271,54 @@ impl Max for StudentsT { } } -impl Distribution for StudentsT { - /// Returns the mean of the student's t-distribution - /// - /// # None - /// - /// If `freedom <= 1.0` - /// - /// # Formula - /// - /// ```text - /// μ - /// ``` - /// - /// where `μ` is the location - fn mean(&self) -> Option { +/// Returns the mean of the student's t-distribution +/// +/// # None +/// +/// If `freedom <= 1.0` +/// +/// # Formula +/// +/// ```text +/// μ +/// ``` +/// +/// where `μ` is the location + +impl Mean for StudentsT { + type Mu = Option; + fn mean(&self) -> Self::Mu { if self.freedom <= 1.0 { None } else { Some(self.location) } } +} - /// Returns the variance of the student's t-distribution - /// - /// # None - /// - /// If `freedom <= 2.0` - /// - /// # Formula - /// - /// ```text - /// if v == f64::INFINITY { - /// Some(σ^2) - /// } else if freedom > 2.0 { - /// Some(v * σ^2 / (v - 2)) - /// } else { - /// None - /// } - /// ``` - /// - /// where `σ` is the scale and `v` is the freedom - fn variance(&self) -> Option { +/// Returns the variance of the student's t-distribution +/// +/// # None +/// +/// If `freedom <= 2.0` +/// +/// # Formula +/// +/// ```text +/// if v == f64::INFINITY { +/// Some(σ^2) +/// } else if freedom > 2.0 { +/// Some(v * σ^2 / (v - 2)) +/// } else { +/// None +/// } +/// ``` +/// +/// where `σ` is the scale and `v` is the freedom + +impl Variance for StudentsT { + type Var = Option; + fn variance(&self) -> Self::Var { if self.freedom.is_infinite() { Some(self.scale * self.scale) } else if self.freedom > 2.0 { @@ -321,7 +327,48 @@ impl Distribution for StudentsT { None } } +} +/// Returns the skewness of the student's t-distribution +/// +/// # None +/// +/// If `x <= 3.0` +/// +/// # Formula +/// +/// ```text +/// 0 +/// ``` + +impl Skewness for StudentsT { + type Skew = Option; + fn skewness(&self) -> Self::Skew { + if self.freedom <= 3.0 { + None + } else { + Some(0.0) + } + } +} + +/// docs +/// moar docs + +impl ExcessKurtosis for StudentsT { + type Kurt = Option; + fn excess_kurtosis(&self) -> Self::Kurt { + if self.freedom <= 2.0 { + None + } else if self.freedom <= 4.0 { + Some(f64::INFINITY) + } else { + Some(6.0 / (self.freedom - 4.0)) + } + } +} + +impl Entropy for StudentsT { /// Returns the entropy for the student's t-distribution /// /// # Formula @@ -333,7 +380,7 @@ impl Distribution for StudentsT { /// /// where `σ` is the scale, `v` is the freedom, `ψ` is the digamma function, and `B` is the /// beta function - fn entropy(&self) -> Option { + fn entropy(&self) -> f64 { // generalised Student's T is related to normal Student's T by `Y = μ + σ X` // where `X` is distributed as Student's T, plugging into the definition // of entropy shows scaling affects the entropy by an additive constant `- ln σ` @@ -341,26 +388,7 @@ impl Distribution for StudentsT { let result = (self.freedom + 1.0) / 2.0 * (gamma::digamma((self.freedom + 1.0) / 2.0) - gamma::digamma(self.freedom / 2.0)) + (self.freedom.sqrt() * beta::beta(self.freedom / 2.0, 0.5)).ln(); - Some(result + shift) - } - - /// Returns the skewness of the student's t-distribution - /// - /// # None - /// - /// If `x <= 3.0` - /// - /// # Formula - /// - /// ```text - /// 0 - /// ``` - fn skewness(&self) -> Option { - if self.freedom <= 3.0 { - None - } else { - Some(0.0) - } + result + shift } } diff --git a/src/distribution/triangular.rs b/src/distribution/triangular.rs index 6be868bb..268c9969 100644 --- a/src/distribution/triangular.rs +++ b/src/distribution/triangular.rs @@ -10,10 +10,10 @@ use std::f64; /// /// ``` /// use statrs::distribution::{Triangular, Continuous}; -/// use statrs::statistics::Distribution; +/// use statrs::statistics::*; /// /// let n = Triangular::new(0.0, 5.0, 2.5).unwrap(); -/// assert_eq!(n.mean().unwrap(), 7.5 / 3.0); +/// assert_eq!(n.mean(), 7.5 / 3.0); /// assert_eq!(n.pdf(2.5), 5.0 / 12.5); /// ``` #[derive(Copy, Clone, PartialEq, Debug)] @@ -236,60 +236,75 @@ impl Max for Triangular { } } -impl Distribution for Triangular { - /// Returns the mean of the triangular distribution - /// - /// # Formula - /// - /// ```text - /// (min + max + mode) / 3 - /// ``` - fn mean(&self) -> Option { - Some((self.min + self.max + self.mode) / 3.0) +/// Returns the mean of the triangular distribution +/// +/// # Formula +/// +/// ```text +/// (min + max + mode) / 3 +/// ``` +impl Mean for Triangular { + type Mu = f64; + fn mean(&self) -> Self::Mu { + (self.min + self.max + self.mode) / 3.0 } +} - /// Returns the variance of the triangular distribution - /// - /// # Formula - /// - /// ```text - /// (min^2 + max^2 + mode^2 - min * max - min * mode - max * mode) / 18 - /// ``` - fn variance(&self) -> Option { +/// Returns the variance of the triangular distribution +/// +/// # Formula +/// +/// ```text +/// (min^2 + max^2 + mode^2 - min * max - min * mode - max * mode) / 18 +/// ``` +impl Variance for Triangular { + type Var = f64; + fn variance(&self) -> Self::Var { let a = self.min; let b = self.max; let c = self.mode; - Some((a * a + b * b + c * c - a * b - a * c - b * c) / 18.0) + (a * a + b * b + c * c - a * b - a * c - b * c) / 18.0 } +} +/// Returns the skewness of the triangular distribution +/// +/// # Formula +/// +/// ```text +/// (sqrt(2) * (min + max - 2 * mode) * (2 * min - max - mode) * (min - 2 * +/// max + mode)) / +/// ( 5 * (min^2 + max^2 + mode^2 - min * max - min * mode - max * mode)^(3 +/// / 2)) +/// ``` +impl Skewness for Triangular { + type Skew = f64; + fn skewness(&self) -> Self::Skew { + let a = self.min; + let b = self.max; + let c = self.mode; + let q = f64::consts::SQRT_2 * (a + b - 2.0 * c) * (2.0 * a - b - c) * (a - 2.0 * b + c); + let d = 5.0 * (a * a + b * b + c * c - a * b - a * c - b * c).powf(3.0 / 2.0); + q / d + } +} - /// Returns the entropy of the triangular distribution - /// - /// # Formula - /// - /// ```text - /// 1 / 2 + ln((max - min) / 2) - /// ``` - fn entropy(&self) -> Option { - Some(0.5 + ((self.max - self.min) / 2.0).ln()) +impl ExcessKurtosis for Triangular { + type Kurt = f64; + fn excess_kurtosis(&self) -> Self::Kurt { + -0.6 } +} - /// Returns the skewness of the triangular distribution +impl Entropy for Triangular { + /// Returns the entropy of the triangular distribution /// /// # Formula /// /// ```text - /// (sqrt(2) * (min + max - 2 * mode) * (2 * min - max - mode) * (min - 2 * - /// max + mode)) / - /// ( 5 * (min^2 + max^2 + mode^2 - min * max - min * mode - max * mode)^(3 - /// / 2)) + /// 1 / 2 + ln((max - min) / 2) /// ``` - fn skewness(&self) -> Option { - let a = self.min; - let b = self.max; - let c = self.mode; - let q = f64::consts::SQRT_2 * (a + b - 2.0 * c) * (2.0 * a - b - c) * (a - 2.0 * b + c); - let d = 5.0 * (a * a + b * b + c * c - a * b - a * c - b * c).powf(3.0 / 2.0); - Some(q / d) + fn entropy(&self) -> f64 { + 0.5 + ((self.max - self.min) / 2.0).ln() } } @@ -437,7 +452,7 @@ mod tests { #[test] fn test_variance() { - let variance = |x: Triangular| x.variance().unwrap(); + let variance = |x: Triangular| x.variance(); test_exact(0.0, 1.0, 0.5, 0.75 / 18.0, variance); test_exact(0.0, 1.0, 0.75, 0.8125 / 18.0, variance); test_exact(-5.0, 8.0, -3.5, 151.75 / 18.0, variance); @@ -448,7 +463,7 @@ mod tests { #[test] fn test_entropy() { - let entropy = |x: Triangular| x.entropy().unwrap(); + let entropy = |x: Triangular| x.entropy(); test_absolute(0.0, 1.0, 0.5, -0.1931471805599453094172, 1e-16, entropy); test_absolute(0.0, 1.0, 0.75, -0.1931471805599453094172, 1e-16, entropy); test_exact(-5.0, 8.0, -3.5, 2.371802176901591426636, entropy); @@ -459,7 +474,7 @@ mod tests { #[test] fn test_skewness() { - let skewness = |x: Triangular| x.skewness().unwrap(); + let skewness = |x: Triangular| x.skewness(); test_exact(0.0, 1.0, 0.5, 0.0, skewness); test_exact(0.0, 1.0, 0.75, -0.4224039833745502226059, skewness); test_exact(-5.0, 8.0, -3.5, 0.5375093589712976359809, skewness); diff --git a/src/distribution/uniform.rs b/src/distribution/uniform.rs index 96a2f49c..6289948f 100644 --- a/src/distribution/uniform.rs +++ b/src/distribution/uniform.rs @@ -11,10 +11,10 @@ use std::fmt::Debug; /// /// ``` /// use statrs::distribution::{Uniform, Continuous}; -/// use statrs::statistics::Distribution; +/// use statrs::statistics::*; /// /// let n = Uniform::new(0.0, 1.0).unwrap(); -/// assert_eq!(n.mean().unwrap(), 0.5); +/// assert_eq!(n.mean(), 0.5); /// assert_eq!(n.pdf(0.5), 1.0); /// ``` #[derive(Debug, Copy, Clone, PartialEq)] @@ -105,6 +105,11 @@ impl Uniform { pub fn standard() -> Self { Self { min: 0.0, max: 1.0 } } + + pub fn std_dev(&self) -> f64 { + let dist = self.max - self.min; + (dist * dist / 12.).sqrt() + } } impl Default for Uniform { @@ -192,49 +197,70 @@ impl Max for Uniform { } } -impl Distribution for Uniform { - /// Returns the mean for the continuous uniform distribution - /// - /// # Formula - /// - /// ```text - /// (min + max) / 2 - /// ``` - fn mean(&self) -> Option { - Some((self.min + self.max) / 2.0) +/// Returns the mean for the continuous uniform distribution +/// +/// # Formula +/// +/// ```text +/// (min + max) / 2 +/// ``` + +impl Mean for Uniform { + type Mu = f64; + fn mean(&self) -> Self::Mu { + (self.min + self.max) / 2.0 } +} - /// Returns the variance for the continuous uniform distribution - /// - /// # Formula - /// - /// ```text - /// (max - min)^2 / 12 - /// ``` - fn variance(&self) -> Option { - Some((self.max - self.min) * (self.max - self.min) / 12.0) +/// Returns the variance for the continuous uniform distribution +/// +/// # Formula +/// +/// ```text +/// (max - min)^2 / 12 +/// ``` + +impl Variance for Uniform { + type Var = f64; + fn variance(&self) -> Self::Var { + (self.max - self.min) * (self.max - self.min) / 12.0 } +} - /// Returns the entropy for the continuous uniform distribution - /// - /// # Formula - /// - /// ```text - /// ln(max - min) - /// ``` - fn entropy(&self) -> Option { - Some((self.max - self.min).ln()) +/// Returns the skewness for the continuous uniform distribution +/// +/// # Formula +/// +/// ```text +/// 0 +/// ``` + +impl Skewness for Uniform { + type Skew = f64; + fn skewness(&self) -> Self::Skew { + 0.0 + } +} + +/// docs + +impl ExcessKurtosis for Uniform { + type Kurt = f64; + fn excess_kurtosis(&self) -> Self::Kurt { + -1.2 } +} - /// Returns the skewness for the continuous uniform distribution +impl Entropy for Uniform { + /// Returns the entropy for the continuous uniform distribution /// /// # Formula /// /// ```text - /// 0 + /// ln(max - min) /// ``` - fn skewness(&self) -> Option { - Some(0.0) + fn entropy(&self) -> f64 { + (self.max - self.min).ln() } } @@ -347,7 +373,7 @@ mod tests { #[test] fn test_variance() { - let variance = |x: Uniform| x.variance().unwrap(); + let variance = |x: Uniform| x.variance(); test_exact(-0.0, 2.0, 1.0 / 3.0, variance); test_exact(0.0, 2.0, 1.0 / 3.0, variance); test_absolute(0.1, 4.0, 1.2675, 1e-15, variance); @@ -356,7 +382,7 @@ mod tests { #[test] fn test_entropy() { - let entropy = |x: Uniform| x.entropy().unwrap(); + let entropy = |x: Uniform| x.entropy(); test_exact(-0.0, 2.0, 0.6931471805599453094172, entropy); test_exact(0.0, 2.0, 0.6931471805599453094172, entropy); test_absolute(0.1, 4.0, 1.360976553135600743431, 1e-15, entropy); @@ -366,7 +392,7 @@ mod tests { #[test] fn test_skewness() { - let skewness = |x: Uniform| x.skewness().unwrap(); + let skewness = |x: Uniform| x.skewness(); test_exact(-0.0, 2.0, 0.0, skewness); test_exact(0.0, 2.0, 0.0, skewness); test_exact(0.1, 4.0, 0.0, skewness); @@ -524,8 +550,8 @@ mod tests { fn test_default() { let n = Uniform::default(); - let n_mean = n.mean().unwrap(); - let n_std = n.std_dev().unwrap(); + let n_mean = n.mean(); + let n_std = n.std_dev(); // Check that the mean of the distribution is close to 1 / 2 assert_almost_eq!(n_mean, 0.5, 1e-15); diff --git a/src/distribution/weibull.rs b/src/distribution/weibull.rs index 6447aaa1..ed2f29a6 100644 --- a/src/distribution/weibull.rs +++ b/src/distribution/weibull.rs @@ -11,13 +11,12 @@ use std::f64; /// /// ``` /// use statrs::distribution::{Weibull, Continuous}; -/// use statrs::statistics::Distribution; -/// use statrs::prec; +/// use statrs::statistics::*; +/// use approx::assert_relative_eq; /// /// let n = Weibull::new(10.0, 1.0).unwrap(); -/// assert!(prec::almost_eq(n.mean().unwrap(), -/// 0.95135076986687318362924871772654021925505786260884, 1e-15)); -/// assert_eq!(n.pdf(1.0), 3.6787944117144232159552377016146086744581113103177); +/// assert_relative_eq!(n.mean(), 0.951350769866873183629, epsilon = 1e-14); +/// assert_relative_eq!(n.pdf(1.0), 3.67879441171442321595, epsilon = 1e-14); /// ``` #[derive(Copy, Clone, PartialEq, Debug)] pub struct Weibull { @@ -213,36 +212,89 @@ impl Max for Weibull { } } -impl Distribution for Weibull { - /// Returns the mean of the weibull distribution - /// - /// # Formula - /// - /// ```text - /// λΓ(1 + 1 / k) - /// ``` - /// - /// where `k` is the shape, `λ` is the scale, and `Γ` is - /// the gamma function - fn mean(&self) -> Option { - Some(self.scale * gamma::gamma(1.0 + 1.0 / self.shape)) +/// Returns the mean of the weibull distribution +/// +/// # Formula +/// +/// ```text +/// λΓ(1 + 1 / k) +/// ``` +/// +/// where `k` is the shape, `λ` is the scale, and `Γ` is +/// the gamma function + +impl Mean for Weibull { + type Mu = f64; + fn mean(&self) -> Self::Mu { + self.scale * gamma::gamma(1.0 + 1.0 / self.shape) } +} - /// Returns the variance of the weibull distribution - /// - /// # Formula - /// - /// ```text - /// λ^2 * (Γ(1 + 2 / k) - Γ(1 + 1 / k)^2) - /// ``` - /// - /// where `k` is the shape, `λ` is the scale, and `Γ` is - /// the gamma function - fn variance(&self) -> Option { - let mean = self.mean()?; - Some(self.scale * self.scale * gamma::gamma(1.0 + 2.0 / self.shape) - mean * mean) +/// Returns the variance of the weibull distribution +/// +/// # Formula +/// +/// ```text +/// λ^2 * (Γ(1 + 2 / k) - Γ(1 + 1 / k)^2) +/// ``` +/// +/// where `k` is the shape, `λ` is the scale, and `Γ` is +/// the gamma function + +impl Variance for Weibull { + type Var = f64; + fn variance(&self) -> Self::Var { + let mean = self.mean(); + self.scale * self.scale * gamma::gamma(1.0 + 2.0 / self.shape) - mean * mean + } +} + +/// Returns the skewness of the weibull distribution +/// +/// # Formula +/// +/// ```text +/// (Γ(1 + 3 / k) * λ^3 - 3μσ^2 - μ^3) / σ^3 +/// ``` +/// +/// where `k` is the shape, `λ` is the scale, and `Γ` is +/// the gamma function, `μ` is the mean of the distribution. +/// and `σ` the standard deviation of the distribution + +impl Skewness for Weibull { + type Skew = f64; + fn skewness(&self) -> Self::Skew { + let mu = self.mean(); + let sigma2 = self.variance(); + let sigma3 = self.variance().powf(1.5); + + (self.scale * self.scale * self.scale * gamma::gamma(1.0 + 3.0 / self.shape) + - 3.0 * sigma2 * mu + - (mu * mu * mu)) + / sigma3 + } +} + +impl ExcessKurtosis for Weibull { + type Kurt = f64; + fn excess_kurtosis(&self) -> Self::Kurt { + let gam_i = |i: f64| gamma::gamma(1. + i / self.shape); + let gam_1 = gam_i(1.); + let gam_2 = gam_i(2.); + let gam_3 = gam_i(3.); + let gam_4 = gam_i(4.); + + let numer = -6. * gam_1.powi(4) + 12. * gam_1.powi(2) * gam_2 + - 3. * gam_2.powi(2) + - 4. * gam_1 * gam_3 + + gam_4; + + let denom = (gam_2 - gam_1.powi(2)).powi(2); + numer / denom } +} +impl Entropy for Weibull { /// Returns the entropy of the weibull distribution /// /// # Formula @@ -253,34 +305,8 @@ impl Distribution for Weibull { /// /// where `k` is the shape, `λ` is the scale, and `γ` is /// the Euler-Mascheroni constant - fn entropy(&self) -> Option { - let entr = consts::EULER_MASCHERONI * (1.0 - 1.0 / self.shape) - + (self.scale / self.shape).ln() - + 1.0; - Some(entr) - } - - /// Returns the skewness of the weibull distribution - /// - /// # Formula - /// - /// ```text - /// (Γ(1 + 3 / k) * λ^3 - 3μσ^2 - μ^3) / σ^3 - /// ``` - /// - /// where `k` is the shape, `λ` is the scale, and `Γ` is - /// the gamma function, `μ` is the mean of the distribution. - /// and `σ` the standard deviation of the distribution - fn skewness(&self) -> Option { - let mu = self.mean()?; - let sigma = self.std_dev()?; - let sigma2 = sigma * sigma; - let sigma3 = sigma2 * sigma; - let skew = (self.scale * self.scale * self.scale * gamma::gamma(1.0 + 3.0 / self.shape) - - 3.0 * sigma2 * mu - - (mu * mu * mu)) - / sigma3; - Some(skew) + fn entropy(&self) -> f64 { + consts::EULER_MASCHERONI * (1.0 - 1.0 / self.shape) + (self.scale / self.shape).ln() + 1.0 } } @@ -406,7 +432,7 @@ mod tests { #[test] fn test_mean() { - let mean = |x: Weibull| x.mean().unwrap(); + let mean = |x: Weibull| x.mean(); test_exact(1.0, 0.1, 0.1, mean); test_exact(1.0, 1.0, 1.0, mean); test_absolute(10.0, 10.0, 9.5135076986687318362924871772654021925505786260884, 1e-14, mean); @@ -415,7 +441,7 @@ mod tests { #[test] fn test_variance() { - let variance = |x: Weibull| x.variance().unwrap(); + let variance = |x: Weibull| x.variance(); test_absolute(1.0, 0.1, 0.01, 1e-16, variance); test_absolute(1.0, 1.0, 1.0, 1e-14, variance); test_absolute(10.0, 10.0, 1.3100455073468309147154581687505295026863354547057, 1e-12, variance); @@ -424,7 +450,7 @@ mod tests { #[test] fn test_entropy() { - let entropy = |x: Weibull| x.entropy().unwrap(); + let entropy = |x: Weibull| x.entropy(); test_absolute(1.0, 0.1, -1.302585092994045684018, 1e-15, entropy); test_exact(1.0, 1.0, 1.0, entropy); test_exact(10.0, 10.0, 1.519494098411379574546, entropy); @@ -433,7 +459,7 @@ mod tests { #[test] fn test_skewnewss() { - let skewness = |x: Weibull| x.skewness().unwrap(); + let skewness = |x: Weibull| x.skewness(); test_absolute(1.0, 0.1, 2.0, 1e-13, skewness); test_absolute(1.0, 1.0, 2.0, 1e-13, skewness); test_absolute(10.0, 10.0, -0.63763713390314440916597757156663888653981696212127, 1e-11, skewness); diff --git a/src/lib.rs b/src/lib.rs index d4d80ebc..59942837 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -20,13 +20,13 @@ //! Statrs also comes with a number of useful utility traits for more detailed introspection of distributions. //! ``` //! use statrs::distribution::{Exp, Continuous, ContinuousCDF}; // `cdf` and `pdf` -//! use statrs::statistics::Distribution; // statistical moments and entropy +//! use statrs::statistics::*; // statistical moments and entropy //! //! let n = Exp::new(1.0).unwrap(); -//! assert_eq!(n.mean(), Some(1.0)); -//! assert_eq!(n.variance(), Some(1.0)); -//! assert_eq!(n.entropy(), Some(1.0)); -//! assert_eq!(n.skewness(), Some(2.0)); +//! assert_eq!(n.mean(), 1.0); +//! assert_eq!(n.variance(), 1.0); +//! assert_eq!(n.entropy(), 1.0); +//! assert_eq!(n.skewness(), 2.0); //! assert_eq!(n.cdf(1.0), 0.6321205588285576784045); //! assert_eq!(n.pdf(1.0), 0.3678794411714423215955); //! ``` @@ -36,7 +36,7 @@ //! //! ``` //! use statrs::distribution::FisherSnedecor; -//! use statrs::statistics::Distribution; +//! use statrs::statistics::*; //! //! let n = FisherSnedecor::new(1.0, 1.0).unwrap(); //! assert!(n.variance().is_none()); diff --git a/src/statistics/slice_statistics.rs b/src/statistics/slice_statistics.rs index 824346d8..5e548bcd 100644 --- a/src/statistics/slice_statistics.rs +++ b/src/statistics/slice_statistics.rs @@ -298,71 +298,89 @@ impl + AsRef<[f64]>> Max for Data { } } -impl + AsRef<[f64]>> Distribution for Data { - /// Evaluates the sample mean, an estimate of the population - /// mean. - /// - /// # Remarks - /// - /// Returns `f64::NAN` if data is empty or an entry is `f64::NAN` - /// - /// # Examples - /// - /// ``` - /// #[macro_use] - /// extern crate statrs; - /// - /// use statrs::statistics::Distribution; - /// use statrs::statistics::Data; - /// - /// # fn main() { - /// let x = []; - /// let x = Data::new(x); - /// assert!(x.mean().unwrap().is_nan()); - /// - /// let y = [0.0, f64::NAN, 3.0, -2.0]; - /// let y = Data::new(y); - /// assert!(y.mean().unwrap().is_nan()); - /// - /// let z = [0.0, 3.0, -2.0]; - /// let z = Data::new(z); - /// assert_almost_eq!(z.mean().unwrap(), 1.0 / 3.0, 1e-15); - /// # } - /// ``` - fn mean(&self) -> Option { - Some(Statistics::mean(self.iter())) +/// Evaluates the sample mean, an estimate of the population +/// mean. +/// +/// # Remarks +/// +/// Returns `None` if data is empty or `f64::NAN` if entry is `f64::NAN` +/// +/// # Examples +/// +/// ``` +/// #[macro_use] +/// extern crate statrs; +/// +/// use statrs::statistics::Mean; +/// use statrs::statistics::Data; +/// +/// # fn main() { +/// let x = []; +/// let x = Data::new(x); +/// assert!(x.mean().is_none()); +/// +/// let y = [0.0, f64::NAN, 3.0, -2.0]; +/// let y = Data::new(y); +/// assert!(y.mean().unwrap().is_nan()); +/// +/// let z = [0.0, 3.0, -2.0]; +/// let z = Data::new(z); +/// assert_almost_eq!(z.mean().unwrap(), 1.0 / 3.0, 1e-15); +/// # } +/// ``` +impl Mean for Data +where + D: AsMut<[f64]> + AsRef<[f64]>, +{ + type Mu = Option; + fn mean(&self) -> Self::Mu { + if self.0.as_ref().is_empty() { + return None; + } + + Some(self.0.as_ref().iter().mean()) } +} - /// Estimates the unbiased population variance from the provided samples - /// - /// # Remarks - /// - /// On a dataset of size `N`, `N-1` is used as a normalizer (Bessel's - /// correction). - /// - /// Returns `f64::NAN` if data has less than two entries or if any entry is - /// `f64::NAN` - /// - /// # Examples - /// - /// ``` - /// use statrs::statistics::Distribution; - /// use statrs::statistics::Data; - /// - /// let x = []; - /// let x = Data::new(x); - /// assert!(x.variance().unwrap().is_nan()); - /// - /// let y = [0.0, f64::NAN, 3.0, -2.0]; - /// let y = Data::new(y); - /// assert!(y.variance().unwrap().is_nan()); - /// - /// let z = [0.0, 3.0, -2.0]; - /// let z = Data::new(z); - /// assert_eq!(z.variance().unwrap(), 19.0 / 3.0); - /// ``` - fn variance(&self) -> Option { - Some(Statistics::variance(self.iter())) +/// Estimates the unbiased population variance from the provided samples +/// +/// # Remarks +/// +/// On a dataset of size `N`, `N-1` is used as a normalizer (Bessel's +/// correction). +/// +/// Returns `None` if data has less than two entries or `f64::NAN` if any entry is +/// `f64::NAN` +/// +/// # Examples +/// +/// ``` +/// use statrs::statistics::Variance; +/// use statrs::statistics::Data; +/// +/// let x = []; +/// let x = Data::new(x); +/// assert!(x.variance().is_none()); +/// +/// let y = [0.0, f64::NAN, 3.0, -2.0]; +/// let y = Data::new(y); +/// assert!(y.variance().unwrap().is_nan()); +/// +/// let z = [0.0, 3.0, -2.0]; +/// let z = Data::new(z); +/// assert_eq!(z.variance().unwrap(), 19.0 / 3.0); +/// ``` +impl Variance for Data +where + D: AsMut<[f64]> + AsRef<[f64]>, +{ + type Var = Option; + fn variance(&self) -> Self::Var { + // TODO: rewrite iterator statistics and rely on them instead + if self.0.as_ref().len() < 2 { + return None; + } + Some(self.0.as_ref().iter().variance()) } } diff --git a/src/statistics/traits.rs b/src/statistics/traits.rs index 9140eab4..474dc00d 100644 --- a/src/statistics/traits.rs +++ b/src/statistics/traits.rs @@ -1,4 +1,4 @@ -use ::num_traits::float::Float; +use num_traits::float::Float; /// The `Min` trait specifies than an object has a minimum value pub trait Min { @@ -9,10 +9,11 @@ pub trait Min { /// /// ``` /// use statrs::statistics::Min; - /// use statrs::distribution::Uniform; + /// use statrs::distribution::{Uniform, UniformError}; /// - /// let n = Uniform::new(0.0, 1.0).unwrap(); + /// let n = Uniform::new(0.0, 1.0)?; /// assert_eq!(0.0, n.min()); + /// # Ok::<(), UniformError>(()) /// ``` fn min(&self) -> T; } @@ -26,118 +27,173 @@ pub trait Max { /// /// ``` /// use statrs::statistics::Max; - /// use statrs::distribution::Uniform; + /// use statrs::distribution::{Uniform, UniformError}; /// - /// let n = Uniform::new(0.0, 1.0).unwrap(); + /// let n = Uniform::new(0.0, 1.0)?; /// assert_eq!(1.0, n.max()); + /// # Ok::<(), UniformError>(()) /// ``` fn max(&self) -> T; } -pub trait DiscreteDistribution { - /// Returns the mean, if it exists. - fn mean(&self) -> Option { - None - } - /// Returns the variance, if it exists. - fn variance(&self) -> Option { - None - } - /// Returns the standard deviation, if it exists. - fn std_dev(&self) -> Option { - self.variance().map(|var| var.sqrt()) - } - /// Returns the entropy, if it exists. - fn entropy(&self) -> Option { - None + +/// Exposes an entropy method, uses base e, (not Shannon, which base 2). +pub trait Entropy { + fn entropy(&self) -> T; +} + +/// Trait to express covariance operations as if the implementing type were an operator. +/// +/// For scalars this is variance and scalar +/// +/// ```text +/// Sigma_ij = Cov[X_i, X_j] +/// ``` +pub trait CovarianceMatrix { + /// Type that acts as matrix-like operator + type M; + /// Type that acts as vector-like state + type V; + + /// returns a covariance matrix, M_ij = Sigma_ij + fn dense(&self) -> Self::M; + + /// transforms vector of uncorrelated samples to have correlations from self, (Sigma)^1/2 * vec{v} + #[doc(alias = "matvec")] + fn colorize(&self, other: Self::V) -> Self::V; + + /// transforms vector sampled with correlations from self to uncorrelated, (Sigma)^1/2 \ vec{v} + #[doc(alias = "solve")] + fn whiten(&self, other: Self::V) -> Self::V; + + /// returns the determinant of the covariance matrix, det(Sigma) + fn determinant(&self) -> T; +} + +#[cfg(feature = "nalgebra")] +mod multivariate { + use nalgebra::{Cholesky, Dim, OMatrix, OVector, Scalar}; + + impl super::CovarianceMatrix for OVector + where + T: Scalar + nalgebra::RealField, + D: Dim, + nalgebra::DefaultAllocator: + nalgebra::allocator::Allocator + nalgebra::allocator::Allocator, + { + type M = OMatrix; + type V = Self; + + fn dense(&self) -> Self::M { + OMatrix::from_diagonal(self) + } + + fn colorize(&self, other: Self::V) -> Self::V { + self.clone().map(|x| x.sqrt()).component_mul(&other) + } + + fn whiten(&self, other: Self::V) -> Self::V { + other.component_div(&self.clone().map(|x| x.sqrt())) + } + + fn determinant(&self) -> T { + self.product() + } } - /// Returns the skewness, if it exists. - fn skewness(&self) -> Option { - None + + impl super::CovarianceMatrix for Cholesky + where + T: Scalar + nalgebra::RealField, + D: Dim, + nalgebra::DefaultAllocator: + nalgebra::allocator::Allocator + nalgebra::allocator::Allocator, + { + type M = OMatrix; + type V = OVector; + + fn dense(&self) -> Self::M { + self.l() * self.l().transpose() + } + + fn colorize(&self, other: Self::V) -> Self::V { + self.l().transpose() * other + } + + fn whiten(&self, other: Self::V) -> Self::V { + self.l_dirty().solve_lower_triangular_unchecked(&other) + } + + fn determinant(&self) -> T { + self.determinant() + } } } -pub trait Distribution { - /// Returns the mean, if it exists. - /// - /// # Examples - /// - /// ``` - /// use statrs::statistics::Distribution; - /// use statrs::distribution::Uniform; - /// - /// let n = Uniform::new(0.0, 1.0).unwrap(); - /// assert_eq!(0.5, n.mean().unwrap()); - /// ``` - fn mean(&self) -> Option { - None +impl CovarianceMatrix for T { + type M = T; + type V = T; + fn dense(&self) -> Self::M { + *self } - /// Returns the variance, if it exists. - /// - /// # Examples - /// - /// ``` - /// use statrs::statistics::Distribution; - /// use statrs::distribution::Uniform; - /// - /// let n = Uniform::new(0.0, 1.0).unwrap(); - /// assert_eq!(1.0 / 12.0, n.variance().unwrap()); - /// ``` - fn variance(&self) -> Option { - None + fn colorize(&self, other: Self::V) -> Self::V { + self.sqrt() * other } - /// Returns the standard deviation, if it exists. - /// - /// # Examples - /// - /// ``` - /// use statrs::statistics::Distribution; - /// use statrs::distribution::Uniform; - /// - /// let n = Uniform::new(0.0, 1.0).unwrap(); - /// assert_eq!((1f64 / 12f64).sqrt(), n.std_dev().unwrap()); - /// ``` - fn std_dev(&self) -> Option { - self.variance().map(|var| var.sqrt()) + fn whiten(&self, other: Self::V) -> Self::V { + other / self.sqrt() } - /// Returns the entropy, if it exists. - /// - /// # Examples - /// - /// ``` - /// use statrs::statistics::Distribution; - /// use statrs::distribution::Uniform; - /// - /// let n = Uniform::new(0.0, 1.0).unwrap(); - /// assert_eq!(0.0, n.entropy().unwrap()); - /// ``` - fn entropy(&self) -> Option { - None - } - /// Returns the skewness, if it exists. - /// - /// # Examples - /// - /// ``` - /// use statrs::statistics::Distribution; - /// use statrs::distribution::Uniform; - /// - /// let n = Uniform::new(0.0, 1.0).unwrap(); - /// assert_eq!(0.0, n.skewness().unwrap()); - /// ``` - fn skewness(&self) -> Option { - None + fn determinant(&self) -> T { + *self } } -/// The `Mean` trait implements the calculation of a mean. -// TODO: Clarify the traits of multidimensional distributions -pub trait MeanN { - fn mean(&self) -> Option; +/// Trait to express the mean, i.e. the lowest non-trivial [standardized moment](https://en.wikipedia.org/wiki/Moment_(mathematics)#Standardized_moments) +/// of a distribution's realization implemented on its distribution's type +/// +/// # Note for implementors +/// Associated types should capture semantics of the distribution, e.g. +/// [`Option`] should be used where moments is defined or undefined based on parameter values. +/// [`()`] is used for moments that are never defined +pub trait Mean { + type Mu; + fn mean(&self) -> Self::Mu; +} + +/// Trait to express the variance, the second [standardized moment](https://en.wikipedia.org/wiki/Moment_(mathematics)#Standardized_moments) +/// of a distribution's realization implemented on its distribution's type. +/// +/// Multivariates generalize this to covariance. +/// # Note for implementors +/// Associated types should capture semantics of the distribution, e.g. +/// [`Option`] should be used where moments is defined or undefined based on parameter values. +/// [`()`] is used for moments that are never defined +pub trait Variance: Mean { + type Var; + fn variance(&self) -> Self::Var; } -// TODO: Clarify the traits of multidimensional distributions -pub trait VarianceN { - fn variance(&self) -> Option; +/// Trait to express the skewness, the third [standardized moment](https://en.wikipedia.org/wiki/Moment_(mathematics)#Standardized_moments) +/// of a distribution's realization implemented on its distribution's type. +/// +/// This is the third central moment normalized by the variance to 3/2 power. +/// # Note for implementors +/// Associated types should capture semantics of the distribution, e.g. +/// [`Option`] should be used where moments is defined or undefined based on parameter values. +/// [`()`] is used for moments that are never defined +pub trait Skewness: Variance { + type Skew; + fn skewness(&self) -> Self::Skew; +} + +/// Trait to express the excess kurtosis, the fourth [standardized moment](https://en.wikipedia.org/wiki/Moment_(mathematics)#Standardized_moments) +/// of a distribution's realization implemented on its distribution's type. +/// +/// This is the third central moment normalized by the variance's second power. +/// # Note for implementors +/// Associated types should capture semantics of the distribution, e.g. +/// [`Option`] should be used where moments is defined or undefined based on parameter values. +/// [`()`] is used for moments that are never defined +pub trait ExcessKurtosis: Variance { + type Kurt; + fn excess_kurtosis(&self) -> Self::Kurt; } /// The `Median` trait returns the median of the distribution. @@ -148,10 +204,11 @@ pub trait Median { /// /// ``` /// use statrs::statistics::Median; - /// use statrs::distribution::Uniform; + /// use statrs::distribution::{Uniform, UniformError}; /// - /// let n = Uniform::new(0.0, 1.0).unwrap(); + /// let n = Uniform::new(0.0, 1.0)?; /// assert_eq!(0.5, n.median()); + /// # Ok::<(), UniformError>(()) /// ``` fn median(&self) -> T; } @@ -165,10 +222,59 @@ pub trait Mode { /// /// ``` /// use statrs::statistics::Mode; - /// use statrs::distribution::Uniform; + /// use statrs::distribution::{Uniform, UniformError}; /// - /// let n = Uniform::new(0.0, 1.0).unwrap(); + /// let n = Uniform::new(0.0, 1.0)?; /// assert_eq!(Some(0.5), n.mode()); + /// # Ok::<(), UniformError>(()) /// ``` fn mode(&self) -> T; } + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn covariance_scalar() { + let x = 4.0_f64; + assert_eq!(x.colorize(1.0), 2.0); + assert_eq!(x.whiten(2.0), 1.0); + assert_eq!(x.determinant(), 4.0); + } + + #[cfg(feature = "nalgebra")] + #[test] + fn covariance_vector() { + use approx::assert_relative_eq; + use nalgebra::{vector, OMatrix}; + let v = vector![1.0_f64, 4.0, 9.0]; + assert_relative_eq!(v.dense(), OMatrix::from_diagonal(&v)); + assert_relative_eq!(v.colorize(vector![1.0, 1.0, 1.0]), vector![1.0, 2.0, 3.0]); + assert_relative_eq!(v.whiten(vector![1.0, 2.0, 3.0]), vector![1.0, 1.0, 1.0]); + } + + #[cfg(feature = "nalgebra")] + #[test] + fn covariance_matrix() { + use approx::assert_relative_eq; + use nalgebra::{matrix, vector, Cholesky}; + let m = matrix![5.0_f64, 8.0; 8.0, 13.0]; + let c = Cholesky::new(m).unwrap(); + assert_relative_eq!(c.dense(), m); + + for v in [vector![1.0, 0.0], vector![0.0, 1.0], vector![1.0, 1.0]] { + assert_relative_eq!( + c.colorize(v).norm_squared(), + (v.transpose() * m * v)[(0, 0)] + ); + } + + for v in [vector![1.0, 0.0], vector![0.0, 1.0], vector![1.0, 1.0]] { + assert_relative_eq!( + c.whiten(v).norm_squared(), + (v.transpose() * (m.solve_lower_triangular_unchecked(&v)))[(0, 0)] + ); + } + } +}