diff --git a/src/distribution/levy.rs b/src/distribution/levy.rs index 087ad350..88208e8b 100644 --- a/src/distribution/levy.rs +++ b/src/distribution/levy.rs @@ -4,12 +4,24 @@ use crate::function::erf::{erf, erfc, erfc_inv}; use crate::statistics::*; use core::f64; +/// Implements the [Levy](https://en.wikipedia.org/wiki/L%C3%A9vy_distribution) distribution. +/// +/// # Example +/// +/// ``` +/// use statrs::distribution::{Levy, Continuous}; +/// use statrs::statistics::Distribution; +/// +/// let n = Levy::new(1.0, 1.0).unwrap(); +/// assert_eq!(n.pdf(0.0), 0.0); +/// ``` #[derive(Copy, Clone, PartialEq, Debug)] pub struct Levy { mu: f64, c: f64, } +/// Represents the errors that can occur when creating a [`Levy`]. #[derive(Copy, Clone, PartialEq, Eq, Debug, Hash)] #[non_exhaustive] pub enum LevyError { @@ -32,6 +44,23 @@ impl std::fmt::Display for LevyError { impl std::error::Error for LevyError {} impl Levy { + /// Constructs a new Levy distribution with a location (μ) and dispersion (c) + /// + /// # Errors + /// + /// Returns and error if `mu` is NaN or infinite or if `c` is NaN, infinite or nonpositive + /// + /// # Example + /// + /// ``` + /// use statrs::distribution::Levy; + /// + /// let mut result = Levy::new(0.0, 1.0); + /// assert!(result.is_ok()); + /// + /// result = Levy::new(0.0, 0.0); + /// assert!(result.is_err()); + /// ``` pub fn new(mu: f64, c: f64) -> Result { if mu.is_nan() || mu.is_infinite() { return Err(LevyError::LocationInvalid); @@ -42,10 +71,30 @@ impl Levy { Ok(Levy { mu, c }) } + /// Returns the location (μ) of the Levy distribution + /// + /// # Example + /// + /// ``` + /// use statrs::distribution::Levy; + /// + /// let n = Levy::new(1.0, 1.0).unwrap(); + /// assert_eq!(n.mu(), 1.0); + /// ``` pub fn mu(&self) -> f64 { self.mu } + /// Returns the dispersion (c) of the Levy distribution + /// + /// # Example + /// + /// ``` + /// use statrs::distribution::Levy; + /// + /// let n = Levy::new(1.0, 1.0).unwrap(); + /// assert_eq!(n.c(), 1.0); + /// ``` pub fn c(&self) -> f64 { self.c } @@ -70,6 +119,17 @@ impl ::rand::distributions::Distribution for Levy { } impl ContinuousCDF for Levy { + /// Calculates the cumulative distribution function for the Levy distribution at `x` + /// + /// # Formula + /// + /// ```text + /// 0 if x <= μ + /// erfc(sqrt(c / (2 * (x - μ)))) if x > μ + /// ``` + /// + /// where `μ` is the location, `c` is the dispersion, and `erfc` is the + /// complementary error function. fn cdf(&self, x: f64) -> f64 { if x <= self.mu { 0.0 @@ -80,6 +140,17 @@ impl ContinuousCDF for Levy { } } + /// Calculates the survival function for the Levy distribution at `x` + /// + /// # Formula + /// + /// ```text + /// 1 if x <= μ + /// erf(sqrt(c / (2 * (x - μ)))) if x > μ + /// ``` + /// + /// where `μ` is the location, `c` is the dispersion, and `erf` is the error + /// function. fn sf(&self, x: f64) -> f64 { if x <= self.mu { 1.0 @@ -90,6 +161,21 @@ impl ContinuousCDF for Levy { } } + /// Calculates the inverse cumulative distribution function for the + /// normal distribution at `x`. + /// + /// # Panics + /// + /// If `x < 0.0` or `x > 1.0` + /// + /// # Formula + /// + /// ```text + /// μ + c * (erfc_inv(x)^2)/2 + /// ``` + /// + /// where `μ` is the mean, `σ` is the standard deviation and `erfc_inv` is + /// the inverse of the complementary error function fn inverse_cdf(&self, x: f64) -> f64 { if !(0.0..=1.0).contains(&x) { panic!("x must be in [0, 1]"); @@ -100,26 +186,74 @@ impl ContinuousCDF for Levy { } impl Min for Levy { + /// Returns the minimum value in the domain of the + /// Levy distribution representable by a double precision float + /// + /// # Formula + /// + /// ```text + /// μ + /// ``` fn min(&self) -> f64 { self.mu } } impl Max for Levy { + /// Returns the maximum value in the domain of the + /// Levy distribution representable by a double precision float + /// + /// # Formula + /// + /// ```text + /// f64::INFINITY + /// ``` fn max(&self) -> f64 { f64::INFINITY } } impl Distribution for Levy { + /// Returns the mean of the Levy distribution + /// + /// # Formula + /// + /// ```text + /// f64::INFINITY + /// ``` fn mean(&self) -> Option { Some(f64::INFINITY) } + /// Returns the variance of the Levy distribution + /// + /// # Formula + /// + /// ```text + /// f64::INFINITY + /// ``` fn variance(&self) -> Option { Some(f64::INFINITY) } + /// Returns the standard deviation of the Levy distribution + /// + /// # Formula + /// + /// ```text + /// f64::INFINITY + /// ``` + fn std_dev(&self) -> Option { + Some(f64::INFINITY) + } + + /// Returns the entropy of the Levy distribution + /// + /// # Formula + /// + /// ```text + /// (1 + 3γ + ln(16πc^2))/2 + /// ``` fn entropy(&self) -> Option { /// CONSTANT_PART = 1.5 * EULER_MASCHERONI + 0.5 * (1.0 + LN_PI + 16.0_f64.ln()) const CONSTANT_PART: f64 = 3.32448280139688989720525569282472133636474609375; @@ -128,18 +262,46 @@ impl Distribution for Levy { } impl Median for Levy { + /// Returns the median of the Levy distribution + /// + /// # Formula + /// + /// ```text + /// μ + c/(2 * erfc_inv(0.5)^2) + /// ``` + /// + /// where `μ` is the mean, `c` is the dispersion and `erfc_inv` is + /// the inverse of the complementary error function. fn median(&self) -> f64 { self.mu + self.c * 0.5 * erfc_inv(0.5).powf(-2.0) } } impl Mode> for Levy { + /// Returns the mode of the Levy distribution + /// + /// # Formula + /// + /// ```text + /// μ + c/3 + /// ``` + /// + /// where `μ` is the mean and `c` is the dispersion. fn mode(&self) -> Option { Some(self.mu + self.c / 3.0) } } impl Continuous for Levy { + /// Calculates the probability density function for the Levy distribution at `x` + /// + /// # Formula + /// + /// ```text + /// (sqrt(c / 2 * π)) * e^(-1/2 * (c / (x - μ))) * (1 / (x - μ)^(3/2)) + /// ``` + /// + /// where `μ` is the mean and `c` is the dispersion. fn pdf(&self, x: f64) -> f64 { if x <= self.mu { 0.0 @@ -149,6 +311,15 @@ impl Continuous for Levy { } } + /// Calculates the log probability density function for the Levy distribution at `x` + /// + /// # Formula + /// + /// ```text + /// 1/2 * (ln(c) - ln(2 * π) - (c / (x - μ))) - 3/2 * ln(x - μ) + /// ``` + /// + /// where `μ` is the mean and `σ` is the standard deviation fn ln_pdf(&self, x: f64) -> f64 { use crate::consts::LN_SQRT_2PI;