1use num_traits::Float;
12use crate::{Distribution, OpenClosed01};
13use rand::Rng;
14use core::fmt;
15
16#[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#[derive(Clone, Copy, Debug, PartialEq, Eq)]
37pub enum Error {
38 ScaleTooSmall,
40 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 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}