12. Random number generation

---

import numpy as np
import scipy.stats

import iqplot

import bokeh.io
bokeh.io.output_notebook()

Loading BokehJS ...

---

Random number generation (RNG) is the process by which a sequence of random numbers may be drawn. Of course, the numbers are not completely random for several reasons.

1. They are drawn from a probability distribution. The most common one is the uniform distribution on the domain \(0 \le x < 1\), i.e., random numbers between zero and one. (“Completely random” does not make sense because of the infinite magnitude of numbers.)

2. In most computer applications, including the ones we’ll use in bootcamp, the random numbers are actually pseudorandom. They depend entirely on an input seed and are then generated by a deterministic algorithm from that seed.

This is a bit academic. Let’s jump right in generating random numbers.

Uniform random numbers

Much of the random number generation functionality you will need is in the np.random module. Let’s start by generating random numbers from a Uniform distribution.

np.random.uniform(low=0, high=1, size=10)

array([0.61224477, 0.405525  , 0.89427367, 0.07970693, 0.35622669,
       0.57078386, 0.9480125 , 0.50490913, 0.03408578, 0.45752362])

The function uniform() in the np.random module generates random numbers on the interval [low, high) from a Uniform distribution. The size kwarg is how many random numbers you wish to generate, and is a kwarg in all of Numpy’s random number generators. The random numbers are returned as a NumPy array.

We can check to make sure it is appropriately drawing random numbers out of the uniform distribution by plotting the cumulative distribution function, just like we did last time. We’ll generate 1,000 random numbers and plot them along with the CDF of a Uniform distribution.

# Generate random numbers
x = np.random.uniform(low=0, high=1, size=1000)

# Plot the ECDF of randomly generated numbers
p = iqplot.ecdf(x, marker_kwargs={"fill_color": None},)

p.line(
    x=[0, 1], y=[0, 1], line_width=2, line_color="orange",
)

bokeh.io.show(p)

So, it looks like our random number generator is doing a good job.

Generating random numbers on the uniform interval is one of the most commonly used RNG applications. In fact, many of the other contexts of RNG are derived from draws from the uniform distribution. For example, you can simulate flipping a biased (unfair) coin by drawing from a Uniform distribution and then asking if the random number if less than the bias.

# Generate 20 random numbers on uniform interval
x = np.random.uniform(low=0, high=1, size=20)

# Make the coin flips (< 0.7 means we have a 70% chance of heads)
heads = x < 0.7

# Show which were heads, and count the number of heads
print(heads)
print("\nThere were", np.sum(heads), " heads.")

[ True  True  True False  True  True  True  True  True  True  True  True
  True  True  True  True False  True  True False]

There were 17  heads.

Choice of generator

As of version 1.23 of Numpy, the algorithm under the hood of calls to functions like np.random.uniform() is the Mersenne Twister Algorithm for generating random numbers. It is a very widely used and reliable method for generating random numbers. However, starting with version 1.17, the numpy.random module offers random number generators with better speed and statistical performance, including the PCG64 generator. Going forward, the preferred approach to doing random number generation is to first instantiate a generator of your choice, and then use its methods to generate numbers out of probability distributions.

Let’s set up a PCG64 generator, which is Numpy’s default.

rng = np.random.default_rng()

Now that we have the generator, we can use it to draw numbers out of distributions. The syntax is the same as before, except rng replaces np.random.

rng.uniform(low=0, high=1, size=20)

array([0.06104799, 0.77031758, 0.33252575, 0.6420329 , 0.32012886,
       0.571194  , 0.29060156, 0.09443274, 0.5389949 , 0.0704456 ,
       0.2604035 , 0.29588342, 0.65112538, 0.33054576, 0.36692791,
       0.97881176, 0.19486477, 0.69013615, 0.15856495, 0.316688  ])

Seeding random number generators

Now, just to demonstrate that random number generation is deterministic, we will explicitly seed the random number generator (which is usually seeded with a number representing the date/time to avoid repeats) to show that we get the same random numbers.

# Instantiate generator with a seed
rng = np.random.default_rng(seed=3252)

# Draw random numbers
rng.uniform(size=10)

array([0.18866535, 0.04418857, 0.02961285, 0.22083971, 0.43341773,
       0.13166813, 0.42112164, 0.43507845, 0.61380912, 0.30627603])

If we reinstantiate with the same seed, we get the same sequence of random numbers.

# Re-seed the RNG
rng = np.random.default_rng(seed=3252)

# Draw random numbers
rng.uniform(size=10)

array([0.18866535, 0.04418857, 0.02961285, 0.22083971, 0.43341773,
       0.13166813, 0.42112164, 0.43507845, 0.61380912, 0.30627603])

The random number sequence is exactly the same. If we choose a different seed, we get totally different random numbers.

rng = np.random.default_rng(seed=3253)
rng.uniform(size=10)

array([0.31390226, 0.73012457, 0.05800998, 0.01557021, 0.29825701,
       0.10106784, 0.06329107, 0.58614237, 0.52023168, 0.52779988])

If you are writing tests, it is often useful to seed the random number generator to get reproducible results. Otherwise, it is best to use the default seed, based on the date and time, so that you get a new set of random numbers in your applications each time do computations.

Drawing random numbers out of other distributions

We can also draw random numbers from other probability distributions besides Uniform. For example, say we wanted to draw random samples from a Normal distribution with mean μ and standard deviation σ.

# Set parameters
mu = 10
sigma = 1

# Draw 100000 random samples
x = rng.normal(mu, sigma, size=100000)

# Plot the histogram
p = iqplot.histogram(x, rug=False, density=True, y_axis_label="approximate PDF",)

bokeh.io.show(p)

It looks Normal, but, again, comparing the resulting ECDF is a better way to look at this. We’ll check out the ECDF with 1000 samples so as not to choke the browser with the display. I will also make use of the theoretical CDF for the Normal distribution available from the scipy.stats module.

# Compute theoretical CDF
x_theor = np.linspace(6, 14, 400)
y_theor = scipy.stats.norm.cdf(x_theor, mu, sigma)

# Plot the ECDF of randomly generated numbers
p = iqplot.ecdf(x, marker_kwargs={"fill_color": None},)

p.line(
    x=x_theor, y=y_theor, line_width=2, line_color="orange",
)

bokeh.io.show(p)

Yup, right on!

Selections from discrete distributions

The random numbers we have generated so far come from continuous probability distributions. We can also draw random numbers from discrete distributions. We already showed how we can do this for “coin flips,” but we can do it for other distributions as well. Say we wanted to draw from a Binomial distribution. We can use rng.binomial().

# Draw how many coin flips land heads in 10 files
rng.binomial(10, 0.5)

7

There are other discrete distributions we can draw from, such as Binomial, Geometric, Poisson, etc., and the documentation describes how to use them.

Choosing elements from an array

It is often useful to randomly choose elements from an existing array. (Actually, this is probably the functionality we will use the most, since it is used in bootstrapping.) The rng.choice() function does this. You equivalently could do this using rng.integers(), where the integers represent indices in the array, except rng.choice() has a great keyword argument, replace, which allows random draws with or without replacement. For example, say you had 50 samples that you wanted to send to a facility for analysis, but you can only afford to send 20. If we used rng.integers(), we might have a problem.

rng = np.random.default_rng(seed=126969234)
rng.integers(0, 51, size=20)

array([12, 31, 47, 26,  3,  5, 46, 49, 26, 38, 24, 17, 46, 26,  6, 17, 35,
        4, 13, 29])

Sample 17 was selected twice and sample 26 was selected thrice. This is not unexpected. We can use rng.choice() instead.

rng.choice(np.arange(51), size=20, replace=False)

array([27, 34,  0, 46,  2, 48, 35, 50, 40, 12, 28, 19, 37, 38, 11, 23, 45,
       15, 29, 32])

Now, because we chose replace=False, we do not get any repeats.

Generating random sequences

Because it works with selecting characters as well as numbers, we can use the rng.choice() function to generate random DNA sequences.

"".join(rng.choice(list("ATGC"), replace=True, size=70))

'CTTGTAAAGGATCTGCCCCGTGTCCCGCGCTAGGTCCGCACTCAAAGTATCTATGAACCGTCCAACCGAG'

Shuffling an array

Similarly, the rng.permutation() function is useful. It takes the entries in an array and shuffles them! Let’s shuffle a deck of cards.

rng.permutation(np.arange(53))

array([ 4,  1, 43, 46, 11, 21, 44, 18, 13, 49, 32, 47,  2, 28, 40, 52, 17,
       29, 27,  0, 14, 48, 37, 31, 20, 30, 36,  3, 45, 24, 50,  9, 22, 25,
        5, 12, 42, 19, 16, 39,  7, 34, 35, 10, 33, 41, 51, 26,  6,  8, 15,
       38, 23])

When do we need RNG?

Answer: VERY OFTEN! We will use random number generator extensively for the rest of the course as we explore probability distributions and perform hacker stats calculations. In particular, we will use rng.choice() extensively when doing bootstrap calculations and rng.permutation() when doing null significance hypothesis testing.

In many ways, probability is the language of biology. Molecular processes have energetics that are comparable to the thermal energy, which means they are always influenced by random thermal forces. The processes of the central dogma, including DNA replication, are no exceptions. This gives rise to random mutations, which are central to understanding how evolution works. If we want to understand them, it is often useful to use random number generators to model the processes.

Computing environment

%load_ext watermark
%watermark -v -p numpy,scipy,bokeh,iqplot,jupyterlab

Python implementation: CPython
Python version       : 3.12.4
IPython version      : 8.25.0

numpy     : 1.26.4
scipy     : 1.13.1
bokeh     : 3.4.1
iqplot    : 0.3.7
jupyterlab: 4.0.13