diff --git a/src/libstd/rand/distributions/exponential.rs b/src/libstd/rand/distributions/exponential.rs new file mode 100644 index 00000000000..4244e6bacdb --- /dev/null +++ b/src/libstd/rand/distributions/exponential.rs @@ -0,0 +1,141 @@ +// Copyright 2013 The Rust Project Developers. See the COPYRIGHT +// file at the top-level directory of this distribution and at +// http://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 exponential distribution. + +use rand::{Rng, Rand}; +use rand::distributions::{ziggurat, ziggurat_tables, Sample, IndependentSample}; + +/// A wrapper around an `f64` to generate Exp(1) random numbers. +/// +/// See `Exp` for the general exponential distribution.Note that this + // has to be unwrapped before use as an `f64` (using either +/// `*` or `cast::transmute` is safe). +/// +/// Implemented via the ZIGNOR variant[1] of the Ziggurat method. The +/// exact description in the paper was adjusted to use tables for the +/// exponential distribution rather than normal. +/// +/// [1]: Jurgen A. Doornik (2005). [*An Improved Ziggurat Method to +/// Generate Normal Random +/// Samples*](http://www.doornik.com/research/ziggurat.pdf). Nuffield +/// College, Oxford +pub struct Exp1(f64); + +// This could be done via `-rng.gen::().ln()` but that is slower. +impl Rand for Exp1 { + #[inline] + fn rand(rng: &mut R) -> Exp1 { + #[inline] + fn pdf(x: f64) -> f64 { + (-x).exp() + } + #[inline] + fn zero_case(rng: &mut R, _u: f64) -> f64 { + ziggurat_tables::ZIG_EXP_R - rng.gen::().ln() + } + + Exp1(ziggurat(rng, false, + &ziggurat_tables::ZIG_EXP_X, + &ziggurat_tables::ZIG_EXP_F, + pdf, zero_case)) + } +} + +/// The exponential distribution `Exp(lambda)`. +/// +/// This distribution has density function: `f(x) = lambda * +/// exp(-lambda * x)` for `x > 0`. +/// +/// # Example +/// +/// ```rust +/// use std::rand; +/// use std::rand::distributions::{Exp, IndependentSample}; +/// +/// fn main() { +/// let exp = Exp::new(2.0); +/// let v = exp.ind_sample(&mut rand::task_rng()); +/// println!("{} is from a Exp(2) distribution", v); +/// } +/// ``` +pub struct Exp { + /// `lambda` stored as `1/lambda`, since this is what we scale by. + priv lambda_inverse: f64 +} + +impl Exp { + /// Construct a new `Exp` with the given shape parameter + /// `lambda`. Fails if `lambda <= 0`. + pub fn new(lambda: f64) -> Exp { + assert!(lambda > 0.0, "Exp::new called with `lambda` <= 0"); + Exp { lambda_inverse: 1.0 / lambda } + } +} + +impl Sample for Exp { + fn sample(&mut self, rng: &mut R) -> f64 { self.ind_sample(rng) } +} +impl IndependentSample for Exp { + fn ind_sample(&self, rng: &mut R) -> f64 { + (*rng.gen::()) * self.lambda_inverse + } +} + +#[cfg(test)] +mod test { + use rand::*; + use super::*; + use iter::range; + use option::{Some, None}; + + #[test] + fn test_exp() { + let mut exp = Exp::new(10.0); + let mut rng = task_rng(); + for _ in range(0, 1000) { + assert!(exp.sample(&mut rng) >= 0.0); + assert!(exp.ind_sample(&mut rng) >= 0.0); + } + } + #[test] + #[should_fail] + fn test_exp_invalid_lambda_zero() { + Exp::new(0.0); + } + #[test] + #[should_fail] + fn test_exp_invalid_lambda_neg() { + Exp::new(-10.0); + } +} + +#[cfg(test)] +mod bench { + use extra::test::BenchHarness; + use rand::{XorShiftRng, RAND_BENCH_N}; + use super::*; + use iter::range; + use option::{Some, None}; + use mem::size_of; + + #[bench] + fn rand_exp(bh: &mut BenchHarness) { + let mut rng = XorShiftRng::new(); + let mut exp = Exp::new(2.71828 * 3.14159); + + bh.iter(|| { + for _ in range(0, RAND_BENCH_N) { + exp.sample(&mut rng); + } + }); + bh.bytes = size_of::() as u64 * RAND_BENCH_N; + } +} diff --git a/src/libstd/rand/distributions/gamma.rs b/src/libstd/rand/distributions/gamma.rs index e0428742459..61929a7c479 100644 --- a/src/libstd/rand/distributions/gamma.rs +++ b/src/libstd/rand/distributions/gamma.rs @@ -11,7 +11,8 @@ //! The Gamma distribution. use rand::{Rng, Open01}; -use super::{IndependentSample, Sample, StandardNormal, Exp}; +use super::{IndependentSample, Sample, Exp}; +use super::normal::StandardNormal; use num; /// The Gamma distribution `Gamma(shape, scale)` distribution. diff --git a/src/libstd/rand/distributions/mod.rs b/src/libstd/rand/distributions/mod.rs index 4778e81f951..2a2932bdb4c 100644 --- a/src/libstd/rand/distributions/mod.rs +++ b/src/libstd/rand/distributions/mod.rs @@ -23,14 +23,18 @@ use iter::range; use option::{Some, None}; use num; -use rand::{Rng, Rand, Open01}; +use rand::{Rng, Rand}; use clone::Clone; pub use self::range::Range; pub use self::gamma::Gamma; +pub use self::normal::Normal; +pub use self::exponential::Exp; pub mod range; pub mod gamma; +pub mod normal; +pub mod exponential; /// Types that can be used to create a random instance of `Support`. pub trait Sample { @@ -246,181 +250,10 @@ fn ziggurat( } } -/// A wrapper around an `f64` to generate N(0, 1) random numbers -/// (a.k.a. a standard normal, or Gaussian). -/// -/// See `Normal` for the general normal distribution. That this has to -/// be unwrapped before use as an `f64` (using either `*` or -/// `cast::transmute` is safe). -/// -/// Implemented via the ZIGNOR variant[1] of the Ziggurat method. -/// -/// [1]: Jurgen A. Doornik (2005). [*An Improved Ziggurat Method to -/// Generate Normal Random -/// Samples*](http://www.doornik.com/research/ziggurat.pdf). Nuffield -/// College, Oxford -pub struct StandardNormal(f64); - -impl Rand for StandardNormal { - fn rand(rng: &mut R) -> StandardNormal { - #[inline] - fn pdf(x: f64) -> f64 { - ((-x*x/2.0) as f64).exp() - } - #[inline] - fn zero_case(rng: &mut R, u: f64) -> f64 { - // compute a random number in the tail by hand - - // strange initial conditions, because the loop is not - // do-while, so the condition should be true on the first - // run, they get overwritten anyway (0 < 1, so these are - // good). - let mut x = 1.0f64; - let mut y = 0.0f64; - - while -2.0 * y < x * x { - let x_ = *rng.gen::>(); - let y_ = *rng.gen::>(); - - x = x_.ln() / ziggurat_tables::ZIG_NORM_R; - y = y_.ln(); - } - - if u < 0.0 { x - ziggurat_tables::ZIG_NORM_R } else { ziggurat_tables::ZIG_NORM_R - x } - } - - StandardNormal(ziggurat( - rng, - true, // this is symmetric - &ziggurat_tables::ZIG_NORM_X, - &ziggurat_tables::ZIG_NORM_F, - pdf, zero_case)) - } -} - -/// The normal distribution `N(mean, std_dev**2)`. -/// -/// This uses the ZIGNOR variant of the Ziggurat method, see -/// `StandardNormal` for more details. -/// -/// # Example -/// -/// ``` -/// use std::rand; -/// use std::rand::distributions::{Normal, IndependentSample}; -/// -/// fn main() { -/// let normal = Normal::new(2.0, 3.0); -/// let v = normal.ind_sample(rand::task_rng()); -/// println!("{} is from a N(2, 9) distribution", v) -/// } -/// ``` -pub struct Normal { - priv mean: f64, - priv std_dev: f64 -} - -impl Normal { - /// Construct a new `Normal` distribution with the given mean and - /// standard deviation. Fails if `std_dev < 0`. - pub fn new(mean: f64, std_dev: f64) -> Normal { - assert!(std_dev >= 0.0, "Normal::new called with `std_dev` < 0"); - Normal { - mean: mean, - std_dev: std_dev - } - } -} -impl Sample for Normal { - fn sample(&mut self, rng: &mut R) -> f64 { self.ind_sample(rng) } -} -impl IndependentSample for Normal { - fn ind_sample(&self, rng: &mut R) -> f64 { - self.mean + self.std_dev * (*rng.gen::()) - } -} - -/// A wrapper around an `f64` to generate Exp(1) random numbers. -/// -/// See `Exp` for the general exponential distribution.Note that this - // has to be unwrapped before use as an `f64` (using either -/// `*` or `cast::transmute` is safe). -/// -/// Implemented via the ZIGNOR variant[1] of the Ziggurat method. The -/// exact description in the paper was adjusted to use tables for the -/// exponential distribution rather than normal. -/// -/// [1]: Jurgen A. Doornik (2005). [*An Improved Ziggurat Method to -/// Generate Normal Random -/// Samples*](http://www.doornik.com/research/ziggurat.pdf). Nuffield -/// College, Oxford -pub struct Exp1(f64); - -// This could be done via `-rng.gen::().ln()` but that is slower. -impl Rand for Exp1 { - #[inline] - fn rand(rng: &mut R) -> Exp1 { - #[inline] - fn pdf(x: f64) -> f64 { - (-x).exp() - } - #[inline] - fn zero_case(rng: &mut R, _u: f64) -> f64 { - ziggurat_tables::ZIG_EXP_R - rng.gen::().ln() - } - - Exp1(ziggurat(rng, false, - &ziggurat_tables::ZIG_EXP_X, - &ziggurat_tables::ZIG_EXP_F, - pdf, zero_case)) - } -} - -/// The exponential distribution `Exp(lambda)`. -/// -/// This distribution has density function: `f(x) = lambda * -/// exp(-lambda * x)` for `x > 0`. -/// -/// # Example -/// -/// ``` -/// use std::rand; -/// use std::rand::distributions::{Exp, IndependentSample}; -/// -/// fn main() { -/// let exp = Exp::new(2.0); -/// let v = exp.ind_sample(rand::task_rng()); -/// println!("{} is from a Exp(2) distribution", v); -/// } -/// ``` -pub struct Exp { - /// `lambda` stored as `1/lambda`, since this is what we scale by. - priv lambda_inverse: f64 -} - -impl Exp { - /// Construct a new `Exp` with the given shape parameter - /// `lambda`. Fails if `lambda <= 0`. - pub fn new(lambda: f64) -> Exp { - assert!(lambda > 0.0, "Exp::new called with `lambda` <= 0"); - Exp { lambda_inverse: 1.0 / lambda } - } -} - -impl Sample for Exp { - fn sample(&mut self, rng: &mut R) -> f64 { self.ind_sample(rng) } -} -impl IndependentSample for Exp { - fn ind_sample(&self, rng: &mut R) -> f64 { - (*rng.gen::()) * self.lambda_inverse - } -} - #[cfg(test)] mod tests { use rand::*; use super::*; - use iter::range; use option::{Some, None}; struct ConstRand(uint); @@ -449,42 +282,6 @@ fn test_rand_sample() { assert_eq!(*rand_sample.sample(&mut task_rng()), 0); assert_eq!(*rand_sample.ind_sample(&mut task_rng()), 0); } - - #[test] - fn test_normal() { - let mut norm = Normal::new(10.0, 10.0); - let mut rng = task_rng(); - for _ in range(0, 1000) { - norm.sample(&mut rng); - norm.ind_sample(&mut rng); - } - } - #[test] - #[should_fail] - fn test_normal_invalid_sd() { - Normal::new(10.0, -1.0); - } - - #[test] - fn test_exp() { - let mut exp = Exp::new(10.0); - let mut rng = task_rng(); - for _ in range(0, 1000) { - assert!(exp.sample(&mut rng) >= 0.0); - assert!(exp.ind_sample(&mut rng) >= 0.0); - } - } - #[test] - #[should_fail] - fn test_exp_invalid_lambda_zero() { - Exp::new(0.0); - } - #[test] - #[should_fail] - fn test_exp_invalid_lambda_neg() { - Exp::new(-10.0); - } - #[test] fn test_weighted_choice() { // this makes assumptions about the internal implementation of @@ -556,38 +353,3 @@ fn test_weighted_choice_weight_overflows() { Weighted { weight: 1, item: 3 }]); } } - -#[cfg(test)] -mod bench { - use extra::test::BenchHarness; - use rand::{XorShiftRng, RAND_BENCH_N}; - use super::*; - use iter::range; - use option::{Some, None}; - use mem::size_of; - - #[bench] - fn rand_normal(bh: &mut BenchHarness) { - let mut rng = XorShiftRng::new(); - let mut normal = Normal::new(-2.71828, 3.14159); - - bh.iter(|| { - for _ in range(0, RAND_BENCH_N) { - normal.sample(&mut rng); - } - }); - bh.bytes = size_of::() as u64 * RAND_BENCH_N; - } - #[bench] - fn rand_exp(bh: &mut BenchHarness) { - let mut rng = XorShiftRng::new(); - let mut exp = Exp::new(2.71828 * 3.14159); - - bh.iter(|| { - for _ in range(0, RAND_BENCH_N) { - exp.sample(&mut rng); - } - }); - bh.bytes = size_of::() as u64 * RAND_BENCH_N; - } -} diff --git a/src/libstd/rand/distributions/normal.rs b/src/libstd/rand/distributions/normal.rs new file mode 100644 index 00000000000..19ee9006cbb --- /dev/null +++ b/src/libstd/rand/distributions/normal.rs @@ -0,0 +1,155 @@ +// Copyright 2013 The Rust Project Developers. See the COPYRIGHT +// file at the top-level directory of this distribution and at +// http://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 normal distribution. + +use rand::{Rng, Rand, Open01}; +use rand::distributions::{ziggurat, ziggurat_tables, Sample, IndependentSample}; + +/// A wrapper around an `f64` to generate N(0, 1) random numbers +/// (a.k.a. a standard normal, or Gaussian). +/// +/// See `Normal` for the general normal distribution. That this has to +/// be unwrapped before use as an `f64` (using either `*` or +/// `cast::transmute` is safe). +/// +/// Implemented via the ZIGNOR variant[1] of the Ziggurat method. +/// +/// [1]: Jurgen A. Doornik (2005). [*An Improved Ziggurat Method to +/// Generate Normal Random +/// Samples*](http://www.doornik.com/research/ziggurat.pdf). Nuffield +/// College, Oxford +pub struct StandardNormal(f64); + +impl Rand for StandardNormal { + fn rand(rng: &mut R) -> StandardNormal { + #[inline] + fn pdf(x: f64) -> f64 { + ((-x*x/2.0) as f64).exp() + } + #[inline] + fn zero_case(rng: &mut R, u: f64) -> f64 { + // compute a random number in the tail by hand + + // strange initial conditions, because the loop is not + // do-while, so the condition should be true on the first + // run, they get overwritten anyway (0 < 1, so these are + // good). + let mut x = 1.0f64; + let mut y = 0.0f64; + + while -2.0 * y < x * x { + let x_ = *rng.gen::>(); + let y_ = *rng.gen::>(); + + x = x_.ln() / ziggurat_tables::ZIG_NORM_R; + y = y_.ln(); + } + + if u < 0.0 { x - ziggurat_tables::ZIG_NORM_R } else { ziggurat_tables::ZIG_NORM_R - x } + } + + StandardNormal(ziggurat( + rng, + true, // this is symmetric + &ziggurat_tables::ZIG_NORM_X, + &ziggurat_tables::ZIG_NORM_F, + pdf, zero_case)) + } +} + +/// The normal distribution `N(mean, std_dev**2)`. +/// +/// This uses the ZIGNOR variant of the Ziggurat method, see +/// `StandardNormal` for more details. +/// +/// # Example +/// +/// ```rust +/// use std::rand; +/// use std::rand::distributions::{Normal, IndependentSample}; +/// +/// fn main() { +/// // mean 2, standard deviation 3 +/// let normal = Normal::new(2.0, 3.0); +/// let v = normal.ind_sample(&mut rand::task_rng()); +/// println!("{} is from a N(2, 9) distribution", v) +/// } +/// ``` +pub struct Normal { + priv mean: f64, + priv std_dev: f64 +} + +impl Normal { + /// Construct a new `Normal` distribution with the given mean and + /// standard deviation. Fails if `std_dev < 0`. + pub fn new(mean: f64, std_dev: f64) -> Normal { + assert!(std_dev >= 0.0, "Normal::new called with `std_dev` < 0"); + Normal { + mean: mean, + std_dev: std_dev + } + } +} +impl Sample for Normal { + fn sample(&mut self, rng: &mut R) -> f64 { self.ind_sample(rng) } +} +impl IndependentSample for Normal { + fn ind_sample(&self, rng: &mut R) -> f64 { + self.mean + self.std_dev * (*rng.gen::()) + } +} + +#[cfg(test)] +mod tests { + use rand::*; + use super::*; + use iter::range; + use option::{Some, None}; + + #[test] + fn test_normal() { + let mut norm = Normal::new(10.0, 10.0); + let mut rng = task_rng(); + for _ in range(0, 1000) { + norm.sample(&mut rng); + norm.ind_sample(&mut rng); + } + } + #[test] + #[should_fail] + fn test_normal_invalid_sd() { + Normal::new(10.0, -1.0); + } +} + +#[cfg(test)] +mod bench { + use extra::test::BenchHarness; + use rand::{XorShiftRng, RAND_BENCH_N}; + use super::*; + use iter::range; + use option::{Some, None}; + use mem::size_of; + + #[bench] + fn rand_normal(bh: &mut BenchHarness) { + let mut rng = XorShiftRng::new(); + let mut normal = Normal::new(-2.71828, 3.14159); + + bh.iter(|| { + for _ in range(0, RAND_BENCH_N) { + normal.sample(&mut rng); + } + }); + bh.bytes = size_of::() as u64 * RAND_BENCH_N; + } +}