]>
Commit | Line | Data |
---|---|---|
b7449926 XL |
1 | // Copyright 2016-2017 The Rust Project Developers. See the COPYRIGHT |
2 | // file at the top-level directory of this distribution and at | |
3 | // https://rust-lang.org/COPYRIGHT. | |
4 | // | |
5 | // Licensed under the Apache License, Version 2.0 <LICENSE-APACHE or | |
6 | // https://www.apache.org/licenses/LICENSE-2.0> or the MIT license | |
7 | // <LICENSE-MIT or https://opensource.org/licenses/MIT>, at your | |
8 | // option. This file may not be copied, modified, or distributed | |
9 | // except according to those terms. | |
10 | ||
11 | //! The Cauchy distribution. | |
12 | ||
13 | use Rng; | |
14 | use distributions::Distribution; | |
15 | use std::f64::consts::PI; | |
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 | /// # Example | |
23 | /// | |
24 | /// ``` | |
25 | /// use rand::distributions::{Cauchy, Distribution}; | |
26 | /// | |
27 | /// let cau = Cauchy::new(2.0, 5.0); | |
28 | /// let v = cau.sample(&mut rand::thread_rng()); | |
29 | /// println!("{} is from a Cauchy(2, 5) distribution", v); | |
30 | /// ``` | |
31 | #[derive(Clone, Copy, Debug)] | |
32 | pub struct Cauchy { | |
33 | median: f64, | |
34 | scale: f64 | |
35 | } | |
36 | ||
37 | impl Cauchy { | |
38 | /// Construct a new `Cauchy` with the given shape parameters | |
39 | /// `median` the peak location and `scale` the scale factor. | |
40 | /// Panics if `scale <= 0`. | |
41 | pub fn new(median: f64, scale: f64) -> Cauchy { | |
42 | assert!(scale > 0.0, "Cauchy::new called with scale factor <= 0"); | |
43 | Cauchy { | |
44 | median, | |
45 | scale | |
46 | } | |
47 | } | |
48 | } | |
49 | ||
50 | impl Distribution<f64> for Cauchy { | |
51 | fn sample<R: Rng + ?Sized>(&self, rng: &mut R) -> f64 { | |
52 | // sample from [0, 1) | |
53 | let x = rng.gen::<f64>(); | |
54 | // get standard cauchy random number | |
55 | // note that π/2 is not exactly representable, even if x=0.5 the result is finite | |
56 | let comp_dev = (PI * x).tan(); | |
57 | // shift and scale according to parameters | |
58 | let result = self.median + self.scale * comp_dev; | |
59 | result | |
60 | } | |
61 | } | |
62 | ||
63 | #[cfg(test)] | |
64 | mod test { | |
65 | use distributions::Distribution; | |
66 | use super::Cauchy; | |
67 | ||
68 | fn median(mut numbers: &mut [f64]) -> f64 { | |
69 | sort(&mut numbers); | |
70 | let mid = numbers.len() / 2; | |
71 | numbers[mid] | |
72 | } | |
73 | ||
74 | fn sort(numbers: &mut [f64]) { | |
75 | numbers.sort_by(|a, b| a.partial_cmp(b).unwrap()); | |
76 | } | |
77 | ||
78 | #[test] | |
79 | fn test_cauchy_median() { | |
80 | let cauchy = Cauchy::new(10.0, 5.0); | |
81 | let mut rng = ::test::rng(123); | |
82 | let mut numbers: [f64; 1000] = [0.0; 1000]; | |
83 | for i in 0..1000 { | |
84 | numbers[i] = cauchy.sample(&mut rng); | |
85 | } | |
86 | let median = median(&mut numbers); | |
87 | println!("Cauchy median: {}", median); | |
88 | assert!((median - 10.0).abs() < 0.5); // not 100% certain, but probable enough | |
89 | } | |
90 | ||
91 | #[test] | |
92 | fn test_cauchy_mean() { | |
93 | let cauchy = Cauchy::new(10.0, 5.0); | |
94 | let mut rng = ::test::rng(123); | |
95 | let mut sum = 0.0; | |
96 | for _ in 0..1000 { | |
97 | sum += cauchy.sample(&mut rng); | |
98 | } | |
99 | let mean = sum / 1000.0; | |
100 | println!("Cauchy mean: {}", mean); | |
101 | // for a Cauchy distribution the mean should not converge | |
102 | assert!((mean - 10.0).abs() > 0.5); // not 100% certain, but probable enough | |
103 | } | |
104 | ||
105 | #[test] | |
106 | #[should_panic] | |
107 | fn test_cauchy_invalid_scale_zero() { | |
108 | Cauchy::new(0.0, 0.0); | |
109 | } | |
110 | ||
111 | #[test] | |
112 | #[should_panic] | |
113 | fn test_cauchy_invalid_scale_neg() { | |
114 | Cauchy::new(0.0, -10.0); | |
115 | } | |
116 | } |