Lesson 26: Random number generation¶
(c) 2016 Justin Bois. This work is licensed under a Creative Commons Attribution License CC-BY 4.0. All code contained herein is licensed under an MIT license.
This tutorial was generated from a Jupyter notebook. You can download the notebook here.
import random
import numpy as np
# This is how we import the module of Matplotlib we'll be using
import matplotlib.pyplot as plt
# Some pretty Seaborn settings
import seaborn as sns
rc={'lines.linewidth': 2, 'axes.labelsize': 18, 'axes.titlesize': 18}
sns.set(rc=rc)
# The following is specific Jupyter notebooks
%matplotlib inline
%config InlineBackend.figure_formats = {'png', 'retina'}
Random number generation (RNG), besides being a song in the original off-Broadway run of Hedwig and the Angry Inch, is the process by which a string of random numbers may be drawn. Of course, the numbers are not completely random for several reasons.
- 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.
- 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. Much of the random number generation functionality you will need is in the np.random module. Let's start by generating random numbers on a uniform interval
np.random.random(size=10)
The function random() in the np.random module generates random numbers on the interval $[0,1)$. The size kwarg is how many random numbers you wish to generate. They 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 functions, just like we did last time. We'll generate 100,000 random numbers and plot them along with the CDF.
# Generate sorted random numbers
x = np.sort(np.random.random(size=100000))
# Generate y-axis for CDF
y = np.arange(1, len(x)+1) / len(x)
# Plot CDF from random numbers (for plotting purposes, only plot 100 points)
plt.plot(x[::1000], y[::1000], marker='.', linestyle='none', markersize=10)
# Plot expected CDF (just a straight line from (0,0) to (1,1)
plt.plot([0, 1], [0, 1], 'k-')
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 do "coin flips," that is random draws that are either zero or one, like this:
# Generate 20 random numbers on uniform interval
x = np.random.random(size=20)
# Make them coin flips
heads = x > 0.5
# Show which were heads, and count the number of heads
print(heads)
print('\nThere were', np.sum(heads), ' heads.')
Seeding random number generators¶
Now, just to demonstrate that random number generation is deterministic, we will explicitly seed the random number generator (usually seeded with a number representing the date/time to avoid repeats) to show that we get the same random numbers.
# Seed the RNG
np.random.seed(42)
# Generate random numbers
np.random.random(size=10)
# Re-seed the RNG
np.random.seed(42)
# Generate random numbers
np.random.random(size=10)
The random numbers are exactly the same. If we choose a different seed, we get totally different random numbers.
# Seed with a number that is close to the answer to everything
np.random.seed(43)
np.random.random(size=10)
If you are writing unit tests (which we will describe tomorrow when we do test driven development), 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. For example, say we wanted to draw random samples from a Normal distribution with mean $\mu$ and standard deviation $\sigma$.
# Set parameters
mu = 10
sigma = 1
# Draw 10000 random samples
x = np.random.normal(mu, sigma, size=10000)
# Plot a histogram of our draws
_ = plt.hist(x, bins=100)
It looks Normal, but let's test the mean and standard deviation.
np.mean(x), np.std(x)
Yup, mean of 10 and standard deviation of 1!
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 that we can do this for "coin flips," but we can do it for other distributions as well. Perhaps the most common example is a choice of a random integer from a set of integers. For example, if we wanted to draw one of four integers in {0, 1, 2, 3}, we can use the np.random.randint() function.
# Draw random integers on [0, 4), i.e., exclusive of last one.
np.random.randint(0, 4, 20)
I bet you can guess where this might be useful! We can generate random DNA sequences.
# Key of bases
bases = 'ATGC'
# Draw random numbers for sequence
x = np.random.randint(0, 4, 50)
# Make sequence
seq_list = [None]*50
for i, b in enumerate(x):
seq_list[i] = bases[b]
# Join the sequence
''.join(seq_list)
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 np.random.choice() function does this. You equivalently could do this using np.random.randint(), where the integers represent indices in the array, except np.random.choice() has a great keyword argument, replace, which allows random draws with or without replacement. For example, say you had 100 samples that you wanted to send to a facility for analysis, but you can only afford to send 20. If we used np.random.randint(), we might have a problem.
np.random.randint(0, 51, 20)
Sample 18 was selected thrice and samples 11 and 23 were selected twice! We can use np.random.choice() instead.
np.random.choice(np.arange(51), 20, replace=False)
Shuffling an array¶
Similarly, the np.random.permutation() function is useful. It takes the entries in an array and shuffles them! Let's shuffle a deck of cards.
np.random.permutation(np.arange(53))
The random module¶
Though it has far less functionality than np.random, the random module in the standard library has some useful functionality. In particular, the random.choice() method can take a string as an input and choose characters out of that string. For example, to generate random DNA sequences, we could just do the following.
# Make sequence
seq_list = [None]*50
for i in range(len(seq_list)):
seq_list[i] = random.choice('ATGC')
# Join the sequence
''.join(seq_list)
When do we need RNG?¶
Answer: VERY OFTEN! We will see many examples in the next lessons and in the exercises.
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.
RNG also comes up A LOT in data analysis. We'll close this lesson with a brief example.