Definition 18.1. Sample Space.
The sample space of an experiment is the set of all possible outcomes.
Chapter 18 Probability Distributions and Pseudorandom Numbers
One big part of scientific computation is the subject of Monte Carlo simulations in which random numbers are used to model some situation. We’ll spend the next part of this course covering this subject.
I’m sure you have some sense of what probability is, that is a way to measure the randomness of something or how sure you are that something will happen. This chapter starts with basic definitions of probability and randomness and then covers both discrete and continuous probability distributions and then continue with a firm mathematical description. We will take a pause from coding to build up a bit of probability background.
Section 18.1 Basic Defintions
An experiment is something that produces outcomes. The terminology arises from scientific experiments and may produce images of chemistry labs, however, it is more general than that. An experiment may be a chemistry lab with measurements, but could be a drug-effectiveness study with tens of thousands of individuals, and also may include flipping a coin.
In the example drug study, it may be the possible measurable levels in a blood sample of whatever is being analyzed. In the coin flip it is the set \(\{H,T\}\text{.}\)
Example 18.2. Sample Space.
The sample space of flipping 3 coins is
\begin{equation*}
\{HHH,HHT,HTH,HTT,THH,THT,TTH, TTT\}
\end{equation*}
Definition 18.3. Event.
An event is a subset of the sample space.
Example 18.4. Event.
What is the event that you flip 2 heads out of 3 coins.
Solution.
If we examine the set that has 2 heads out of the set in Example 18.2, then it is the set
\begin{equation*}
A = \{HHT,HTH, THH\}
\end{equation*}
It will be important to understand events as sets and will help finding other information.
If the number of elements in a sample space is finite and denoted \(n(S)\text{,}\) then we can find the probability that an event \(A\) occurs with
\begin{equation*}
P(A) = \frac{n(A)}{n(S)}
\end{equation*}
Section 18.2 Random Variables
We often will need to know what the probability that an event occurs. Because events correspond to elements of the space, we often examine how every element in the sample space behaves in some way.
Definition 18.5. .
A random variable is a function from the sample space \(S\) to some real number.
Example 18.6. .
Let’s consider a six-sided die. The sample space is
\begin{equation*}
S=\{\text{\drawdie{1}},\text{\drawdie{2}},\text{\drawdie{3}},\text{\drawdie{4}},\text{\drawdie{5}},\text{\drawdie{6}}\}.
\end{equation*}
Let \(X\) be the number of dots on each die. Since \(X\) is a function from the sample space to the integers (subset of the reals), \(X\) is a random variable.
Example 18.7.
Let the sample space be that of the 3 coin flips from Example 18.2. Let \(X\) be the number of \(H\) in each coin flip.
Example 18.8.
Consider an experiment in which the amount of rain (in either inches or cm) is measured in your hometown on a given day. The sample space is \(\{ x \geq 0 \; | \; x \in \mathbb{R}\}\) or nonnegative real numbers. A random variable for this could simply be the amount of rain in the day.
Section 18.3 Discrete Random Variables
A set is discrete if it either contains finite number of elements or if the number of elements are infinite, they can be ordered (technically mapped to the natural numbers). This is important because if the sample space is discrete or not will determine how we treat the probabilities. A random variable is discrete if the sample space is discrete. The examples above with the coin flip (either 1 or 3 coins) and the dice are examples of discrete random variables. The random variable in Example \ref{ex:rain} is not.
Definition 18.9.
If a random variable is discrete, that is, the underlying sample space is discrete, the probability distribution or probability density function or pdf is a function \(f(x)\) defined by \(p(x) = P(X=x)\) for all \(x \in S\text{,}\) the sample space.
More colloquially, the probability density function gives the probability of every outcome in the sample space.
Example 18.10.
Consider the six-sided die. If the die is fair then the probability of each die coming up are the same.
\begin{equation*}
p(x) = 1/6 \qquad \text{for all $x \in \{1,2,3,4,5,6\}$.}
\end{equation*}
We can plot the distribution \(p(x)\) with
using Plots bar(1:6, [1/6 for i=1:6], legend=false)
which results in the plot
Example 18.12.
Find the probability distribution of a 3-coin flip experiment, where \(X\) is the number of heads. Assume the coin is fair.
Subsection 18.3.1 Properties of a discrete probability distribution
A function let \(x_i\) be the domain of a discrete random variable. The function \(p(x)\) is a probability density function if
-
\(\displaystyle \displaystyle p(x_1)+p(x_2) + \cdots + p(x_n) = 1\)
Example 18.14.
Consider rolling two fair dice and let \(X\text{,}\) the random variable be the sum of the number of dots on the dice. Find the probability density function.
Solution.
Think about the two dice have different colors, a red die and a green die. Since each die have 6 outcomes, there are \(6^2=36\) different outcomes and here are all the possible rolls (where the first number is the roll on the red die and the second number is the number on the green die).
\begin{align*}
(1,1) \amp\amp (1,2) \amp\amp (1,3) \amp\amp (1,4) \amp\amp (1,5) \amp\amp (1,6) \\
(2,1) \amp\amp (2,2) \amp\amp (2,3) \amp\amp (2,4) \amp\amp (2,5) \amp\amp (2,6) \\
(3,1) \amp\amp (3,2) \amp\amp (3,3) \amp\amp (3,4) \amp\amp (3,5) \amp\amp (3,6) \\
(4,1) \amp\amp (4,2) \amp\amp (4,3) \amp\amp (4,4) \amp\amp (4,5) \amp\amp (4,6) \\
(5,1) \amp\amp (5,2) \amp\amp (5,3) \amp\amp (5,4) \amp\amp (5,5) \amp\amp (5,6) \\
(6,1) \amp\amp (6,2) \amp\amp (6,3) \amp\amp (6,4) \amp\amp (6,5) \amp\amp (6,6)
\end{align*}
Now to find \(p(x)\text{,}\) we need to determine the probability of each of sum of the dice. The values of \(x\) can be any integer 2 through 12 and then the probability is the number of the sum over 36. For example, \(p(2)=1/36\) is because the only way to get a 2 is with the results (1,1). The value of \(p(3)=2/36\) because one can get 3 with a (1,2) or a (2,1). This logic continues for all values from 2 to 12. The values of \(p(x)\) are
\begin{align*}
p(2) \amp = 1/36 \amp p(3) \amp = 2/36 \amp p(4) \amp = 3/36 \amp p(5) \amp = 4/36 \amp p(6) \amp = 5/36 \amp p(7) \amp = 6/36\\
p(8) \amp = 5/36 \amp p(9) \amp = 4/36 \amp p(10) \amp = 3/36 \amp p(11) \amp = 2/36 \amp p(12) \amp = 1/36
\end{align*}
A plot of this distribution is
bar(2:12,[(6-abs(i-7))/36 for i=2:12], xticks=2:12, legend=false)
Section 18.4 Continuous Probability
If a random variable is not discrete it is called a continuous random variable and generally occurs if the values it takes on are subsets of the real line. Let’s consider a set \(A=\{x \; | \; 0 \leq x \leq 1\}\) or all real numbers between 0 and 1. Events are still subsets of the set \(A\text{,}\) however the probability that events occur is the fraction of the set.
Subsection 18.4.1 The Calculus of Continuous Probabilities
The probability distribution function or pdf is a function, \(f\) defined on a set \(S\) with the following properties:
\begin{align}
f(x) \amp \geq 0 \qquad \text{for all $x \in S$}\tag{18.1}\\
\int_S f(x) \,dx \amp = 1 \tag{18.2}
\end{align}
where the set \(S\) is generally a real interval or the infinite or semi-infinite line. The first property ensures that probability will always be positive and the second property ensures that probabilities are at most 1 and only 1 with the entire set. Also, if \(f\) is a continuous pdf, then the probability that a random variable \(X\) takes on a value in the set \(A\) is
\begin{equation*}
P(A) = \int_A f(x) \, dx
\end{equation*}
Subsection 18.4.2 Cumulative Distribution Functions
Another helpful function of random variables is that of the cumulative distribution function or cdf for a given pdf, \(f(x)\text{.}\) Denote \(F(x)\text{,}\) the cdf, which is defined as
\begin{equation*}
F(x) = \int_{-\infty}^x f(t) \, dt
\end{equation*}
or simply the antiderivative. Graphically, if we have a function \(f\) as shown below, then \(F(x)\) represents the area under the curve to the right of the value \(x\text{.}\)
Also, since \(F\) is the antiderivative, then the cdf can be used to find the probability.
Remark 18.18.
Let \(F(x)\) be a cdf where \(f(x)\) is the given pdf of the random variable \(X\text{.}\) The probability can be found using
\begin{equation*}
P(a \leq X \leq b) = F(b) - F(a)
\end{equation*}
Subsection 18.4.3 Inverse Cumulative Distribution Functions
There are many nice properties of \(F(x)\text{,}\) the cdf. We know:
-
\(\displaystyle \lim_{x \rightarrow -\infty} F(x) = 0 \)
-
\(\lim_{x \rightarrow \infty} F(x) = 1\text{.}\)
-
\(F(x)\) is continuous.
-
\(F(x)\) is nondecreasing.
where some of these properties are explored in the exercises. The last 2 properties above result give that the function \(F(x)\) is one-to-one and therefore has an inverse. The function \(F^{-1}(x)\) is called the inverse cumulative distribution function. The section above, showed the importance of evaluating a cdf for probabilities. That is given
\begin{equation*}
F(x) = \int_{-\infty}^x f(t) \, dt
\end{equation*}
evaluating \(F(a)\) for a given value of \(a\text{.}\) We will see that it is just as important to solve \(F(x)=p\) for \(x\text{.}\) This corresponds to the same plot as above:
where given a value of \(p\text{,}\) we are seeking \(x\text{.}\) Since \(F(x)=p\text{,}\) we also write \(x=F^{-1}(p)\text{.}\) The function \(F^{-1}(x)\) have the properties:
and in short, this function finds \(x\) values for given probabilities values.
Section 18.5 Standard Continuous Distributions
Subsection 18.5.1 Uniform Distribution
The example above at the top of the previous section was an example of a uniform distribution with pdf:
\begin{equation*}
u(x) = \begin{cases}
1 \amp 0 \leq x \leq 1 \\
0 \amp \text{otherwise}
\end{cases}
\end{equation*}
It is clear that the property in (\ref{eq:nonneg}) is satisfied and note that
\begin{equation*}
\int_{\mathbb{R}} u(x) = \int_0^1 1 \, dx = 1
\end{equation*}
\begin{equation*}
X=\left\{x \; \biggr| \; \frac{1}{3} \leq x \leq \frac{1}{2}\right\},
\end{equation*}
then
\begin{equation*}
P(X) = \int_{1/3}^{1/2} 1 \,dx = \frac{1}{2} - \frac{1}{3} = \frac{1}{6}\text{.}
\end{equation*}
Also, a plot of the distribution is
Check Your Understanding 18.21.
Find the pdf of a uniform distribution that is a positive constant \(c\) on \(a \leq x \leq b\) and 0 elsewhere. Note, you must find the value of \(c\) so the pdf satisfies (18.2).
Subsection 18.5.2 The Normal Distribution
Probably the most important continuous distribution is the normal distribution (also known as the bell curve) which is symmetric about some center value and tapers off to zero in both directions. An example is:
1
No pun intended.
The shape of the normal distribution has two important parameters, the mean, \(\mu\text{,}\) which is the line of symmetry and the standard deviation, \(\sigma\) which determines how spread out it is. The functional form of the pdf is
\begin{equation*}
f(x \, | \, \mu, \sigma) = \frac{1}{\sqrt{2\pi \sigma^2}} e^{-(x-\mu)^2/(2\sigma^2)}
\end{equation*}
Section 18.6 Pseudo Random Number Generator
In the next few chapters of the text, we will solving problems with random numbers. Thus far, we talked about flipping coins and rolling dice and we will also see about drawing cards from a shuffled deck. Since each of these is a physical process, we can’t directly do them, however, we will simulate them with random numbers, and technically on a computer we will use a pseudorandom number generator.
A pseudo-random number generator is a function or procedure that produces a sequence of numbers that behave like random numbers. Let’s examine what this means in terms of the discrete probability with discrete probability space \(S=\{1,2,3,4,5,6\}\) which would generate a simulation of rolling a single fair die. If a pseudo-random number generator produces a sequence from this set then the sequence should have the following property:
-
If the event \(X=\{i\}\) for any \(i\) between 1 and 6, then \(P(X)\approx 1/6\text{.}\) And by approximately, as the sequence gets larger, the approximation becomes closer to 1/6.Is this enough? No, the sequence\begin{equation*} \{1,2,3,4,5,6,1,2,3,4,5,6,\ldots\} \end{equation*}satisfies the above property, but I don’t think anyone would consider this random. Another property would be:
-
If we know the sequence \(\{a_1,a_2,a_3,\ldots,a_n\}\) then we can’t predict the next number \(a_{n+1}\text{.}\)
This is obviously violated in the sequence above. A little more technical definition of a sequence of pseudo-random numbers Let\\ \((a_1,a_2,a_3, \ldots)\) be a random sequence. Typically we mean the following properties need to hold:
-
any number in the range 1 to 6 is equally likely to occur.
-
Take \(N\) random numbers and let \(s_n\) be the number of times the number \(n\) occurs. The fraction \(s_n/N\) should go to 1/6 in the limit as \(n\rightarrow \infty\text{.}\)
Subsection 18.6.1 Seeds and Pseudorandom Numbers
For pseudorandom numbers, the third condition above is not true. In fact, there is a function that generates the next number in the sequence or there is a function \(F\) such that \(s_{n+1}=F(s_n)\text{,}\) however, it is not evident that there is such a function. Because pseudorandom numbers are generated with a function, it needs a number to determine the next one. This starting number is often chosen using the current date/time, however, you can set this, called the seed. In julia, we can set the seed with the
Random.seed! function.
using Random Random.seed!(1234)
and any positive integer can be set here. Note: you don’t need to set the seed for random numbers. It is useful in the case of writing a book in which you hope that the readers will follow along. In this case, everything should be randomly generated in an identical way.
Subsection 18.6.2 Using Julia to simulate the rolling of a die
First, the main commands that are built-in to Julia are listed in the Julia Manual for Random Numbers. We can generate 10 random numbers between 1 and 6 using
S=rand(1:6,10)
which returns the array
[2, 4, 2, 6, 3, 3, 6, 5, 3, 5] .
Notice that if you rerun the command, you’ll get a different sequence of random numbers. We can check that this is doing what we expect by checking the probability that we get a 1 (or any other number).
Subsection 18.6.3 Coin flips in Julia
We can flip coins by selecting a random boolean. For example, the following will simulate 10 coin flips.
coins=rand(Bool,10)
which returns the array
[0, 0, 0, 1, 1, 0, 1, 0, 1, 0] This is a boolean Vector, note the Bool at the beginning. The 1s represent true and 0, false for compactness. Generally, as we did above, we will flip multiple coins and ask questions about the number of heads. We will examine this in the next chapter. If we are interested in unfair coin flips, we can choose randomly from an array where the unfair balance is known. For example,
rand([false,false,false,true,true],10)
gives 10 coin flips or the array
[0 ,0 ,0 ,1 ,1 ,1 ,1 ,1 ,0 ,0] as a weighted coin where \(P(H)=2/5\) and \(P(T)=3/5\text{.}\)
Subsection 18.6.4 Floating Point random numbers
For continuous random variable, like the uniform distribution or the Normal distribution, there are some different properties of the pseudo random number generate. Assume that the floating point number is in the range \(0 < x < 1\text{.}\) Then the sequence \((a_1,a_2,a_3, \ldots)\) is random if
To generate an array of 20 random numbers distributed uniformly on \([0,1)\text{,}\) type
S=rand(20)
which returns
20-element Vector{Float64}:
0.3020595409518265
0.05842018425830109
0.20552381705814504
0.7079246996200363
0.41775698944882156
0.22721762712581473
0.7836312465908486
0.20071380866783706
0.6020380841683062
0.5653680870733934
0.7008893474454004
0.4527965951678993
0.13067526601187207
0.4578870482605456
0.6665192594508963
0.5008741327606322
0.8106990681415074
0.46484850303763725
0.9770123967639948
0.5001612179097075
Section 18.7 The Normal Distribution
The normal distribution is the most important continuous distributions. We can investigate it with the
Distributions package in julia:
using Distributions
To produce a standard normal (mean of 0, standard deviation of 1),
std_normal = Normal()
which returns
Normal{Float64}(μ=0.0, σ=1.0) and notice the mean and standard deviation are what we expect. Other means and standard deviations can be found by adding arguments. We can plot the normal distribution with
using StatsPlots plot(std_normal)
To generate an array of normally distributed values,
random_normals = rand(std_normal,10)
returns
10-element Vector{Float64}:
2.121560382124172
-0.26224077844309135
0.3455010838287982
0.44394906831563774
0.44133351371484625
-1.3397752732187238
0.21507132788419303
0.28927403687088354
-0.6365687810422291
-0.8196796137020412
Section 18.8 Other Distributions
Julia has many other distributions available in the
Distributions package and the best place to start is at the Distributions package website
