compute/distributions/
pareto.rs

1use crate::distributions::*;
2
3/// Implements the [Pareto](https://en.wikipedia.org/wiki/Pareto_distribution) distribution.
4#[derive(Debug, Clone, Copy)]
5pub struct Pareto {
6    /// Shape parameter α.
7    alpha: f64,
8    /// Parameter which controls the minimum value of the distribution.
9    minval: f64,
10}
11
12impl Pareto {
13    /// Create a new Pareto distribution with shape `alpha` and minimum value `minval`.
14    ///
15    /// # Errors
16    /// Panics if `alpha <= 0` or `minval <= 0`.
17    pub fn new(alpha: f64, minval: f64) -> Self {
18        if alpha <= 0. || minval <= 0. {
19            panic!("Both alpha and beta must be positive.");
20        }
21        Pareto { alpha, minval }
22    }
23    pub fn set_alpha(&mut self, alpha: f64) -> &mut Self {
24        if alpha <= 0. {
25            panic!("Alpha must be positive.");
26        }
27        self.alpha = alpha;
28        self
29    }
30    pub fn set_minval(&mut self, minval: f64) -> &mut Self {
31        if minval <= 0. {
32            panic!("minval must be positive.");
33        }
34        self.minval = minval;
35        self
36    }
37}
38
39impl Default for Pareto {
40    fn default() -> Self {
41        Self::new(1., 1.)
42    }
43}
44
45impl Distribution for Pareto {
46    type Output = f64;
47    /// Samples from the given Pareto distribution using inverse transform sampling.
48    fn sample(&self) -> f64 {
49        let u = alea::f64();
50        self.minval / u.powf(1. / self.alpha)
51    }
52}
53
54impl Distribution1D for Pareto {
55    fn update(&mut self, params: &[f64]) {
56        assert!(params.len() == 2);
57        self.set_alpha(params[0]).set_minval(params[1]);
58    }
59}
60
61impl Continuous for Pareto {
62    type PDFType = f64;
63    /// Calculates the probability density function for the given Pareto function at `x`.
64    ///
65    /// # Remarks
66    /// This returns 0 if `x < minval`
67    fn pdf(&self, x: f64) -> f64 {
68        if x < self.minval {
69            return 0.;
70        }
71        self.alpha * self.minval.powf(self.alpha) / x.powf(self.alpha - 1.)
72    }
73}
74
75impl Mean for Pareto {
76    type MeanType = f64;
77    /// Calculates the mean of the Pareto distribution.
78    fn mean(&self) -> f64 {
79        if self.alpha <= 1. {
80            f64::INFINITY
81        } else {
82            self.alpha * self.minval / (self.alpha - 1.)
83        }
84    }
85}
86
87impl Variance for Pareto {
88    type VarianceType = f64;
89    /// Calculates the variance of the Pareto distribution.
90    fn var(&self) -> f64 {
91        if self.alpha <= 2. {
92            f64::INFINITY
93        } else {
94            self.minval.powi(2) * self.alpha / ((self.alpha - 1.).powi(2) * (self.alpha - 2.))
95        }
96    }
97}
98
99#[cfg(test)]
100mod tests {
101    use super::*;
102    use crate::statistics::{mean, var};
103    use approx_eq::assert_approx_eq;
104
105    #[test]
106    fn test_moments() {
107        let dist = Pareto::new(4., 4.);
108        let data = dist.sample_n(1e6 as usize);
109        assert_approx_eq!(dist.mean(), mean(&data), 0.05);
110        assert_approx_eq!(dist.var(), var(&data), 0.05);
111    }
112}
compute/distributions/pareto.rs

compute/distributions/
pareto.rs
