rand_distr/
cauchy.rs

1// Copyright 2018 Developers of the Rand project.
2// Copyright 2016-2017 The Rust Project Developers.
3//
4// Licensed under the Apache License, Version 2.0 <LICENSE-APACHE or
5// https://www.apache.org/licenses/LICENSE-2.0> or the MIT license
6// <LICENSE-MIT or https://opensource.org/licenses/MIT>, at your
7// option. This file may not be copied, modified, or distributed
8// except according to those terms.
9
10//! The Cauchy distribution.
11
12use num_traits::{Float, FloatConst};
13use crate::{Distribution, Standard};
14use rand::Rng;
15use core::fmt;
16
17/// The Cauchy distribution `Cauchy(median, scale)`.
18///
19/// This distribution has a density function:
20/// `f(x) = 1 / (pi * scale * (1 + ((x - median) / scale)^2))`
21///
22/// Note that at least for `f32`, results are not fully portable due to minor
23/// differences in the target system's *tan* implementation, `tanf`.
24///
25/// # Example
26///
27/// ```
28/// use rand_distr::{Cauchy, Distribution};
29///
30/// let cau = Cauchy::new(2.0, 5.0).unwrap();
31/// let v = cau.sample(&mut rand::thread_rng());
32/// println!("{} is from a Cauchy(2, 5) distribution", v);
33/// ```
34#[derive(Clone, Copy, Debug)]
35#[cfg_attr(feature = "serde1", derive(serde::Serialize, serde::Deserialize))]
36pub struct Cauchy<F>
37where F: Float + FloatConst, Standard: Distribution<F>
38{
39    median: F,
40    scale: F,
41}
42
43/// Error type returned from `Cauchy::new`.
44#[derive(Clone, Copy, Debug, PartialEq, Eq)]
45pub enum Error {
46    /// `scale <= 0` or `nan`.
47    ScaleTooSmall,
48}
49
50impl fmt::Display for Error {
51    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
52        f.write_str(match self {
53            Error::ScaleTooSmall => "scale is not positive in Cauchy distribution",
54        })
55    }
56}
57
58#[cfg(feature = "std")]
59#[cfg_attr(doc_cfg, doc(cfg(feature = "std")))]
60impl std::error::Error for Error {}
61
62impl<F> Cauchy<F>
63where F: Float + FloatConst, Standard: Distribution<F>
64{
65    /// Construct a new `Cauchy` with the given shape parameters
66    /// `median` the peak location and `scale` the scale factor.
67    pub fn new(median: F, scale: F) -> Result<Cauchy<F>, Error> {
68        if !(scale > F::zero()) {
69            return Err(Error::ScaleTooSmall);
70        }
71        Ok(Cauchy { median, scale })
72    }
73}
74
75impl<F> Distribution<F> for Cauchy<F>
76where F: Float + FloatConst, Standard: Distribution<F>
77{
78    fn sample<R: Rng + ?Sized>(&self, rng: &mut R) -> F {
79        // sample from [0, 1)
80        let x = Standard.sample(rng);
81        // get standard cauchy random number
82        // note that π/2 is not exactly representable, even if x=0.5 the result is finite
83        let comp_dev = (F::PI() * x).tan();
84        // shift and scale according to parameters
85        self.median + self.scale * comp_dev
86    }
87}
88
89#[cfg(test)]
90mod test {
91    use super::*;
92
93    fn median(numbers: &mut [f64]) -> f64 {
94        sort(numbers);
95        let mid = numbers.len() / 2;
96        numbers[mid]
97    }
98
99    fn sort(numbers: &mut [f64]) {
100        numbers.sort_by(|a, b| a.partial_cmp(b).unwrap());
101    }
102
103    #[test]
104    fn test_cauchy_averages() {
105        // NOTE: given that the variance and mean are undefined,
106        // this test does not have any rigorous statistical meaning.
107        let cauchy = Cauchy::new(10.0, 5.0).unwrap();
108        let mut rng = crate::test::rng(123);
109        let mut numbers: [f64; 1000] = [0.0; 1000];
110        let mut sum = 0.0;
111        for number in &mut numbers[..] {
112            *number = cauchy.sample(&mut rng);
113            sum += *number;
114        }
115        let median = median(&mut numbers);
116        #[cfg(feature = "std")]
117        std::println!("Cauchy median: {}", median);
118        assert!((median - 10.0).abs() < 0.4); // not 100% certain, but probable enough
119        let mean = sum / 1000.0;
120        #[cfg(feature = "std")]
121        std::println!("Cauchy mean: {}", mean);
122        // for a Cauchy distribution the mean should not converge
123        assert!((mean - 10.0).abs() > 0.4); // not 100% certain, but probable enough
124    }
125
126    #[test]
127    #[should_panic]
128    fn test_cauchy_invalid_scale_zero() {
129        Cauchy::new(0.0, 0.0).unwrap();
130    }
131
132    #[test]
133    #[should_panic]
134    fn test_cauchy_invalid_scale_neg() {
135        Cauchy::new(0.0, -10.0).unwrap();
136    }
137
138    #[test]
139    fn value_stability() {
140        fn gen_samples<F: Float + FloatConst + core::fmt::Debug>(m: F, s: F, buf: &mut [F])
141        where Standard: Distribution<F> {
142            let distr = Cauchy::new(m, s).unwrap();
143            let mut rng = crate::test::rng(353);
144            for x in buf {
145                *x = rng.sample(&distr);
146            }
147        }
148
149        let mut buf = [0.0; 4];
150        gen_samples(100f64, 10.0, &mut buf);
151        assert_eq!(&buf, &[
152            77.93369152808678,
153            90.1606912098641,
154            125.31516221323625,
155            86.10217834773925
156        ]);
157
158        // Unfortunately this test is not fully portable due to reliance on the
159        // system's implementation of tanf (see doc on Cauchy struct).
160        let mut buf = [0.0; 4];
161        gen_samples(10f32, 7.0, &mut buf);
162        let expected = [15.023088, -5.446413, 3.7092876, 3.112482];
163        for (a, b) in buf.iter().zip(expected.iter()) {
164            assert_almost_eq!(*a, *b, 1e-5);
165        }
166    }
167}
rand_distr/cauchy.rs

rand_distr/
cauchy.rs
