13. Random number generation

---

# Colab setup ------------------
import os, sys, subprocess
if "google.colab" in sys.modules:
    cmd = "pip install --upgrade iqplot colorcet watermark"
    process = subprocess.Popen(cmd.split(), stdout=subprocess.PIPE, stderr=subprocess.PIPE)
    stdout, stderr = process.communicate()
# ------------------------------

import numpy as np
import pandas as pd
import scipy.stats

import iqplot

import bokeh.io
bokeh.io.output_notebook()

Loading BokehJS ...

---

Random number generation (RNG) is the process by which a string 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.

Uniform random numbers

This is a bit academic. Let’s jump right in generating 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.06549033, 0.559875  , 0.80693176, 0.42704107, 0.24706077,
       0.38824185, 0.78960865, 0.05427908, 0.45800661, 0.07814837])

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 techniques. 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.

# 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.")

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

There were 13  heads.

Choice of generator

As of version 1.17.3 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 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.

rg = 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 rg replaces np.random.

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

array([0.08497523, 0.68510652, 0.99816581, 0.0161308 , 0.93697552,
       0.01346615, 0.78713019, 0.60621575, 0.70869909, 0.79716342,
       0.72021754, 0.21321055, 0.02691178, 0.90362922, 0.175835  ,
       0.5290432 , 0.24906699, 0.54187212, 0.37458403, 0.31563962])

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
rg = np.random.default_rng(seed=3252)

# Draw random numbers
rg.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
rg = np.random.default_rng(seed=3252)

# Draw random numbers
rg.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.

rg = np.random.default_rng(seed=3253)
rg.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.

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 = rg.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. 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 from 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 rg.binomial().

# Draw how many coin flips land heads in 10 files
rg.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. The rg.choice() function does this. You equivalently could do this using rg.integers(), where the integers represent indices in the array, except rg.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 rg.integers(), we might have a problem.

rg = np.random.default_rng(seed=126969234)
rg.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 rg.choice() instead.

rg.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 rg.choice() function to generate random DNA sequences.

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

'CTTGTAAAGGATCTGCCCCGTGTCCCGCGCTAGGTCCGCACTCAAAGTATCTATGAACCGTCCAACCGAG'

Shuffling an array

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

rg.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 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,pandas,bokeh,iqplot,jupyterlab

CPython 3.8.5
IPython 7.18.1

numpy 1.19.1
scipy 1.5.0
pandas 1.1.3
bokeh 2.2.3
iqplot 0.1.6
jupyterlab 2.2.6