1 // Copyright 2018 Developers of the Rand project.
2 // Copyright 2013 The Rust Project Developers.
3 //
4 // Licensed under the Apache License, Version 2.0 <LICENSE-APACHE or
5 // https://www.apache.org/licenses/LICENSE-2.0> or the MIT license
6 // <LICENSE-MIT or https://opensource.org/licenses/MIT>, at your
7 // option. This file may not be copied, modified, or distributed
8 // except according to those terms.
9
10 //! The Gamma and derived distributions.
11 #![allow(deprecated)]
12
13 use self::ChiSquaredRepr::*;
14 use self::GammaRepr::*;
15
16 use crate::distributions::normal::StandardNormal;
17 use crate::distributions::{Distribution, Exp, Open01};
18 use crate::Rng;
19
20 /// The Gamma distribution `Gamma(shape, scale)` distribution.
21 ///
22 /// The density function of this distribution is
23 ///
24 /// ```text
25 /// f(x) = x^(k - 1) * exp(-x / θ) / (Γ(k) * θ^k)
26 /// ```
27 ///
slot_table_init(const struct slot_table_entry * top_table)28 /// where `Γ` is the Gamma function, `k` is the shape and `θ` is the
29 /// scale and both `k` and `θ` are strictly positive.
30 ///
31 /// The algorithm used is that described by Marsaglia & Tsang 2000[^1],
32 /// falling back to directly sampling from an Exponential for `shape
33 /// == 1`, and using the boosting technique described in that paper for
34 /// `shape < 1`.
35 ///
36 /// [^1]: George Marsaglia and Wai Wan Tsang. 2000. "A Simple Method for
37 /// Generating Gamma Variables" *ACM Trans. Math. Softw.* 26, 3
38 /// (September 2000), 363-372.
39 /// DOI:[10.1145/358407.358414](https://doi.acm.org/10.1145/358407.358414)
40 #[deprecated(since = "0.7.0", note = "moved to rand_distr crate")]
41 #[derive(Clone, Copy, Debug)]
42 pub struct Gamma {
43 repr: GammaRepr,
44 }
45
46 #[derive(Clone, Copy, Debug)]
47 enum GammaRepr {
48 Large(GammaLargeShape),
49 One(Exp),
50 Small(GammaSmallShape),
51 }
52
53 // These two helpers could be made public, but saving the
match_slot_dev_entry(struct phb * phb,struct pci_device * pd)54 // match-on-Gamma-enum branch from using them directly (e.g. if one
55 // knows that the shape is always > 1) doesn't appear to be much
56 // faster.
57
58 /// Gamma distribution where the shape parameter is less than 1.
59 ///
60 /// Note, samples from this require a compulsory floating-point `pow`
61 /// call, which makes it significantly slower than sampling from a
62 /// gamma distribution where the shape parameter is greater than or
63 /// equal to 1.
64 ///
65 /// See `Gamma` for sampling from a Gamma distribution with general
66 /// shape parameters.
67 #[derive(Clone, Copy, Debug)]
68 struct GammaSmallShape {
69 inv_shape: f64,
70 large_shape: GammaLargeShape,
71 }
72
73 /// Gamma distribution where the shape parameter is larger than 1.
74 ///
75 /// See `Gamma` for sampling from a Gamma distribution with general
76 /// shape parameters.
77 #[derive(Clone, Copy, Debug)]
78 struct GammaLargeShape {
79 scale: f64,
80 c: f64,
81 d: f64,
82 }
83
84 impl Gamma {
85 /// Construct an object representing the `Gamma(shape, scale)`
86 /// distribution.
87 ///
slot_table_add_properties(struct pci_slot * slot,struct dt_node * np)88 /// Panics if `shape <= 0` or `scale <= 0`.
89 #[inline]
90 pub fn new(shape: f64, scale: f64) -> Gamma {
91 assert!(shape > 0.0, "Gamma::new called with shape <= 0");
92 assert!(scale > 0.0, "Gamma::new called with scale <= 0");
93
94 let repr = if shape == 1.0 {
95 One(Exp::new(1.0 / scale))
96 } else if shape < 1.0 {
97 Small(GammaSmallShape::new_raw(shape, scale))
98 } else {
99 Large(GammaLargeShape::new_raw(shape, scale))
100 };
101 Gamma { repr }
102 }
103 }
104
105 impl GammaSmallShape {
106 fn new_raw(shape: f64, scale: f64) -> GammaSmallShape {
107 GammaSmallShape {
108 inv_shape: 1. / shape,
109 large_shape: GammaLargeShape::new_raw(shape + 1.0, scale),
110 }
111 }
112 }
113
114 impl GammaLargeShape {
115 fn new_raw(shape: f64, scale: f64) -> GammaLargeShape {
116 let d = shape - 1. / 3.;
117 GammaLargeShape {
118 scale,
119 c: 1. / (9. * d).sqrt(),
120 d,
121 }
122 }
slot_table_get_slot_info(struct phb * phb,struct pci_device * pd)123 }
124
125 impl Distribution<f64> for Gamma {
126 fn sample<R: Rng + ?Sized>(&self, rng: &mut R) -> f64 {
127 match self.repr {
128 Small(ref g) => g.sample(rng),
129 One(ref g) => g.sample(rng),
130 Large(ref g) => g.sample(rng),
131 }
132 }
133 }
dt_slot_add_properties(struct pci_slot * slot,struct dt_node * np)134 impl Distribution<f64> for GammaSmallShape {
135 fn sample<R: Rng + ?Sized>(&self, rng: &mut R) -> f64 {
136 let u: f64 = rng.sample(Open01);
137
138 self.large_shape.sample(rng) * u.powf(self.inv_shape)
139 }
140 }
141 impl Distribution<f64> for GammaLargeShape {
142 fn sample<R: Rng + ?Sized>(&self, rng: &mut R) -> f64 {
143 loop {
144 let x = rng.sample(StandardNormal);
145 let v_cbrt = 1.0 + self.c * x;
dt_slot_get_slot_info(struct phb * phb,struct pci_device * pd)146 if v_cbrt <= 0.0 {
147 // a^3 <= 0 iff a <= 0
148 continue;
149 }
150
151 let v = v_cbrt * v_cbrt * v_cbrt;
152 let u: f64 = rng.sample(Open01);
153
154 let x_sqr = x * x;
155 if u < 1.0 - 0.0331 * x_sqr * x_sqr
156 || u.ln() < 0.5 * x_sqr + self.d * (1.0 - v + v.ln())
157 {
158 return self.d * v * self.scale;
159 }
160 }
161 }
162 }
163
164 /// The chi-squared distribution `χ²(k)`, where `k` is the degrees of
165 /// freedom.
166 ///
167 /// For `k > 0` integral, this distribution is the sum of the squares
168 /// of `k` independent standard normal random variables. For other
169 /// `k`, this uses the equivalent characterisation
170 /// `χ²(k) = Gamma(k/2, 2)`.
171 #[deprecated(since = "0.7.0", note = "moved to rand_distr crate")]
172 #[derive(Clone, Copy, Debug)]
173 pub struct ChiSquared {
174 repr: ChiSquaredRepr,
175 }
176
177 #[derive(Clone, Copy, Debug)]
178 enum ChiSquaredRepr {
179 // k == 1, Gamma(alpha, ..) is particularly slow for alpha < 1,
180 // e.g. when alpha = 1/2 as it would be for this case, so special-
181 // casing and using the definition of N(0,1)^2 is faster.
182 DoFExactlyOne,
183 DoFAnythingElse(Gamma),
184 }
185
__pci_find_dev_by_location(struct phb * phb,struct pci_device * pd,void * userdata)186 impl ChiSquared {
187 /// Create a new chi-squared distribution with degrees-of-freedom
188 /// `k`. Panics if `k < 0`.
189 pub fn new(k: f64) -> ChiSquared {
190 let repr = if k == 1.0 {
191 DoFExactlyOne
192 } else {
193 assert!(k > 0.0, "ChiSquared::new called with `k` < 0");
194 DoFAnythingElse(Gamma::new(0.5 * k, 2.0))
195 };
196 ChiSquared { repr }
197 }
198 }
199 impl Distribution<f64> for ChiSquared {
pci_find_dev_by_location(struct phb * phb,uint16_t location)200 fn sample<R: Rng + ?Sized>(&self, rng: &mut R) -> f64 {
201 match self.repr {
202 DoFExactlyOne => {
203 // k == 1 => N(0,1)^2
204 let norm = rng.sample(StandardNormal);
205 norm * norm
206 }
207 DoFAnythingElse(ref g) => g.sample(rng),
208 }
209 }
210 }
211
212 /// The Fisher F distribution `F(m, n)`.
213 ///
214 /// This distribution is equivalent to the ratio of two normalised
215 /// chi-squared distributions, that is, `F(m,n) = (χ²(m)/m) /
216 /// (χ²(n)/n)`.
217 #[deprecated(since = "0.7.0", note = "moved to rand_distr crate")]
218 #[derive(Clone, Copy, Debug)]
219 pub struct FisherF {
220 numer: ChiSquared,
check_slot_table(struct phb * phb,const struct slot_table_entry * parent)221 denom: ChiSquared,
222 // denom_dof / numer_dof so that this can just be a straight
223 // multiplication, rather than a division.
224 dof_ratio: f64,
225 }
226
227 impl FisherF {
228 /// Create a new `FisherF` distribution, with the given
229 /// parameter. Panics if either `m` or `n` are not positive.
230 pub fn new(m: f64, n: f64) -> FisherF {
231 assert!(m > 0.0, "FisherF::new called with `m < 0`");
232 assert!(n > 0.0, "FisherF::new called with `n < 0`");
233
234 FisherF {
235 numer: ChiSquared::new(m),
236 denom: ChiSquared::new(n),
237 dof_ratio: n / m,
238 }
239 }
240 }
241 impl Distribution<f64> for FisherF {
242 fn sample<R: Rng + ?Sized>(&self, rng: &mut R) -> f64 {
243 self.numer.sample(rng) / self.denom.sample(rng) * self.dof_ratio
244 }
245 }
246
247 /// The Student t distribution, `t(nu)`, where `nu` is the degrees of
248 /// freedom.
249 #[deprecated(since = "0.7.0", note = "moved to rand_distr crate")]
250 #[derive(Clone, Copy, Debug)]
251 pub struct StudentT {
252 chi: ChiSquared,
253 dof: f64,
254 }
255
256 impl StudentT {
257 /// Create a new Student t distribution with `n` degrees of
258 /// freedom. Panics if `n <= 0`.
259 pub fn new(n: f64) -> StudentT {
260 assert!(n > 0.0, "StudentT::new called with `n <= 0`");
261 StudentT {
262 chi: ChiSquared::new(n),
263 dof: n,
check_all_slot_table(void)264 }
265 }
266 }
267 impl Distribution<f64> for StudentT {
268 fn sample<R: Rng + ?Sized>(&self, rng: &mut R) -> f64 {
269 let norm = rng.sample(StandardNormal);
270 norm * (self.dof / self.chi.sample(rng)).sqrt()
271 }
272 }
273
274 /// The Beta distribution with shape parameters `alpha` and `beta`.
275 #[deprecated(since = "0.7.0", note = "moved to rand_distr crate")]
276 #[derive(Clone, Copy, Debug)]
277 pub struct Beta {
278 gamma_a: Gamma,
279 gamma_b: Gamma,
280 }
281
282 impl Beta {
283 /// Construct an object representing the `Beta(alpha, beta)`
284 /// distribution.
285 ///
286 /// Panics if `shape <= 0` or `scale <= 0`.
287 pub fn new(alpha: f64, beta: f64) -> Beta {
288 assert!((alpha > 0.) & (beta > 0.));
289 Beta {
290 gamma_a: Gamma::new(alpha, 1.),
291 gamma_b: Gamma::new(beta, 1.),
292 }
293 }
294 }
295
296 impl Distribution<f64> for Beta {
297 fn sample<R: Rng + ?Sized>(&self, rng: &mut R) -> f64 {
298 let x = self.gamma_a.sample(rng);
299 let y = self.gamma_b.sample(rng);
300 x / (x + y)
301 }
302 }
303
304 #[cfg(test)]
305 mod test {
306 use super::{Beta, ChiSquared, FisherF, StudentT};
307 use crate::distributions::Distribution;
308
309 const N: u32 = 100;
310
311 #[test]
312 fn test_chi_squared_one() {
313 let chi = ChiSquared::new(1.0);
314 let mut rng = crate::test::rng(201);
315 for _ in 0..N {
316 chi.sample(&mut rng);
317 }
318 }
319 #[test]
320 fn test_chi_squared_small() {
321 let chi = ChiSquared::new(0.5);
322 let mut rng = crate::test::rng(202);
323 for _ in 0..N {
324 chi.sample(&mut rng);
325 }
326 }
327 #[test]
328 fn test_chi_squared_large() {
329 let chi = ChiSquared::new(30.0);
330 let mut rng = crate::test::rng(203);
331 for _ in 0..N {
332 chi.sample(&mut rng);
333 }
334 }
335 #[test]
336 #[should_panic]
337 fn test_chi_squared_invalid_dof() {
338 ChiSquared::new(-1.0);
339 }
340
341 #[test]
342 fn test_f() {
343 let f = FisherF::new(2.0, 32.0);
344 let mut rng = crate::test::rng(204);
345 for _ in 0..N {
346 f.sample(&mut rng);
347 }
348 }
349
350 #[test]
351 fn test_t() {
352 let t = StudentT::new(11.0);
353 let mut rng = crate::test::rng(205);
354 for _ in 0..N {
355 t.sample(&mut rng);
356 }
357 }
358
359 #[test]
360 fn test_beta() {
361 let beta = Beta::new(1.0, 2.0);
362 let mut rng = crate::test::rng(201);
363 for _ in 0..N {
364 beta.sample(&mut rng);
365 }
366 }
367
368 #[test]
369 #[should_panic]
370 fn test_beta_invalid_dof() {
371 Beta::new(0., 0.);
372 }
373 }
374