src/librand/distributions/exponential.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 exponential distribution.
  12
  13 use std::num::Float;
  14 use {Rng, Rand};
  15 use distributions::{ziggurat, ziggurat_tables, Sample, IndependentSample};
  16
  17 /// A wrapper around an `f64` to generate Exp(1) random numbers.
  18 ///
  19 /// See `Exp` for the general exponential distribution.Note that this
  20  // has to be unwrapped before use as an `f64` (using either
  21 /// `*` or `cast::transmute` is safe).
  22 ///
  23 /// Implemented via the ZIGNOR variant[1] of the Ziggurat method. The
  24 /// exact description in the paper was adjusted to use tables for the
  25 /// exponential distribution rather than normal.
  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 Exp1(pub f64);
  32
  33 // This could be done via `-rng.gen::<f64>().ln()` but that is slower.
  34 impl Rand for Exp1 {
  35     #[inline]
  36     fn rand<R:Rng>(rng: &mut R) -> Exp1 {
  37         #[inline]
  38         fn pdf(x: f64) -> f64 {
  39             (-x).exp()
  40         }
  41         #[inline]
  42         fn zero_case<R:Rng>(rng: &mut R, _u: f64) -> f64 {
  43             ziggurat_tables::ZIG_EXP_R - rng.gen::<f64>().ln()
  44         }
  45
  46         Exp1(ziggurat(rng, false,
  47                       &ziggurat_tables::ZIG_EXP_X,
  48                       &ziggurat_tables::ZIG_EXP_F,
  49                       pdf, zero_case))
  50     }
  51 }
  52
  53 /// The exponential distribution `Exp(lambda)`.
  54 ///
  55 /// This distribution has density function: `f(x) = lambda *
  56 /// exp(-lambda * x)` for `x > 0`.
  57 ///
  58 /// # Example
  59 ///
  60 /// ```rust
  61 /// use rand::distributions::{Exp, IndependentSample};
  62 ///
  63 /// let exp = Exp::new(2.0);
  64 /// let v = exp.ind_sample(&mut rand::task_rng());
  65 /// println!("{} is from a Exp(2) distribution", v);
  66 /// ```
  67 pub struct Exp {
  68     /// `lambda` stored as `1/lambda`, since this is what we scale by.
  69     lambda_inverse: f64
  70 }
  71
  72 impl Exp {
  73     /// Construct a new `Exp` with the given shape parameter
  74     /// `lambda`. Fails if `lambda <= 0`.
  75     pub fn new(lambda: f64) -> Exp {
  76         assert!(lambda > 0.0, "Exp::new called with `lambda` <= 0");
  77         Exp { lambda_inverse: 1.0 / lambda }
  78     }
  79 }
  80
  81 impl Sample<f64> for Exp {
  82     fn sample<R: Rng>(&mut self, rng: &mut R) -> f64 { self.ind_sample(rng) }
  83 }
  84 impl IndependentSample<f64> for Exp {
  85     fn ind_sample<R: Rng>(&self, rng: &mut R) -> f64 {
  86         let Exp1(n) = rng.gen::<Exp1>();
  87         n * self.lambda_inverse
  88     }
  89 }
  90
  91 #[cfg(test)]
  92 mod test {
  93     use distributions::{Sample, IndependentSample};
  94     use {Rng, task_rng};
  95     use super::Exp;
  96
  97     #[test]
  98     fn test_exp() {
  99         let mut exp = Exp::new(10.0);
 100         let mut rng = task_rng();
 101         for _ in range(0, 1000) {
 102             assert!(exp.sample(&mut rng) >= 0.0);
 103             assert!(exp.ind_sample(&mut rng) >= 0.0);
 104         }
 105     }
 106     #[test]
 107     #[should_fail]
 108     fn test_exp_invalid_lambda_zero() {
 109         Exp::new(0.0);
 110     }
 111     #[test]
 112     #[should_fail]
 113     fn test_exp_invalid_lambda_neg() {
 114         Exp::new(-10.0);
 115     }
 116 }
 117
 118 #[cfg(test)]
 119 mod bench {
 120     extern crate test;
 121     use self::test::BenchHarness;
 122     use std::mem::size_of;
 123     use {XorShiftRng, RAND_BENCH_N};
 124     use super::Exp;
 125     use distributions::Sample;
 126
 127     #[bench]
 128     fn rand_exp(bh: &mut BenchHarness) {
 129         let mut rng = XorShiftRng::new().unwrap();
 130         let mut exp = Exp::new(2.71828 * 3.14159);
 131
 132         bh.iter(|| {
 133             for _ in range(0, RAND_BENCH_N) {
 134                 exp.sample(&mut rng);
 135             }
 136         });
 137         bh.bytes = size_of::<f64>() as u64 * RAND_BENCH_N;
 138     }
 139 }