std::rand: implement the Gamma distribution #10223
Conversation
|
The Gamma distribution is pretty fundamental, and it actually opens the way for extremely short implementations of a variety of other reasonably common distributions (for most, the actual sampling implementation is a line or two of arithmetic, so the largest part of the code will be the documentation & tests). Such distributions include:
I'm not sure how many of these are desired to be in libstd, but as some points of reference: C++11 (and Boost) has all of the above but Beta, and Python's stdlib has only Beta. (Personally, I'd say all, since they are literally some simple arithmetic done to one or two of the other distributions (once this lands) in |
Implements the [Gamma distribution](https://en.wikipedia.org/wiki/Gamma_distribution), using the algorithm described by Marsaglia & Tsang 2000[1]. I added tests checking that the mean and variance of this implementation is as expected for a range of values of the parameters in huonw/random-tests@5d87c00 (they pass locally, but obviously won't even build on Travis until this is merged). Also, moves `std::rand::distributions` to a subfolder, and performs a minor clean-up of the benchmarking (makes the number of iterations shared by the whole `std::rand` subtree). [1]: George Marsaglia and Wai Wan Tsang. 2000. "A Simple Method for Generating Gamma Variables" *ACM Trans. Math. Softw.* 26, 3 (September 2000), 363-372. DOI:[10.1145/358407.358414](http://doi.acm.org/10.1145/358407.358414).
| assert!(scale > 0.0, "Gamma::new called with scale <= 0"); | ||
|
|
||
| match shape { | ||
| 1.0 => One(Exp::new(1.0 / scale)), |
There was a problem hiding this comment.
Can you do equality testing on a floating point value like this?
There was a problem hiding this comment.
Yes, floats can be compared for equality, it just doesn't give exactly the answer one would expect in every situation, but it still gives sane & deterministic results, e.g.:
>>> 1 + 0.1 - 1 == 0.1
False
>>> 1.0 == 1.0
TrueHowever, in this case this optimisation is only valid when shape is exactly 1, i.e. it will be triggered by a call like Gamma::new(1.0, 2.0), so this is fine (doing an approximation like if (shape - 1.0).abs() < 1e-5 would mean that Gamma::new(1.000001, 1.0) does not have the Gamma(1.000001, 1) distribution, i.e. it is incorrect).
Implements the Gamma distribution, using the algorithm described by Marsaglia & Tsang 2000[1]. I added tests checking that the mean and variance of this implementation is as expected for a range of values of the parameters in huonw/random-tests@5d87c00 (they pass locally, but obviously won't even build on Travis until this is merged).
Also, moves
std::rand::distributionsto a subfolder, and performs a minor clean-up of the benchmarking (makes the number of iterations shared by the wholestd::randsubtree).[1]: George Marsaglia and Wai Wan Tsang. 2000. "A Simple Method for Generating Gamma Variables" ACM Trans. Math. Softw. 26, 3 (September 2000), 363-372. DOI:10.1145/358407.358414.