Access method for distribution-sampled random numbers generation in Numpy

Question

I am interested in knowing which method Numpy uses in numpy.random to generate random numbers distributed as beta, lognormal, Weibull distributions, among others. In the documentation they often discuss the shape of the distribution, but how can I access to the method? For example, in numpy.random.chisquare() the method is expained, but not in others.

Answer 1

Pseudorandom numbers start their life as uniformly distributed integers over the whole number representation space and then transformed into the desired shapes and ranges using https://en.wikipedia.org/wiki/Inverse_transform_sampling. While the descriptions usually suppose a uniform distribution of real numbers on [0,1], it is straightforward to imagine the transformation of getting the uniform distribution on [0,1] from one on the range [0,2^{64}]. In numpy, it is done like this:

def uniform_from_uint64(x):
    return (x >> np.uint64(11)) * (1.0 / 9007199254740992.0)


def uniform_from_uint32(x):
    out = np.empty(len(x) // 2)
    for i in range(0, len(x), 2):
        a = x[i] >> 5
        b = x[i + 1] >> 6
        out[i // 2] = (a * 67108864.0 + b) / 9007199254740992.0
    return out

Things get non-trivial in higher dimensions (as early as in 2D), especially when it comes to efficiency.