|
| 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 | + let mut x = rng.gen::<f64>(); |
| 53 | + // guard against the extremely unlikely event we get 0 or 1 from the generator |
| 54 | + while x <= 0.0 || x >= 1.0 { |
| 55 | + x = rng.gen::<f64>(); |
| 56 | + } |
| 57 | + // get standard cauchy random number |
| 58 | + let comp_dev = (2.0 * PI * x).tan(); |
| 59 | + // shift and scale according to parameters |
| 60 | + let result = self.median + self.scale * comp_dev; |
| 61 | + result |
| 62 | + } |
| 63 | +} |
| 64 | + |
| 65 | +#[cfg(test)] |
| 66 | +mod test { |
| 67 | + use distributions::Distribution; |
| 68 | + use super::Cauchy; |
| 69 | + |
| 70 | + fn median(mut numbers: &mut [f64]) -> f64 { |
| 71 | + sort(&mut numbers); |
| 72 | + let mid = numbers.len() / 2; |
| 73 | + numbers[mid] |
| 74 | + } |
| 75 | + |
| 76 | + fn sort(numbers: &mut [f64]) { |
| 77 | + numbers.sort_by(|a, b| a.partial_cmp(b).unwrap()); |
| 78 | + } |
| 79 | + |
| 80 | + #[test] |
| 81 | + fn test_cauchy_median() { |
| 82 | + let cauchy = Cauchy::new(10.0, 5.0); |
| 83 | + let mut rng = ::test::rng(123); |
| 84 | + let mut numbers: [f64; 1000] = [0.0; 1000]; |
| 85 | + for i in 0..1000 { |
| 86 | + numbers[i] = cauchy.sample(&mut rng); |
| 87 | + } |
| 88 | + let median = median(&mut numbers); |
| 89 | + println!("Cauchy median: {}", median); |
| 90 | + assert!((median - 10.0).abs() < 0.5); // not 100% certain, but probable enough |
| 91 | + } |
| 92 | + |
| 93 | + #[test] |
| 94 | + fn test_cauchy_mean() { |
| 95 | + let cauchy = Cauchy::new(10.0, 5.0); |
| 96 | + let mut rng = ::test::rng(123); |
| 97 | + let mut sum = 0.0; |
| 98 | + for _ in 0..1000 { |
| 99 | + sum += cauchy.sample(&mut rng); |
| 100 | + } |
| 101 | + let mean = sum / 1000.0; |
| 102 | + println!("Cauchy mean: {}", mean); |
| 103 | + // for a Cauchy distribution the mean should not converge |
| 104 | + assert!((mean - 10.0).abs() > 0.5); // not 100% certain, but probable enough |
| 105 | + } |
| 106 | + |
| 107 | + #[test] |
| 108 | + #[should_panic] |
| 109 | + fn test_cauchy_invalid_scale_zero() { |
| 110 | + Cauchy::new(0.0, 0.0); |
| 111 | + } |
| 112 | + |
| 113 | + #[test] |
| 114 | + #[should_panic] |
| 115 | + fn test_cauchy_invalid_scale_neg() { |
| 116 | + Cauchy::new(0.0, -10.0); |
| 117 | + } |
| 118 | +} |
0 commit comments