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