src/librand/distributions/normal.rs

   1 // Copyright 2013 The Rust Project Developers. See the COPYRIGHT
   2 // file at the top-level directory of this distribution and at
   3 // http://rust-lang.org/COPYRIGHT.
   4 //
   5 // Licensed under the Apache License, Version 2.0 <LICENSE-APACHE or
   6 // http://www.apache.org/licenses/LICENSE-2.0> or the MIT license
   7 // <LICENSE-MIT or http://opensource.org/licenses/MIT>, at your
   8 // option. This file may not be copied, modified, or distributed
   9 // except according to those terms.
  10
  11 //! The normal and derived distributions.
  12
  13 use core::num::Float;
  14
  15 use {Rng, Rand, Open01};
  16 use distributions::{ziggurat, ziggurat_tables, Sample, IndependentSample};
  17
  18 /// A wrapper around an `f64` to generate N(0, 1) random numbers
  19 /// (a.k.a.  a standard normal, or Gaussian).
  20 ///
  21 /// See `Normal` for the general normal distribution. That this has to
  22 /// be unwrapped before use as an `f64` (using either `*` or
  23 /// `mem::transmute` is safe).
  24 ///
  25 /// Implemented via the ZIGNOR variant[1] of the Ziggurat method.
  26 ///
  27 /// [1]: Jurgen A. Doornik (2005). [*An Improved Ziggurat Method to
  28 /// Generate Normal Random
  29 /// Samples*](http://www.doornik.com/research/ziggurat.pdf). Nuffield
  30 /// College, Oxford
  31 pub struct StandardNormal(pub f64);
  32
  33 impl Rand for StandardNormal {
  34     fn rand<R:Rng>(rng: &mut R) -> StandardNormal {
  35         #[inline]
  36         fn pdf(x: f64) -> f64 {
  37             (-x*x/2.0).exp()
  38         }
  39         #[inline]
  40         fn zero_case<R:Rng>(rng: &mut R, u: f64) -> f64 {
  41             // compute a random number in the tail by hand
  42
  43             // strange initial conditions, because the loop is not
  44             // do-while, so the condition should be true on the first
  45             // run, they get overwritten anyway (0 < 1, so these are
  46             // good).
  47             let mut x = 1.0f64;
  48             let mut y = 0.0f64;
  49
  50             while -2.0 * y < x * x {
  51                 let Open01(x_) = rng.gen::<Open01<f64>>();
  52                 let Open01(y_) = rng.gen::<Open01<f64>>();
  53
  54                 x = x_.ln() / ziggurat_tables::ZIG_NORM_R;
  55                 y = y_.ln();
  56             }
  57
  58             if u < 0.0 { x - ziggurat_tables::ZIG_NORM_R } else { ziggurat_tables::ZIG_NORM_R - x }
  59         }
  60
  61         StandardNormal(ziggurat(
  62             rng,
  63             true, // this is symmetric
  64             &ziggurat_tables::ZIG_NORM_X,
  65             &ziggurat_tables::ZIG_NORM_F,
  66             pdf, zero_case))
  67     }
  68 }
  69
  70 /// The normal distribution `N(mean, std_dev**2)`.
  71 ///
  72 /// This uses the ZIGNOR variant of the Ziggurat method, see
  73 /// `StandardNormal` for more details.
  74 ///
  75 /// # Example
  76 ///
  77 /// ```rust
  78 /// use std::rand;
  79 /// use std::rand::distributions::{Normal, IndependentSample};
  80 ///
  81 /// // mean 2, standard deviation 3
  82 /// let normal = Normal::new(2.0, 3.0);
  83 /// let v = normal.ind_sample(&mut rand::task_rng());
  84 /// println!("{} is from a N(2, 9) distribution", v)
  85 /// ```
  86 pub struct Normal {
  87     mean: f64,
  88     std_dev: f64,
  89 }
  90
  91 impl Normal {
  92     /// Construct a new `Normal` distribution with the given mean and
  93     /// standard deviation.
  94     ///
  95     /// # Panics
  96     ///
  97     /// Panics if `std_dev < 0`.
  98     pub fn new(mean: f64, std_dev: f64) -> Normal {
  99         assert!(std_dev >= 0.0, "Normal::new called with `std_dev` < 0");
 100         Normal {
 101             mean: mean,
 102             std_dev: std_dev
 103         }
 104     }
 105 }
 106 impl Sample<f64> for Normal {
 107     fn sample<R: Rng>(&mut self, rng: &mut R) -> f64 { self.ind_sample(rng) }
 108 }
 109 impl IndependentSample<f64> for Normal {
 110     fn ind_sample<R: Rng>(&self, rng: &mut R) -> f64 {
 111         let StandardNormal(n) = rng.gen::<StandardNormal>();
 112         self.mean + self.std_dev * n
 113     }
 114 }
 115
 116
 117 /// The log-normal distribution `ln N(mean, std_dev**2)`.
 118 ///
 119 /// If `X` is log-normal distributed, then `ln(X)` is `N(mean,
 120 /// std_dev**2)` distributed.
 121 ///
 122 /// # Example
 123 ///
 124 /// ```rust
 125 /// use std::rand;
 126 /// use std::rand::distributions::{LogNormal, IndependentSample};
 127 ///
 128 /// // mean 2, standard deviation 3
 129 /// let log_normal = LogNormal::new(2.0, 3.0);
 130 /// let v = log_normal.ind_sample(&mut rand::task_rng());
 131 /// println!("{} is from an ln N(2, 9) distribution", v)
 132 /// ```
 133 pub struct LogNormal {
 134     norm: Normal
 135 }
 136
 137 impl LogNormal {
 138     /// Construct a new `LogNormal` distribution with the given mean
 139     /// and standard deviation.
 140     ///
 141     /// # Panics
 142     ///
 143     /// Panics if `std_dev < 0`.
 144     pub fn new(mean: f64, std_dev: f64) -> LogNormal {
 145         assert!(std_dev >= 0.0, "LogNormal::new called with `std_dev` < 0");
 146         LogNormal { norm: Normal::new(mean, std_dev) }
 147     }
 148 }
 149 impl Sample<f64> for LogNormal {
 150     fn sample<R: Rng>(&mut self, rng: &mut R) -> f64 { self.ind_sample(rng) }
 151 }
 152 impl IndependentSample<f64> for LogNormal {
 153     fn ind_sample<R: Rng>(&self, rng: &mut R) -> f64 {
 154         self.norm.ind_sample(rng).exp()
 155     }
 156 }
 157
 158 #[cfg(test)]
 159 mod tests {
 160     use std::prelude::*;
 161
 162     use distributions::{Sample, IndependentSample};
 163     use super::{Normal, LogNormal};
 164
 165     #[test]
 166     fn test_normal() {
 167         let mut norm = Normal::new(10.0, 10.0);
 168         let mut rng = ::test::rng();
 169         for _ in range(0u, 1000) {
 170             norm.sample(&mut rng);
 171             norm.ind_sample(&mut rng);
 172         }
 173     }
 174     #[test]
 175     #[should_fail]
 176     fn test_normal_invalid_sd() {
 177         Normal::new(10.0, -1.0);
 178     }
 179
 180
 181     #[test]
 182     fn test_log_normal() {
 183         let mut lnorm = LogNormal::new(10.0, 10.0);
 184         let mut rng = ::test::rng();
 185         for _ in range(0u, 1000) {
 186             lnorm.sample(&mut rng);
 187             lnorm.ind_sample(&mut rng);
 188         }
 189     }
 190     #[test]
 191     #[should_fail]
 192     fn test_log_normal_invalid_sd() {
 193         LogNormal::new(10.0, -1.0);
 194     }
 195 }
 196
 197 #[cfg(test)]
 198 mod bench {
 199     extern crate test;
 200     use std::prelude::*;
 201     use self::test::Bencher;
 202     use std::mem::size_of;
 203     use distributions::{Sample};
 204     use super::Normal;
 205
 206     #[bench]
 207     fn rand_normal(b: &mut Bencher) {
 208         let mut rng = ::test::weak_rng();
 209         let mut normal = Normal::new(-2.71828, 3.14159);
 210
 211         b.iter(|| {
 212             for _ in range(0, ::RAND_BENCH_N) {
 213                 normal.sample(&mut rng);
 214             }
 215         });
 216         b.bytes = size_of::<f64>() as u64 * ::RAND_BENCH_N;
 217     }
 218 }