From 9d5639d11ba5d3ef4ca6be323d6604a5b9e0caaa Mon Sep 17 00:00:00 2001 From: Huon Wilson Date: Sat, 7 Dec 2013 22:14:20 +1100 Subject: [PATCH 1/5] std::rand: move normal and exponential to their own file. --- src/libstd/rand/distributions/exponential.rs | 141 +++++++++++ src/libstd/rand/distributions/gamma.rs | 3 +- src/libstd/rand/distributions/mod.rs | 248 +------------------ src/libstd/rand/distributions/normal.rs | 155 ++++++++++++ 4 files changed, 303 insertions(+), 244 deletions(-) create mode 100644 src/libstd/rand/distributions/exponential.rs create mode 100644 src/libstd/rand/distributions/normal.rs diff --git a/src/libstd/rand/distributions/exponential.rs b/src/libstd/rand/distributions/exponential.rs new file mode 100644 index 0000000000000..4244e6bacdbbd --- /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 e042874245967..61929a7c479c1 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 4778e81f95169..2a2932bdb4cf9 100644 --- a/src/libstd/rand/distributions/mod.rs +++ b/src/libstd/rand/distributions/mod.rs @@ -23,14 +23,18 @@ that do not need to record state. 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 @@ mod tests { 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 @@ mod tests { 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 0000000000000..19ee9006cbb84 --- /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; + } +} From 1d986de2488249763f730ba83fd3d8235391e74d Mon Sep 17 00:00:00 2001 From: Huon Wilson Date: Sat, 7 Dec 2013 22:22:55 +1100 Subject: [PATCH 2/5] std::rand: implement the log-normal distribution. --- src/libstd/rand/distributions/mod.rs | 2 +- src/libstd/rand/distributions/normal.rs | 58 ++++++++++++++++++++++++- 2 files changed, 58 insertions(+), 2 deletions(-) diff --git a/src/libstd/rand/distributions/mod.rs b/src/libstd/rand/distributions/mod.rs index 2a2932bdb4cf9..ee24d724731bc 100644 --- a/src/libstd/rand/distributions/mod.rs +++ b/src/libstd/rand/distributions/mod.rs @@ -28,7 +28,7 @@ use clone::Clone; pub use self::range::Range; pub use self::gamma::Gamma; -pub use self::normal::Normal; +pub use self::normal::{Normal, LogNormal}; pub use self::exponential::Exp; pub mod range; diff --git a/src/libstd/rand/distributions/normal.rs b/src/libstd/rand/distributions/normal.rs index 19ee9006cbb84..a779f4f60660e 100644 --- a/src/libstd/rand/distributions/normal.rs +++ b/src/libstd/rand/distributions/normal.rs @@ -8,7 +8,7 @@ // option. This file may not be copied, modified, or distributed // except according to those terms. -//! The normal distribution. +//! The normal and derived distributions. use rand::{Rng, Rand, Open01}; use rand::distributions::{ziggurat, ziggurat_tables, Sample, IndependentSample}; @@ -108,6 +108,46 @@ impl IndependentSample for Normal { } } + +/// The log-normal distribution `ln N(mean, std_dev**2)`. +/// +/// If `X` is log-normal distributed, then `ln(X)` is `N(mean, +/// std_dev**2)` distributed. +/// +/// # Example +/// +/// ```rust +/// use std::rand; +/// use std::rand::distributions::{LogNormal, IndependentSample}; +/// +/// fn main() { +/// // mean 2, standard deviation 3 +/// let log_normal = LogNormal::new(2.0, 3.0); +/// let v = normal.ind_sample(&mut rand::task_rng()); +/// println!("{} is from an ln N(2, 9) distribution", v) +/// } +/// ``` +pub struct LogNormal { + priv norm: Normal +} + +impl LogNormal { + /// Construct a new `LogNormal` distribution with the given mean + /// and standard deviation. Fails if `std_dev < 0`. + pub fn new(mean: f64, std_dev: f64) -> LogNormal { + assert!(std_dev >= 0.0, "LogNormal::new called with `std_dev` < 0"); + LogNormal { norm: Normal::new(mean, std_dev) } + } +} +impl Sample for LogNormal { + fn sample(&mut self, rng: &mut R) -> f64 { self.ind_sample(rng) } +} +impl IndependentSample for LogNormal { + fn ind_sample(&self, rng: &mut R) -> f64 { + self.norm.ind_sample(rng).exp() + } +} + #[cfg(test)] mod tests { use rand::*; @@ -129,6 +169,22 @@ mod tests { fn test_normal_invalid_sd() { Normal::new(10.0, -1.0); } + + + #[test] + fn test_log_normal() { + let mut lnorm = LogNormal::new(10.0, 10.0); + let mut rng = task_rng(); + for _ in range(0, 1000) { + lnorm.sample(&mut rng); + lnorm.ind_sample(&mut rng); + } + } + #[test] + #[should_fail] + fn test_log_normal_invalid_sd() { + LogNormal::new(10.0, -1.0); + } } #[cfg(test)] From 1ee42912e1c83106082012e9c23090934be4ae69 Mon Sep 17 00:00:00 2001 From: Huon Wilson Date: Sat, 7 Dec 2013 22:51:08 +1100 Subject: [PATCH 3/5] std::rand: implement the chi-squared distribution. --- src/libstd/rand/distributions/gamma.rs | 99 +++++++++++++++++++++++++- src/libstd/rand/distributions/mod.rs | 2 +- 2 files changed, 99 insertions(+), 2 deletions(-) diff --git a/src/libstd/rand/distributions/gamma.rs b/src/libstd/rand/distributions/gamma.rs index 61929a7c479c1..7d583771230a9 100644 --- a/src/libstd/rand/distributions/gamma.rs +++ b/src/libstd/rand/distributions/gamma.rs @@ -8,7 +8,7 @@ // option. This file may not be copied, modified, or distributed // except according to those terms. -//! The Gamma distribution. +//! The Gamma and derived distributions. use rand::{Rng, Open01}; use super::{IndependentSample, Sample, Exp}; @@ -169,6 +169,103 @@ impl IndependentSample for GammaLargeShape { } } +/// The chi-squared distribution `χ²(k)`, where `k` is the degrees of +/// freedom. +/// +/// For `k > 0` integral, this distribution is the sum of the squares +/// of `k` independent standard normal random variables. For other +/// `k`, this uses the equivalent characterisation `χ²(k) = Gamma(k/2, +/// 2)`. +/// +/// # Example +/// +/// ```rust +/// use std::rand; +/// use std::rand::distributions::{ChiSquared, IndependentSample}; +/// +/// fn main() { +/// let chi = ChiSquared::new(11.0); +/// let v = chi.ind_sample(&mut rand::task_rng()); +/// println!("{} is from a χ²(11) distribution", v) +/// } +/// ``` +pub enum ChiSquared { + // k == 1, Gamma(alpha, ..) is particularly slow for alpha < 1, + // e.g. when alpha = 1/2 as it would be for this case, so special- + // casing and using the definition of N(0,1)^2 is faster. + priv DoFExactlyOne, + priv DoFAnythingElse(Gamma) +} + +impl ChiSquared { + /// Create a new chi-squared distribution with degrees-of-freedom + /// `k`. Fails if `k < 0`. + pub fn new(k: f64) -> ChiSquared { + if k == 1.0 { + DoFExactlyOne + } else { + assert!(k > 0.0, "ChiSquared::new called with `k` < 0"); + DoFAnythingElse(Gamma::new(0.5 * k, 2.0)) + } + } +} +impl Sample for ChiSquared { + fn sample(&mut self, rng: &mut R) -> f64 { self.ind_sample(rng) } +} +impl IndependentSample for ChiSquared { + fn ind_sample(&self, rng: &mut R) -> f64 { + match *self { + DoFExactlyOne => { + // k == 1 => N(0,1)^2 + let norm = *rng.gen::(); + norm * norm + } + DoFAnythingElse(ref g) => g.ind_sample(rng) + } + } +} + +#[cfg(test)] +mod test { + use rand::*; + use super::*; + use iter::range; + use option::{Some, None}; + + #[test] + fn test_chi_squared_one() { + let mut chi = ChiSquared::new(1.0); + let mut rng = task_rng(); + for _ in range(0, 1000) { + chi.sample(&mut rng); + chi.ind_sample(&mut rng); + } + } + #[test] + fn test_chi_squared_small() { + let mut chi = ChiSquared::new(0.5); + let mut rng = task_rng(); + for _ in range(0, 1000) { + chi.sample(&mut rng); + chi.ind_sample(&mut rng); + } + } + #[test] + fn test_chi_squared_large() { + let mut chi = ChiSquared::new(30.0); + let mut rng = task_rng(); + for _ in range(0, 1000) { + chi.sample(&mut rng); + chi.ind_sample(&mut rng); + } + } + #[test] + #[should_fail] + fn test_log_normal_invalid_dof() { + ChiSquared::new(-1.0); + } +} + #[cfg(test)] mod bench { use super::*; diff --git a/src/libstd/rand/distributions/mod.rs b/src/libstd/rand/distributions/mod.rs index ee24d724731bc..4e789a28ad275 100644 --- a/src/libstd/rand/distributions/mod.rs +++ b/src/libstd/rand/distributions/mod.rs @@ -27,7 +27,7 @@ use rand::{Rng, Rand}; use clone::Clone; pub use self::range::Range; -pub use self::gamma::Gamma; +pub use self::gamma::{Gamma, ChiSquared}; pub use self::normal::{Normal, LogNormal}; pub use self::exponential::Exp; From 6155a1c98012f9565744dad3832660accd135b83 Mon Sep 17 00:00:00 2001 From: Huon Wilson Date: Sat, 7 Dec 2013 23:13:37 +1100 Subject: [PATCH 4/5] std::rand: implement the F distribution. --- src/libstd/rand/distributions/gamma.rs | 59 ++++++++++++++++++++++++++ src/libstd/rand/distributions/mod.rs | 2 +- 2 files changed, 60 insertions(+), 1 deletion(-) diff --git a/src/libstd/rand/distributions/gamma.rs b/src/libstd/rand/distributions/gamma.rs index 7d583771230a9..dec0d3922d326 100644 --- a/src/libstd/rand/distributions/gamma.rs +++ b/src/libstd/rand/distributions/gamma.rs @@ -225,6 +225,55 @@ impl IndependentSample for ChiSquared { } } +/// The Fisher F distribution `F(m, n)`. +/// +/// This distribution is equivalent to the ratio of two normalised +/// chi-squared distributions, that is, `F(m,n) = (χ²(m)/m) / +/// (χ²(n)/n)`. +/// +/// # Example +/// +/// ```rust +/// use std::rand; +/// use std::rand::distributions::{FisherF, IndependentSample}; +/// +/// fn main() { +/// let f = FisherF::new(2.0, 32.0); +/// let v = f.ind_sample(&mut rand::task_rng()); +/// println!("{} is from an F(2, 32) distribution", v) +/// } +/// ``` +pub struct FisherF { + priv numer: ChiSquared, + priv denom: ChiSquared, + // denom_dof / numer_dof so that this can just be a straight + // multiplication, rather than a division. + priv dof_ratio: f64, +} + +impl FisherF { + /// Create a new `FisherF` distribution, with the given + /// parameter. Fails if either `m` or `n` are not positive. + pub fn new(m: f64, n: f64) -> FisherF { + assert!(m > 0.0, "FisherF::new called with `m < 0`"); + assert!(n > 0.0, "FisherF::new called with `n < 0`"); + + FisherF { + numer: ChiSquared::new(m), + denom: ChiSquared::new(n), + dof_ratio: n / m + } + } +} +impl Sample for FisherF { + fn sample(&mut self, rng: &mut R) -> f64 { self.ind_sample(rng) } +} +impl IndependentSample for FisherF { + fn ind_sample(&self, rng: &mut R) -> f64 { + self.numer.ind_sample(rng) / self.denom.ind_sample(rng) * self.dof_ratio + } +} + #[cfg(test)] mod test { use rand::*; @@ -264,6 +313,16 @@ mod test { fn test_log_normal_invalid_dof() { ChiSquared::new(-1.0); } + + #[test] + fn test_f() { + let mut f = FisherF::new(2.0, 32.0); + let mut rng = task_rng(); + for _ in range(0, 1000) { + f.sample(&mut rng); + f.ind_sample(&mut rng); + } + } } #[cfg(test)] diff --git a/src/libstd/rand/distributions/mod.rs b/src/libstd/rand/distributions/mod.rs index 4e789a28ad275..261b77fb24562 100644 --- a/src/libstd/rand/distributions/mod.rs +++ b/src/libstd/rand/distributions/mod.rs @@ -27,7 +27,7 @@ use rand::{Rng, Rand}; use clone::Clone; pub use self::range::Range; -pub use self::gamma::{Gamma, ChiSquared}; +pub use self::gamma::{Gamma, ChiSquared, FisherF}; pub use self::normal::{Normal, LogNormal}; pub use self::exponential::Exp; From 705b705ba5182745e12652dc7eeeb10828f520d0 Mon Sep 17 00:00:00 2001 From: Huon Wilson Date: Sat, 7 Dec 2013 23:39:34 +1100 Subject: [PATCH 5/5] std::rand: implement the student t distribution. --- src/libstd/rand/distributions/gamma.rs | 51 ++++++++++++++++++++++++++ src/libstd/rand/distributions/mod.rs | 2 +- 2 files changed, 52 insertions(+), 1 deletion(-) diff --git a/src/libstd/rand/distributions/gamma.rs b/src/libstd/rand/distributions/gamma.rs index dec0d3922d326..ae7ff99af9242 100644 --- a/src/libstd/rand/distributions/gamma.rs +++ b/src/libstd/rand/distributions/gamma.rs @@ -274,6 +274,47 @@ impl IndependentSample for FisherF { } } +/// The Student t distribution, `t(nu)`, where `nu` is the degrees of +/// freedom. +/// +/// # Example +/// +/// ```rust +/// use std::rand; +/// use std::rand::distributions::{StudentT, IndependentSample}; +/// +/// fn main() { +/// let t = StudentT::new(11.0); +/// let v = t.ind_sample(&mut rand::task_rng()); +/// println!("{} is from a t(11) distribution", v) +/// } +/// ``` +pub struct StudentT { + priv chi: ChiSquared, + priv dof: f64 +} + +impl StudentT { + /// Create a new Student t distribution with `n` degrees of + /// freedom. Fails if `n <= 0`. + pub fn new(n: f64) -> StudentT { + assert!(n > 0.0, "StudentT::new called with `n <= 0`"); + StudentT { + chi: ChiSquared::new(n), + dof: n + } + } +} +impl Sample for StudentT { + fn sample(&mut self, rng: &mut R) -> f64 { self.ind_sample(rng) } +} +impl IndependentSample for StudentT { + fn ind_sample(&self, rng: &mut R) -> f64 { + let norm = *rng.gen::(); + norm * (self.dof / self.chi.ind_sample(rng)).sqrt() + } +} + #[cfg(test)] mod test { use rand::*; @@ -323,6 +364,16 @@ mod test { f.ind_sample(&mut rng); } } + + #[test] + fn test_t() { + let mut t = StudentT::new(11.0); + let mut rng = task_rng(); + for _ in range(0, 1000) { + t.sample(&mut rng); + t.ind_sample(&mut rng); + } + } } #[cfg(test)] diff --git a/src/libstd/rand/distributions/mod.rs b/src/libstd/rand/distributions/mod.rs index 261b77fb24562..a381ac35d30c3 100644 --- a/src/libstd/rand/distributions/mod.rs +++ b/src/libstd/rand/distributions/mod.rs @@ -27,7 +27,7 @@ use rand::{Rng, Rand}; use clone::Clone; pub use self::range::Range; -pub use self::gamma::{Gamma, ChiSquared, FisherF}; +pub use self::gamma::{Gamma, ChiSquared, FisherF, StudentT}; pub use self::normal::{Normal, LogNormal}; pub use self::exponential::Exp;