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