numpy lets you generate random samples from a beta distribution (or any other arbitrary distribution) with this API:
samples = np.random.beta(a,b, size=1000)
What is this doing beneath the hood? Inverse Transform Sampling?
$\begingroup$
$\endgroup$
1
-
$\begingroup$ I think you can find your answer here, github.com/numpy/numpy/blob/…. $\endgroup$– GijsDec 23, 2020 at 20:26
Add a comment
|
1 Answer
$\begingroup$
$\endgroup$
5
The code for numpy.random.beta is found at legacy-distributions.c at the time of this writing.
- When $a$ and $b$ are both 1 or less, then Jöhnk's beta generator is used (see page 418 of Non-Uniform Random Variate Generation), with a modification to avoid divisions by zero.
- Otherwise, it uses the formula $X/(X+Y)$ where $X$ and $Y$ are gamma($a$) and gamma($b$), respectively, generated by the legacy gamma generator. (In general, the gamma generator uses Marsaglia and Tsang's algorithm [2000] unless the parameter is less than 1.)
Notice that numpy.random.* functions (including numpy.random.beta) became legacy functions since NumPy 1.17 introduced a new system for pseudorandom number generation (see the NumPy RNG Policy). Thus, the implementation of numpy.random.beta is not expected to change for as long as numpy.random.* functions are still present in NumPy, and the beta generator used in the new RNG system may differ from the one presented here.
In fact, at the time of this writing, the new generator is slightly different from the legacy generator (see distributions.c at the time of this writing), and one of the reasons for introducing the new RNG system is to allow the non-uniform random generators to be improved without having to maintain backward compatibility.
REFERENCES:
- Marsaglia, G., Tsang, W.W., "A simple method for generating gamma variables", ACM Transactions on Mathematical Software 26(3), 2000.
-
$\begingroup$ These are some interesting facts about
numpy.random.binomialandnumpy.random.*overall, but I don't think the answers the question. What method doesnumpy.random.betause to generate random values? $\endgroup$– Sycorax ♦Dec 23, 2020 at 21:07 -
3$\begingroup$ @Sycorax: I already answered it: It uses Jöhnk's algorithm for a and b both 1 or less, and otherwise generates a beta random number with two gamma random numbers. $\endgroup$– Peter O.Dec 23, 2020 at 21:10
-
$\begingroup$ Thank you for editing your post to clarify how
numpy.random.betaworks. $\endgroup$– Sycorax ♦Dec 23, 2020 at 21:52 -
$\begingroup$ ah yes, I saw that Gamma trick on the bottom Beta's wikipedia page and was able to sample from a Beta using the ratio of Gammas. Is there a name for this technique? Then, the next logical question is how to generate samples from a Gamma. You answered that w/ the Marsalglia and Tsang's algorithm called Ziggurat. $\endgroup$– eggie5Dec 24, 2020 at 3:28
-
1$\begingroup$ No, it's not the Ziggurat algorithm, but an algorithm from 2000 called "A simple method for generating gamma variables". $\endgroup$– Peter O.Dec 24, 2020 at 4:15