# How does numpy generate samples from a beta distribution?

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?

The code for numpy.random.beta is found at legacy-distributions.c at the time of this writing.

1. 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.
2. 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.
• These are some interesting facts about numpy.random.binomial and numpy.random.* overall, but I don't think the answers the question. What method does numpy.random.beta use to generate random values? – Sycorax Dec 23 '20 at 21:07
• @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. – Peter O. Dec 23 '20 at 21:10
• Thank you for editing your post to clarify how numpy.random.beta works. – Sycorax Dec 23 '20 at 21:52
• I meant to say "both 1 or less" rather than "both 0". – Peter O. Dec 23 '20 at 23:17
• No, it's not the Ziggurat algorithm, but an algorithm from 2000 called "A simple method for generating gamma variables". – Peter O. Dec 24 '20 at 4:15