(index<- ) ./libstd/rand/distributions.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 //! Sampling from random distributions
12
13 // Some implementations use the Ziggurat method
14 // https://en.wikipedia.org/wiki/Ziggurat_algorithm
15 //
16 // The version used here is ZIGNOR [Doornik 2005, "An Improved
17 // Ziggurat Method to Generate Normal Random Samples"] which is slower
18 // (about double, it generates an extra random number) than the
19 // canonical version [Marsaglia & Tsang 2000, "The Ziggurat Method for
20 // Generating Random Variables"], but more robust. If one wanted, one
21 // could implement VIZIGNOR the ZIGNOR paper for more speed.
22
23 use num;
24 use rand::{Rng,Rand};
25
26 mod ziggurat_tables;
27
28 // inlining should mean there is no performance penalty for this
29 #[inline]
30 fn ziggurat<R:Rng>(rng: &mut R,
31 center_u: bool,
32 X: ziggurat_tables::ZigTable,
33 F: ziggurat_tables::ZigTable,
34 F_DIFF: ziggurat_tables::ZigTable,
35 pdf: &'static fn(f64) -> f64, // probability density function
36 zero_case: &'static fn(&mut R, f64) -> f64) -> f64 {
37 loop {
38 let u = if center_u {2.0 * rng.gen() - 1.0} else {rng.gen()};
39 let i: uint = rng.gen::<uint>() & 0xff;
40 let x = u * X[i];
41
42 let test_x = if center_u {num::abs(x)} else {x};
43
44 // algebraically equivalent to |u| < X[i+1]/X[i] (or u < X[i+1]/X[i])
45 if test_x < X[i + 1] {
46 return x;
47 }
48 if i == 0 {
49 return zero_case(rng, u);
50 }
51 // algebraically equivalent to f1 + DRanU()*(f0 - f1) < 1
52 if F[i+1] + F_DIFF[i+1] * rng.gen() < pdf(x) {
53 return x;
54 }
55 }
56 }
57
58 /// A wrapper around an `f64` to generate N(0, 1) random numbers (a.k.a. a
59 /// standard normal, or Gaussian). Multiplying the generated values by the
60 /// desired standard deviation `sigma` then adding the desired mean `mu` will
61 /// give N(mu, sigma^2) distributed random numbers.
62 ///
63 /// Note that this has to be unwrapped before use as an `f64` (using either
64 /// `*` or `cast::transmute` is safe).
65 ///
66 /// # Example
67 ///
68 /// ```
69 /// use std::rand::distributions::StandardNormal;
70 ///
71 /// fn main() {
72 /// let normal = 2.0 + (*rand::random::<StandardNormal>()) * 3.0;
73 /// println!("{} is from a N(2, 9) distribution", normal)
74 /// }
75 /// ```
76 pub struct StandardNormal(f64);
77
78 impl Rand for StandardNormal {
79 fn rand<R:Rng>(rng: &mut R) -> StandardNormal {
80 #[inline]
81 fn pdf(x: f64) -> f64 {
82 ((-x*x/2.0) as f64).exp()
83 }
84 #[inline]
85 fn zero_case<R:Rng>(rng: &mut R, u: f64) -> f64 {
86 // compute a random number in the tail by hand
87
88 // strange initial conditions, because the loop is not
89 // do-while, so the condition should be true on the first
90 // run, they get overwritten anyway (0 < 1, so these are
91 // good).
92 let mut x = 1.0f64;
93 let mut y = 0.0f64;
94
95 // FIXME #7755: infinities?
96 while -2.0 * y < x * x {
97 x = rng.gen::<f64>().ln() / ziggurat_tables::ZIG_NORM_R;
98 y = rng.gen::<f64>().ln();
99 }
100
101 if u < 0.0 { x - ziggurat_tables::ZIG_NORM_R } else { ziggurat_tables::ZIG_NORM_R - x }
102 }
103
104 StandardNormal(ziggurat(
105 rng,
106 true, // this is symmetric
107 &ziggurat_tables::ZIG_NORM_X,
108 &ziggurat_tables::ZIG_NORM_F, &ziggurat_tables::ZIG_NORM_F_DIFF,
109 pdf, zero_case))
110 }
111 }
112
113 /// A wrapper around an `f64` to generate Exp(1) random numbers. Dividing by
114 /// the desired rate `lambda` will give Exp(lambda) distributed random
115 /// numbers.
116 ///
117 /// Note that this has to be unwrapped before use as an `f64` (using either
118 /// `*` or `cast::transmute` is safe).
119 ///
120 /// # Example
121 ///
122 /// ```
123 /// use std::rand::distributions::Exp1;
124 ///
125 /// fn main() {
126 /// let exp2 = (*rand::random::<Exp1>()) * 0.5;
127 /// println!("{} is from a Exp(2) distribution", exp2);
128 /// }
129 /// ```
130 pub struct Exp1(f64);
131
132 // This could be done via `-rng.gen::<f64>().ln()` but that is slower.
133 impl Rand for Exp1 {
134 #[inline]
135 fn rand<R:Rng>(rng: &mut R) -> Exp1 {
136 #[inline]
137 fn pdf(x: f64) -> f64 {
138 (-x).exp()
139 }
140 #[inline]
141 fn zero_case<R:Rng>(rng: &mut R, _u: f64) -> f64 {
142 ziggurat_tables::ZIG_EXP_R - rng.gen::<f64>().ln()
143 }
144
145 Exp1(ziggurat(rng, false,
146 &ziggurat_tables::ZIG_EXP_X,
147 &ziggurat_tables::ZIG_EXP_F, &ziggurat_tables::ZIG_EXP_F_DIFF,
148 pdf, zero_case))
149 }
150 }