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 /// 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 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 /// 88 /// Panics if `shape <= 0` or `scale <= 0`. 89 #[inline] new(shape: f64, scale: f64) -> Gamma90 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 { new_raw(shape: f64, scale: f64) -> GammaSmallShape106 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 { new_raw(shape: f64, scale: f64) -> GammaLargeShape115 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 } 123 } 124 125 impl Distribution<f64> for Gamma { sample<R: Rng + ?Sized>(&self, rng: &mut R) -> f64126 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 } 134 impl Distribution<f64> for GammaSmallShape { sample<R: Rng + ?Sized>(&self, rng: &mut R) -> f64135 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 { sample<R: Rng + ?Sized>(&self, rng: &mut R) -> f64142 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; 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 186 impl ChiSquared { 187 /// Create a new chi-squared distribution with degrees-of-freedom 188 /// `k`. Panics if `k < 0`. new(k: f64) -> ChiSquared189 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 { sample<R: Rng + ?Sized>(&self, rng: &mut R) -> f64200 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, 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. new(m: f64, n: f64) -> FisherF230 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 { sample<R: Rng + ?Sized>(&self, rng: &mut R) -> f64242 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`. new(n: f64) -> StudentT259 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, 264 } 265 } 266 } 267 impl Distribution<f64> for StudentT { sample<R: Rng + ?Sized>(&self, rng: &mut R) -> f64268 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`. new(alpha: f64, beta: f64) -> Beta287 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 { sample<R: Rng + ?Sized>(&self, rng: &mut R) -> f64297 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] test_chi_squared_one()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] test_chi_squared_small()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] test_chi_squared_large()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] test_chi_squared_invalid_dof()337 fn test_chi_squared_invalid_dof() { 338 ChiSquared::new(-1.0); 339 } 340 341 #[test] test_f()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] test_t()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] test_beta()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] test_beta_invalid_dof()370 fn test_beta_invalid_dof() { 371 Beta::new(0., 0.); 372 } 373 } 374