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