rand_distr/
pareto.rs

1// Copyright 2018 Developers of the Rand project.
2//
3// Licensed under the Apache License, Version 2.0 <LICENSE-APACHE or
4// https://www.apache.org/licenses/LICENSE-2.0> or the MIT license
5// <LICENSE-MIT or https://opensource.org/licenses/MIT>, at your
6// option. This file may not be copied, modified, or distributed
7// except according to those terms.
8
9//! The Pareto distribution `Pareto(xₘ, α)`.
10
11use crate::{Distribution, OpenClosed01};
12use core::fmt;
13use num_traits::Float;
14use rand::Rng;
15
16/// The [Pareto distribution](https://en.wikipedia.org/wiki/Pareto_distribution) `Pareto(xₘ, α)`.
17///
18/// The Pareto distribution is a continuous probability distribution with
19/// scale parameter `xₘ` ( or `k`) and shape parameter `α`.
20///
21/// # Plot
22///
23/// The following plot shows the Pareto distribution with various values of
24/// `xₘ` and `α`.
25/// Note how the shape parameter `α` corresponds to the height of the jump
26/// in density at `x = xₘ`, and to the rate of decay in the tail.
27///
28/// ![Pareto distribution](https://raw.githubusercontent.com/rust-random/charts/main/charts/pareto.svg)
29///
30/// # Example
31/// ```
32/// use rand::prelude::*;
33/// use rand_distr::Pareto;
34///
35/// let val: f64 = rand::rng().sample(Pareto::new(1., 2.).unwrap());
36/// println!("{}", val);
37/// ```
38#[derive(Clone, Copy, Debug, PartialEq)]
39#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
40pub struct Pareto<F>
41where
42    F: Float,
43    OpenClosed01: Distribution<F>,
44{
45    scale: F,
46    inv_neg_shape: F,
47}
48
49/// Error type returned from [`Pareto::new`].
50#[derive(Clone, Copy, Debug, PartialEq, Eq)]
51pub enum Error {
52    /// `scale <= 0` or `nan`.
53    ScaleTooSmall,
54    /// `shape <= 0` or `nan`.
55    ShapeTooSmall,
56}
57
58impl fmt::Display for Error {
59    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
60        f.write_str(match self {
61            Error::ScaleTooSmall => "scale is not positive in Pareto distribution",
62            Error::ShapeTooSmall => "shape is not positive in Pareto distribution",
63        })
64    }
65}
66
67#[cfg(feature = "std")]
68impl std::error::Error for Error {}
69
70impl<F> Pareto<F>
71where
72    F: Float,
73    OpenClosed01: Distribution<F>,
74{
75    /// Construct a new Pareto distribution with given `scale` and `shape`.
76    ///
77    /// In the literature, `scale` is commonly written as x<sub>m</sub> or k and
78    /// `shape` is often written as α.
79    pub fn new(scale: F, shape: F) -> Result<Pareto<F>, Error> {
80        let zero = F::zero();
81
82        if !(scale > zero) {
83            return Err(Error::ScaleTooSmall);
84        }
85        if !(shape > zero) {
86            return Err(Error::ShapeTooSmall);
87        }
88        Ok(Pareto {
89            scale,
90            inv_neg_shape: F::from(-1.0).unwrap() / shape,
91        })
92    }
93}
94
95impl<F> Distribution<F> for Pareto<F>
96where
97    F: Float,
98    OpenClosed01: Distribution<F>,
99{
100    fn sample<R: Rng + ?Sized>(&self, rng: &mut R) -> F {
101        let u: F = OpenClosed01.sample(rng);
102        self.scale * u.powf(self.inv_neg_shape)
103    }
104}
105
106#[cfg(test)]
107mod tests {
108    use super::*;
109    use core::fmt::{Debug, Display, LowerExp};
110    use rand::RngExt;
111
112    #[test]
113    #[should_panic]
114    fn invalid() {
115        Pareto::new(0., 0.).unwrap();
116    }
117
118    #[test]
119    fn sample() {
120        let scale = 1.0;
121        let shape = 2.0;
122        let d = Pareto::new(scale, shape).unwrap();
123        let mut rng = crate::test::rng(1);
124        for _ in 0..1000 {
125            let r = d.sample(&mut rng);
126            assert!(r >= scale);
127        }
128    }
129
130    #[test]
131    fn value_stability() {
132        fn test_samples<F: Float + Debug + Display + LowerExp, D: Distribution<F>>(
133            distr: D,
134            thresh: F,
135            expected: &[F],
136        ) {
137            let mut rng = crate::test::rng(213);
138            for v in expected {
139                let x = rng.sample(&distr);
140                assert_almost_eq!(x, *v, thresh);
141            }
142        }
143
144        test_samples(
145            Pareto::new(1f32, 1.0).unwrap(),
146            1e-6,
147            &[1.0423688, 2.1235929, 4.132709, 1.4679428],
148        );
149        test_samples(
150            Pareto::new(2.0, 0.5).unwrap(),
151            1e-14,
152            &[
153                9.019295276219136,
154                4.3097126018270595,
155                6.837815045397157,
156                105.8826669383772,
157            ],
158        );
159    }
160
161    #[test]
162    fn pareto_distributions_can_be_compared() {
163        assert_eq!(Pareto::new(1.0, 2.0), Pareto::new(1.0, 2.0));
164    }
165}
rand_distr/pareto.rs

rand_distr/
pareto.rs
