]>
Commit | Line | Data |
---|---|---|
1a4d82fc JJ |
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 | ||
32a655c1 SL |
13 | use core::fmt; |
14 | ||
a7813a04 | 15 | #[cfg(not(test))] // only necessary for no_std |
e9174d1e | 16 | use FloatMath; |
1a4d82fc | 17 | |
3157f602 XL |
18 | use {Open01, Rand, Rng}; |
19 | use distributions::{IndependentSample, Sample, ziggurat, ziggurat_tables}; | |
1a4d82fc JJ |
20 | |
21 | /// A wrapper around an `f64` to generate N(0, 1) random numbers | |
22 | /// (a.k.a. a standard normal, or Gaussian). | |
23 | /// | |
24 | /// See `Normal` for the general normal distribution. That this has to | |
25 | /// be unwrapped before use as an `f64` (using either `*` or | |
26 | /// `mem::transmute` is safe). | |
27 | /// | |
28 | /// Implemented via the ZIGNOR variant[1] of the Ziggurat method. | |
29 | /// | |
30 | /// [1]: Jurgen A. Doornik (2005). [*An Improved Ziggurat Method to | |
31 | /// Generate Normal Random | |
32 | /// Samples*](http://www.doornik.com/research/ziggurat.pdf). Nuffield | |
33 | /// College, Oxford | |
c34b1796 | 34 | #[derive(Copy, Clone)] |
1a4d82fc JJ |
35 | pub struct StandardNormal(pub f64); |
36 | ||
37 | impl Rand for StandardNormal { | |
b039eaaf | 38 | fn rand<R: Rng>(rng: &mut R) -> StandardNormal { |
1a4d82fc JJ |
39 | #[inline] |
40 | fn pdf(x: f64) -> f64 { | |
b039eaaf | 41 | (-x * x / 2.0).exp() |
1a4d82fc JJ |
42 | } |
43 | #[inline] | |
b039eaaf | 44 | fn zero_case<R: Rng>(rng: &mut R, u: f64) -> f64 { |
1a4d82fc JJ |
45 | // compute a random number in the tail by hand |
46 | ||
47 | // strange initial conditions, because the loop is not | |
48 | // do-while, so the condition should be true on the first | |
49 | // run, they get overwritten anyway (0 < 1, so these are | |
50 | // good). | |
51 | let mut x = 1.0f64; | |
52 | let mut y = 0.0f64; | |
53 | ||
54 | while -2.0 * y < x * x { | |
55 | let Open01(x_) = rng.gen::<Open01<f64>>(); | |
56 | let Open01(y_) = rng.gen::<Open01<f64>>(); | |
57 | ||
58 | x = x_.ln() / ziggurat_tables::ZIG_NORM_R; | |
59 | y = y_.ln(); | |
60 | } | |
61 | ||
b039eaaf SL |
62 | if u < 0.0 { |
63 | x - ziggurat_tables::ZIG_NORM_R | |
64 | } else { | |
65 | ziggurat_tables::ZIG_NORM_R - x | |
66 | } | |
1a4d82fc JJ |
67 | } |
68 | ||
b039eaaf SL |
69 | StandardNormal(ziggurat(rng, |
70 | true, // this is symmetric | |
71 | &ziggurat_tables::ZIG_NORM_X, | |
72 | &ziggurat_tables::ZIG_NORM_F, | |
73 | pdf, | |
74 | zero_case)) | |
1a4d82fc JJ |
75 | } |
76 | } | |
77 | ||
32a655c1 SL |
78 | impl fmt::Debug for StandardNormal { |
79 | fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result { | |
80 | f.debug_tuple("StandardNormal") | |
81 | .field(&self.0) | |
82 | .finish() | |
83 | } | |
84 | } | |
85 | ||
1a4d82fc JJ |
86 | /// The normal distribution `N(mean, std_dev**2)`. |
87 | /// | |
88 | /// This uses the ZIGNOR variant of the Ziggurat method, see | |
89 | /// `StandardNormal` for more details. | |
c34b1796 | 90 | #[derive(Copy, Clone)] |
1a4d82fc JJ |
91 | pub struct Normal { |
92 | mean: f64, | |
93 | std_dev: f64, | |
94 | } | |
95 | ||
96 | impl Normal { | |
97 | /// Construct a new `Normal` distribution with the given mean and | |
98 | /// standard deviation. | |
99 | /// | |
100 | /// # Panics | |
101 | /// | |
102 | /// Panics if `std_dev < 0`. | |
103 | pub fn new(mean: f64, std_dev: f64) -> Normal { | |
104 | assert!(std_dev >= 0.0, "Normal::new called with `std_dev` < 0"); | |
105 | Normal { | |
3b2f2976 XL |
106 | mean, |
107 | std_dev, | |
1a4d82fc JJ |
108 | } |
109 | } | |
110 | } | |
32a655c1 | 111 | |
1a4d82fc | 112 | impl Sample<f64> for Normal { |
b039eaaf SL |
113 | fn sample<R: Rng>(&mut self, rng: &mut R) -> f64 { |
114 | self.ind_sample(rng) | |
115 | } | |
1a4d82fc | 116 | } |
32a655c1 | 117 | |
1a4d82fc JJ |
118 | impl IndependentSample<f64> for Normal { |
119 | fn ind_sample<R: Rng>(&self, rng: &mut R) -> f64 { | |
120 | let StandardNormal(n) = rng.gen::<StandardNormal>(); | |
121 | self.mean + self.std_dev * n | |
122 | } | |
123 | } | |
124 | ||
32a655c1 SL |
125 | impl fmt::Debug for Normal { |
126 | fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result { | |
127 | f.debug_struct("Normal") | |
128 | .field("mean", &self.mean) | |
129 | .field("std_dev", &self.std_dev) | |
130 | .finish() | |
131 | } | |
132 | } | |
133 | ||
1a4d82fc JJ |
134 | |
135 | /// The log-normal distribution `ln N(mean, std_dev**2)`. | |
136 | /// | |
137 | /// If `X` is log-normal distributed, then `ln(X)` is `N(mean, | |
138 | /// std_dev**2)` distributed. | |
c34b1796 | 139 | #[derive(Copy, Clone)] |
1a4d82fc | 140 | pub struct LogNormal { |
b039eaaf | 141 | norm: Normal, |
1a4d82fc JJ |
142 | } |
143 | ||
144 | impl LogNormal { | |
145 | /// Construct a new `LogNormal` distribution with the given mean | |
146 | /// and standard deviation. | |
147 | /// | |
148 | /// # Panics | |
149 | /// | |
150 | /// Panics if `std_dev < 0`. | |
151 | pub fn new(mean: f64, std_dev: f64) -> LogNormal { | |
152 | assert!(std_dev >= 0.0, "LogNormal::new called with `std_dev` < 0"); | |
153 | LogNormal { norm: Normal::new(mean, std_dev) } | |
154 | } | |
155 | } | |
32a655c1 | 156 | |
1a4d82fc | 157 | impl Sample<f64> for LogNormal { |
b039eaaf SL |
158 | fn sample<R: Rng>(&mut self, rng: &mut R) -> f64 { |
159 | self.ind_sample(rng) | |
160 | } | |
1a4d82fc | 161 | } |
32a655c1 | 162 | |
1a4d82fc JJ |
163 | impl IndependentSample<f64> for LogNormal { |
164 | fn ind_sample<R: Rng>(&self, rng: &mut R) -> f64 { | |
165 | self.norm.ind_sample(rng).exp() | |
166 | } | |
167 | } | |
168 | ||
32a655c1 SL |
169 | impl fmt::Debug for LogNormal { |
170 | fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result { | |
171 | f.debug_struct("LogNormal") | |
172 | .field("norm", &self.norm) | |
173 | .finish() | |
174 | } | |
175 | } | |
176 | ||
1a4d82fc JJ |
177 | #[cfg(test)] |
178 | mod tests { | |
3157f602 XL |
179 | use distributions::{IndependentSample, Sample}; |
180 | use super::{LogNormal, Normal}; | |
1a4d82fc JJ |
181 | |
182 | #[test] | |
183 | fn test_normal() { | |
184 | let mut norm = Normal::new(10.0, 10.0); | |
185 | let mut rng = ::test::rng(); | |
85aaf69f | 186 | for _ in 0..1000 { |
1a4d82fc JJ |
187 | norm.sample(&mut rng); |
188 | norm.ind_sample(&mut rng); | |
189 | } | |
190 | } | |
191 | #[test] | |
c34b1796 | 192 | #[should_panic] |
1a4d82fc JJ |
193 | fn test_normal_invalid_sd() { |
194 | Normal::new(10.0, -1.0); | |
195 | } | |
196 | ||
197 | ||
198 | #[test] | |
199 | fn test_log_normal() { | |
200 | let mut lnorm = LogNormal::new(10.0, 10.0); | |
201 | let mut rng = ::test::rng(); | |
85aaf69f | 202 | for _ in 0..1000 { |
1a4d82fc JJ |
203 | lnorm.sample(&mut rng); |
204 | lnorm.ind_sample(&mut rng); | |
205 | } | |
206 | } | |
207 | #[test] | |
c34b1796 | 208 | #[should_panic] |
1a4d82fc JJ |
209 | fn test_log_normal_invalid_sd() { |
210 | LogNormal::new(10.0, -1.0); | |
211 | } | |
212 | } | |
213 | ||
214 | #[cfg(test)] | |
215 | mod bench { | |
216 | extern crate test; | |
1a4d82fc JJ |
217 | use self::test::Bencher; |
218 | use std::mem::size_of; | |
b039eaaf | 219 | use distributions::Sample; |
1a4d82fc JJ |
220 | use super::Normal; |
221 | ||
222 | #[bench] | |
223 | fn rand_normal(b: &mut Bencher) { | |
224 | let mut rng = ::test::weak_rng(); | |
225 | let mut normal = Normal::new(-2.71828, 3.14159); | |
226 | ||
227 | b.iter(|| { | |
85aaf69f | 228 | for _ in 0..::RAND_BENCH_N { |
1a4d82fc JJ |
229 | normal.sample(&mut rng); |
230 | } | |
231 | }); | |
232 | b.bytes = size_of::<f64>() as u64 * ::RAND_BENCH_N; | |
233 | } | |
234 | } |