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 Gamma and derived distributions.
12
13 use std::num::Float;
14 use {Rng, Open01};
15 use super::normal::StandardNormal;
16 use super::{IndependentSample, Sample, Exp};
17
18 /// The Gamma distribution `Gamma(shape, scale)` distribution.
19 ///
20 /// The density function of this distribution is
21 ///
22 /// ```notrust
23 /// f(x) = x^(k - 1) * exp(-x / θ) / (Î(k) * θ^k)
24 /// ```
25 ///
26 /// where `Î` is the Gamma function, `k` is the shape and `θ` is the
27 /// scale and both `k` and `θ` are strictly positive.
28 ///
29 /// The algorithm used is that described by Marsaglia & Tsang 2000[1],
30 /// falling back to directly sampling from an Exponential for `shape
31 /// == 1`, and using the boosting technique described in [1] for
32 /// `shape < 1`.
33 ///
34 /// # Example
35 ///
36 /// ```rust
37 /// use rand::distributions::{IndependentSample, Gamma};
38 ///
39 /// let gamma = Gamma::new(2.0, 5.0);
40 /// let v = gamma.ind_sample(&mut rand::task_rng());
41 /// println!("{} is from a Gamma(2, 5) distribution", v);
42 /// ```
43 ///
44 /// [1]: George Marsaglia and Wai Wan Tsang. 2000. "A Simple Method
45 /// for Generating Gamma Variables" *ACM Trans. Math. Softw.* 26, 3
46 /// (September 2000),
47 /// 363-372. DOI:[10.1145/358407.358414](http://doi.acm.org/10.1145/358407.358414)
48 pub struct Gamma {
49 repr: GammaRepr,
50 }
51
52 enum GammaRepr {
53 Large(GammaLargeShape),
54 One(Exp),
55 Small(GammaSmallShape)
56 }
57
58 // These two helpers could be made public, but saving the
59 // match-on-Gamma-enum branch from using them directly (e.g. if one
60 // knows that the shape is always > 1) doesn't appear to be much
61 // faster.
62
63 /// Gamma distribution where the shape parameter is less than 1.
64 ///
65 /// Note, samples from this require a compulsory floating-point `pow`
66 /// call, which makes it significantly slower than sampling from a
67 /// gamma distribution where the shape parameter is greater than or
68 /// equal to 1.
69 ///
70 /// See `Gamma` for sampling from a Gamma distribution with general
71 /// shape parameters.
72 struct GammaSmallShape {
73 inv_shape: f64,
74 large_shape: GammaLargeShape
75 }
76
77 /// Gamma distribution where the shape parameter is larger than 1.
78 ///
79 /// See `Gamma` for sampling from a Gamma distribution with general
80 /// shape parameters.
81 struct GammaLargeShape {
82 shape: f64,
83 scale: f64,
84 c: f64,
85 d: f64
86 }
87
88 impl Gamma {
89 /// Construct an object representing the `Gamma(shape, scale)`
90 /// distribution.
91 ///
92 /// Fails if `shape <= 0` or `scale <= 0`.
93 pub fn new(shape: f64, scale: f64) -> Gamma {
94 assert!(shape > 0.0, "Gamma::new called with shape <= 0");
95 assert!(scale > 0.0, "Gamma::new called with scale <= 0");
96
97 let repr = match shape {
98 1.0 => One(Exp::new(1.0 / scale)),
99 0.0 .. 1.0 => Small(GammaSmallShape::new_raw(shape, scale)),
100 _ => Large(GammaLargeShape::new_raw(shape, scale))
101 };
102 Gamma { repr: repr }
103 }
104 }
105
106 impl GammaSmallShape {
107 fn new_raw(shape: f64, scale: f64) -> GammaSmallShape {
108 GammaSmallShape {
109 inv_shape: 1. / shape,
110 large_shape: GammaLargeShape::new_raw(shape + 1.0, scale)
111 }
112 }
113 }
114
115 impl GammaLargeShape {
116 fn new_raw(shape: f64, scale: f64) -> GammaLargeShape {
117 let d = shape - 1. / 3.;
118 GammaLargeShape {
119 shape: shape,
120 scale: scale,
121 c: 1. / (9. * d).sqrt(),
122 d: d
123 }
124 }
125 }
126
127 impl Sample<f64> for Gamma {
128 fn sample<R: Rng>(&mut self, rng: &mut R) -> f64 { self.ind_sample(rng) }
129 }
130 impl Sample<f64> for GammaSmallShape {
131 fn sample<R: Rng>(&mut self, rng: &mut R) -> f64 { self.ind_sample(rng) }
132 }
133 impl Sample<f64> for GammaLargeShape {
134 fn sample<R: Rng>(&mut self, rng: &mut R) -> f64 { self.ind_sample(rng) }
135 }
136
137 impl IndependentSample<f64> for Gamma {
138 fn ind_sample<R: Rng>(&self, rng: &mut R) -> f64 {
139 match self.repr {
140 Small(ref g) => g.ind_sample(rng),
141 One(ref g) => g.ind_sample(rng),
142 Large(ref g) => g.ind_sample(rng),
143 }
144 }
145 }
146 impl IndependentSample<f64> for GammaSmallShape {
147 fn ind_sample<R: Rng>(&self, rng: &mut R) -> f64 {
148 let Open01(u) = rng.gen::<Open01<f64>>();
149
150 self.large_shape.ind_sample(rng) * u.powf(self.inv_shape)
151 }
152 }
153 impl IndependentSample<f64> for GammaLargeShape {
154 fn ind_sample<R: Rng>(&self, rng: &mut R) -> f64 {
155 loop {
156 let StandardNormal(x) = rng.gen::<StandardNormal>();
157 let v_cbrt = 1.0 + self.c * x;
158 if v_cbrt <= 0.0 { // a^3 <= 0 iff a <= 0
159 continue
160 }
161
162 let v = v_cbrt * v_cbrt * v_cbrt;
163 let Open01(u) = rng.gen::<Open01<f64>>();
164
165 let x_sqr = x * x;
166 if u < 1.0 - 0.0331 * x_sqr * x_sqr ||
167 u.ln() < 0.5 * x_sqr + self.d * (1.0 - v + v.ln()) {
168 return self.d * v * self.scale
169 }
170 }
171 }
172 }
173
174 /// The chi-squared distribution `ϲ(k)`, where `k` is the degrees of
175 /// freedom.
176 ///
177 /// For `k > 0` integral, this distribution is the sum of the squares
178 /// of `k` independent standard normal random variables. For other
179 /// `k`, this uses the equivalent characterisation `ϲ(k) = Gamma(k/2,
180 /// 2)`.
181 ///
182 /// # Example
183 ///
184 /// ```rust
185 /// use rand::distributions::{ChiSquared, IndependentSample};
186 ///
187 /// let chi = ChiSquared::new(11.0);
188 /// let v = chi.ind_sample(&mut rand::task_rng());
189 /// println!("{} is from a ϲ(11) distribution", v)
190 /// ```
191 pub struct ChiSquared {
192 repr: ChiSquaredRepr,
193 }
194
195 enum ChiSquaredRepr {
196 // k == 1, Gamma(alpha, ..) is particularly slow for alpha < 1,
197 // e.g. when alpha = 1/2 as it would be for this case, so special-
198 // casing and using the definition of N(0,1)^2 is faster.
199 DoFExactlyOne,
200 DoFAnythingElse(Gamma),
201 }
202
203 impl ChiSquared {
204 /// Create a new chi-squared distribution with degrees-of-freedom
205 /// `k`. Fails if `k < 0`.
206 pub fn new(k: f64) -> ChiSquared {
207 let repr = if k == 1.0 {
208 DoFExactlyOne
209 } else {
210 assert!(k > 0.0, "ChiSquared::new called with `k` < 0");
211 DoFAnythingElse(Gamma::new(0.5 * k, 2.0))
212 };
213 ChiSquared { repr: repr }
214 }
215 }
216 impl Sample<f64> for ChiSquared {
217 fn sample<R: Rng>(&mut self, rng: &mut R) -> f64 { self.ind_sample(rng) }
218 }
219 impl IndependentSample<f64> for ChiSquared {
220 fn ind_sample<R: Rng>(&self, rng: &mut R) -> f64 {
221 match self.repr {
222 DoFExactlyOne => {
223 // k == 1 => N(0,1)^2
224 let StandardNormal(norm) = rng.gen::<StandardNormal>();
225 norm * norm
226 }
227 DoFAnythingElse(ref g) => g.ind_sample(rng)
228 }
229 }
230 }
231
232 /// The Fisher F distribution `F(m, n)`.
233 ///
234 /// This distribution is equivalent to the ratio of two normalised
235 /// chi-squared distributions, that is, `F(m,n) = (ϲ(m)/m) /
236 /// (ϲ(n)/n)`.
237 ///
238 /// # Example
239 ///
240 /// ```rust
241 /// use rand::distributions::{FisherF, IndependentSample};
242 ///
243 /// let f = FisherF::new(2.0, 32.0);
244 /// let v = f.ind_sample(&mut rand::task_rng());
245 /// println!("{} is from an F(2, 32) distribution", v)
246 /// ```
247 pub struct FisherF {
248 numer: ChiSquared,
249 denom: ChiSquared,
250 // denom_dof / numer_dof so that this can just be a straight
251 // multiplication, rather than a division.
252 dof_ratio: f64,
253 }
254
255 impl FisherF {
256 /// Create a new `FisherF` distribution, with the given
257 /// parameter. Fails if either `m` or `n` are not positive.
258 pub fn new(m: f64, n: f64) -> FisherF {
259 assert!(m > 0.0, "FisherF::new called with `m < 0`");
260 assert!(n > 0.0, "FisherF::new called with `n < 0`");
261
262 FisherF {
263 numer: ChiSquared::new(m),
264 denom: ChiSquared::new(n),
265 dof_ratio: n / m
266 }
267 }
268 }
269 impl Sample<f64> for FisherF {
270 fn sample<R: Rng>(&mut self, rng: &mut R) -> f64 { self.ind_sample(rng) }
271 }
272 impl IndependentSample<f64> for FisherF {
273 fn ind_sample<R: Rng>(&self, rng: &mut R) -> f64 {
274 self.numer.ind_sample(rng) / self.denom.ind_sample(rng) * self.dof_ratio
275 }
276 }
277
278 /// The Student t distribution, `t(nu)`, where `nu` is the degrees of
279 /// freedom.
280 ///
281 /// # Example
282 ///
283 /// ```rust
284 /// use rand::distributions::{StudentT, IndependentSample};
285 ///
286 /// let t = StudentT::new(11.0);
287 /// let v = t.ind_sample(&mut rand::task_rng());
288 /// println!("{} is from a t(11) distribution", v)
289 /// ```
290 pub struct StudentT {
291 chi: ChiSquared,
292 dof: f64
293 }
294
295 impl StudentT {
296 /// Create a new Student t distribution with `n` degrees of
297 /// freedom. Fails if `n <= 0`.
298 pub fn new(n: f64) -> StudentT {
299 assert!(n > 0.0, "StudentT::new called with `n <= 0`");
300 StudentT {
301 chi: ChiSquared::new(n),
302 dof: n
303 }
304 }
305 }
306 impl Sample<f64> for StudentT {
307 fn sample<R: Rng>(&mut self, rng: &mut R) -> f64 { self.ind_sample(rng) }
308 }
309 impl IndependentSample<f64> for StudentT {
310 fn ind_sample<R: Rng>(&self, rng: &mut R) -> f64 {
311 let StandardNormal(norm) = rng.gen::<StandardNormal>();
312 norm * (self.dof / self.chi.ind_sample(rng)).sqrt()
313 }
314 }
315
316 #[cfg(test)]
317 mod test {
318 use distributions::{Sample, IndependentSample};
319 use {Rng, task_rng};
320 use super::{ChiSquared, StudentT, FisherF};
321
322 #[test]
323 fn test_chi_squared_one() {
324 let mut chi = ChiSquared::new(1.0);
325 let mut rng = task_rng();
326 for _ in range(0, 1000) {
327 chi.sample(&mut rng);
328 chi.ind_sample(&mut rng);
329 }
330 }
331 #[test]
332 fn test_chi_squared_small() {
333 let mut chi = ChiSquared::new(0.5);
334 let mut rng = task_rng();
335 for _ in range(0, 1000) {
336 chi.sample(&mut rng);
337 chi.ind_sample(&mut rng);
338 }
339 }
340 #[test]
341 fn test_chi_squared_large() {
342 let mut chi = ChiSquared::new(30.0);
343 let mut rng = task_rng();
344 for _ in range(0, 1000) {
345 chi.sample(&mut rng);
346 chi.ind_sample(&mut rng);
347 }
348 }
349 #[test]
350 #[should_fail]
351 fn test_chi_squared_invalid_dof() {
352 ChiSquared::new(-1.0);
353 }
354
355 #[test]
356 fn test_f() {
357 let mut f = FisherF::new(2.0, 32.0);
358 let mut rng = task_rng();
359 for _ in range(0, 1000) {
360 f.sample(&mut rng);
361 f.ind_sample(&mut rng);
362 }
363 }
364
365 #[test]
366 fn test_t() {
367 let mut t = StudentT::new(11.0);
368 let mut rng = task_rng();
369 for _ in range(0, 1000) {
370 t.sample(&mut rng);
371 t.ind_sample(&mut rng);
372 }
373 }
374 }
375
376 #[cfg(test)]
377 mod bench {
378 extern crate test;
379 use self::test::Bencher;
380 use std::mem::size_of;
381 use distributions::IndependentSample;
382 use {XorShiftRng, RAND_BENCH_N};
383 use super::Gamma;
384
385
386 #[bench]
387 fn bench_gamma_large_shape(b: &mut Bencher) {
388 let gamma = Gamma::new(10., 1.0);
389 let mut rng = XorShiftRng::new().unwrap();
390
391 b.iter(|| {
392 for _ in range(0, RAND_BENCH_N) {
393 gamma.ind_sample(&mut rng);
394 }
395 });
396 b.bytes = size_of::<f64>() as u64 * RAND_BENCH_N;
397 }
398
399 #[bench]
400 fn bench_gamma_small_shape(b: &mut Bencher) {
401 let gamma = Gamma::new(0.1, 1.0);
402 let mut rng = XorShiftRng::new().unwrap();
403
404 b.iter(|| {
405 for _ in range(0, RAND_BENCH_N) {
406 gamma.ind_sample(&mut rng);
407 }
408 });
409 b.bytes = size_of::<f64>() as u64 * RAND_BENCH_N;
410 }
411 }
librand/distributions/gamma.rs:289:8-289:8 -struct- definition:
/// ```
pub struct StudentT {
chi: ChiSquared,
references:- 5299: assert!(n > 0.0, "StudentT::new called with `n <= 0`");
300: StudentT {
301: chi: ChiSquared::new(n),
--
305: }
306: impl Sample<f64> for StudentT {
307: fn sample<R: Rng>(&mut self, rng: &mut R) -> f64 { self.ind_sample(rng) }
308: }
309: impl IndependentSample<f64> for StudentT {
310: fn ind_sample<R: Rng>(&self, rng: &mut R) -> f64 {
librand/distributions/gamma.rs:190:8-190:8 -struct- definition:
/// ```
pub struct ChiSquared {
repr: ChiSquaredRepr,
references:- 8212: };
213: ChiSquared { repr: repr }
214: }
215: }
216: impl Sample<f64> for ChiSquared {
217: fn sample<R: Rng>(&mut self, rng: &mut R) -> f64 { self.ind_sample(rng) }
--
247: pub struct FisherF {
248: numer: ChiSquared,
249: denom: ChiSquared,
250: // denom_dof / numer_dof so that this can just be a straight
--
290: pub struct StudentT {
291: chi: ChiSquared,
292: dof: f64
librand/distributions/gamma.rs:246:8-246:8 -struct- definition:
/// ```
pub struct FisherF {
numer: ChiSquared,
references:- 5262: FisherF {
263: numer: ChiSquared::new(m),
--
268: }
269: impl Sample<f64> for FisherF {
270: fn sample<R: Rng>(&mut self, rng: &mut R) -> f64 { self.ind_sample(rng) }
271: }
272: impl IndependentSample<f64> for FisherF {
273: fn ind_sample<R: Rng>(&self, rng: &mut R) -> f64 {
librand/distributions/gamma.rs:80:22-80:22 -struct- definition:
/// shape parameters.
struct GammaLargeShape {
shape: f64,
references:- 7117: let d = shape - 1. / 3.;
118: GammaLargeShape {
119: shape: shape,
--
132: }
133: impl Sample<f64> for GammaLargeShape {
134: fn sample<R: Rng>(&mut self, rng: &mut R) -> f64 { self.ind_sample(rng) }
--
152: }
153: impl IndependentSample<f64> for GammaLargeShape {
154: fn ind_sample<R: Rng>(&self, rng: &mut R) -> f64 {
librand/distributions/gamma.rs:71:22-71:22 -struct- definition:
/// shape parameters.
struct GammaSmallShape {
inv_shape: f64,
references:- 6107: fn new_raw(shape: f64, scale: f64) -> GammaSmallShape {
108: GammaSmallShape {
109: inv_shape: 1. / shape,
--
129: }
130: impl Sample<f64> for GammaSmallShape {
131: fn sample<R: Rng>(&mut self, rng: &mut R) -> f64 { self.ind_sample(rng) }
--
145: }
146: impl IndependentSample<f64> for GammaSmallShape {
147: fn ind_sample<R: Rng>(&self, rng: &mut R) -> f64 {
librand/distributions/gamma.rs:47:83-47:83 -struct- definition:
/// 363-372. DOI:[10.1145/358407.358414](http://doi.acm.org/10.1145/358407.358414)
pub struct Gamma {
repr: GammaRepr,
references:- 6137: impl IndependentSample<f64> for Gamma {
138: fn ind_sample<R: Rng>(&self, rng: &mut R) -> f64 {
--
199: DoFExactlyOne,
200: DoFAnythingElse(Gamma),
201: }