Ideal and robust soliton distribution

src/distribution/ideal_soliton.rs

@@ -0,0 +1,248 @@
+use crate::distribution::Soliton;
+use crate::statistics::*;
+use crate::{Result, StatsError};
+use rand::distributions::Distribution;
+use rand::Rng;
+use std::f64;
+
+/// Implements the [Discrete
+/// Uniform](https://en.wikipedia.org/wiki/Discrete_uniform_distribution)
+/// distribution and an related ideal soliton of the discrete uniform distribuiton
+///
+/// # Examples
+///
+/// ```
+/// use statrs::distribution::{DiscreteUniform, Discrete};
+/// use statrs::statistics::Mean;
+/// use statrs::distribution::IdealSoliton;
+///
+/// let sol = IdealSoliton::new(5).unwrap();
+/// let n = DiscreteUniform::new(0, 5).unwrap();
+/// assert_eq!(n.mean(), 2.5);
+/// assert_eq!(n.pmf(3), 1.0 / 6.0);
+/// assert_eq!(sol.soliton(1), 1.0/5.0);
+/// ```
+
+#[derive(Debug, Clone, PartialEq)]
+pub struct IdealSoliton {
+    min: i64,
+    max: i64,
+}
+
+impl IdealSoliton {
+    /// Constructs a new discrete uniform distribution with a minimum value
+    /// of `min` and a maximum value of `max`.
+    ///
+    /// Additionally construct the ideal soliton of the same max value of `max`
+    /// falling between (1, max)
+    ///
+    /// # Errors
+    ///
+    /// Returns an error if `max < min`
+    ///
+    /// # Examples
+    ///
+    /// ```
+    /// use statrs::distribution::DiscreteUniform;
+    /// use statrs::distribution::IdealSoliton;
+    ///
+    /// let mut sol = IdealSoliton::new(5);
+    /// let mut result = DiscreteUniform::new(0, 5);
+    /// assert!(result.is_ok());
+    /// assert!(sol.is_ok());
+    ///
+    /// result = DiscreteUniform::new(5, 0);
+    /// sol = IdealSoliton::new(-1);
+    /// assert!(result.is_err());
+    /// assert!(sol.is_err());
+    /// ```
+    pub fn new(max: i64) -> Result<IdealSoliton> {
+        if max < 1 {
+            Err(StatsError::BadParams)
+        } else {
+            Ok(IdealSoliton { min: 1, max })
+        }
+    }
+}
+
+impl Distribution<f64> for IdealSoliton {
+    fn sample<R: Rng + ?Sized>(&self, r: &mut R) -> f64 {
+        r.gen_range(0, 1) as f64
+    }
+}
+
+impl Min<i64> for IdealSoliton {
+    /// Returns the minimum value in the domain of the discrete uniform
+    /// distribution
+    ///
+    /// # Remarks
+    ///
+    /// This is the same value as the minimum passed into the constructor
+    fn min(&self) -> i64 {
+        self.min
+    }
+}
+
+impl Max<i64> for IdealSoliton {
+    /// Returns the maximum value in the domain of the discrete uniform
+    /// distribution
+    ///
+    /// # Remarks
+    ///
+    /// This is the same value as the maximum passed into the constructor
+    fn max(&self) -> i64 {
+        self.max
+    }
+}
+
+impl Mean<f64> for IdealSoliton {
+    /// Returns the mean of the discrete uniform distribution
+    ///
+    /// # Formula
+    ///
+    /// ```ignore
+    /// (min + max) / 2
+    /// ```
+    fn mean(&self) -> f64 {
+        (self.min + self.max) as f64 / 2.0
+    }
+}
+
+impl Variance<f64> for IdealSoliton {
+    /// Returns the variance of the discrete uniform distribution
+    ///
+    /// # Formula
+    ///
+    /// ```ignore
+    /// ((max - min + 1)^2 - 1) / 12
+    /// ```
+    fn variance(&self) -> f64 {
+        let diff = (self.max - self.min) as f64;
+        ((diff + 1.0) * (diff + 1.0) - 1.0) / 12.0
+    }
+
+    /// Returns the standard deviation of the discrete uniform distribution
+    ///
+    /// # Formula
+    ///
+    /// ```ignore
+    /// sqrt(((max - min + 1)^2 - 1) / 12)
+    /// ```
+    fn std_dev(&self) -> f64 {
+        self.variance().sqrt()
+    }
+}
+
+impl Soliton<i64, f64> for IdealSoliton {
+    /// Calculates the ideal soliton for the
+    /// discrete uniform distribution at `x`
+    ///
+    /// # Remarks
+    ///
+    /// Returns `0.0` if `x` is not in `[min, max]`
+    ///
+    /// # Formula
+    ///
+    /// ```ignore
+    /// p(1) = 1 / (max)
+    /// p(x) = 1/(x(x-1))
+    /// ```
+    fn soliton(&self, x: i64) -> f64 {
+        if x > 1 && x < self.max {
+            1.0 / ((x as f64) * (x as f64 - 1.0))
+        } else if x == 1 {
+            1.0 / self.max as f64
+        } else {
+            // Point must be in range (0, limit]
+            0.0
+        }
+    }
+
+    fn normalization_factor(&self) -> f64 {
+        0.0
+    }
+
+    fn additive_probability(&self, _x: i64) -> f64 {
+        0.0
+    }
+}
+
+#[cfg_attr(rustfmt, rustfmt_skip)]
+#[cfg(test)]
+mod test {
+    use std::fmt::Debug;
+    use std::f64;
+    use crate::statistics::*;
+    use crate::distribution::IdealSoliton;
+
+    fn try_create(max: i64) -> IdealSoliton {
+        let n = IdealSoliton::new(max);
+        assert!(n.is_ok());
+        n.unwrap()
+    }
+
+    fn create_case(max: i64) {
+        let n = try_create(max);
+        assert_eq!(1, n.min());
+        assert_eq!(max, n.max());
+    }
+
+    fn bad_create_case(max: i64) {
+        let n = IdealSoliton::new(max);
+        assert!(n.is_err());
+    }
+
+    fn get_value<T, F>(max: i64, eval: F) -> T
+        where T: PartialEq + Debug,
+              F: Fn(IdealSoliton) -> T
+    {
+        let n = try_create(max);
+        eval(n)
+    }
+
+    fn test_case<T, F>(max: i64, expected: T, eval: F)
+        where T: PartialEq + Debug,
+              F: Fn(IdealSoliton) -> T
+    {
+        let x = get_value(max, eval);
+        assert_eq!(expected, x);
+    }
+
+    fn test_case_greater<T, F>(max: i64, expected: T, eval: F)
+        where T: PartialEq + Debug + Into<f64>,
+              F: Fn(IdealSoliton) -> T
+    {
+        let sol = get_value(max, eval);
+        let a: f64 = sol.into();
+        let b = expected.into();
+        assert!(a > b, "{} greater than {}", a, b);
+    }
+
+    #[test]
+    fn test_create() {
+        create_case(10);
+        create_case(4);
+        create_case(20);
+    }
+
+    #[test]
+    fn test_bad_create() {
+        bad_create_case(-2);
+        bad_create_case(0);
+    }
+
+    #[test]
+    fn test_mean() {
+        test_case_greater(10, 0.9, |x| x.mean());
+    }
+
+    #[test]
+    fn test_variance() {
+        test_case(10, 8.25, |x| x.variance());
+    }
+
+    #[test]
+    fn test_std_dev() {
+        test_case(10, (8.25f64).sqrt(), |x| x.std_dev());
+    }
+}

src/distribution/mod.rs

@@ -30,6 +30,8 @@ pub use self::students_t::StudentsT;
 pub use self::triangular::Triangular;
 pub use self::uniform::Uniform;
 pub use self::weibull::Weibull;
+pub use self::robust_soliton::RobustSoliton;
+pub use self::ideal_soliton::IdealSoliton;
 use crate::statistics::{Max, Min};
 
 mod bernoulli;
@@ -48,6 +50,7 @@ mod fisher_snedecor;
 mod gamma;
 mod geometric;
 mod hypergeometric;
+mod ideal_soliton;
 mod internal;
 mod inverse_gamma;
 mod log_normal;
@@ -57,6 +60,7 @@ mod negative_binomial;
 mod normal;
 mod pareto;
 mod poisson;
+mod robust_soliton;
 mod students_t;
 mod triangular;
 mod uniform;
@@ -269,3 +273,24 @@ pub trait CheckedDiscrete<T, K> {
     /// ```
     fn checked_ln_pmf(&self, x: T) -> Result<K>;
 }
+
+/// The 'IdealSoliton' trait provides an interface for interacting
+/// with discrete statistical distributions from integers 1..N with
+/// N as the single parameter for the distribution
+pub trait Soliton<T, K> {
+    /// Returns the probability mass function calculated at `x` for
+    /// a given distribution
+    ///
+    ///  # Examples
+    ///
+    ///  ```
+    ///  use statrs::distribution::{IdealSoliton, Soliton};
+    ///  use statrs::prec;
+    ///
+    ///  let n = IdealSoliton::new(5).unwrap();
+    ///  assert_eq!(n.soliton(1), 0.2);
+    ///  ```
+    fn soliton(&self, x: T) -> K;
+    fn normalization_factor(&self) -> K;
+    fn additive_probability(&self, x: T) -> K;
+}