Merge pull request #474 from MaximoB/add_cauchy_distribution

Added Cauchy distribution

Author: dhardy

benches/distributions.rs

@@ -109,6 +109,7 @@ distr_float!(distr_normal, f64, Normal::new(-1.23, 4.56));
 distr_float!(distr_log_normal, f64, LogNormal::new(-1.23, 4.56));
 distr_float!(distr_gamma_large_shape, f64, Gamma::new(10., 1.0));
 distr_float!(distr_gamma_small_shape, f64, Gamma::new(0.1, 1.0));
+distr_float!(distr_cauchy, f64, Cauchy::new(4.2, 6.9));
 distr_int!(distr_binomial, u64, Binomial::new(20, 0.7));
 distr_int!(distr_poisson, u64, Poisson::new(4.0));
 distr!(distr_bernoulli, bool, Bernoulli::new(0.18));

src/distributions/binomial.rs

@@ -11,9 +11,8 @@
 //! The binomial distribution.
 
 use Rng;
-use distributions::{Distribution, Bernoulli};
+use distributions::{Distribution, Bernoulli, Cauchy};
 use distributions::log_gamma::log_gamma;
-use std::f64::consts::PI;
 
 /// The binomial distribution `Binomial(n, p)`.
 ///
@@ -90,13 +89,14 @@ impl Distribution<u64> for Binomial {
 
         let mut lresult;
 
+        // we use the Cauchy distribution as the comparison distribution
+        // f(x) ~ 1/(1+x^2)
+        let cauchy = Cauchy::new(0.0, 1.0);
         loop {
             let mut comp_dev: f64;
-            // we use the lorentzian distribution as the comparison distribution
-            // f(x) ~ 1/(1+x/^2)
             loop {
-                // draw from the lorentzian distribution
-                comp_dev = (PI*rng.gen::<f64>()).tan();
+                // draw from the Cauchy distribution
+                comp_dev = rng.sample(cauchy);
                 // shift the peak of the comparison ditribution
                 lresult = expected + sq * comp_dev;
                 // repeat the drawing until we are in the range of possible values

src/distributions/cauchy.rs

@@ -0,0 +1,119 @@
+// Copyright 2016-2017 The Rust Project Developers. See the COPYRIGHT
+// file at the top-level directory of this distribution and at
+// https://rust-lang.org/COPYRIGHT.
+//
+// Licensed under the Apache License, Version 2.0 <LICENSE-APACHE or
+// https://www.apache.org/licenses/LICENSE-2.0> or the MIT license
+// <LICENSE-MIT or https://opensource.org/licenses/MIT>, at your
+// option. This file may not be copied, modified, or distributed
+// except according to those terms.
+
+//! The Cauchy distribution.
+
+use Rng;
+use distributions::Distribution;
+use std::f64::consts::PI;
+
+/// The Cauchy distribution `Cauchy(median, scale)`.
+///
+/// This distribution has a density function:
+/// `f(x) = 1 / (pi * scale * (1 + ((x - median) / scale)^2))`
+///
+/// # Example
+///
+/// ```
+/// use rand::distributions::{Cauchy, Distribution};
+///
+/// let cau = Cauchy::new(2.0, 5.0);
+/// let v = cau.sample(&mut rand::thread_rng());
+/// println!("{} is from a Cauchy(2, 5) distribution", v);
+/// ```
+#[derive(Clone, Copy, Debug)]
+pub struct Cauchy {
+    median: f64,
+    scale: f64
+}
+
+impl Cauchy {
+    /// Construct a new `Cauchy` with the given shape parameters
+    /// `median` the peak location and `scale` the scale factor.
+    /// Panics if `scale <= 0`.
+    pub fn new(median: f64, scale: f64) -> Cauchy {
+        assert!(scale > 0.0, "Cauchy::new called with scale factor <= 0");
+        Cauchy {
+            median,
+            scale
+        }
+    }
+}
+
+impl Distribution<f64> for Cauchy {
+    fn sample<R: Rng + ?Sized>(&self, rng: &mut R) -> f64 {
+        // sample from [0, 1)
+        let mut x = rng.gen::<f64>();
+        // guard against the extremely unlikely case we get the invalid 0.5
+        while x == 0.5 {
+            x = rng.gen::<f64>();
+        }
+        // get standard cauchy random number
+        let comp_dev = (PI * x).tan();
+        // shift and scale according to parameters
+        let result = self.median + self.scale * comp_dev;
+        result
+    }
+}
+
+#[cfg(test)]
+mod test {
+    use distributions::Distribution;
+    use super::Cauchy;
+
+    fn median(mut numbers: &mut [f64]) -> f64 {
+        sort(&mut numbers);
+        let mid = numbers.len() / 2;
+        numbers[mid]
+    }
+
+    fn sort(numbers: &mut [f64]) {
+        numbers.sort_by(|a, b| a.partial_cmp(b).unwrap());
+    }
+
+    #[test]
+    fn test_cauchy_median() {
+        let cauchy = Cauchy::new(10.0, 5.0);
+        let mut rng = ::test::rng(123);
+        let mut numbers: [f64; 1000] = [0.0; 1000];
+        for i in 0..1000 {
+            numbers[i] = cauchy.sample(&mut rng);
+        }
+        let median = median(&mut numbers);
+        println!("Cauchy median: {}", median);
+        assert!((median - 10.0).abs() < 0.5); // not 100% certain, but probable enough
+    }
+
+    #[test]
+    fn test_cauchy_mean() {
+        let cauchy = Cauchy::new(10.0, 5.0);
+        let mut rng = ::test::rng(123);
+        let mut sum = 0.0;
+        for _ in 0..1000 {
+            sum += cauchy.sample(&mut rng);
+        }
+        let mean = sum / 1000.0;
+        println!("Cauchy mean: {}", mean);
+        // for a Cauchy distribution the mean should not converge
+        assert!((mean - 10.0).abs() > 0.5); // not 100% certain, but probable enough
+    }
+
+    #[test]
+    #[should_panic]
+    fn test_cauchy_invalid_scale_zero() {
+        Cauchy::new(0.0, 0.0);
+    }
+
+    #[test]
+    #[should_panic]
+    fn test_cauchy_invalid_scale_neg() {
+        Cauchy::new(0.0, -10.0);
+    }
+}

src/distributions/mod.rs

@@ -81,6 +81,7 @@
 //! - Related to real-valued quantities that grow linearly
 //!   (e.g. errors, offsets):
 //!   - [`Normal`] distribution, and [`StandardNormal`] as a primitive
+//!   - [`Cauchy`] distribution
 //! - Related to Bernoulli trials (yes/no events, with a given probability):
 //!   - [`Binomial`] distribution
 //!   - [`Bernoulli`] distribution, similar to [`Rng::gen_bool`].
@@ -147,6 +148,7 @@
 //! [`Alphanumeric`]: struct.Alphanumeric.html
 //! [`Bernoulli`]: struct.Bernoulli.html
 //! [`Binomial`]: struct.Binomial.html
+//! [`Cauchy`]: struct.Cauchy.html
 //! [`ChiSquared`]: struct.ChiSquared.html
 //! [`Exp`]: struct.Exp.html
 //! [`Exp1`]: struct.Exp1.html
@@ -180,6 +182,8 @@ pub use self::uniform::Uniform as Range;
 #[cfg(feature = "std")]
 #[doc(inline)] pub use self::binomial::Binomial;
 #[doc(inline)] pub use self::bernoulli::Bernoulli;
+#[cfg(feature = "std")]
+#[doc(inline)] pub use self::cauchy::Cauchy;
 
 pub mod uniform;
 #[cfg(feature="std")]
@@ -193,6 +197,8 @@ pub mod uniform;
 #[cfg(feature = "std")]
 #[doc(hidden)] pub mod binomial;
 #[doc(hidden)] pub mod bernoulli;
+#[cfg(feature = "std")]
+#[doc(hidden)] pub mod cauchy;
 
 mod float;
 mod integer;

src/distributions/poisson.rs

@@ -11,9 +11,8 @@
 //! The Poisson distribution.
 
 use Rng;
-use distributions::Distribution;
+use distributions::{Distribution, Cauchy};
 use distributions::log_gamma::log_gamma;
-use std::f64::consts::PI;
 
 /// The Poisson distribution `Poisson(lambda)`.
 ///
@@ -73,23 +72,25 @@ impl Distribution<u64> for Poisson {
         else {
             let mut int_result: u64;
 
+            // we use the Cauchy distribution as the comparison distribution
+            // f(x) ~ 1/(1+x^2)
+            let cauchy = Cauchy::new(0.0, 1.0);
+
             loop {
                 let mut result;
                 let mut comp_dev;
 
-                // we use the lorentzian distribution as the comparison distribution
-                // f(x) ~ 1/(1+x/^2)
                 loop {
-                    // draw from the lorentzian distribution
-                    comp_dev = (PI * rng.gen::<f64>()).tan();
+                    // draw from the Cauchy distribution
+                    comp_dev = rng.sample(cauchy);
                     // shift the peak of the comparison ditribution
                     result = self.sqrt_2lambda * comp_dev + self.lambda;
                     // repeat the drawing until we are in the range of possible values
                     if result >= 0.0 {
                         break;
                     }
                 }
-                // now the result is a random variable greater than 0 with Lorentzian distribution
+                // now the result is a random variable greater than 0 with Cauchy distribution
                 // the result should be an integer value
                 result = result.floor();
                 int_result = result as u64;