statrs/distribution/
pareto.rs

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

statrs/distribution/
pareto.rs
