From a391f2c4afa5c93ee2437114106d8371124ac154 Mon Sep 17 00:00:00 2001 From: mberocay Date: Fri, 18 May 2018 16:57:11 -0500 Subject: [PATCH 1/3] Added Cauchy distribution --- benches/distributions.rs | 1 + src/distributions/binomial.rs | 23 +++---- src/distributions/cauchy.rs | 118 ++++++++++++++++++++++++++++++++++ src/distributions/mod.rs | 6 ++ src/distributions/poisson.rs | 26 +++----- 5 files changed, 142 insertions(+), 32 deletions(-) create mode 100644 src/distributions/cauchy.rs diff --git a/benches/distributions.rs b/benches/distributions.rs index b2d21529866..fd6b5ae0304 100644 --- a/benches/distributions.rs +++ b/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)); diff --git a/src/distributions/binomial.rs b/src/distributions/binomial.rs index 85506275c6a..c4fba1b2eb7 100644 --- a/src/distributions/binomial.rs +++ b/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,20 +89,14 @@ impl Distribution 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::()).tan(); - // shift the peak of the comparison ditribution - lresult = expected + sq * comp_dev; - // repeat the drawing until we are in the range of possible values - if lresult >= 0.0 && lresult < float_n + 1.0 { - break; - } - } + // draw from the Cauchy distribution + let comp_dev = rng.sample(cauchy); + // shift the peak of the comparison distribution + lresult = expected + sq * comp_dev; // the result should be discrete lresult = lresult.floor(); diff --git a/src/distributions/cauchy.rs b/src/distributions/cauchy.rs new file mode 100644 index 00000000000..cbbd219ab3b --- /dev/null +++ b/src/distributions/cauchy.rs @@ -0,0 +1,118 @@ +// 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 or the MIT license +// , 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 for Cauchy { + fn sample(&self, rng: &mut R) -> f64 { + let mut x = rng.gen::(); + // guard against the extremely unlikely event we get 0 or 1 from the generator + while x <= 0.0 || x >= 1.0 { + x = rng.gen::(); + } + // get standard cauchy random number + let comp_dev = (2.0 * 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); + } +} diff --git a/src/distributions/mod.rs b/src/distributions/mod.rs index e73fec9fd39..f6dccc1bbf7 100644 --- a/src/distributions/mod.rs +++ b/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; diff --git a/src/distributions/poisson.rs b/src/distributions/poisson.rs index 8fbf031d0c7..a3bf8e858fa 100644 --- a/src/distributions/poisson.rs +++ b/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,16 @@ impl Distribution 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; + // draw from the Cauchy distribution + let comp_dev = rng.sample(cauchy); + // shift the peak of the comparison ditribution + let mut result = self.sqrt_2lambda * comp_dev + self.lambda; - // 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::()).tan(); - // 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; From e956a4823b004ffb50684501bd4a79ab155ada31 Mon Sep 17 00:00:00 2001 From: mberocay Date: Wed, 23 May 2018 13:47:37 -0500 Subject: [PATCH 2/3] Fixed mistake in the domain getting passed into the tangent function --- src/distributions/cauchy.rs | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/src/distributions/cauchy.rs b/src/distributions/cauchy.rs index cbbd219ab3b..d7084b6ccdb 100644 --- a/src/distributions/cauchy.rs +++ b/src/distributions/cauchy.rs @@ -49,13 +49,14 @@ impl Cauchy { impl Distribution for Cauchy { fn sample(&self, rng: &mut R) -> f64 { - let mut x = rng.gen::(); - // guard against the extremely unlikely event we get 0 or 1 from the generator - while x <= 0.0 || x >= 1.0 { + // sample from [0, 1) + let mut x: f64 = rng.gen::(); + // guard against the extremely unlikely case we get the invalid 0.5 + while x == 0.5 { x = rng.gen::(); } // get standard cauchy random number - let comp_dev = (2.0 * PI * x).tan(); + let comp_dev = (PI * x).tan(); // shift and scale according to parameters let result = self.median + self.scale * comp_dev; result From cc377b2280ed313556df92876e88dcedc5dfb783 Mon Sep 17 00:00:00 2001 From: mberocay Date: Thu, 24 May 2018 15:56:50 -0500 Subject: [PATCH 3/3] Adding guard loops back to binomial and poission --- src/distributions/binomial.rs | 15 +++++++++++---- src/distributions/cauchy.rs | 2 +- src/distributions/poisson.rs | 17 +++++++++++++---- 3 files changed, 25 insertions(+), 9 deletions(-) diff --git a/src/distributions/binomial.rs b/src/distributions/binomial.rs index c4fba1b2eb7..5eb03e38850 100644 --- a/src/distributions/binomial.rs +++ b/src/distributions/binomial.rs @@ -93,10 +93,17 @@ impl Distribution for Binomial { // f(x) ~ 1/(1+x^2) let cauchy = Cauchy::new(0.0, 1.0); loop { - // draw from the Cauchy distribution - let comp_dev = rng.sample(cauchy); - // shift the peak of the comparison distribution - lresult = expected + sq * comp_dev; + let mut comp_dev: f64; + loop { + // 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 + if lresult >= 0.0 && lresult < float_n + 1.0 { + break; + } + } // the result should be discrete lresult = lresult.floor(); diff --git a/src/distributions/cauchy.rs b/src/distributions/cauchy.rs index d7084b6ccdb..4b6c7e69624 100644 --- a/src/distributions/cauchy.rs +++ b/src/distributions/cauchy.rs @@ -50,7 +50,7 @@ impl Cauchy { impl Distribution for Cauchy { fn sample(&self, rng: &mut R) -> f64 { // sample from [0, 1) - let mut x: f64 = rng.gen::(); + let mut x = rng.gen::(); // guard against the extremely unlikely case we get the invalid 0.5 while x == 0.5 { x = rng.gen::(); diff --git a/src/distributions/poisson.rs b/src/distributions/poisson.rs index a3bf8e858fa..bdecb771536 100644 --- a/src/distributions/poisson.rs +++ b/src/distributions/poisson.rs @@ -75,12 +75,21 @@ impl Distribution for Poisson { // we use the Cauchy distribution as the comparison distribution // f(x) ~ 1/(1+x^2) let cauchy = Cauchy::new(0.0, 1.0); + loop { - // draw from the Cauchy distribution - let comp_dev = rng.sample(cauchy); - // shift the peak of the comparison ditribution - let mut result = self.sqrt_2lambda * comp_dev + self.lambda; + let mut result; + let mut comp_dev; + loop { + // 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 Cauchy distribution // the result should be an integer value result = result.floor();