import numpy as np8 Probability theory
8.1 Randomness
There are various ways to generate random variables and sets in Python. Throughout this chapter we will rely on basic concepts from probability theory such as the probability density function (pdf) and cumulative distribution function (cdf).
Furthermore, we will be working with the normal, uniform and other distributions. We assume you are familiar with these distributions. Whenever you run into a mathematical or probabilistic aspect you do not understand, please look it up and otherwise ask the teacher about it.
We will use functionality from the random subpackage of NumPy. This subpackage is not to be confused with Python’s standard-library random module, which has a lot of similar functionality.
8.1.1 Samples
The random subpackage has functions that can generate so-called pseudo-random numbers. Informally speaking, these numbers behave like random numbers when performing numerical experiments with them, but they are not truly random (which can be an issue, e.g., if you used these numbers for cryptographic purposes).
We next discuss some functions to generate random numbers from well-known distributions.
np.random.rand(m,n): Generate an m \times n two-dimensional array with numbers uniformly drawn from the interval [0,1). That is, subintervals of equal length have equal probability, which intuitively means that every number is equally likely to be sampled.
m,n = 2, 4
M = np.random.rand(m,n)
print(M)[[0.75848035 0.33499741 0.75427562 0.8172453 ]
[0.91739344 0.58462821 0.68086611 0.80690302]]When no dimensions m and n are supplied, a single random number is returned.
x = np.random.rand()
print(x)0.5792754564208292np.random.randn(m,n): Generates an m \times n two-dimensional array with numbers drawn from the normal distribution with mean \mu = 0 and standard deviation \sigma = 1.
m,n = 3, 5
M = np.random.randn(m,n)
print(M)[[-0.243273 0.41639416 -0.65640647 0.1024049 0.37905796]
[ 0.65492487 1.97533233 -1.29199186 1.22705906 0.67947861]
[ 0.16859646 1.10668774 0.21278151 0.18181822 -0.22282751]]We conclude with sampling numbers from the discrete uniform distribution over a set \{a,a+1,\dots,b-1\}. Here every number i in this discrete interval is generated with equal probability 1/(b - a).
np.random.randint(a,b,(m,n)): Generates an m \times n two-dimensional array with numbers drawn from the discrete uniform distribution on \{a,a+1,\dots,b-1\}.
a, b = 4, 10
m,n = 3, 15
M = np.random.randint(a,b,(m,n))
print(M)[[7 8 5 5 9 4 5 6 4 8 5 7 4 8 6]
[4 9 9 5 4 4 8 7 4 5 7 9 6 4 6]
[9 5 7 4 4 8 9 5 8 9 6 9 7 4 5]]NumPy provides sampling methods for many common probability distributions. See the documentation for a list of all the distributions.
Finally, we remark that it is also possible to generate samples from different distributions from the same family. As an example, consider a normal distribution with mean \mu and standard deviation \sigma. We can generate samples from this distribution with np.random.normal().
As input it takes two keyword arguments: loc for the mean \mu and scale for the standard deviation \sigma. You can figure this out by inspecting the documentation.
Most distributions have optional loc and scale arguments, which have default values if none are specified; loc=0 and scale=1 are the defaults for np.random.normal().
mu = 1
sigma = 2
samples = np.random.normal(loc=mu, scale=sigma)
print(samples)3.803412802646984This function also allows vectorized inputs: If we input an array of means and standard deviations, then a sample for all corresponding mean–standard-deviation pairs.
mu = np.array([0,1,2,2])
sigma = np.array([1,2,4,2])
samples = np.random.normal(loc=mu, scale=sigma)
print(samples)[1.70428604 1.5102823 2.39072117 0.19740327]Finally, if you want to have m samples from all corresponding mean–standard-deviation pairs of location and scale parameters, you can specify this in the size keyword argument by setting it to (m,n) where n is the common length of the location and scale array. For every scale-location pair, the m samples of this pair can be found in a column of the output.
mu = np.array([0,1,2,2])
sigma = np.array([1,2,4,2])
m = 3
n = np.size(mu)
samples = np.random.normal(loc=mu, scale=sigma, size=(m,n))
print(samples)[[-0.98797274 -2.20406818 8.75405895 3.17724463]
[ 1.68112359 2.86267349 2.92929963 0.69686498]
[-0.37543368 0.92896628 -1.97601306 0.67003226]]8.1.2 Subsets
In addition to generating random numbers, it is also possible to generate a random subset (with duplicates allowed) of elements of a given array using the choice() function. It takes as input an array from which we want to obtain a subset and we can set the size of the subset that we want to have using the size keyword argument.
k = 15
x = np.arange(0,11,1)
subset = np.random.choice(x,size=k)
print(subset)[ 6 1 2 10 7 9 7 4 5 4 8 0 0 2 1]As can be seen from the output above, some numbers appear twice in the subset, meaning that choice() samples a subset with replacement. If you want to sample without replacement, you can set the keyword replace to False.
k = 8
x = np.arange(0,11,1)
subset = np.random.choice(x,size=k,replace=False)
print(subset)[ 3 7 2 9 0 5 4 10]Finally, you can also specify the probability with which each element is sampled using the keyword argument p.
k = 8
# Set {0,1,2,...,9}
x = np.arange(0,10,1)
# Probabilities for elements in set
prob = np.array([1/3,1/3,0,0,0,0,0,0,0,1/3])
# Generating random subset
subset = np.random.choice(x,size=k,p=prob)
# Only 0, 1, 9 have positive probabilities
print(subset)[0 1 0 1 1 9 9 0]There is also the permutation() function that returns the elements in an array in a random order. That is, it creates a so-called random permutation of the elements in the array.
x_perm = np.random.permutation(x)
print(x_perm)[7 6 9 8 0 1 5 3 2 4]If you apply this function to a two-dimensional array, it returns the same array in which the rows are randomly permuted (i.e., the inner lists are randomly permuted).
X = np.arange(0,18,1).reshape(3,6)
x_perm = np.random.permutation(X)
print(x_perm)[[ 0 1 2 3 4 5]
[12 13 14 15 16 17]
[ 6 7 8 9 10 11]]8.1.3 Seed
When writing code that involves random numbers or objects, it can sometimes be useful to “fix” the randomness in the script, e.g., when debugging. This can be done by setting a so-called random seed using the seed() function.
For the sake of comparison, let us first generate two random numbers from [0,1].
a = np.random.rand()
b = np.random.rand()
print(a,b)0.8465842091410098 0.12091672686991706If you rerun the code above it will give different outputs every time. Try this yourself by copying the code into, e.g., Spyder.
We next do the same, but with a fixed seed s using seed(s). Different choices of s fix the randomness of a Python script in a different way. If you store the code below into a Python script and rerun it a couple of times, the output will always be the same. We emphasize that the command np.random.rand() gives different outputs for a and b within the Python script.
#Set seed to be s = 3
np.random.seed(3)
a = np.random.rand()
b = np.random.rand()
print(a,b)0.5507979025745755 0.70814782261810488.2 Probability distributions
The stats module of SciPy has many built-in probability distributions. Each distribution can be seen as an object on which various methods can be performed (such as accessing its probability density function or summary statistics like the mean and median). You should think of an ‘object’ in the context of object-oriented programming, see, e.g., here to recall the basics of this programming paradigm.
import scipyIn this section we will focus on continuous probability distributions. SciPy also has many built-in discrete probability distributions.
A list of all continuous distributions that are present in the stats module can be found here; they are so-called stats.rv_continuous objects. We can instantiate a distributional object by using scipy.stats.dist_name where dist_name is the name of a built-in (continuous) probability distribution in the mentioned list, together with values for its required parameters.
Many distributions have required parameters scale and loc that model the scale and location of the distribution, respectively. Depending on the distribution that is considered, these parameters have different meanings.
As an example, the normal distribution has probability density function f(x) = \frac{1}{\sqrt{2\pi\sigma^2} } e^{-\frac{(x-\mu)^2}{2\sigma^2}} which is parameterized by \mu and \sigma > 0.
In Python \mu is the loc parameter, and \sigma the scale parameter. To figure out the function of the scale and loc parameter, you can check the documentation, which can be found here for the normal distribution.
All distributions have default values for the loc and/or scale parameters, which are loc=0 and scale = 1 for the normal distribution.
# Instantiate normal distribution object with mu=0, sigma=1
dist_norm = scipy.stats.norm(loc=0, scale=1)Once a distribution object has been instantiated, we can use methods (i.e., functions) to obtain various properties of the distribution, such as its probability density function (pdf), cumulative distribution function (cdf) and summary statistics such as the mean, variance and median (or, more generally, quantiles).
We give a list of some common methods for a distribution object named dist_name. We start with common functions associated with a probability distribution.
dist_name.pdf(x): Value f(x) where f is the pdf of the distribution.dist_name.cdf(x): Value F(x) where F is the cdf of the distribution.dist_name.sf(x): Value S(x) where S is the survival function (1 - F) of the distribution.dist_name.ppf(alpha): Returns x so that F(x) = \alpha where \alpha \in (0,1).
x = 1
print(dist_norm.pdf(x))0.24197072451914337alpha = 1/2
print(dist_norm.ppf(alpha))0.0All the above functions are vectorized, in the sense that they can also handle higher-dimensional arrays as input. This is convenient, e.g., for visualizing these functions as the example below illustrates.
Show code generating the plot below
import matplotlib.pyplot as plt
# Define the x range of x-values
x = np.linspace(-10,10,600)
# Function values of pdf
y = dist_norm.pdf(x) # Computes pdf values of all elements in x
#Create the figure
plt.figure()
# Create the plot
plt.plot(x, y)
# Create labels for axes
plt.xlabel('x')
plt.ylabel('Function value')
# Fix the range of the axes
plt.xlim(-10,10)
plt.ylim(0,1)
# Add title to the plot
plt.title('PDF of the normal distribution with $\mu=0$ and $\sigma=1$')
# Add grid to the background
plt.grid()
# Show the plot
plt.show()
We can also access various summary statistics:
dist_name.mean(): Returns mean of the distributiondist_name.var(): Returns variance of the distributiondist_name.median(): Returns median of the distribution
dist_norm = scipy.stats.norm(loc=0,scale=2)
mean = dist_norm.mean()
variance = dist_norm.var()
median = dist_norm.median()
print("Mean of the distribution is", mean)
print("Variance of the distribution is", variance)
print("Median of the distribution is", median)Mean of the distribution is 0.0
Variance of the distribution is 4.0
Median of the distribution is 0.0Note that in the example above the standard deviation equals \sigma = 2; the variance is then \sigma^2 = 4.
Finally, you can obtain the lower and upper bounds of the support.
dist_name.support(): Returns lower bound a and upper bound b of (smallest) interval [a,b] for which all probability mass is contained in it.
Let us consider the uniform distribution.
Here the scale and loc parameters result in a uniform distribution on the interval [loc, loc + scale].
#Uniform distribution on [3,3+4] = [3,7]
dist_unif = scipy.stats.uniform(3,4)
a, b = dist_unif.support()
print(f"The distribution is supported on the interval [{a},{b}]")The distribution is supported on the interval [3.0,7.0]Some distributions have an unbounded support. In the case of the normal distribution, we have (a,b) = (-\infty,\infty). For the exponential Distribution with loc=0, we have [a,b) = [0,\infty).
dist_norm = scipy.stats.norm(loc=0, scale=1)
a, b = dist_norm.support()
print(f"The normal distribution is \
supported on the interval [{a},{b}]")The normal distribution is supported on the interval [-inf,inf]dist_exp = scipy.stats.expon(loc=0,scale=1)
a, b = dist_exp.support()
print(f"The exponential distribution is \
supported on the interval [{a},{b}]")The exponential distribution is supported on the interval [0.0,inf]The value inf that b has in the example above is in fact a number within NumPy, namely np.inf that is a floating-point number representing positive infinity. Let us check that the upper bound of the exponential distribution is indeed np.inf.
statement = (b == np.inf)
print(statement)True8.2.1 Distributions as input arguments
Distributional objects can also serve as input arguments of a function. In that case, you can access the methods of the object inside the function.
Suppose we want to write a function that outputs a message saying whether the mean is smaller than, equal to, or larger than the median. We can do this as follows for an arbitrary distribution that is assumed to have finite mean and median.
def mean_median(dist):
if dist.mean() < dist.median():
return "The mean is smaller than the median."
elif np.isclose(dist.mean(), dist.median()): # Checks dist.mean() == dist.median()
return "The mean and median are equal." # but allows for some numerical imprecision
elif dist.mean() > dist.median():
return "The mean is larger than the median."In the function above, the input argument dist is a distribution object whose methods mean() and median() we access within the function.
dist_norm = scipy.stats.norm(loc=1,scale=2)
comparison = mean_median(dist_norm)
print(comparison)The mean and median are equal.