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 
12 use Rng;
13 use distributions::Distribution;
14 use std::f64::consts::PI;
15 
16 /// The Cauchy distribution `Cauchy(median, scale)`.
17 ///
18 /// This distribution has a density function:
19 /// `f(x) = 1 / (pi * scale * (1 + ((x - median) / scale)^2))`
20 ///
21 /// # Example
22 ///
23 /// ```
24 /// use rand::distributions::{Cauchy, Distribution};
25 ///
26 /// let cau = Cauchy::new(2.0, 5.0);
27 /// let v = cau.sample(&mut rand::thread_rng());
28 /// println!("{} is from a Cauchy(2, 5) distribution", v);
29 /// ```
30 #[derive(Clone, Copy, Debug)]
31 pub struct Cauchy {
32     median: f64,
33     scale: f64
34 }
35 
36 impl Cauchy {
37     /// Construct a new `Cauchy` with the given shape parameters
38     /// `median` the peak location and `scale` the scale factor.
39     /// Panics if `scale <= 0`.
new(median: f64, scale: f64) -> Cauchy40     pub fn new(median: f64, scale: f64) -> Cauchy {
41         assert!(scale > 0.0, "Cauchy::new called with scale factor <= 0");
42         Cauchy {
43             median,
44             scale
45         }
46     }
47 }
48 
49 impl Distribution<f64> for Cauchy {
sample<R: Rng + ?Sized>(&self, rng: &mut R) -> f6450     fn sample<R: Rng + ?Sized>(&self, rng: &mut R) -> f64 {
51         // sample from [0, 1)
52         let x = rng.gen::<f64>();
53         // get standard cauchy random number
54         // note that π/2 is not exactly representable, even if x=0.5 the result is finite
55         let comp_dev = (PI * x).tan();
56         // shift and scale according to parameters
57         let result = self.median + self.scale * comp_dev;
58         result
59     }
60 }
61 
62 #[cfg(test)]
63 mod test {
64     use distributions::Distribution;
65     use super::Cauchy;
66 
median(mut numbers: &mut [f64]) -> f6467     fn median(mut numbers: &mut [f64]) -> f64 {
68         sort(&mut numbers);
69         let mid = numbers.len() / 2;
70         numbers[mid]
71     }
72 
sort(numbers: &mut [f64])73     fn sort(numbers: &mut [f64]) {
74         numbers.sort_by(|a, b| a.partial_cmp(b).unwrap());
75     }
76 
77     #[test]
test_cauchy_median()78     fn test_cauchy_median() {
79         let cauchy = Cauchy::new(10.0, 5.0);
80         let mut rng = ::test::rng(123);
81         let mut numbers: [f64; 1000] = [0.0; 1000];
82         for i in 0..1000 {
83             numbers[i] = cauchy.sample(&mut rng);
84         }
85         let median = median(&mut numbers);
86         println!("Cauchy median: {}", median);
87         assert!((median - 10.0).abs() < 0.5); // not 100% certain, but probable enough
88     }
89 
90     #[test]
test_cauchy_mean()91     fn test_cauchy_mean() {
92         let cauchy = Cauchy::new(10.0, 5.0);
93         let mut rng = ::test::rng(123);
94         let mut sum = 0.0;
95         for _ in 0..1000 {
96             sum += cauchy.sample(&mut rng);
97         }
98         let mean = sum / 1000.0;
99         println!("Cauchy mean: {}", mean);
100         // for a Cauchy distribution the mean should not converge
101         assert!((mean - 10.0).abs() > 0.5); // not 100% certain, but probable enough
102     }
103 
104     #[test]
105     #[should_panic]
test_cauchy_invalid_scale_zero()106     fn test_cauchy_invalid_scale_zero() {
107         Cauchy::new(0.0, 0.0);
108     }
109 
110     #[test]
111     #[should_panic]
test_cauchy_invalid_scale_neg()112     fn test_cauchy_invalid_scale_neg() {
113         Cauchy::new(0.0, -10.0);
114     }
115 }
116