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..5eb03e38850 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,13 +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(); + // 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 diff --git a/src/distributions/cauchy.rs b/src/distributions/cauchy.rs new file mode 100644 index 00000000000..4b6c7e69624 --- /dev/null +++ b/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 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 { + // sample from [0, 1) + let mut x = 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 = (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..bdecb771536 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,15 +72,17 @@ 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; - // 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(); + // 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 @@ -89,7 +90,7 @@ impl Distribution for Poisson { 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;