12  Random number generation

| Download notebook

---

Random number generation (RNG) is the process by which a sequence of random numbers may be drawn. When using a computer to draw the random numbers, the numbers are not completely random. The notion of “completely random” is nonsensical because of the infinitude of numbers. Random numbers must be drawn from some probability distribution. Furthermore, in most computer applications, the random numbers are actually pseudorandom. They depend entirely on an input seed and are then generated by a deterministic algorithm from that seed.

So what do (pseudo)random number generators do? RNGs are capable of (approximately) drawing integers from a Discrete Uniform distribution. For example, NumPy’s default built-in generator, the PCG64 generator, generates 128 bit numbers, allowing for \(2^{128}\), or about \(10^{38}\), possible integers. Importantly, each draw of of a random integer is (approximately) independent of all others.

In practice, the drawn integers are converted to floating-point numbers (since a double floating-point number has far less than 128 bits) on the interval [0, 1) by dividing a generated random integer by \(2^{138}\). Effectively, then, the random number generators provide draws out of a Uniform distribution on the interval [0, 1).

To convert from random numbers on a Uniform distribution to random numbers from a nonuniform distribution, we need a transform. For many named distributions convenient transforms exist. For example, the Box-Muller transform is often used to get random draws from a Normal distribution. In the absence of a clever transform, we can use a distribution’s quantile function, also known as a percent-point function, which is the inverse CDF, \(F^{-1}(y)\). For example, the quantile function for the Exponential distribution is

\[ \begin{aligned} F^{-1}(p) = -\beta^{-1}\,\ln(1-p), \end{aligned} \tag{12.1}\]

where \(p\) is the value of the CDF, ranging from zero to one. We first draw \(p\) out of a Uniform distribution on [0, 1), and then compute \(F^{-1}(p)\) to get a draw from an Exponential distribution. A graphical illustration of using a quantile function to draw 50 random numbers out of Gamma(5, 2) is shown below.

import numpy as np
import scipy.stats as st

import bokeh.io
import bokeh.plotting
bokeh.io.output_notebook(hide_banner=True)

rng = np.random.default_rng(seed=12341234)
alpha = 5
beta = 2
y = np.linspace(0, 8, 400)
cdf = st.gamma.cdf(y, alpha, loc=0, scale=1 / beta)

udraws = rng.uniform(size=50)

p = bokeh.plotting.figure(
    width=300,
    height=200,
    x_axis_label="y",
    y_axis_label="F(y; 2, 5)",
    x_range=[0, 8],
    y_range=[0, 1],
)
p.xgrid.grid_line_color = None
p.ygrid.grid_line_color = None

p.line(y, cdf, line_width=2)

for u in udraws:
    x_vals = [0] + [st.gamma.ppf(u, alpha, loc=0, scale=1 / beta)] * 2
    y_vals = [u, u, 0]
    p.line(x_vals, y_vals, color="gray", line_width=0.5)

p.scatter(np.zeros_like(udraws), udraws, marker="x", color="black", line_width=0.5, size=7)
p.scatter(
    st.gamma.ppf(udraws, alpha, loc=0, scale=1 / beta),
    np.zeros_like(udraws),
    line_width=0.5,
    fill_color=None,
    line_color="black",
)

p.renderers[0].level = "overlay"
p.renderers[-1].level = "overlay"
p.renderers[-2].level = "overlay"

bokeh.io.show(p)

Each draw from a Uniform distribution, marked by an × on the vertical axis, is converted to a draw from a Gamma distribution, marked by ○ on the horizontal axis, by computing the inverse CDF. Because the draws of the random integers from which the draws from a Uniform distribution on \([0, 1)\) are independent, we get independent draws from the Gamma distribution.