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.
10
11use num_traits::Float;
12use crate::{Distribution, OpenClosed01};
13use rand::Rng;
14use core::fmt;
15
16/// Samples floating-point numbers according to the Pareto distribution
17///
18/// # Example
19/// ```
20/// use rand::prelude::*;
21/// use rand_distr::Pareto;
22///
23/// let val: f64 = thread_rng().sample(Pareto::new(1., 2.).unwrap());
24/// println!("{}", val);
25/// ```
26#[derive(Clone, Copy, Debug)]
27#[cfg_attr(feature = "serde1", derive(serde::Serialize, serde::Deserialize))]
28pub struct Pareto<F>
29where F: Float, OpenClosed01: Distribution<F>
30{
31    scale: F,
32    inv_neg_shape: F,
33}
34
35/// Error type returned from `Pareto::new`.
36#[derive(Clone, Copy, Debug, PartialEq, Eq)]
37pub enum Error {
38    /// `scale <= 0` or `nan`.
39    ScaleTooSmall,
40    /// `shape <= 0` or `nan`.
41    ShapeTooSmall,
42}
43
44impl fmt::Display for Error {
45    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
46        f.write_str(match self {
47            Error::ScaleTooSmall => "scale is not positive in Pareto distribution",
48            Error::ShapeTooSmall => "shape is not positive in Pareto distribution",
49        })
50    }
51}
52
53#[cfg(feature = "std")]
54#[cfg_attr(doc_cfg, doc(cfg(feature = "std")))]
55impl std::error::Error for Error {}
56
57impl<F> Pareto<F>
58where F: Float, OpenClosed01: Distribution<F>
59{
60    /// Construct a new Pareto distribution with given `scale` and `shape`.
61    ///
62    /// In the literature, `scale` is commonly written as x<sub>m</sub> or k and
63    /// `shape` is often written as α.
64    pub fn new(scale: F, shape: F) -> Result<Pareto<F>, Error> {
65        let zero = F::zero();
66
67        if !(scale > zero) {
68            return Err(Error::ScaleTooSmall);
69        }
70        if !(shape > zero) {
71            return Err(Error::ShapeTooSmall);
72        }
73        Ok(Pareto {
74            scale,
75            inv_neg_shape: F::from(-1.0).unwrap() / shape,
76        })
77    }
78}
79
80impl<F> Distribution<F> for Pareto<F>
81where F: Float, OpenClosed01: Distribution<F>
82{
83    fn sample<R: Rng + ?Sized>(&self, rng: &mut R) -> F {
84        let u: F = OpenClosed01.sample(rng);
85        self.scale * u.powf(self.inv_neg_shape)
86    }
87}
88
89#[cfg(test)]
90mod tests {
91    use super::*;
92    use core::fmt::{Debug, Display, LowerExp};
93
94    #[test]
95    #[should_panic]
96    fn invalid() {
97        Pareto::new(0., 0.).unwrap();
98    }
99
100    #[test]
101    fn sample() {
102        let scale = 1.0;
103        let shape = 2.0;
104        let d = Pareto::new(scale, shape).unwrap();
105        let mut rng = crate::test::rng(1);
106        for _ in 0..1000 {
107            let r = d.sample(&mut rng);
108            assert!(r >= scale);
109        }
110    }
111
112    #[test]
113    fn value_stability() {
114        fn test_samples<F: Float + Debug + Display + LowerExp, D: Distribution<F>>(
115            distr: D, thresh: F, expected: &[F],
116        ) {
117            let mut rng = crate::test::rng(213);
118            for v in expected {
119                let x = rng.sample(&distr);
120                assert_almost_eq!(x, *v, thresh);
121            }
122        }
123
124        test_samples(Pareto::new(1f32, 1.0).unwrap(), 1e-6, &[
125            1.0423688, 2.1235929, 4.132709, 1.4679428,
126        ]);
127        test_samples(Pareto::new(2.0, 0.5).unwrap(), 1e-14, &[
128            9.019295276219136,
129            4.3097126018270595,
130            6.837815045397157,
131            105.8826669383772,
132        ]);
133    }
134}