Lesson 28: Random number generation¶
(c) 2019 Justin Bois. With the exception of pasted graphics, where the source is noted, 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 document was prepared at Caltech with financial support from the Donna and Benjamin M. Rosen Bioengineering Center.
This tutorial was generated from a Jupyter notebook. You can download the notebook here.
import numpy as np
import pandas as pd
import scipy.stats
import bokeh_catplot
import bokeh.io
import bokeh.plotting
bokeh.io.output_notebook()
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. ("Completely random" does not make sense because of the infinite magnitude of numbers.)
- 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 from a Uniform distribution.
np.random.uniform(low=0, high=1, size=10)
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 = bokeh_catplot.ecdf(
data=pd.DataFrame({'x': x}),
cats=None,
val='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.')
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.
# Seed the RNG
np.random.seed(42)
# Generate random numbers
np.random.uniform(size=10)
# Re-seed the RNG
np.random.seed(42)
# Generate random numbers
np.random.uniform(size=10)
The random number sequence is 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.uniform(size=10)
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. For example, say we wanted to draw random samples from a Normal distribution with mean μ and standard deviation σ. (We already saw this example when we were looking at histograms, but we repeat it here.)
# Set parameters
mu = 10
sigma = 1
# Draw 100000 random samples
x = np.random.normal(mu, sigma, size=100000)
# Plot the histogram
p = bokeh_catplot.histogram(
data=pd.DataFrame(data={'x': x}),
cats=None,
val='x',
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.
# 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 = bokeh_catplot.ecdf(
data=pd.DataFrame({'x': x[:1000]}),
cats=None,
val='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 np.random.binomial().
# Draw how many coin flips land heads in 10 files
np.random.binomial(10, 0.5)
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.seed(42)
np.random.randint(0, 51, size=20)
Sample 10 was selected twice and sample 23 was selected thrice! We can use np.random.choice() instead.
np.random.choice(np.arange(51), size=20, replace=False)
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 np.random.choice() function to generate random DNA sequences.
''.join(np.random.choice(list('ATGC'), replace=True, size=70))
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))
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, which we will see in the lesson on hacker stats.
Computing environment¶
%load_ext watermark
%watermark -v -p numpy,scipy,pandas,bokeh,bokeh_catplot,jupyterlab