statrs/distribution/
pareto.rs

1use crate::distribution::{Continuous, ContinuousCDF};
2use crate::statistics::*;
3use std::f64;
4
5/// Implements the [Pareto](https://en.wikipedia.org/wiki/Pareto_distribution)
6/// distribution
7///
8/// # Examples
9///
10/// ```
11/// use statrs::distribution::{Pareto, Continuous};
12/// use statrs::statistics::Distribution;
13/// use statrs::prec;
14///
15/// let p = Pareto::new(1.0, 2.0).unwrap();
16/// assert_eq!(p.mean().unwrap(), 2.0);
17/// assert!(prec::almost_eq(p.pdf(2.0), 0.25, 1e-15));
18/// ```
19#[derive(Copy, Clone, PartialEq, Debug)]
20pub struct Pareto {
21    scale: f64,
22    shape: f64,
23}
24
25/// Represents the errors that can occur when creating a [`Pareto`].
26#[derive(Copy, Clone, PartialEq, Eq, Debug, Hash)]
27#[non_exhaustive]
28pub enum ParetoError {
29    /// The scale is NaN, zero or less than zero.
30    ScaleInvalid,
31
32    /// The shape is NaN, zero or less than zero.
33    ShapeInvalid,
34}
35
36impl std::fmt::Display for ParetoError {
37    #[cfg_attr(coverage_nightly, coverage(off))]
38    fn fmt(&self, f: &mut std::fmt::Formatter) -> std::fmt::Result {
39        match self {
40            ParetoError::ScaleInvalid => write!(f, "Scale is NaN, zero, or less than zero"),
41            ParetoError::ShapeInvalid => write!(f, "Shape is NaN, zero, or less than zero"),
42        }
43    }
44}
45
46impl std::error::Error for ParetoError {}
47
48impl Pareto {
49    /// Constructs a new Pareto distribution with scale `scale`, and `shape`
50    /// shape.
51    ///
52    /// # Errors
53    ///
54    /// Returns an error if any of `scale` or `shape` are `NaN`.
55    /// Returns an error if `scale <= 0.0` or `shape <= 0.0`
56    ///
57    /// # Examples
58    ///
59    /// ```
60    /// use statrs::distribution::Pareto;
61    ///
62    /// let mut result = Pareto::new(1.0, 2.0);
63    /// assert!(result.is_ok());
64    ///
65    /// result = Pareto::new(0.0, 0.0);
66    /// assert!(result.is_err());
67    /// ```
68    pub fn new(scale: f64, shape: f64) -> Result<Pareto, ParetoError> {
69        if scale.is_nan() || scale <= 0.0 {
70            return Err(ParetoError::ScaleInvalid);
71        }
72
73        if shape.is_nan() || shape <= 0.0 {
74            return Err(ParetoError::ShapeInvalid);
75        }
76
77        Ok(Pareto { scale, shape })
78    }
79
80    /// Returns the scale of the Pareto distribution
81    ///
82    /// # Examples
83    ///
84    /// ```
85    /// use statrs::distribution::Pareto;
86    ///
87    /// let n = Pareto::new(1.0, 2.0).unwrap();
88    /// assert_eq!(n.scale(), 1.0);
89    /// ```
90    pub fn scale(&self) -> f64 {
91        self.scale
92    }
93
94    /// Returns the shape of the Pareto distribution
95    ///
96    /// # Examples
97    ///
98    /// ```
99    /// use statrs::distribution::Pareto;
100    ///
101    /// let n = Pareto::new(1.0, 2.0).unwrap();
102    /// assert_eq!(n.shape(), 2.0);
103    /// ```
104    pub fn shape(&self) -> f64 {
105        self.shape
106    }
107}
108
109impl std::fmt::Display for Pareto {
110    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
111        write!(f, "Pareto({},{})", self.scale, self.shape)
112    }
113}
114
115#[cfg(feature = "rand")]
116#[cfg_attr(docsrs, doc(cfg(feature = "rand")))]
117impl ::rand::distributions::Distribution<f64> for Pareto {
118    fn sample<R: ::rand::Rng + ?Sized>(&self, rng: &mut R) -> f64 {
119        use rand::distributions::OpenClosed01;
120
121        // Inverse transform sampling
122        let u: f64 = rng.sample(OpenClosed01);
123        self.scale * u.powf(-1.0 / self.shape)
124    }
125}
126
127impl ContinuousCDF<f64, f64> for Pareto {
128    /// Calculates the cumulative distribution function for the Pareto
129    /// distribution at `x`
130    ///
131    /// # Formula
132    ///
133    /// ```text
134    /// if x < x_m {
135    ///     0
136    /// } else {
137    ///     1 - (x_m/x)^α
138    /// }
139    /// ```
140    ///
141    /// where `x_m` is the scale and `α` is the shape
142    fn cdf(&self, x: f64) -> f64 {
143        if x < self.scale {
144            0.0
145        } else {
146            1.0 - (self.scale / x).powf(self.shape)
147        }
148    }
149
150    /// Calculates the survival function for the Pareto
151    /// distribution at `x`
152    ///
153    /// # Formula
154    ///
155    /// ```text
156    /// if x < x_m {
157    ///     1
158    /// } else {
159    ///     (x_m/x)^α
160    /// }
161    /// ```
162    ///
163    /// where `x_m` is the scale and `α` is the shape
164    fn sf(&self, x: f64) -> f64 {
165        if x < self.scale {
166            1.0
167        } else {
168            (self.scale / x).powf(self.shape)
169        }
170    }
171
172    /// Calculates the inverse cumulative distribution function for the Pareto
173    /// distribution at `x`
174    ///
175    /// # Formula
176    ///
177    /// ```text
178    /// x_m / (1 - x)^(1 / α)
179    /// ```
180    ///
181    /// where `x_m` is the scale and `α` is the shape
182    fn inverse_cdf(&self, p: f64) -> f64 {
183        if !(0.0..=1.0).contains(&p) {
184            panic!("x must be in [0, 1]");
185        } else {
186            self.scale * (1.0 - p).powf(-1.0 / self.shape)
187        }
188    }
189}
190
191impl Min<f64> for Pareto {
192    /// Returns the minimum value in the domain of the Pareto distribution
193    /// representable by a double precision float
194    ///
195    /// # Formula
196    ///
197    /// ```text
198    /// x_m
199    /// ```
200    ///
201    /// where `x_m` is the scale
202    fn min(&self) -> f64 {
203        self.scale
204    }
205}
206
207impl Max<f64> for Pareto {
208    /// Returns the maximum value in the domain of the Pareto distribution
209    /// representable by a double precision float
210    ///
211    /// # Formula
212    ///
213    /// ```text
214    /// f64::INFINITY
215    /// ```
216    fn max(&self) -> f64 {
217        f64::INFINITY
218    }
219}
220
221impl Distribution<f64> for Pareto {
222    /// Returns the mean of the Pareto distribution
223    ///
224    /// # Formula
225    ///
226    /// ```text
227    /// if α <= 1 {
228    ///     f64::INFINITY
229    /// } else {
230    ///     (α * x_m)/(α - 1)
231    /// }
232    /// ```
233    ///
234    /// where `x_m` is the scale and `α` is the shape
235    fn mean(&self) -> Option<f64> {
236        if self.shape <= 1.0 {
237            None
238        } else {
239            Some((self.shape * self.scale) / (self.shape - 1.0))
240        }
241    }
242
243    /// Returns the variance of the Pareto distribution
244    ///
245    /// # Formula
246    ///
247    /// ```text
248    /// if α <= 2 {
249    ///     f64::INFINITY
250    /// } else {
251    ///     (x_m/(α - 1))^2 * (α/(α - 2))
252    /// }
253    /// ```
254    ///
255    /// where `x_m` is the scale and `α` is the shape
256    fn variance(&self) -> Option<f64> {
257        if self.shape <= 2.0 {
258            None
259        } else {
260            let a = self.scale / (self.shape - 1.0); // just a temporary variable
261            Some(a * a * self.shape / (self.shape - 2.0))
262        }
263    }
264
265    /// Returns the entropy for the Pareto distribution
266    ///
267    /// # Formula
268    ///
269    /// ```text
270    /// ln(α/x_m) - 1/α - 1
271    /// ```
272    ///
273    /// where `x_m` is the scale and `α` is the shape
274    fn entropy(&self) -> Option<f64> {
275        Some(self.shape.ln() - self.scale.ln() - (1.0 / self.shape) - 1.0)
276    }
277
278    /// Returns the skewness of the Pareto distribution
279    ///
280    /// # Panics
281    ///
282    /// If `α <= 3.0`
283    ///
284    /// where `α` is the shape
285    ///
286    /// # Formula
287    ///
288    /// ```text
289    ///     (2*(α + 1)/(α - 3))*sqrt((α - 2)/α)
290    /// ```
291    ///
292    /// where `α` is the shape
293    fn skewness(&self) -> Option<f64> {
294        if self.shape <= 3.0 {
295            None
296        } else {
297            Some(
298                (2.0 * (self.shape + 1.0) / (self.shape - 3.0))
299                    * ((self.shape - 2.0) / self.shape).sqrt(),
300            )
301        }
302    }
303}
304
305impl Median<f64> for Pareto {
306    /// Returns the median of the Pareto distribution
307    ///
308    /// # Formula
309    ///
310    /// ```text
311    /// x_m*2^(1/α)
312    /// ```
313    ///
314    /// where `x_m` is the scale and `α` is the shape
315    fn median(&self) -> f64 {
316        self.scale * (2f64.powf(1.0 / self.shape))
317    }
318}
319
320impl Mode<Option<f64>> for Pareto {
321    /// Returns the mode of the Pareto distribution
322    ///
323    /// # Formula
324    ///
325    /// ```text
326    /// x_m
327    /// ```
328    ///
329    /// where `x_m` is the scale
330    fn mode(&self) -> Option<f64> {
331        Some(self.scale)
332    }
333}
334
335impl Continuous<f64, f64> for Pareto {
336    /// Calculates the probability density function for the Pareto distribution
337    /// at `x`
338    ///
339    /// # Formula
340    ///
341    /// ```text
342    /// if x < x_m {
343    ///     0
344    /// } else {
345    ///     (α * x_m^α)/(x^(α + 1))
346    /// }
347    /// ```
348    ///
349    /// where `x_m` is the scale and `α` is the shape
350    fn pdf(&self, x: f64) -> f64 {
351        if x < self.scale {
352            0.0
353        } else {
354            (self.shape * self.scale.powf(self.shape)) / x.powf(self.shape + 1.0)
355        }
356    }
357
358    /// Calculates the log probability density function for the Pareto
359    /// distribution at `x`
360    ///
361    /// # Formula
362    ///
363    /// ```text
364    /// if x < x_m {
365    ///     f64::NEG_INFINITY
366    /// } else {
367    ///     ln(α) + α*ln(x_m) - (α + 1)*ln(x)
368    /// }
369    /// ```
370    ///
371    /// where `x_m` is the scale and `α` is the shape
372    fn ln_pdf(&self, x: f64) -> f64 {
373        if x < self.scale {
374            f64::NEG_INFINITY
375        } else {
376            self.shape.ln() + self.shape * self.scale.ln() - (self.shape + 1.0) * x.ln()
377        }
378    }
379}
380
381#[rustfmt::skip]
382#[cfg(test)]
383mod tests {
384    use super::*;
385    use crate::distribution::internal::*;
386    use crate::testing_boiler;
387
388    testing_boiler!(scale: f64, shape: f64; Pareto; ParetoError);
389
390    #[test]
391    fn test_create() {
392        create_ok(10.0, 0.1);
393        create_ok(5.0, 1.0);
394        create_ok(0.1, 10.0);
395        create_ok(10.0, 100.0);
396        create_ok(1.0, f64::INFINITY);
397        create_ok(f64::INFINITY, f64::INFINITY);
398    }
399
400    #[test]
401    fn test_bad_create() {
402        test_create_err(1.0, -1.0, ParetoError::ShapeInvalid);
403        test_create_err(-1.0, 1.0, ParetoError::ScaleInvalid);
404        create_err(0.0, 0.0);
405        create_err(-1.0, -1.0);
406        create_err(f64::NAN, 1.0);
407        create_err(1.0, f64::NAN);
408        create_err(f64::NAN, f64::NAN);
409    }
410
411    #[test]
412    fn test_variance() {
413        let variance = |x: Pareto| x.variance().unwrap();
414        test_exact(1.0, 3.0, 0.75, variance);
415        test_absolute(10.0, 10.0, 125.0 / 81.0, 1e-13, variance);
416    }
417
418    #[test]
419    fn test_variance_degen() {
420        test_none(1.0, 1.0, |dist| dist.variance()); // shape <= 2.0
421    }
422
423    #[test]
424    fn test_entropy() {
425        let entropy = |x: Pareto| x.entropy().unwrap();
426        test_exact(0.1, 0.1, -11.0, entropy);
427        test_exact(1.0, 1.0, -2.0, entropy);
428        test_exact(10.0, 10.0, -1.1, entropy);
429        test_exact(3.0, 1.0, -2.0 - 3f64.ln(), entropy);
430        test_exact(1.0, 3.0, -4.0/3.0 + 3f64.ln(), entropy);
431    }
432
433    #[test]
434    fn test_skewness() {
435        let skewness = |x: Pareto| x.skewness().unwrap();
436        test_exact(1.0, 4.0, 5.0*2f64.sqrt(), skewness);
437        test_exact(1.0, 100.0, (707.0/485.0)*2f64.sqrt(), skewness);
438    }
439
440    #[test]
441    fn test_skewness_invalid_shape() {
442        test_none(1.0, 3.0, |dist| dist.skewness());
443    }
444
445    #[test]
446    fn test_mode() {
447        let mode = |x: Pareto| x.mode().unwrap();
448        test_exact(0.1, 1.0, 0.1, mode);
449        test_exact(2.0, 1.0, 2.0, mode);
450        test_exact(10.0, f64::INFINITY, 10.0, mode);
451        test_exact(f64::INFINITY, 1.0, f64::INFINITY, mode);
452    }
453
454    #[test]
455    fn test_median() {
456        let median = |x: Pareto| x.median();
457        test_exact(0.1, 0.1, 102.4, median);
458        test_exact(1.0, 1.0, 2.0, median);
459        test_exact(10.0, 10.0, 10.0*2f64.powf(0.1), median);
460        test_exact(3.0, 0.5, 12.0, median);
461        test_exact(10.0, f64::INFINITY, 10.0, median);
462    }
463
464    #[test]
465    fn test_min_max() {
466        let min = |x: Pareto| x.min();
467        let max = |x: Pareto| x.max();
468        test_exact(0.2, f64::INFINITY, 0.2, min);
469        test_exact(10.0, f64::INFINITY, 10.0, min);
470        test_exact(f64::INFINITY, 1.0, f64::INFINITY, min);
471        test_exact(1.0, 0.1, f64::INFINITY, max);
472        test_exact(3.0, 10.0, f64::INFINITY, max);
473    }
474
475    #[test]
476    fn test_pdf() {
477        let pdf = |arg: f64| move |x: Pareto| x.pdf(arg);
478        test_exact(1.0, 1.0, 0.0, pdf(0.1));
479        test_exact(1.0, 1.0, 1.0, pdf(1.0));
480        test_exact(1.0, 1.0, 4.0/9.0, pdf(1.5));
481        test_exact(1.0, 1.0, 1.0/25.0, pdf(5.0));
482        test_exact(1.0, 1.0, 1.0/2500.0, pdf(50.0));
483        test_exact(1.0, 4.0, 4.0, pdf(1.0));
484        test_exact(1.0, 4.0, 128.0/243.0, pdf(1.5));
485        test_exact(1.0, 4.0, 1.0/78125000.0, pdf(50.0));
486        test_exact(3.0, 2.0, 2.0/3.0, pdf(3.0));
487        test_exact(3.0, 2.0, 18.0/125.0, pdf(5.0));
488        test_absolute(25.0, 100.0, 1.5777218104420236e-30, 1e-50, pdf(50.0));
489        test_absolute(100.0, 25.0, 6.6003546737276816e-6, 1e-16, pdf(150.0));
490        test_exact(1.0, 2.0, 0.0, pdf(f64::INFINITY));
491    }
492
493    #[test]
494    fn test_ln_pdf() {
495        let ln_pdf = |arg: f64| move |x: Pareto| x.ln_pdf(arg);
496        test_exact(1.0, 1.0, f64::NEG_INFINITY, ln_pdf(0.1));
497        test_exact(1.0, 1.0, 0.0, ln_pdf(1.0));
498        test_absolute(1.0, 1.0, 4f64.ln() - 9f64.ln(), 1e-14, ln_pdf(1.5));
499        test_absolute(1.0, 1.0, -(25f64.ln()), 1e-14, ln_pdf(5.0));
500        test_absolute(1.0, 1.0, -(2500f64.ln()), 1e-14, ln_pdf(50.0));
501        test_absolute(1.0, 4.0, 4f64.ln(), 1e-14, ln_pdf(1.0));
502        test_absolute(1.0, 4.0, 128f64.ln() - 243f64.ln(), 1e-14, ln_pdf(1.5));
503        test_absolute(1.0, 4.0, -(78125000f64.ln()), 1e-14, ln_pdf(50.0));
504        test_absolute(3.0, 2.0, 2f64.ln() - 3f64.ln(), 1e-14, ln_pdf(3.0));
505        test_absolute(3.0, 2.0, 18f64.ln() - 125f64.ln(), 1e-14, ln_pdf(5.0));
506        test_absolute(25.0, 100.0, 1.5777218104420236e-30f64.ln(), 1e-12, ln_pdf(50.0));
507        test_absolute(100.0, 25.0, 6.6003546737276816e-6f64.ln(), 1e-12, ln_pdf(150.0));
508        test_exact(1.0, 2.0, f64::NEG_INFINITY, ln_pdf(f64::INFINITY));
509    }
510
511    #[test]
512    fn test_cdf() {
513        let cdf = |arg: f64| move |x: Pareto| x.cdf(arg);
514        test_exact(0.1, 0.1, 0.0, cdf(0.1));
515        test_exact(1.0, 1.0, 0.0, cdf(1.0));
516        test_exact(5.0, 5.0, 0.0, cdf(2.0));
517        test_exact(7.0, 7.0, 0.9176457, cdf(10.0));
518        test_exact(10.0, 10.0, 50700551.0/60466176.0, cdf(12.0));
519        test_exact(5.0, 1.0, 0.5, cdf(10.0));
520        test_exact(3.0, 10.0, 1023.0/1024.0, cdf(6.0));
521        test_exact(1.0, 1.0, 1.0, cdf(f64::INFINITY));
522    }
523
524    #[test]
525    fn test_sf() {
526        let sf = |arg: f64| move |x: Pareto| x.sf(arg);
527        test_exact(0.1, 0.1, 1.0, sf(0.1));
528        test_exact(1.0, 1.0, 1.0, sf(1.0));
529        test_exact(5.0, 5.0, 1.0, sf(2.0));
530        test_absolute(7.0, 7.0, 0.08235429999999999, 1e-14, sf(10.0));
531        test_absolute(10.0, 10.0, 0.16150558288984573, 1e-14, sf(12.0));
532        test_exact(5.0, 1.0, 0.5, sf(10.0));
533        test_absolute(3.0, 10.0, 0.0009765625, 1e-14, sf(6.0));
534        test_exact(1.0, 1.0, 0.0, sf(f64::INFINITY));
535    }
536
537    #[test]
538    fn test_inverse_cdf() {
539        let func = |arg: f64| move |x: Pareto| x.inverse_cdf(x.cdf(arg));
540        test_exact(0.1, 0.1, 0.1, func(0.1));
541        test_exact(1.0, 1.0, 1.0, func(1.0));
542        test_exact(7.0, 7.0, 10.0, func(10.0));
543        test_exact(10.0, 10.0, 12.0, func(12.0));
544        test_exact(5.0, 1.0, 10.0, func(10.0));
545        test_exact(3.0, 10.0, 6.0, func(6.0));
546    }
547
548    #[test]
549    fn test_continuous() {
550        test::check_continuous_distribution(&create_ok(1.0, 10.0), 1.0, 10.0);
551        test::check_continuous_distribution(&create_ok(0.1, 2.0), 0.1, 100.0);
552    }
553}
statrs/distribution/pareto.rs

statrs/distribution/
pareto.rs
