Probability CheatSheet

Overview

Teaching: min
Exercises: min

Questions

Objectives

Counting

Multiplication Rule

Let’s say we have a compound experiment (an experiment with multiple components). If the 1st component has \(n_1\) possible outcomes, the 2nd component has \(n_2\) possible outcomes, …, and the \(r\)th component has \(n_r\) possible outcomes, then overall there are \(n_1n_2 \dots n_r\) possibilities for the whole experiment.

Sampling Table

The sampling table gives the number of possible samples of size \(k\) out of a population of size \(n\), under various assumptions about how the sample is collected.

Naive Definition of Probability

If all outcomes are equally likely, the probability of an event \(A\) happening is:

[P_{\textrm{naive}}(A) = \frac{\text{number of outcomes favorable to $A$}}{\text{number of outcomes}}]

Independence

Independent Events:

\(A\) and \(B\) are independent if knowing whether \(A\) occurred gives no information about whether \(B\) occurred. More formally, \(A\) and \(B\) (which have nonzero probability) are independent if and only if one of the following equivalent statements holds:
 ֍ \(P(A \cap B) = P(A)P(B)\)
 ֍ \(P(A|B) = P(A)\)
 ֍ \(P(B|A) = P(B)\)

Conditional Independence:

\(A\) and \(B\) are conditionally independent given \(C\) if
\(P(A \cap B\|C) = P(A\|C)P(B\|C)\).
Conditional independence does not imply independence, and independence does not imply conditional independence.

Unions, Intersections, and Complements

 ֍ De Morgan’s Laws: A useful identity that can make calculating probabilities of unions easier by relating them to intersections, and vice versa. Analogous results hold with more than two sets.
 ֍ \((A \cup B)^c = A^c \cap B^c\)
 ֍ \((A \cap B)^c = A^c \cup B^c\)

Joint, Marginal, and Conditional

 ֍ Joint Probability: \(P(A \cap B)\) or \(P(A, B)\) – Probability of \(A\) and \(B\).
 ֍ Marginal (Unconditional) Probability: \(P(A)\) – Probability of \(A\).
 ֍ Conditional Probability: \(P(A\|B) = \frac{P(A, B)}{P(B)}\) – Probability of \(A\), given that \(B\) occurred.
 ֍ Conditional Probability is Probability: \(P(A\|B)\) is a probability function for any fixed \(B\). Any theorem that holds for probability also holds for conditional probability.

Probability of an Intersection or Union

Intersections via Conditioning

[\begin{align} P(A,B) &= P(A)P(B|A)
P(A,B,C) &= P(A)P(B|A)P(C|A,B) \end{align}]

Unions via Inclusion-Exclusion

[\begin{align} P(A \cup B) &= P(A) + P(B) - P(A \cap B)
P(A \cup B \cup C) &= P(A) + P(B) + P(C)
&\quad - P(A \cap B) - P(A \cap C) - P(B \cap C)
&\quad + P(A \cap B \cap C) \end{align}]

Simpson’s Paradox

It is possible to have \(P(A\mid B,C) < P(A\mid B^c, C) \text{ and } P(A\mid B, C^c) < P(A \mid B^c, C^c)\) \(\text{yet also } P(A\mid B) > P(A \mid B^c).\)

Law of Total Probability (LOTP)

Let \({ B}_1, { B}_2, { B}_3, ... { B}_n\) be a partition of the sample space (i.e., they are disjoint and their union is the entire sample space).

[\begin{align} P({ A}) &= P({ A} | { B}_1)P({ B}_1) + P({ A} | { B}_2)P({ B}_2) + \dots + P({ A} | { B}_n)P({ B}_n)
P({ A}) &= P({ A} \cap { B}_1)+ P({ A} \cap { B}_2)+ \dots + P({ A} \cap { B}_n) \end{align}]

For LOTP with extra conditioning, just add in another event \(C\)!

[\begin{align} P({ A}| { C}) &= P({ A} | { B}_1, { C})P({ B}_1 | { C}) + \dots + P({ A} | { B}_n, { C})P({ B}_n | { C})
P({ A}| { C}) &= P({ A} \cap { B}_1 | { C})+ P({ A} \cap { B}_2 | { C})+ \dots + P({ A} \cap { B}_n | { C}) \end{align}]

Special case of LOTP with \({ B}\) and \({ B^c}\) as partition:

[\begin{align} P({ A}) &= P({ A} | { B})P({ B}) + P({ A} | { B^c})P({ B^c})
P({ A}) &= P({ A} \cap { B})+ P({ A} \cap { B^c}) \end{align}]

Bayes’ Rule

Bayes’ Rule, with extra conditioning (just add in \(C\)!)

We can also write

Odds Form of Bayes’ Rule
\(\frac{P({ A}\| { B})}{P({ A^c}\| { B})} = \frac{P({ B}\|{ A})}{P({ B}| { A^c})}\frac{P({ A})}{P({ A^c})}\)

The posterior odds of \(A\) are the likelihood ratio times the prior odds.

Random Variables and their Distributions

PMF, CDF, and Independence

Probability Mass Function (PMF)
Gives the probability that a discrete random variable takes on the value x.

[p_X(x) = P(X=x)]

The PMF satisfies:
\(p_X(x) \geq 0 \quad \textrm{and} \quad \sum_x p_X(x) = 1\)

Cumulative Distribution Function (CDF)
Gives the probability that a random variable is less than or equal to x.

[F_X(x) = P(X \leq x)]

The CDF is an increasing, right-continuous function with:
\(F_X(x) \to 0 \quad \textrm{as} \quad x \to -\infty \quad \textrm{and} \quad F_X(x) \to 1 \quad \textrm{as} \quad x \to \infty\)

Independence
Intuitively, two random variables are independent if knowing the value of one gives no information about the other. Discrete random variables X and Y are independent if for all values of x and y:

[P(X=x, Y=y) = P(X = x)P(Y = y)]

Expected Value and Indicators

Expected Value and Linearity

Expected Value
(a.k.a. mean, expectation, or average) is a weighted average of the possible outcomes of our random variable.
Mathematically, if x₁, x₂, x₃, … are all of the distinct possible values that X can take, the expected value of X is:

\(E(X) = \sum\limits_{i}x_iP(X=x_i)\) Linearity
For any random variables X and Y, and constants a, b, c:

[E(aX + bY + c) = aE(X) + bE(Y) + c]

Same distribution implies same mean
If X and Y have the same distribution, then E(X) = E(Y) and, more generally:

[E(g(X)) = E(g(Y))]

Conditional Expected Value
Conditional expected value is defined like expectation, only conditioned on any event A:

Indicator Random Variables

Indicator Random Variable An indicator random variable is a random variable that takes on the value 1 or 0. It is always an indicator of some event: if the event occurs, the indicator is 1; otherwise, it is 0. They are useful for many problems about counting how many events of some kind occur. Write:

[I_A = \begin{cases} 1 & \text{if $A$ occurs},
0 & \text{if $A$ does not occur}. \end{cases}]

Note that: \(I_A^2 = I_A\),  \(I_A I_B = I_{A \cap B}\),  \(I_{A \cup B} = I_A + I_B - I_A I_B\)

Distribution
\(I_A \sim Bern(p) where p = P(A)\).

Fundamental Bridge
The expectation of the indicator for event A is the probability of event A: \(E(I_A) = P(A)\).

Variance and Standard Deviation

[var(X) = E \left(X - E(X)\right)^2 = E(X^2) - (E(X))^2]

[\textrm{SD}(X) = \sqrt{var(X)}]

Continuous RVs, LOTUS, UoU

Continuous Random Variables (CRVs)

What is a Continuous Random Variable (CRV)?

A continuous random variable can take on any possible value within a certain interval (for example, [0, 1]), whereas a discrete random variable can only take on variables in a list of countable values (for example, all the integers, or the values 1, 1/2, 1/4, 1/8, etc.)

Do Continuous Random Variables have PMFs?

No. The probability that a continuous random variable takes on any specific value is 0.

What’s the probability that a CRV is in an interval?

Take the difference in CDF values (or use the PDF as described later).

[P(a \leq X \leq b) = P(X \leq b) - P(X \leq a) = F_X(b) - F_X(a)]

For X ~ N(μ, σ^2), this becomes

[P(a \leq X \leq b) = \Phi\left(\frac{b-\mu}{\sigma}\right) - \Phi\left(\frac{a-\mu}{\sigma}\right)]

What is the Probability Density Function (PDF)?

The PDF f is the derivative of the CDF F.

[F’(x) = f(x)]

A PDF is nonnegative and integrates to 1. By the fundamental theorem of calculus, to get from PDF back to CDF we can integrate:

\(F(x) = \int_{-\infty}^x f(t)dt\)

To find the probability that a CRV takes on a value in an interval, integrate the PDF over that interval.

[F(b) - F(a) = \int_a^b f(x)dx]

How do I find the expected value of a CRV?

Analogous to the discrete case, where you sum x times the PMF, for CRVs you integrate x times the PDF.

[E(X) = \int_{-\infty}^\infty xf(x)dx]

LOTUS

Expected value of a function of an r.v.

The expected value of X is defined this way:

\(E(X) = \sum_x xP(X=x) \text{ (for discrete X)}\) \(E(X) = \int_{-\infty}^\infty xf(x)dx \text{ (for continuous X)}\)

The Law of the Unconscious Statistician (LOTUS) states that you can find the expected value of a function of a random variable, g(X), in a similar way, by replacing the x in front of the PMF/PDF by g(x) but still working with the PMF/PDF of X:

\(E(g(X)) = \sum_x g(x)P(X=x) \text{ (for discrete X)}\) \(E(g(X)) = \int_{-\infty}^\infty g(x)f(x)dx \text{ (for continuous X)}\)

What’s a function of a random variable?

A function of a random variable is also a random variable. For example, if X is the number of bikes you see in an hour, then g(X) = 2X is the number of bike wheels you see in that hour, and h(X) = \({X \choose 2} = \frac{X(X-1)}{2}\) is the number of pairs of bikes such that you see both of those bikes in that hour.

What’s the point?

You don’t need to know the PMF/PDF of g(X) to find its expected value. All you need is the PMF/PDF of X.

Universality of Uniform (UoU)

When you plug any CRV into its own CDF, you get a Uniform(0,1) random variable. When you plug a Uniform(0,1) r.v. into an inverse CDF, you get an r.v. with that CDF. For example, let’s say that a random variable X has CDF

[F(x) = 1 - e^{-x}, \textrm{ for } x>0]

By UoU, if we plug X into this function then we get a uniformly distributed random variable.

[F(X) = 1 - e^{-X} \sim \textrm{Unif}(0,1)]

Similarly, if \(U ~ Unif(0,1)\) then \(F^{-1}(U)\) has CDF \(F\). The key point is that for any continuous random variable \(X\), we can transform it into a Uniform random variable and back by using its CDF.

Moments

Moments describe the shape of a distribution. Let X have mean μ and standard deviation σ, and Z=(X-μ)/σ be the standardized version of X. The kth moment of X is μₖ = E(Xᵏ), and the kth standardized moment of X is mₖ = E(Zᵏ). The mean, variance, skewness, and kurtosis are important summaries of the shape of a distribution.

 ֍ Mean: \(E(X) = μ₁\)
 ֍ Variance: \(var(X) = μ₂ - μ₁²\)
 ֍ Skewness: \(skew(X) = m₃\)
 ֍ Kurtosis: \(kurt(X) = m₄ - 3\)

Moment Generating Functions

• MGF (Moment Generating Function): For any random variable X, the function
  \(Mₓ(t) = E(e^{tX})\)
  is the moment generating function (MGF) of X, if it exists for all t in some open interval containing 0. The variable t could just as well have been called u or v. It’s a bookkeeping device that lets us work with the function Mₓ rather than the sequence of moments.
• Why is it called the Moment Generating Function? Because the kth derivative of the moment generating function, evaluated at 0, is the kth moment of X.
  \(μₖ = E(Xᵏ) = Mₓ⁽ᵏ⁾(0)\)
  This is true by Taylor expansion of \(e^{tX}\) since
  \(Mₓ(t) = E(e^{tX}) = ∑_(k=0)^(∞) [\frac{E(Xᵏ)tᵏ}{k!}] = ∑_(k=0)^(∞) [\frac{μₖtᵏ}{k!}]\)
• MGF of linear functions: If \(Y = aX + b\), then \(Mₙ(t) = E(e^(t(aX + b))) = e^(bt)Mₓ(at)\)
• Uniqueness: If it exists, the MGF uniquely determines the distribution. This means that for any two random variables X and Y, they are distributed the same (their PMFs/PDFs are equal) if and only if their MGFs are equal.
• Summing Independent RVs by Multiplying MGFs: If X and Y are independent, then
  \(Mₓ₊ᵧ(t) = E(e^(t(X + Y))) = E(e^{tX})E(e^(tY)) = Mₓ(t) ⋅ Mᵧ(t)\)
  The MGF of the sum of two random variables is the product of the MGFs of those two random variables.

Joint PDFs and CDFs

Joint Distributions

• Joint CDF: The joint cumulative distribution function (CDF) of X and Y is \(F(x,y) = P(X ≤ x, Y ≤ y)\).
• Joint PMF: In the discrete case, X and Y have a joint probability mass function \((PMF)\)
  \(pₓᵧ(x,y) = P(X=x, Y=y)\).
• Joint PDF: In the continuous case, X and Y have a joint probability density function \((PDF)\)
  \(fₓᵧ(x,y) = (∂²/∂x∂y)Fₓᵧ(x,y)\).
  The joint PMF/PDF must be nonnegative and sum/integrate to 1.

Conditional Distributions

 ֍ Conditioning and Bayes’ rule for discrete random variables:
\(P(Y=y|X=x) = P(X=x,Y=y) / P(X=x)\)
\(= P(X=x|Y=y)P(Y=y) / ∑ᵧ P(X=x|Y=ᵧ)P(Y=ᵧ)\)
 ֍ Conditioning and Bayes’ rule for continuous random variables:
\(fᵧ|ₓ(y|x) = fₓᵧ(x, y) / fₓ(x)\)
\(= (fₓ|ᵧ(x|y)fᵧ(y)) / fₓ(x)\)
 ֍ Hybrid Bayes’ rule:
\(fₓ(x|A) = (P(A | X = x)fₓ(x)) / P(A)\)

Marginal Distributions

To find the distribution of one (or more) random variables from a joint PMF/PDF, sum/integrate over the unwanted random variables.\

 ֍ Marginal PMF from joint PMF:
\(P(X = x) = ∑ₓ P(X=x, Y=y)\)
 ֍ Marginal PDF from joint PDF:
\(fₓ(x) = ∫[∞, -∞] fₓᵧ(x, y) dy\)

Independence of Random Variables

Random variables X and Y are independent if and only if any of the following conditions holds:
 ֍ Joint CDF is the product of the marginal CDFs.
 ֍ Joint PMF/PDF is the product of the marginal PMFs/PDFs.
 ֍ Conditional distribution of Y given X is the marginal distribution of Y. Write X ⫫ Y to denote that X and Y are independent.

Multivariate LOTUS

Law of the unconscious statistician (LOTUS) in more than one dimension is analogous to the 1D LOTUS.
For discrete random variables:
\(E(g(X, Y)) = ∑ₓ∑y g(x, y)P(X=x, Y=y)\)
For continuous random variables:
\(E(g(X, Y)) = ∫[-∞, ∞]∫[-∞, ∞] g(x, y)fₓᵧ(x, y)dxdy\)

Covariance and Transformations

Covariance and Correlation

Covariance is the analog of variance for two random variables. \(\text{cov}(X, Y) = E\left((X - E(X))(Y - E(Y))\right) = E(XY) - E(X)E(Y)\) Note that \(\text{cov}(X, X) = E(X^2) - (E(X))^2 = \text{var}(X)\)

Correlation is a standardized version of covariance that is always between \(-1\) and \(1\). \(\text{corr}(X, Y) = \frac{\text{cov}(X, Y)}{\sqrt{\text{var}(X)\text{var}(Y)}}\)

Covariance and Independence If two random variables are independent, then they are uncorrelated. The converse is not necessarily true (e.g., consider \(X \sim \mathcal{N}(0,1)\) and \(Y=X^2\)).
\(X \perp Y \longrightarrow \text{cov}(X, Y) = 0 \longrightarrow E(XY) = E(X)E(Y)\)

Covariance and Variance The variance of a sum can be found by \(var(X + Y) = var(X) + var(Y) + 2\text{cov}(X, Y)\) \(var(X_1 + X_2 + \dots + X_n ) = \sum_{i = 1}^{n}var(X_i) + 2\sum_{i < j} \text{cov}(X_i, X_j)\) If \(X\) and \(Y\) are independent, then they have covariance \(0\), so \(X \perp Y \Longrightarrow var(X + Y) = var(X) + var(Y)\) If \(X_1, X_2, \dots, X_n\) are identically distributed and have the same covariance relationships (often by symmetry), then \(var(X_1 + X_2 + \dots + X_n ) = nvar(X_1) + 2{n \choose 2}\text{cov}(X_1, X_2)\)

Covariance Properties For random variables \(W, X, Y, Z\) and constants \(a, b\): \(\text{cov}(X, Y) = \text{cov}(Y, X)\) \(\text{cov}(X + a, Y + b) = \text{cov}(X, Y)\) \(\text{cov}(aX, bY) = ab\text{cov}(X, Y)\) \(\text{cov}(W + X, Y + Z) = \text{cov}(W, Y) + \text{cov}(W, Z) + \text{cov}(X, Y) + \text{cov}(X, Z)\)

Correlation is location-invariant and scale-invariant For any constants \(a,b,c,d\) with \(a\) and \(c\) nonzero, \(\text{corr}(aX + b, cY + d) = \text{corr}(X, Y)\)

Transformations

One Variable Transformations Let’s say that we have a random variable \(X\) with PDF \(f_X(x)\), but we are also interested in some function of \(X\). We call this function \(Y = g(X)\). Also let \(y=g(x)\). If \(g\) is differentiable and strictly increasing (or strictly decreasing), then the PDF of \(Y\) is \(f_Y(y) = f_X(x)\left|\frac{dx}{dy}\right| = f_X(g^{-1}(y))\left|\frac{d}{dy}g^{-1}(y)\right|\) The derivative of the inverse transformation is called the Jacobian.

Two Variable Transformations Similarly, let’s say we know the joint PDF of \(U\) and \(V\) but are also interested in the random vector \((X, Y)\) defined by \((X, Y) = g(U, V)\). Let \(\frac{\partial (u,v)}{\partial (x,y)} = \begin{pmatrix} \frac{\partial u}{\partial x} & \frac{\partial u}{\partial y} \\ \frac{\partial v}{\partial x} & \frac{\partial v}{\partial y} \\ \end{pmatrix}\) be the Jacobian matrix. If the entries in this matrix exist and are continuous, and the determinant of the matrix is never \(0\), then \(f_{X,Y}(x, y) = f_{U,V}(u,v) \left|\left| \frac{\partial (u,v)}{\partial (x,y)}\right| \right|\) The inner bars tell us to take the matrix’s determinant, and the outer bars tell us to take the absolute value. In a \(2 \times 2\) matrix, \(\left| \left| \begin{array}{ccc} a & b \\ c & d \end{array} \right| \right| = |ad - bc|\)

Convolutions

Convolution Integral If you want to find the PDF of the sum of two independent CRVs \(X\) and \(Y\), you can do the following integral: \(f_{X+Y}(t)=\int_{-\infty}^\infty f_X(x)f_Y(t-x)dx\)

Example Let \(X,Y \sim \mathcal{N}(0,1)\) be i.i.d. Then for each fixed \(t\), \(f_{X+Y}(t)=\int_{-\infty}^\infty \frac{1}{\sqrt{2\pi}}e^{-x^2/2} \frac{1}{\sqrt{2\pi}}e^{-(t-x)^2/2} dx\) By completing the square and using the fact that a Normal PDF integrates to \(1\), this works out to \(f_{X+Y}(t)\) being the \(\mathcal{N}(0,2)\) PDF.

Poisson Process

Definition

We have a Poisson process of rate \(\lambda\) arrivals per unit time if the following conditions hold:

1. The number of arrivals in a time interval of length \(t\) is \(\text{Pois}(\lambda t)\).
2. Numbers of arrivals in disjoint time intervals are independent.

For example, the numbers of arrivals in the time intervals \([0,5]\), \((5,12),\) and \([13,23)\) are independent with \(\text{Pois}(5\lambda)\), \(\text{Pois}(7\lambda)\), and \(\text{Pois}(10\lambda)\) distributions, respectively.

Count-Time Duality

Consider a Poisson process of emails arriving in an inbox at rate \(\lambda\) emails per hour. Let \(T_n\) be the time of arrival of the \(n\)th email (relative to some starting time \(0\)) and \(N_t\) be the number of emails that arrive in \([0,t]\).

Let’s find the distribution of \(T_1\). The event \(T_1 > t\), the event that you have to wait more than \(t\) hours to get the first email, is the same as the event \(N_t = 0\), which is the event that there are no emails in the first \(t\) hours. So,

[P(T_1 > t) = P(N_t = 0) = e^{-\lambda t}]

Therefore, \(P(T_1 \leq t) = 1 - e^{-\lambda t}\), and \(T_1\) follows an exponential distribution with parameter \(\lambda\).

By the memoryless property and similar reasoning, the interarrival times between emails are i.i.d. exponential random variables with parameter \(\lambda\), i.e., the differences \(T_n - T_{n-1}\) are i.i.d. exponential random variables with parameter \(\lambda\).

Order Statistics

Definition

Let’s say you have \(n\) i.i.d. random variables \(X_1, X_2, \dots, X_n\). If you arrange them from smallest to largest, the \(i\)th element in that list is the \(i\)th order statistic, denoted \(X_{(i)}\). So \(X_{(1)}\) is the smallest in the list and \(X_{(n)}\) is the largest in the list.

Note that the order statistics are dependent, e.g., learning \(X_{(4)} = 42\) gives us the information that \(X_{(1)},X_{(2)},X_{(3)}\) are \(\leq 42\) and \(X_{(5)},X_{(6)},\dots,X_{(n)}\) are \(\geq 42\).

Distribution

Taking \(n\) i.i.d. random variables \(X_1, X_2, \dots, X_n\) with CDF \(F(x)\) and PDF \(f(x)\), the CDF and PDF of \(X_{(i)}\) are: \(F_{X_{(i)}}(x) = P (X_{(i)} \leq x) = \sum_{k=i}^n {n \choose k} F(x)^k(1 - F(x))^{n - k}\) \(f_{X_{(i)}}(x) = n{n - 1 \choose i - 1}F(x)^{i-1}(1 - F(x))^{n-i}f(x)\)

Uniform Order Statistics

The \(j\)th order statistic of i.i.d. \(U_1,\dots,U_n \sim \text{Unif}(0,1)\) is \(U_{(j)} \sim \text{Beta}(j, n - j + 1)\).

Conditional Expectation

Conditioning on an Event

We can find \(E(Y|A)\), the expected value of \(Y\) given that event \(A\) occurred. A very important case is when \(A\) is the event \(X=x\). Note that \(E(Y|A)\) is a number.

For example:

• The expected value of a fair die roll, given that it is prime, is \(\frac{1}{3} \cdot 2 + \frac{1}{3} \cdot 3 + \frac{1}{3} \cdot 5 = \frac{10}{3}\).
• Let \(Y\) be the number of successes in \(10\) independent Bernoulli trials with probability \(p\) of success. Let \(A\) be the event that the first \(3\) trials are all successes. Then \(E(Y\|A) = 3 + 7p\) since the number of successes among the last \(7\) trials is \(\text{Bin}(7,p)\).
• Let \(T \sim \text{Expo}(1/10)\) be how long you have to wait until the shuttle comes. Given that you have already waited \(t\) minutes, the expected additional waiting time is \(10\) more minutes, by the memoryless property. That is, \(E(T\|T>t) = t + 10\).

Conditioning on a Random Variable

We can also find \(E(Y|X)\), the expected value of \(Y\) given the random variable \(X\). This is a function of the random variable \(X\). It is not a number except in certain special cases such as if \(X \perp Y\). To find \(E(Y|X)\), find \(E(Y|X = x)\) and then plug in \(X\) for \(x\).

For example:

• If \(E(Y\|X=x) = x^3+5x\), then \(E(Y\|X) = X^3 + 5X\).
• Let \(Y\) be the number of successes in \(10\) independent Bernoulli trials with probability \(p\) of success and \(X\) be the number of successes among the first \(3\) trials. Then \(E(Y\|X)=X+7p\).
• Let \(X \sim \mathcal{N}(0,1)\) and \(Y=X^2\). Then \(E(Y\|X=x) = x^2\) since if we know \(X=x\) then we know \(Y=x^2\). And \(E(X\|Y=y) = 0\) since if we know \(Y=y\) then we know \(X = \pm \sqrt{y}\), with equal probabilities (by symmetry). So \(E(Y\|X)=X^2\), \(E(X\|Y)=0\).

Properties of Conditional Expectation

1. \(E(Y\|X) = E(Y)\) if \(X \perp Y\)
2. \(E(h(X)W\|X) = h(X)E(W\|X)\) (taking out what’s known)
   In particular, \(E(h(X)\|X) = h(X)\).
3. \(E(E(Y\|X)) = E(Y)\) (Adam’s Law, a.k.a. Law of Total Expectation)

Adam’s Law (a.k.a. Law of Total Expectation)

For any events \(A_1, A_2, \dots, A_n\) that partition the sample space: \(E(Y) = E(Y|A_1)P(A_1) + \dots + E(Y|A_n)P(A_n)\)

For the special case where the partition is \(A, A^c\), this says: \(E(Y) = E(Y|A)P(A) + E(Y|A^c)P(A^c)\)

Eve’s Law (a.k.a. Law of Total Variance)

\(\text{Var}(Y) = E(\text{Var}(Y|X)) + \text{Var}(E(Y|X))\)

MVN, LLN, CLT

Law of Large Numbers (LLN)

Let \(X_1, X_2, X_3, \dots\) be i.i.d. with mean \(\mu\). The sample mean is \(\bar{X}_n = \frac{X_1 + X_2 + X_3 + \dots + X_n}{n}\). The Law of Large Numbers states that as \(n \to \infty\), \(\bar{X}_n \to \mu\) with probability \(1\). For example, in flips of a coin with probability \(p\) of Heads, let \(X_j\) be the indicator of the \(j\)th flip being Heads. Then LLN says the proportion of Heads converges to \(p\) (with probability \(1\)).

Central Limit Theorem (CLT)

Approximation using CLT

We use \(\dot{\,\sim\,}\) to denote “is approximately distributed.” We can use the Central Limit Theorem to approximate the distribution of a random variable \(Y = X_1 + X_2 + \dots + X_n\) that is a sum of \(n\) i.i.d. random variables \(X_i\). Let \(E(Y) = \mu_Y\) and \(\text{Var}(Y) = \sigma^2_Y\). The CLT says: \(Y \dot{\,\sim\,} \mathcal{N}(\mu_Y, \sigma^2_Y)\)

If the \(X_i\) are i.i.d. with mean \(\mu_X\) and variance \(\sigma^2_X\), then \(\mu_Y = n \mu_X\) and \(\sigma^2_Y = n \sigma^2_X\). For the sample mean \(\bar{X}_n\), the CLT says: \(\bar{X}_n = \frac{1}{n}(X_1 + X_2 + \dots + X_n) \dot{\,\sim\,} \mathcal{N}(\mu_X, \frac{\sigma^2_X}{n})\)

Asymptotic Distributions using CLT

We use \(\xrightarrow{D}\) to denote “converges in distribution to” as \(n \to \infty\). The CLT says that if we standardize the sum \(X_1 + \dots + X_n\), then the distribution of the sum converges to \(\mathcal{N}(0,1)\) as \(n \to \infty\): \(\frac{1}{\sigma\sqrt{n}}(X_1 + \dots + X_n - n\mu_X) \xrightarrow{D} \mathcal{N}(0, 1)\) In other words, the CDF of the left-hand side goes to the standard Normal CDF, \(\Phi\). In terms of the sample mean, the CLT says: \(\frac{\sqrt{n}(\bar{X}_n - \mu_X)}{\sigma_X} \xrightarrow{D} \mathcal{N}(0, 1)\)

Markov Chains

Definition

A Markov chain is a random walk in a state space, which we will assume is finite, say \(\{1, 2, \dots, M\}\). We let \(X_t\) denote which element of the state space the walk is visiting at time \(t\). The Markov chain is the sequence of random variables tracking where the walk is at all points in time, \(X_0, X_1, X_2, \dots\). By definition, a Markov chain must satisfy the Markov property, which says that if you want to predict where the chain will be at a future time, if we know the present state then the entire past history is irrelevant. Given the present, the past and future are conditionally independent. In symbols: \(P(X_{n+1} = j | X_0 = i_0, X_1 = i_1, \dots, X_n = i) = P(X_{n+1} = j | X_n = i)\)

State Properties

A state is either recurrent or transient.

• If you start at a recurrent state, then you will always return back to that state at some point in the future. \textmusicalnote You can check-out any time you like, but you can never leave. \textmusicalnote
• Otherwise, you are at a transient state. There is some positive probability that once you leave, you will never return. \textmusicalnote You don’t have to go home, but you can’t stay here. \textmusicalnote

A state is either periodic or aperiodic.

• If you start at a periodic state of period \(k\), then the GCD of the possible numbers of steps it would take to return back is \(k > 1\).
• Otherwise, you are at an aperiodic state. The GCD of the possible numbers of steps it would take to return back is 1.

Transition Matrix

Let the state space be \(\{1,2,\dots,M\}\). The transition matrix \(Q\) is the \(M \times M\) matrix where element \(q_{ij}\) is the probability that the chain goes from state \(i\) to state \(j\) in one step:

To find the probability that the chain goes from state \(i\) to state \(j\) in exactly \(m\) steps, take the \((i, j)\) element of \(Q^m\):

If \(X_0\) is distributed according to the row vector PMF \(\vec{p}\), i.e., \(p_j = P(X_0 = j)\), then the PMF of \(X_n\) is \(\vec{p}Q^n\).

Chain Properties

A chain is irreducible if you can get from anywhere to anywhere. If a chain (on a finite state space) is irreducible, then all of its states are recurrent. A chain is periodic if any of its states are periodic, and is aperiodic if none of its states are periodic. In an irreducible chain, all states have the same period.

A chain is reversible with respect to \(\vec{s}\) if \(s_iq_{ij} = s_jq_{ji}\) for all \(i, j\). Examples of reversible chains include any chain with \(q_{ij} = q_{ji}\), with \(\vec{s} = (\frac{1}{M}, \frac{1}{M}, \dots, \frac{1}{M})\), and random walk on an undirected network.

Stationary Distribution

Let \(\vec{s} = (s_1, s_2, \dots, s_M)\) be a PMF (written as a row vector). We will call \(\vec{s}\) the stationary distribution for the chain if \(\vec{s}Q = \vec{s}\). As a consequence, if \(X_t\) has the stationary distribution, then all future \(X_{t+1}, X_{t + 2}, \dots\) also have the stationary distribution.

For irreducible, aperiodic chains, the stationary distribution exists, is unique, and \(s_i\) is the long-run probability of a chain being at state \(i\). The expected number of steps to return to \(i\) starting from \(i\) is \(1/s_i\).

To find the stationary distribution, you can solve the matrix equation \((Q' - I){\vec{s}\,}'= 0\). The stationary distribution is uniform if the columns of \(Q\) sum to 1.

Reversibility Condition Implies Stationarity: If you have a PMF \(\vec{s}\) and a Markov chain with transition matrix \(Q\), then \(s_iq_{ij} = s_jq_{ji}\) for all states \(i, j\) implies that \(\vec{s}\) is stationary.

Random Walk on an Undirected Network

If you have a collection of nodes, pairs of which can be connected by undirected edges, and a Markov chain is run by going from the current node to a uniformly random node that is connected to it by an edge, then this is a random walk on an undirected network. The stationary distribution of this chain is proportional to the degree sequence (this is the sequence of degrees, where the degree of a node is how many edges are attached to it). For example, the stationary distribution of random walk on the network shown above is proportional to \((3,3,2,4,2)\), so it’s \((\frac{3}{14}, \frac{3}{14}, \frac{2}{14}, \frac{4}{14}, \frac{2}{14})\).

Continuous Distributions

Uniform Distribution

Let \(U\) be distributed \(\text{Unif}(a, b)\). We know the following:

• Properties of the Uniform: For a Uniform distribution, the probability of a draw from any interval within the support is proportional to the length of the interval.

• Example: William throws darts really badly, so his darts are uniform over the whole room because they’re equally likely to appear anywhere. William’s darts have a Uniform distribution on the surface of the room. The Uniform is the only distribution where the probability of hitting in any specific region is proportional to the length/area/volume of that region, and where the density of occurrence in any one specific spot is constant throughout the whole support.

• PDF and CDF (top is Unif(0, 1), bottom is Unif(a, b))

For the Uniform distribution:

• Unif(0, 1):
  • PDF: \(f(x) = \left\{ \begin{array}{lr} 1 & x \in [0, 1] \\ 0 & x \notin [0, 1] \end{array} \right.\)
  • CDF: \(F(x) = \left\{ \begin{array}{lr} 0 & x < 0 \\ x & x \in [0, 1] \\ 1 & x > 1 \end{array} \right.\)
• Unif(a, b):
  • PDF: \(f(x) = \left\{ \begin{array}{lr} \frac{1}{b-a} & x \in [a, b] \\ 0 & x \notin [a, b] \end{array} \right.\)
  • CDF: \(F(x) = \left\{ \begin{array}{lr} 0 & x < a \\ \frac{x-a}{b-a} & x \in [a, b] \\ 1 & x > b \end{array} \right.\)

    Normal Distribution

Let us say that X is distributed N(μ, σ^2). We know the following:

Central Limit Theorem The Normal distribution is ubiquitous because of the Central Limit Theorem, which states that the sample mean of i.i.d. r.v.s will approach a Normal distribution as the sample size grows, regardless of the initial distribution.

Location-Scale Transformation Every time we shift a Normal r.v. (by adding a constant) or rescale a Normal (by multiplying by a constant), we change it to another Normal r.v. For any Normal X ~ N(μ, σ^2), we can transform it to the standard N(0, 1) by the following transformation: Z = (X - μ) / σ ~ N(0, 1)

Standard Normal The Standard Normal, Z ~ N(0, 1), has mean 0 and variance 1. Its CDF is denoted by Φ.

Exponential Distribution

Let us say that X is distributed Expo(λ). We know the following:

Story You’re sitting on an open meadow right before the break of dawn, wishing that airplanes in the night sky were shooting stars, because you could really use a wish right now. You know that shooting stars come on average every 15 minutes, but a shooting star is not “due” to come just because you’ve waited so long. Your waiting time is memoryless; the additional time until the next shooting star comes does not depend on how long you’ve waited already.

Example The waiting time until the next shooting star is distributed Expo(4) hours. Here λ=4 is the rate parameter, since shooting stars arrive at a rate of 1 per 1/4 hour on average. The expected time until the next shooting star is 1/λ = 1/4 hour.

Expos as a rescaled Expo(1) Y ~ Expo(λ) ⇒ X = λY ~ Expo(1)

Memorylessness The Exponential Distribution is the only continuous memoryless distribution. The memoryless property says that for X ~ Expo(λ) and any positive numbers s and t: P(X > s + t | X > s) = P(X > t) Equivalently, (X - a | X > a) ~ Expo(λ)

Min of Expos If we have independent Xi ~ Expo(λi), then min(X1, …, Xk) ~ Expo(λ1 + λ2 + … + λk).

Max of Expos If we have i.i.d. Xi ~ Expo(λ), then max(X1, …, Xk) has the same distribution as Y1 + Y2 + … + Yk, where Yj ~ Expo(jλ) and the Yj are independent.

Gamma Distribution

Let us say that X is distributed Gam(a, λ). We know the following:

Story You sit waiting for shooting stars, where the waiting time for a star is distributed Expo(λ). You want to see n shooting stars before you go home. The total waiting time for the nth shooting star is Gam(n, λ).

Example You are at a bank, and there are 3 people ahead of you. The serving time for each person is Exponential with mean 2 minutes. Only one person at a time can be served. The distribution of your waiting time until it’s your turn to be served is Gam(3, 1/2).

Beta Distribution

Conjugate Prior of the Binomial In the Bayesian approach to statistics, parameters are viewed as random variables, to reflect our uncertainty. The prior for a parameter is its distribution before observing data. The posterior is the distribution for the parameter after observing data. Beta is the conjugate prior of the Binomial because if you have a Beta-distributed prior on p in a Binomial, then the posterior distribution on p given the Binomial data is also Beta-distributed. Consider the following two-level model: X | p ~ Bin(n, p) p ~ Beta(a, b) Then after observing X = x, we get the posterior distribution p | (X = x) ~ Beta(a + x, b + n - x)

Order statistics of the Uniform See Order Statistics.

Beta-Gamma relationship If X ~ Gam(a, λ), Y ~ Gam(b, λ), with X independent of Y, then:

• X / (X + Y) ~ Beta(a, b)
• X + Y independent of X / (X + Y)

Chi-Square Distribution

Let us say that X is distributed chi2_n. We know the following:

Story A Chi-Square(n) is the sum of the squares of n independent standard Normal r.v.s.

Properties and Representations

• \(X\) is distributed as \(Z1^2 + Z2^2 + ... + Zn^2\) for i.i.d. \(Z_i ~ N(0,1)\)\
• \[X ~ Gam(n/2, 1/2)\]

Discrete Distributions

Distributions for four sampling schemes

Bernoulli Distribution

The Bernoulli distribution is the simplest case of the Binomial distribution, where we only have one trial (\(n=1\)). Let us say that X is distributed Bern(p). We know the following:

Story A trial is performed with probability p of “success”, and X is the indicator of success: 1 means success, 0 means failure.

Example Let X be the indicator of Heads for a fair coin toss. Then X follows the Bernoulli distribution with parameter p=1/2. Also, 1 - X follows the Bernoulli distribution with parameter p=1/2, representing the indicator of Tails.

Binomial Distribution

Let us say that X is distributed Bin(n, p). We know the following:

Story X is the number of “successes” that we will achieve in n independent trials, where each trial is either a success or a failure, each with the same probability p of success. X can be expressed as the sum of multiple independent Bernoulli random variables with parameter p. If X ~ Bin(n, p) and Xj ~ Bern(p), where all the Bernoullis are independent, then: X = X1 + X2 + X3 + … + Xn

Example If Jeremy Lin makes 10 free throws, and each throw independently has a 3/4 chance of getting in, then the number of successful throws is distributed as Bin(10, 3/4).

Properties

• Redefine success: If X ~ Bin(n, p), then n - X ~ Bin(n, 1 - p)
• Sum: If X ~ Bin(n, p) and Y ~ Bin(m, p) with X and Y being independent, then X + Y ~ Bin(n + m, p)

• Conditional: If X ~ Bin(n, p) and Y ~ Bin(m, p) with X and Y being independent, then X

  (X + Y = r) ~ Hypergeometric(n, m, r)

• Binomial-Poisson Relationship: Bin(n, p) is approximately Pois(np) if p is small.
• Binomial-Normal Relationship: Bin(n, p) is approximately N(np, np(1 - p)) if n is large and p is not near 0 or 1.

Geometric Distribution

Let us say that X is distributed Geom(p). We know the following:

Story X is the number of “failures” that we will achieve before we achieve our first success. The successes have a probability p.

Example If each Pokéball we throw has a 1/10 probability of catching Mew, then the number of failed Pokéballs before catching Mew follows the Geometric distribution with parameter p = 1/10.

The PMF of a Poisson distribution is given by:

[P(X = k) = \frac{e^{-\lambda}\lambda^k}{k!}]

where X is the random variable following a Poisson distribution, and λ is the average rate of events occurring per unit space or time.

Multivariate Distributions

Multinomial Distribution

Let us say that the vector \(\vec{X} = (X_1, X_2, X_3, \dots, X_k) \sim \text{Mult}_k(n, \vec{p})\) where \(\vec{p} = (p_1, p_2, \dots, p_k)\).

• Story: We have \(n\) items, which can fall into any one of the \(k\) buckets independently with the probabilities \(\vec{p} = (p_1, p_2, \dots, p_k)\).
• Example: Let us assume that every year, 100 students in the Harry Potter Universe are randomly and independently sorted into one of four houses with equal probability. The number of people in each of the houses is distributed \(\text{Mult}_4(100, \vec{p})\), where \(\vec{p} = (0.25, 0.25, 0.25, 0.25)\). Note that \(X_1 + X_2 + \dots + X_4 = 100\), and they are dependent.
• Joint PMF: For \(n = n_1 + n_2 + \dots + n_k\), the joint probability mass function is: \(P(\vec{X} = \vec{n}) = \frac{n!}{n_1!n_2!\dots n_k!}p_1^{n_1}p_2^{n_2}\dots p_k^{n_k}\)
• Marginal PMF, Lumping, and Conditionals: Marginally, \(X_i \sim \text{Bin}(n,p_i)\) since we can define “success” to mean category \(i\). If you lump together multiple categories in a Multinomial, then it is still Multinomial. Conditioning on some \(X_j\) also gives a Multinomial.
• Variances and Covariances: We have \(X_i \sim \text{Bin}(n, p_i)\) marginally, so \(\text{Var}(X_i) = np_i(1-p_i)\). Also, \(\text{Cov}(X_i, X_j) = -np_ip_j\) for $$i \neq j.

Multivariate Uniform Distribution

See the univariate Uniform for stories and examples. For the 2D Uniform on some region, probability is proportional to area. Every point in the support has equal density, of value \(1/\text{area of region}\). For the 3D Uniform, probability is proportional to volume.

Multivariate Normal (MVN) Distribution

A vector \(\vec{X} = (X_1, X_2, \dots, X_k)\) is Multivariate Normal if every linear combination is Normally distributed, i.e., \(t_1X_1 + t_2X_2 + \dots + t_kX_k\) is Normal for any constants \(t_1, t_2, \dots, t_k\). The parameters of the Multivariate Normal are the mean vector \(\vec{\mu} = (\mu_1, \mu_2, \dots, \mu_k)\) and the covariance matrix where the \((i, j)\) entry is \(\text{cov}(X_i, X_j)\).

• Properties: The Multivariate Normal has the following properties:
  • Any subvector is also MVN.
  • If any two elements within an MVN are uncorrelated, then they are independent.
  • The joint PDF of a Bivariate Normal \((X,Y)\) with \(\mathcal{N}(0,1)\) marginal distributions and correlation \(\rho \in (-1,1)\) is: \(f_{X,Y}(x,y) = \frac{1}{2 \pi \tau} \exp\left(-\frac{1}{2 \tau^2} (x^2+y^2-2 \rho xy)\right),\) with \(\tau = \sqrt{1-\rho^2}\).

Distribution Properties

Important CDFs

• Standard Normal: \(\Phi\)
• Exponential(\(\lambda\)): \(F(x) = 1 - e^{-\lambda x}\), for \(x \in (0, \infty)\)
• Uniform(0,1): \(F(x) = x\), for \(x \in (0, 1)\)

Convolutions of Random Variables

A convolution of \(n\) random variables is simply their sum. For the following results, let \(X\) and \(Y\) be independent.

1. \(X \sim \text{Pois}(\lambda_1)\), \(Y \sim \text{Pois}(\lambda_2)\) \(\longrightarrow X + Y \sim \text{Pois}(\lambda_1 + \lambda_2)\)
2. \(X \sim \text{Bin}(n_1, p)\), \(Y \sim \text{Bin}(n_2, p)\) \(\longrightarrow X + Y \sim \text{Bin}(n_1 + n_2, p)\). \(\text{Bin}(n,p)\) can be thought of as a sum of i.i.d. \(\text{Bern}(p)\) random variables.
3. \(X \sim \text{Gam}(a_1, \lambda)\), \(Y \sim \text{Gam}(a_2, \lambda)\) \(\longrightarrow X + Y \sim \text{Gam}(a_1 + a_2, \lambda)\). \(\text{Gam}(n,\lambda)\) with \(n\) an integer can be thought of as a sum of i.i.d. \(\text{Expo}(\lambda)\) random variables.
4. \(X \sim \text{NBin}(r_1, p)\), \(Y \sim \text{NBin}(r_2, p)\) \(\longrightarrow X + Y \sim \text{NBin}(r_1 + r_2, p)\). \(\text{NBin}(r,p)\) can be thought of as a sum of i.i.d. \(\text{Geom}(p)\) random variables.
5. \(X \sim \mathcal{N}(\mu_1, \sigma_1^2)\), \(Y \sim \mathcal{N}(\mu_2, \sigma_2^2)\) \(\longrightarrow X + Y \sim \mathcal{N}(\mu_1 + \mu_2, \sigma_1^2 + \sigma_2^2)\)

Special Cases of Distributions

1. \(\text{Bin}(1, p) \sim \text{Bern}(p)\).
2. \(\text{Beta}(1, 1) \sim \text{Unif}(0, 1)\).
3. \(\text{Gam}(1, \lambda) \sim \text{Expo}(\lambda)\).
4. \(\chi^2_n \sim \text{Gam}\left(\frac{n}{2}, \frac{1}{2}\right)\).
5. \(\text{NBin}(1, p) \sim \text{Geom}(p)\).

[\textbf{Inequalities}]

1. Cauchy-Schwarz: \(|E(XY)| \leq \sqrt{E(X^2)E(Y^2)}\)

2. Markov: \(P(X \geq a) \leq \frac{E|X|}{a}\) for \(a>0\)

3. Chebyshev: \(P(|X - \mu| \geq a) \leq \frac{\sigma^2}{a^2}\) for \(E(X)=\mu\), \(\text{Var}(X) = \sigma^2\)

4. Jensen: \(E(g(X)) \geq g(E(X))\) for \(g\) convex; reverse if \(g\) is concave

[\textbf{Formulas}]

[\textbf{Geometric Series}]

[1 + r + r^2 + \dots + r^{n-1} = \sum_{k=0}^{n-1} r^k = \frac{1 - r^n}{1 -r}]

[\textbf{Exponential Function ($e^x$)}]

[e^x = \sum_{n=0}^\infty \frac{x^n}{n!}= 1 + x + \frac{x^2}{2!} + \frac{x^3}{3!} + \dots = \lim_{n \rightarrow \infty} \left( 1 + \frac{x}{n} \right)^n]

[\textbf{Gamma and Beta Integrals}]

You can sometimes solve complicated-looking integrals by pattern-matching to a gamma or beta integral:

[\int_0^\infty x^{t-1}e^{-x}\, dx = \Gamma(t) \hspace{1 cm} \int_0^1 x^{a - 1}(1-x)^{b-1}\, dx = \frac{\Gamma(a)\Gamma(b)}{\Gamma(a + b)}]

Also, \(\Gamma(a+1) = a \Gamma(a)\), and \(\Gamma(n) = (n - 1)!\) if \(n\) is a positive integer.

[\textbf{Euler’s Approximation for Harmonic Sums}]

[1 + \frac{1}{2} + \frac{1}{3} + \dots + \frac{1}{n} \approx \log n + 0.577 \dots]

[\textbf{Stirling’s Approximation for Factorials}]

[n! \approx \sqrt{2\pi n}\left(\frac{n}{e}\right)^n]

Miscellaneous Definitions

Medians and Quantiles Let \(X\) have CDF \(F\). Then \(X\) has median \(m\) if \(F(m) \geq 0.5\) and \(P(X \geq m) \geq 0.5\). For \(X\) continuous, \(m\) satisfies \(F(m) = \frac{1}{2}\). In general, the \(a\)th quantile of \(X\) is \(\min \{x: F(x) \geq a\}\); the median is the case \(a = \frac{1}{2}\).

log Statisticians generally use \(\log\) to refer to natural log (i.e., base \(e\)).

i.i.d r.v.s Independent, identically-distributed random variables.

Example Problems

Calculating Probability

A textbook has \(n\) typos, which are randomly scattered amongst its \(n\) pages, independently. You pick a random page. What is the probability that it has no typos?

Solution

There is a \(\left(1 - \frac{1}{n}\right)\) probability that any specific typo isn’t on your page, and thus a \(\boxed{\left(1 - \frac{1}{n}\right)^n}\) probability that there are no typos on your page. For \(n\) large, this is approximately \(e^{-1} = \frac{1}{e}\).

Linearity and Indicators (1)

In a group of \(n\) people, what is the expected number of distinct birthdays (month and day)? What is the expected number of birthday matches?

Solution

Let \(X\) be the number of distinct birthdays and \(I_j\) be the indicator for the \(j\)th day being represented.
\(E(I_j) = 1 - P(\text{no one born on day }j) = 1 - \left(\frac{364}{365}\right)^n\)
By linearity, \(\boxed{E(X) = 365\left(1-\frac{364}{365}\right)^n}\).\ Now let \(Y\) be the number of birthday matches and \(J_i\) be the indicator that the \(i\)th pair of people have the same birthday. The probability that any two specific people share a birthday is \(\frac{1}{365}\), so \(\boxed{E(Y) = \frac{\binom{n}{2}}{365}}\).

Linearity and Indicators (2)

This problem is commonly known as the hat-matching problem. There are \(n\) people at a party, each with a hat. At the end of the party, they each leave with a random hat. What is the expected number of people who leave with the right hat?

Solution

Each hat has a \(\frac{1}{n}\) chance of going to the right person. By linearity, the average number of hats that go to their owners is \(\boxed{1}\).

Linearity and First Success

This problem is commonly known as the coupon collector problem. There are \(n\) coupon types. At each draw, you get a uniformly random coupon type. What is the expected number of coupons needed until you have a complete set?

Solution

Let \(N\) be the number of coupons needed; we want \(E(N)\). Let \(N = N_1 + \dots + N_n\), where \(N_1\) is the draws to get our first new coupon, \(N_2\) is the additional draws needed to draw our second new coupon, and so on. By the story of the First Success, \(N_2 \sim \text{FS}\left(\frac{n-1}{n}\right)\) (after collecting the first coupon type, there’s \(\frac{n-1}{n}\) chance you’ll get something new). Similarly, \(N_3 \sim \text{FS}\left(\frac{n-2}{n}\right)\), and \(N_j \sim \text{FS}\left(\frac{n-j+1}{n}\right)\). By linearity,
\(E(N) = E(N_1) + \dots + E(N_n) = \frac{n}{n} + \frac{n}{n-1} + \dots + \frac{n}{1} = \boxed{n\sum^n_{j=1} \frac{1}{j}}\) This is approximately \(n (\log(n) + 0.577)\) by Euler’s approximation.

Orderings of i.i.d. random variables

I call 2 UberX’s and 3 Lyfts at the same time. If the time it takes for the rides to reach me are i.i.d., what is the probability that all the Lyfts will arrive first?

Solution

Since the arrival times of the five cars are i.i.d., all \(5!\) orderings of the arrivals are equally likely. There are \(3!2!\) orderings that involve the Lyfts arriving first, so the probability that the Lyfts arrive first is \(\boxed{\frac{3!2!}{5!} = \frac{1}{10}}\). Alternatively, there are \(\binom{5}{3}\) ways to choose 3 of the 5 slots for the Lyfts to occupy, where each of the choices is equally likely. One of these choices has all 3 of the Lyfts arriving first, so the probability is \(\boxed{\frac{1}{5 \choose 3} = \frac{1}{10}}\).

Expectation of Negative Hypergeometric

What is the expected number of cards that you draw before you pick your first Ace in a shuffled deck (not counting the Ace)?

Solution

Consider a non-Ace. Denote this to be card \(j\). Let \(I_j\) be the indicator that card \(j\) will be drawn before the first Ace. Note that \(I_j=1\) says that \(j\) is before all 4 of the Aces in the deck. The probability that this occurs is \(\frac{1}{5}\) by symmetry. Let \(X\) be the number of cards drawn before the first Ace. Then \(X = I_1 + I_2 + \ldots + I_{48}\), where each indicator corresponds to one of the 48 non-Aces. Thus,
\(E(X) = E(I_1) + E(I_2) + \ldots + E(I_{48}) = \frac{48}{5} = \boxed{9.6}\).

Minimum and Maximum of RVs

What is the CDF of the maximum of \(n\) independent Unif(0,1) random variables?

Solution

Note that for r.v.s \(X_1,X_2,\dots,X_n\),
\(P(\min(X_1, X_2, \dots, X_n) \geq a) = P(X_1 \geq a, X_2 \geq a, \dots, X_n \geq a)\)
Similarly,
\(P(\max(X_1, X_2, \dots, X_n) \leq a) = P(X_1 \leq a, X_2 \leq a, \dots, X_n \leq a)\)
We will use this principle to find the CDF of \(U_{(n)}\), where \(U_{(n)} = \max(U_1, U_2, \dots, U_n)\) and \(U_i \sim \text{Unif}(0, 1)\) are i.i.d. \(\ P(\max(U_1, U_2, \dots, U_n) \leq a) = P(U_1 \leq a, U_2 \leq a, \dots, U_n \leq a) = P(U_1 \leq a)P(U_2 \leq a)\dots P(U_n \leq a) = \boxed{a^n}\)
for \(0<a<1\) (and the CDF is \(0\) for \(a \leq 0\) and \(1\) for \(a \geq 1\)).

Pattern-matching with \(e^x\) Taylor series

Solution

For \(X \sim \text{Pois}(\lambda)\), find \(E\left(\frac{1}{X+1}\right)\). By LOTUS,
\(E\left(\frac{1}{X+1}\right) = \sum_{k=0}^\infty \frac{1}{k+1} \frac{e^{-\lambda}\lambda^k}{k!} = \frac{e^{-\lambda}}{\lambda}\sum_{k=0}^\infty \frac{\lambda^{k+1}}{(k+1)!} = \boxed{\frac{e^{-\lambda}}{\lambda}(e^\lambda-1)}\)

Adam’s Law and Eve’s Law

William really likes speedsolving Rubik’s Cubes. But he’s pretty bad at it, so sometimes he fails. On any given day, William will attempt \(N \sim \text{Geom}(s)\) Rubik’s Cubes. Suppose each time, he has probability \(p\) of solving the cube, independently. Let \(T\) be the number of Rubik’s Cubes he solves during a day. Find the mean and variance of \(T\).

Solution

Note that \(T|N \sim \text{Bin}(N,p)\). So by Adam’s Law,
\(E(T) = E(E(T|N)) = E(Np) = \boxed{\frac{p (1-s)}{s}}\)
Similarly, by Eve’s Law, we have
\(\text{Var}(T) = E(\text{Var}(T|N)) + \text{Var}(E(T|N)) = E(Np(1-p)) + \text{Var}(Np) = \frac{p(1-p)(1-s)}{s} + \frac{p^2(1-s)}{s^2} = \boxed{\frac{p(1-s)(p+s(1-p))}{s^2}}\)

MGF – Distribution Matching

(Continuing the Rubik’s Cube question above) Find the MGF of \(T\). What is the name of this distribution and its parameter(s)?

Solution

By Adam’s Law, we have
\(E(e^{tT}) = E(E(e^{tT}|N)) = E((pe^t + q)^N) = s\sum_{n=0}^\infty(pe^t + 1-p)^n(1-s)^n =\frac{s}{1-(1-s)(pe^t+1-p)} =\frac{s}{s+(1-s)p-(1-s)pe^t}\)
Intuitively, we would expect that \(T\) is distributed Geometrically since \(T\) is just a filtered version of \(N\), which itself is Geometrically distributed. The MGF of \(X\sim\text{Geom}(\theta)\) is \(E(e^{tX}) = \frac{\theta}{1-(1-\theta) e^t}\)
So, we would want to try to get our MGF into this form to identify what \(\theta\) is. Taking our original MGF, it would appear that dividing by \(s+(1-s)p\) would allow us to do this. Therefore, we have that
\(E(e^{tT}) = \frac{s}{s+(1-s)p - (1-s)pe^t} = \frac{\frac{s}{s+(1-s)p}}{1-\frac{(1-s)p}{s+(1-s)p}e^t}\)
By pattern-matching, it thus follows that \(\boxed{T \sim \text{Geom}(\theta)}\) where
\(\boxed{\theta = \frac{s}{s+(1-s)p}}\)

MGF – Finding Moments

Find \(E(X^3)\) for \(X \sim \text{Expo}(\lambda)\) using the MGF of \(X\).

Solution

The MGF of an \(\text{Expo}(\lambda)\) is \(M(t) = \frac{\lambda}{\lambda-t}\). To get the third moment, we can take the third derivative of the MGF and evaluate at \(t=0\):
\(\boxed{E(X^3) = \frac{6}{\lambda^3}}\)
But a much nicer way to use the MGF here is via pattern recognition: note that \(M(t)\) looks like it came from a geometric series:
\(\frac{1}{1-\frac{t}{\lambda}} = \sum^{\infty}_{n=0} \left(\frac{t}{\lambda}\right)^n = \sum^{\infty}_{n=0} \frac{n!}{\lambda^n} \frac{t^n}{n!}\)
The coefficient of \(\frac{t^n}{n!}\) here is the \(n\)th moment of \(X\), so we have \(E(X^n) = \frac{n!}{\lambda^n}\) for all nonnegative integers \(n\).

Markov chains (1)

Suppose \(X_n\) is a two-state Markov chain with transition matrix\ \(Q = \begin{bmatrix} 1-\alpha & \alpha \\ \beta & 1-\beta \end{bmatrix}\)
Find the stationary distribution \(\vec{s} = (s_0, s_1)\) of \(X_n\) by solving \(\vec{s} Q = \vec{s}\), and show that the chain is reversible with respect to \(\vec{s}\).\

Solution

The equation \(\vec{s}Q = \vec{s}\) says that
\(s_0 = s_0(1-\alpha) + s_1 \beta \text{ and } s_1 = s_0(\alpha) + s_0(1-\beta)\)
By solving this system of linear equations, we have
\(\boxed{\vec{s} = \left(\frac{\beta}{\alpha+\beta}, \frac{\alpha}{\alpha+\beta}\right)}\)
To show that the chain is reversible with respect to \(\vec{s}\), we must show \(s_i q_{ij} = s_j q_{ji}\) for all \(i, j\). This is done if we can show \(s_0 q_{01} = s_1 q_{10}\). And indeed,
\(s_0 q_{01} = \frac{\alpha\beta}{\alpha+\beta} = s_1 q_{10}\)

Markov chains (2)

William and Sebastian play a modified game of Settlers of Catan, where every turn they randomly move the robber (which starts on the center tile) to one of the adjacent hexagons.

1. Is this Markov chain irreducible? Is it aperiodic?

Solution

Yes to both. The Markov chain is irreducible because it can get from anywhere to anywhere else. The Markov chain is aperiodic because the robber can return back to a square in 2, 3, 4, 5, … moves, and the greatest common divisor (GCD) of those numbers is 1.

2. What is the stationary distribution of this Markov chain?

Solution

Since this is a random walk on an undirected graph, the stationary distribution is proportional to the degree sequence. The degree for the corner pieces is 3, the degree for the edge pieces is 4, and the degree for the center pieces is 6. To normalize this degree sequence, we divide by its sum. The sum of the degrees is 6(3) + 6(4) + 7(6) = 84. Thus, the stationary probability of being on a corner is 3/84 = 1/28, on an edge is 4/84 = 1/21, and in the center is 6/84 = 1/14.

3. What fraction of the time will the robber be in the center tile in this game, in the long run?

Solution

By the above, 1/14.

4. What is the expected amount of moves it will take for the robber to return to the center tile?

Solution

Since this chain is irreducible and aperiodic, to get the expected time to return we can just invert the stationary probability. Thus, on average it will take 14 turns for the robber to return to the center tile.

Problem-Solving Strategies

Contributions from Jessy Hwang, Yuan Jiang, Yuqi Hou

1. Getting started. Start by defining relevant events and random variables. (“Let \(A\) be the event that I pick the fair coin”; “Let \(X\) be the number of successes.”) Clear notion is important for clear thinking! Then decide what it is that you’re supposed to be finding, in terms of your notation (“I want to find \(P(X=3|A)\)”). Think about what type of object your answer should be (a number? A random variable? A PMF? A PDF?) and what it should be in terms of.
   Try simple and extreme cases. To make an abstract experiment more concrete, try drawing a picture or making up numbers that could have happened. Pattern recognition: does the structure of the problem resemble something we’ve seen before?

2. Calculating probability of an event. Use counting principles if the naive definition of probability applies. Is the probability of the complement easier to find? Look for symmetries. Look for something to condition on, then apply Bayes’ Rule or the Law of Total Probability.

3. Finding the distribution of a random variable. First make sure you need the full distribution not just the mean (see next item). Check the support of the random variable: what values can it take on? Use this to rule out distributions that don’t fit. Is there a story for one of the named distributions that fits the problem at hand? Can you write the random variable as a function of an r.v. with a known distribution, say \(Y = g(X)\)?

4. Calculating expectation. If it has a named distribution, check out the table of distributions. If it’s a function of an r.v. with a named distribution, try LOTUS. If it’s a count of something, try breaking it up into indicator r.v.s. If you can condition on something natural, consider using Adam’s law.

5. Calculating variance. Consider independence, named distributions, and LOTUS. If it’s a count of something, break it up into a sum of indicator r.v.s. If it’s a sum, use properties of covariance. If you can condition on something natural, consider using Eve’s Law.

6. Calculating \(E(X^2)\). Do you already know \(E(X)\) or \(\text{Var}(X)\)? Recall that \(\text{Var}(X) = E(X^2) - (E(X))^2\). Otherwise try LOTUS.

7. Calculating covariance. Use the properties of covariance. If you’re trying to find the covariance between two components of a Multinomial distribution, \(X_i, X_j\), then the covariance is \(-np_ip_j\) for \(i \neq j\).

8. Symmetry. If \(X_1,\dots,X_n\) are i.i.d., consider using symmetry.

9. Calculating probabilities of orderings. Remember that all \(n!\) ordering of i.i.d. continuous random variables \(X_1,\dots,X_n\) are equally likely.

10. Determining independence. There are several equivalent definitions. Think about simple and extreme cases to see if you can find a counterexample.

11. Do a painful integral. If your integral looks painful, see if you can write your integral in terms of a known PDF (like Gamma or Beta), and use the fact that PDFs integrate to \(1\)?

12. Before moving on. Check some simple and extreme cases, check whether the answer seems plausible, check for biohazards.

Biohazards

Contributions from Jessy Hwang

1. Don’t misuse the naive definition of probability. When answering “What is the probability that in a group of 3 people, no two have the same birth month?”, it is not correct to treat the people as indistinguishable balls being placed into 12 boxes, since that assumes the list of birth months {January, January, January} is just as likely as the list {January, April, June}, even though the latter is six times more likely.

2. Don’t confuse unconditional, conditional, and joint probabilities. In applying \(P(A|B) = \frac{P(B|A)P(A)}{P(B)}\), it is not correct to say “\(P(B) = 1\) because we know \(B\) happened”; \(P(B)\) is the prior probability of \(B\). Don’t confuse \(P(A|B)\) with \(P(A,B)\).

3. Don’t assume independence without justification. In the matching problem, the probability that card 1 is a match and card 2 is a match is not \(1/n^2\). Binomial and Hypergeometric are often confused; the trials are independent in the Binomial story and dependent in the Hypergeometric story.

4. Don’t forget to do sanity checks. Probabilities must be between \(0\) and \(1\). Variances must be \(\geq 0\). Supports must make sense. PMFs must sum to \(1\). PDFs must integrate to \(1\).

5. Don’t confuse random variables, numbers, and events. Let \(X\) be an r.v. Then \(g(X)\) is an r.v. for any function \(g\). In particular, \(X^2\), \(|X|\), \(F(X)\), and \(I_{X>3}\) are r.v.s. \(P(X^2 < X | X \geq 0)\), \(E(X)\), \(\text{Var}(X)\), and \(g(E(X))\) are numbers. \(X = 2\) and \(F(X) \geq -1\) are events. It does not make sense to write \(\int_{-\infty}^\infty F(X) dx\), because \(F(X)\) is a random variable. It does not make sense to write \(P(X)\), because \(X\) is not an event.

6. Don’t confuse a random variable with its distribution. To get the PDF of \(X^2\), you can’t just square the PDF of \(X\). The right way is to use transformations. To get the PDF of \(X + Y\), you can’t just add the PDF of \(X\) and the PDF of \(Y\). The right way is to compute the convolution.

7. Don’t pull non-linear functions out of expectations. \(E(g(X))\) does not equal \(g(E(X))\) in general. The St. Petersburg paradox is an extreme example. See also Jensen’s inequality. The right way to find \(E(g(X))\) is with LOTUS.

Distributions in R

The table above gives R commands for working with various named distributions. Commands analogous to pbinom, qbinom, and rbinom work for the other distributions in the table. For example, pnorm, qnorm, and rnorm can be used to get the CDF, quantiles, and random generation for the Normal. For the Multinomial, dmultinom can be used for calculating the joint PMF and rmultinom can be used for generating random vectors. For the Multivariate Normal, after installing and loading the mvtnorm package, dmvnorm can be used for calculating the joint PDF and rmvnorm can be used for generating random vectors.

Recommended Resources

• Introduction to Probability Book
• Stat 110 Online
• Stat 110 Quora Blog
• Quora Probability FAQ
• R Studio
• LaTeX File: github.com/wzchen/probability_cheatsheet

Please share this cheatsheet with friends! wzchen.com/probability-cheatsheet

Table of Distributions

Key Points

Permutations and Combinations

Overview

Teaching: min
Exercises: min

Questions

• what are permutations and combinations

Objectives

• Apply Fundamental Counting Principle

• Apply Permutations

• Apply Combinations

Concepts / Definitions

The Fundamental Counting Principle states that if one event has \(m\) possible outcomes and a second independent event has \(n\) possible outcomes, then there are \(m\) x \(n\) total possible outcomes for the two events together.

If you have four flavors of ice cream and two types of cones, then ther are \(4 * 2 = 8\) possible combinations.

In mathematics, the factorial of a non-negative integer \(n\), denoted by \(n!\), is the product of all positive integers less than or equal to \(n\). \(n! = n(n-1)(n-2)(n-3)\ ...\ (2)(1)\) By definition, \(0! = 1\).

Permutations are the number of ways a set of \(n\) distinguishable objects can be arranged in order.
\(4!\) = 24 ways to order four items

The number of permutations on \(n\) objects taken \(r\) at a time is given by \(P\binom{n}{r} = P(n, r) = nPr = \frac{n!}{(n-r)!}\)

The number of ways \(n\) items can be ordered with replacement \(r\) times is \(n^r\)

\(\frac{4!}{(4-3)!}\) = 24 ways of selecting and ordering 3 or 4 letters, but only 4 ways if order does not matter.

Combinations are the number of ways selecting \(r\) items from a group of \(n\) items where order does not matter.
To take out all the ways \(r\) can happen, we divide out all the ways \(r!\) can happen.
The number of combinitions of \(n\) objects taken \(r\) at a time is given by \(C\binom{n}{r} = C(n, r) = nCr = \frac{n!}{r!(n-r)!}\) \(\implies\) Also called \(n\) choose \(k\), noted \(\binom{n}{k} = \frac{n!}{k!(n-k)!}\)

Counting Subsets of an \(n\)-Set
Consider a binomial situation, where there is a yes or no, success or failure, possibility happening \(n\) times. The number of ways this can happen is \(2^n\). There are \(2^n\) subsets of a set with \(n\) objects.

Exercises

1. A four-volume work is placed in random order on a bookshelf. What is the probability of the volumes being in proper order from left to right or from right to left?

Solution

If you have four volumes, and the question is what order to place them in, the question is a simple permutation problem. There are 4 possibilities for the first, 3 for the second, 2 for the third and only one for the last, making the number of permutations of the books 24. Only one order has them in ascending order and only one order will have them in descending order. Thus: \(P=\frac{2}{24}=\frac{1}{12}\label{answer1.1}\)

2.A wooden cube with painted faces is sawed up into 1000 little cubes, all of the same size. The little cubes are then mixed up, and one is chosen at random. What is the probability of its having just 2 painted faces?

Solution

The wooden cube is made of \(10^3\) cubes implying a 10x10x10 cube. The cubes that have two faces painted will be the edges which are not on a corner. Thus, since there are 12 edges of 10 cubes each, 8 of which are not corners, that implies we have 96 edges. 96 thus becomes our number of desireable outcomes and 1000 has always been the total number of outcomes: \(P=\frac{96}{1000} = 0.096\label{answer1.2}\)

3. A batch of n manufactured items contains k defective items. Suppose m items are selected at random from the batch. What is the probability that l of these items are defective?

Solution

If there are \(n\) total items and \(k\) of them are defective. We select \(m\) and want to know the probability of \(l\) of them being defective ones. There are \(\binom{n}{m}\) possible ways to chose \(m\) different items from the population of \(n\) items which will be our denominator. Now we need to know how many of those possibilities have \(l\) bad ones in them for our numerator. If there’s \(k\) total defective ones, then there are \(\binom{k}{l}\) \(P = \frac{\binom{k}{l}}{\binom{n}{m}}\label{answer1.3}\)

4. Ten books are placed in random order on a bookshelf. Find the probability of three given books being side by side.

Solution

There are \(10!\) possible ways to order the ten books. You can imagine that since the three books have to take up three adjacent positions, there are 8 possible locations for the three books to be ordered. Taking just the first position (that being the first three “Slots”), there are 6 ways to order the desired books and then \(7!\) ways to order the remaining books. Thus there are \(8 \cdot{} 6 \cdot{} 7!\) desirable orders giving us a probability of \(P=\frac{8 \cdot{} 6 \cdot{} 7!}{10!}=\frac{1}{15}\label{answer1.4}\)

5. One marksman has an 80% probability of hitting a target, while another has only a 70% probability of hitting the target. What is the probability of the target being hit (at least once) if both marksman fire at it simultaneously?

Solution

There are 4 possibilities to consider here: neither hits, one or the other hits and both hits. While it may be tempting to calculate the probability of each event, since we only care about the probability of at least one hitting the target, we need only calculate the probability that no one hits and subtract that from 1. The first marksman has an \(.8\) probability of hitting meaning he has a missing probability of .2; similarly the second marksman has a .7 chance of hitting and a .3 chance of missing. Thus, the probability of both one and two missing is the product of the two missing probabilities: \(.2*.3=.06\) \(P=1-.06=.94\label{answer1.5}\)

6. Suppose n people sit down at random and independently of each other in an auditorium containing n + k seats. What is the probability that m seats specified in advance (m < n) will be occupied?

Solution

The total number of ways \(n+k\) seats can be occupied by \(n\) people is \(\binom{n+k}{n}=\frac{(n+k)!}{n!k!}\) but once again the difficult part is finding the total number of desirable outcomes. If you were to specify \(m\) seats, you effectively divide the auditorium into two buckets so we need to find the number of ways to partition people into those two buckets which will give us the numerator, \(\binom{n+k}{m}=\frac{(n+k)!}{m!(n+k-m)!}\)
\(P=\frac{(n+k)!}{m!(n+k-m)!}\frac{n!k!}{(n+k)!}\)
\(=\frac{n!k!}{m!(n+k-m)!}\)

7. Three cards are drawn at random from a full deck. What is the probability of getting a three, a seven and an ace?

Solution

Again, the number of ways to get three cards from a 52 card deck is \(\binom{52}{3}\). Since there are 4 each of sevens, threes and aces, there are $4^3$ desirable hands. \(P=\frac{4^3}{\frac{52!}{3!49!}}=\frac{16}{5525}=0.00289593 \label{answer1.7}\) Which it is worth pointing out, is no different than any other 3-card hand of three different cards.

8. What is the probability of being able to form a triangle from three segments chosen at random from five line segments of lengths 1, 3, 5, 7 and 9?

Hint. A triangle cannot be formed if one segment is longer than the sum of the other two.

Solution

If you indiscriminately chose 3 line segments from our bank of 5, you have \(\binom{5}{3}=10\) total possibilities for triangles.
When you look at the bank \((1, 3, 5, 7, 9)\), however, you have to make sure that no two segments are longer than the third segment, otherwise making a triangle is impossible.
The brute-force way to do this is start with 1 and realize there’s no triangle that can be formed with 1.
Then, looking at 3, you realize that you can do \((3, 5, 7)\) and \((3, 7, 9).\ Starting now with 5, you can do only\)(5, 7, 9)\(\ At this point, you realize you're done and that the answer is plain to see\)P=\frac{3}{10}=.3\label{answer1.8}$$

9. Suppose a number from 1 to 1000 is selected at random. What is the proba­bility that the last two digits of its cube are both 1?

Hint There is no need to look through a table of cubes.

Solution

We could do this problem in 5 minutes of programming and an instant of computation but that’s not the point! We need to think our way through this one. How many cubes of the numbers between 1 and 1000, have 11 for the last two digits. Luckily, each cube is unique so there’s no complications there: only 1000 possibilities. Let’s break down the number we’re cubing into two parts, one that encapsulates the part less than 100 and then the rest of the number \(n=a+b=100c+b\) Now, just for fun, let’s cube than number \(n^3=(100c+b)^3=b^3+300 b^2 c+30000 b c^2+1000000 c^3\) Clearly the only term here that will matter to the last two digits of the cube is \(b\), the part less than 100. Now we can reduce our now size 100 subspace by a great deal if you realize that in order for a number’s cube to end in 1, the last number in the cube will need to be 1, leaving us: (1, 11, 21, 31, 41, 51, 61, 71, 81, 91). At this point I recommend just cubing all the numbers and realizing that only \(71^3=357911\) fulfills the requirements giving only one desirable outcome per century(71, 171, 271 etc.) \(P=\frac{10}{1000}=.01\label{answer1.9}\)

10. Find the probability that a randomly selected positive integer will give a number ending in 1 if it is

a) Squared ; b) Raised to the fourth power; c) Multiplied by an arbitrary positive integer. Hint. It is enough to consider one-digit numbers.

Solution

Now we want to know about all positive integers but just the last digit and the probability that the last digit is one. To prove that we can use the hint given, we’ll split the integer into the past less than 10 and the part greater than 10. \(n=10a+b\)

a–Squared

\(n^2=100a^2+10ab+b^2\) OK, so clearly only \(b^2\) contributes to the last digit. At this point... just do the squaring especially since it’s something you can do in your head. (1, 2, 3, 4, 5, 6, 7, 8, 9)\(\rightarrow\)(1, 4, 9, 16, 25, 36, 49, 64, 81) therefore there are two desireable outcomes for every decade of random positive integers.
\(P=\frac{2}{10}=.2\)

b–fourth-powered

\(n^4=100000000 a^4+4000000 a^3 b+60000 a^2 b^2+400 a b^3+b^4\) Great, once again, just \(b^4\). Now, doing the forth power is a little harder so let’s reason through and reduce the subspace. Clearly, only odd numbers will contribute since any power of an even number is again an even number. Now, let’s just do the arithmetic. (1, 3, 5, 7, 9)\(\rightarrow\)(1, 81, 625, 2401, 6561) giving us 4 desirable outcomes per decade of random numbers
\(P=\frac{4}{10}=.4\)

c–multiplied by random positive number

\(n*r=(10a+b)r\) Now let’s let $r$ be a similar integer as n \(n*r=(10a+b)(10c+d)=100ac+10(bc+ad)+bd\) So we just need to consider the first two digits of both the random number and the arbitrary number. This leaves us with $10^2$ possibilities; \(5^2\) possible desirable outcomes once we exclude even numbers: (1, 3, 5, 7, 9). Once again we resort to brute force by just defining the above as a vector and multiplying the transpose of that vector by the vector. The only desirable outcomes turn out to be (1*1, 9*9, 7*3, 3*7), giving four desirable outcomes per century.
\(P=\frac{4}{100}=.04\)

11. One of the numbers 2, 4, 6, 7, 8, 11, 12 and 13 is chosen at random as the numerator of a fraction, and then one of the remaining numbers is chosen at random as the denominator of the fraction. What is the probability of the fraction being in lowest terms?

Solution

Here we are given 8 possible numbers (2, 4, 6, 7, 8, 11, 12, 13) and are told that it making a random fraction out of two of the numbers. There are \(\binom{8}{2}=28\) possible pairs but since \(\frac{a}{b} \neq \frac{b}{a}\), we multiply that by two to get the total number of fractions we can get from these numbers as being 56. Looking at the numbers given, 7, 11 and 13 are all prime and all the others are divisible by two; therefore only the fractions that contain one or more prime numbers will be in lowest terms. The number 7 can be in 7 different pairs with the other numbers which makes for 14 fractions. 11 can be in 6 pairs or 12 fractions (since we already counted 7) and similarly, 13 can be in 10 uncounted fractions giving 36 possible fractions in lowest terms.
\(P=\frac{36}{56}=\frac{9}{14}\label{answer1.11}\)

12. The word “drawer” is spelled with six scrabble tiles. The tiles are then randomly rearranged. What is the probability of the rearranged tiles spelling the word “reward?»

Solution

The word drawer has 6 letters and there are \(6!=720\) possible ways of arranging the letters. It is then natural to say there’s only one way to spell reward correctly and thus the probability of spelling it correctly after a random reordering is \(1/720\) BUT there is a complication in that the two of our letters are the same, meaning there are to distinct arrangement of our distinguishable tiles that will give us the proper spelling. \(P=\frac{2}{720} = \frac{1}{360}\label{answer1.12}\)

13. In throwing 6n dice, what is the probability of getting each face n times? Use Stirling’s formula to estimate this probability for large n.

Solution

For any die, there are 6 different possibilities. Since one dice’s outcome does not depend on another’s, that means that for a roll of \(6n\) die, there are \(6^{6n}\) different possible outcomes for the dice rolls. Now for desirable outcomes, we want each of the 6 faces to show up n times. To accomplish this, we just count the number of ways to apportion \(6n\) things into 6 groups of $n$ each or \(\frac{(6n)!}{(n!)^6}\) which, given Stirling’s approximation \(n! \approx \sqrt{2 \pi n} n^n e^{-n}\) gives us for large \(n\) \(\frac{(6n)!}{(n!)^6} \approx \sqrt{2 \pi 6n} (6n)^{6n} e^{-6n} * \frac{1}{(\sqrt{2 \pi n} n^n e^{-n})^6} = \frac{3 \cdot 6^{6n}}{4 (\pi n)^{5/2}}\) \(P=\frac{(6n)!}{(n!)^6 6^{6n}} \approx \frac{3}{4 (\pi n)^{5/2}}\label{answer1.13}\)

14. A full deck of cards is divided in half at random. Use Stirling’s formula to estimate the probability that each half contains the same number of red and black cards.

Solution

To figure out the total number of possibilities, we must realize that one draw of half a deck, implies the other half of the deck implicitly, therefore there are \(\binom{52}{26}\) total 26 card draws. Then, to get 13 red cards, there are \(\binom{26}{13} = 10400600\) ways to get 13 red cards and the same number for black cards making $10400600^2$ the number of total possible desirable outcomes.
\(P=\frac{16232365000}{74417546961}=0.218126\)
Because we’re cool and modern and have Mathematica, we don’t NEED to do the Stirling’s formula approximation but it’ll be good for us so we shall.
\(P=\frac{\binom{26}{13}^2}{\binom{52}{26}} = \frac{(26!)^4}{(13!)^4 52!}=\frac{((2n)!)^4}{(n!)^4 (4n)!}\) where \(n=13\).
\(P=\frac{(\sqrt{4 \pi n}(2n)^{2n}e^{-2n})^4}{(\sqrt{2 \pi n}n^{n}e^{-n})^4 (\sqrt{8 \pi n}(4n)^{4n}e^{-4n})}\)
\(P=\frac{(4 \pi n)^2}{(2 \pi n)^2 \sqrt{8 \pi n}} \frac{(2n)^{8n}}{n^{4n}(4n)^{4n}} = \frac{4}{\sqrt{8 \pi n}} \frac{2^{8n}}{4^{4n}}\)
\(P=\frac{2}{\sqrt{26 \pi}} = 0.221293\) which when you take the ratio of approximate to exact, you get $1.01452$ so the approximation is less than \(2\%\) off... not bad!

15. Use Stirling’s formula to estimate the probability that all 50 states are represented in a committee of 50 senators chosen at random.

Solution

There are 100 senators at any given time. Much like the previous problem, there are \(\binom{100}{50}=100891344545564193334812497256\) different 50 senator committees. Much like the previous problem, there are \(\binom{2}{1}=2\) ways for California to be represented on the committee making for 50 states \(2^{50}\) possible even committees \(P=\frac{2^{50}}{\binom{100}{50}} = \text{1.115952921347132$\grave{ }$*${}^{\wedge}$-14}\label{answer1.15}\)

16. Suppose 2n customers stand in line at a box office, n with 5-dollar bills and n with IO-dollar bills. Suppose each ticket costs 5 dollars, and the box office has no money initially. What is the probability that none of the customers has to wait for change?

Solution

For \(2n\) people, there are \((2n)!\) different possible lines; we want the number of lines where at any given point in the line, there are more people with 5 dollar bills than 10 dollar bills. The given reference (freely available on Google books) has a fascinating geometrical argument about how \(\binom{2n}{n+1}\) is the number of lines that have one or more too many people with tens in front of people with fives. Also, in this argument, there are \(\binom{2n}{n}\) trajectories instead of \((2n)!\) lines. \(P=\frac{\binom{2n}{n}-\binom{2n}{n+1}}{\binom{2n}{n}} = \frac{1}{n+1}\label{answer1.16}\)

17. Prove that \(\sum_{k=0}^{n} {n \choose k} ^{2} ={n \choose 2} ^{n}\)

Hint. Use the binomial theorem to calculate the coefficient of \(x^{n}\) in the product \((1 + x)^{n} \cdot (1 + x)^{n} = (1 + x)^{2n}\).

Solution

We wish to prove that \(\sum_{k=0}^{n} \binom{n}{k}^2 = \binom{2n}{n}\label{answer1.17}\)
This is easies if we listen to the hint and consider that
\(\binom{2n}{n}\) is the coefficient of \(x^n\) in the polynomial
\((x+1)^{2n}\) which is also equivalent to \((x+1)^n(x+1)^n\)
We make use of: \((x+1)^n = \sum_{k=0}^{n} \binom{n}{k} x^k\)
\((x+1)^n(x+1)^n = \sum_{k=0}^{n}\sum_{j=0}^{n} \binom{n}{k} \binom{n}{j} x^k x^j\)
We want the \(x^n\) term: where \(k+j=n\) or, put another way, where
\(k=n-j\)... \(\sum_{j=0}^{n} \binom{n}{n-j} \binom{n}{j} = \sum_{j=0}^{n} \frac{n!}{j!(n-j)!}\frac{n!}{(n-j)!j!} = \sum_{j=0}^{n} \binom{n}{j}^2\)

Exercises

1. How many ways can 8 runners in a track race finish?
2. Subway is running a special on 6” subs for $3.99 with 6 choices of bread, 3 choices of meat, and 2 choices of cheese. Not including toppings or sause, how many choices are available?
3. How many ways can 2 students of a class of 25 be selected for room reps?
4. How many distinguishable ways can the letters, TOPPINGS be written?
5. A license plate has 3 letters followed by 3 digits. How many possible plates are there?
6. From a group fo 10 people in theater, how many possible ways can a teacher select roles for a show on the 7 dwarfs?
7. In the original version of poker, straight poker, a five-card hand is dealt from a standard deck of 52. How many different hands are possible?
8. A coin is flipped 10 times, and the sequence is recorded.
   1. How many sequences are possible?
   2. How many sequences have exactly 7 heads?
9. A particular subway advertises 256 ways for a sandwich to be fixed with veggies. How many veggie toppings does this subway offer?
10. Calculate \(_{10}C_3\) without a calculator.

Key Points

Combinatorics CheatSheet

Overview

Teaching: min
Exercises: min

Questions

Objectives

Basics of generating functions

• Introduction [Wilf 1–3]:

  • how to define a sequence: exact formula, recurrent relation (Fibonacci), algorithm (the sequence of primes); there are uncomputable sequences (programs that do not stop)

  • a new way: power series (members of the sequence as coefficients in the series)

  • advantages: many advanced tools from analytical theory of functions

  • very powerful: works on many sequences where nothing else is known to work

  • allows to get asymptotic formulas and statistical properties

  • powerful way to prove combinatorial identities

  • “Konečne vidím, že je tá matalýza aj na niečo dobrá. Keby mi to bol niekto predtým povedal…”

• Two examples [Wilf 3–7]:

  • \(a_{n+1} = 2a_n + 1\) for \(n\ge 0\), \(a_0 = 0\)

  • write few members, guess \(a_n = 2^n-1\), provable by induction

  • multiply by \(x^n\), sum over all \(n\), assign gf: \(\displaystyle \qquad{A(x)\over x}=2A(x)+{1\over 1-x}\)

  • partial fraction expansion: \(\displaystyle \qquad A(x)={x\over (1-x)(1-2x)}={1\over 1-2x}-{1\over 1-x}\)

  • the method stays basically the same for harder problems

  • \(a_{n+1}=2a_n+n\) for \(n\ge 0\), \(a_0=1\)

  • exact formula not obvious; no unqualified variables in the recurrence

  • obstacle: \(\sum_{n\ge 0} nx^n = x/(1-x)^2\); solution: differentiation

  • concern: is differentiation allowed? discussed later, but in principle yes: in formal power series (as an algebraic ring) or via convergence (if we care about analytical properties)

  • \[\displaystyle A(x) = {1-2x+2x^2\over (1-x)^2(1-2x)} = {A\over (1-x)^2} + {B\over 1-x} + {C\over 1-2x} = {-1\over (1-x)^2} + {2\over 1-2x}\]

  • \(1/(1-x)^2\) is just \(x/(1-x)^2\) (see above) shifted by \(1\)

  • \[a_n=2^{n+1}-n-1\]

• The method [Wilf 8]:

  • 1. make sure variables in the recurrence are qualified (e.g. range for \(n\))

  • 2. name and define the gf

  • 3. multiply by \(x^n\), sum over all \(n\) in the range

  • 4. express both sides in terms of the gf

  • 5. solve the equation for gf

  • 6. calculate coefficients of gf power series expansion

  • useful notation: \([x^n]f(x)\); e.g. \([x^n]e^x=1/n!\qquad [t^r]{1\over 1-3t}=3^r\qquad [v^m](1+v)^s={s\choose m}\)

• Solve \(a_n=5a_{n-1}-6a_{n-2}\) for \(n\ge 2\), \(a_0 = 0\), \(a_1=1\).

• Fibonacci [Wilf 8–10]:

  • three-term recurrence: \(F_{n+1}=F_n+F_{n-1}\) for \(n\ge 1\), \(F_0=0\), \(F_1=1\).

  • apply the method (\(r_\pm = (1\pm \sqrt 5)/2\)): \(\displaystyle F(x) = {x\over 1-x-x^2} = {x\over (1-xr_{+})(1-xr_{-})}={1\over r_{+}-r_{-}}\left({1\over 1-xr_{+}}-{1\over 1-xr_{-}}\right)\)

  • \[F_n={1\over \sqrt 5}(r_{+}^n-r_{-}^n)\]
  • the second term is \({} < 1\) and goes to zero, so the first term \({1\over \sqrt 5}({1+\sqrt 5\over 2})^n\) gives a good approximation

• Find ogf for the following sequences (always \(n\ge 0\)) [W1.1]:

  ---

  (a) \(a_n = n\)
  (b) \(a_n = \alpha n + \beta\)
  (c) \(a_n = n^2\)
  (d) \(a_n = n^3\)
  (e) \(a_n = P(n)\); \(P\) is a polynomial of degree \(m\)
  (f) \(a_n = 3^n\)
  (g) \(a_n = 5\cdot 7^n-3\cdot 4^n\)
  (h) \(a_n = (-1)^n\)
  ——- ————————————————- –

• Find the following coefficients [W1.5]:

  ---

  (a) \([x^n]\, e^{2x}\)
  (b) \([x^n/n!]\, e^{\alpha x}\)
  (c) \([x^n/n!]\, \sin x\)
  (d) \([x^n]\, 1/(1-ax)(1-bx)\) (\(a\neq b\))
  (e) \([x^n]\, (1+x^2)^m\)
  ——- ————————————– –

• Compute \(\square_n = \sum_{k=1}^n k^2\).

  • assign ogf to the sequence \(1^2, 2^2, \dots, n^2\): \(f(x) = \sum_{k=1}^n{k^2x^k}\)

  • \[(x{\rm D})^2 [(x^{n+1}-1)/(x-1)] = x {-2 n^2 x^{n + 1} + n^2 x^{n + 2} + n^2 x^n - 2 n x^{n + 1} + x^{n + 1} + 2 n x^n + x^n - x - 1)\over (x - 1)^3}\]
  • note that \(\square_n = f(1) = \lim_{x\to 1} (xD)^2 [(x^{n+1}-1)/(x-1)]=n(n+1)(2n+1)/6\)

• Find the sequence with gf \(1/(1-x)^3\).

• Find a linear recurrence going back two sequence members that has a solution that contains \(n\cdot 3^n\) (possibly plus some linear combination of other exponential or polynomial factors).

• Find explicit formulas for the following sequences [W1.6, R2, R3, R7]:

  ---

  (a) \(a_{n+1} = 3a_n+2\) for \(n\ge 0\); \(a_0=0\)
  (b) \(a_{n+1} = \alpha a_n + \beta\) for \(n\ge 0\); \(a_0=0\)
  (c) \(a_{n+1} = a_n/3 +1\) for \(n\ge 0\); \(a_0=1\)
  (d) \(a_{n+2} = 2a_{n+1}-a_n\) for \(n\ge 0\), \(a_0=0\), \(a_1=1\)
  (e) \(a_{n+2} = 3a_{n+1}-2a_n+3\) for \(n\ge 0\); \(a_0=1\), \(a_1=2\)
  (f) \(a_n = 2a_{n-1}-a_{n-2}+(-1)^n\) for \(n>1\); \(a_0=a_1=1\)
  (g) \(a_n = 2a_{n-1}-n\cdot(-1)^n\) for \(n\ge 1\); \(a_0=0\)
  (h) \(a_n = 3a_{n-1} + {n\choose 2}\) for \(n\ge 1\); \(a_0=2\)
  (i) \(a_n = 2a_{n-1}-a_{n-2}-2\) for \(n\ge 2\); \(a_0=a_{10}=0\)
  (j) \(a_n = 4(a_{n-1}-a_{n-2})+(-1)^n\) for \(n\ge 2\); \(a_0=1\), \(a_1=4\)
  (k) \(a_n = -3a_{n-1}+a_{n-2}+3a_{n-3}\) for \(n\ge 3\); \(a_0=20\), \(a_1=-36\), \(a_2=60\)
  ——- ——————————————————————————– –

Ordinary generating functions

• From the homework: solve \(a_n = 2a_{n-1}-a_{n-2}-2\) for \(n\ge 1\); \(a_0=a_{10}=0\).
  Applying the standard method, while keeping \(a_1\) as a parameter, we get \(A(x)={a_1x-a_1x^2-2x^2\over (1-x)^3}={a_1x\over (1-x)^2}+{x(1-x)\over (1-x)^3}-{x^2+x\over (1-x)^3},\) so \(a_n=(a_1+1)n-n^2\). From \(a_{10}=0\) we get \(a_1=9\), thus \(a_n=n(10-n)\).

• Another way for boundary problems (this particular example is motivated by splines, Wilf 10–11):

  • consider \(au_{n+1}+bu_n+cu_{n-1}=d_n\) for \(1\le n\le N-1\); \(u_0=u_N=0\).

  • similar to Fibonacci with two given non-consecutive terms (but more general)

  • define \(U(x)= \sum_{j=0}^N u_jx^j\) (unknown); \(D(x)=\sum_{j=1}^{N-1} d_jx^j\) (known)

  • derive \(\displaystyle a\cdot {U(x)-u_1x\over x}+bU(x)+cx(U(x)-u_{N-1}x^{N-1}) = D(x)\)

  • \((a+bx+cx^2) U(x) = x D(x) +au_1x + cu_{N-1}x^N\) (*)

  • plug in suitable values of \(x\) (roots \(r_{+}\) and \(r_{-}\) of the quadratic polynomial on the LHS)

  • solve the system of two linear equations and two uknowns \(u_1\), \(u_{N-1}\)

  • if the roots are equal, differentiate (*) to obtain the second equation

• Mutually recursive sequences [Knuth 343, Example 3]

  • consider the number \(u_n\) of tilings of \(3\times n\) board with \(2\times 1\) dominoes

  • define \(v_n\) as the number of tilings of \(3\times n\) board without a corner

  • \(u_n = 2v_{n-1} + u_{n-2}\); \(u_0 = 1\); \(u_1 = 0\)

  • \(v_n = v_{n-2} + u_{n-1}\); \(v_0 = 0\); \(v_1 = 1\)

  • derive \(U(x) = {1-x^2\over 1-4x^2+x^4},\qquad V(x) = {x\over 1-4x^2+x^4}\)

  • consider \(W(z) = 1/(1-4z+z^2)\); \(U(x) = (1-x^2)W(x^2)\), so \(u_{2n} = w_n - w_{n-1}\)

  • hence \(u_{2n} = {(2+\sqrt 3)^n\over 3-\sqrt 3} + {(2-\sqrt 3)^n\over 3+\sqrt 3} = \left\lceil {(2+\sqrt 3)^n\over 3-\sqrt 3}\right\rceil\) (derivation as a homework)

• Given \(f(x)\overset{\text{ogf}}{\longleftrightarrow}(a_n)_{n\ge 0}\), express ogf for the following sequences in terms of \(f\) [W1.3]:\

  ---

  (a) \((a_n+c)_{n\ge 0}\)
  (b) \((na_n)_{n\ge 0}\) ; napísať im \((P(n)a_n)_{n\ge 0} \longleftrightarrow P(x{\rm D})f(x)\) (c) \(0, a_1, a_2, a_3, \dots\)
  (d) \(0, 0, 1, a_3, a_4, a_5,\dots\)
  (e) \((a_{n+2}+3a_{n+1}+a_n)_{n\ge 0}\)
  (f) \(a_0, 0, a_2, 0, a_4, 0, a_6, 0\dots\)
  (g) \(a_0, 0, a_1, 0, a_2, 0, a_3, 0\dots\)
  (h) \(a_1, a_2, a_3, a_4,\dots\)
  (i) \(a_0, a_2, a_4, \dots\)
  ——- ————————————— ———————————————————————–

<!-- -->

• introducing a new variable and changing the order of summation can help \(\begin{aligned} \sum_{n\ge 0} {n\choose k}x^n &=& [y^k]\sum_{m\ge 0} \left(\sum_{n\ge 0} {n\choose m}x^n\right)y^m = [y^k]\sum_{n\ge 0} (1+y)^nx^n\nonumber\\ &=& [y^k] {1\over 1-x(1+y)} = {1\over 1-x}[y^k] {1\over 1-{x\over 1-x}y} = {x^k\over (1-x)^{k+1}} \label{binomial} \end{aligned}\)

• alternatively, one can use binomial theorem (Knuth 199/5.56 and 5.57): \(\begin{aligned} {1\over (1-z)^{n+1}} &=& (1-z)^{-n-1} =\sum_{k\ge 0} {-n-1\choose k}(-z)^k\\ &=& \sum_{k\ge 0} {(-n-1)(-n-2)\dots(-n-k)\over k!}(-z)^k = \sum_{k\ge 0} {n+k\choose n}z^k \end{aligned}\)

<!-- -->

• a ring with addition and multiplication \(\sum_n a_nx^n\sum_n b_nx^n = \sum_n \sum_k (a_k b_{n-k})x^n\)

• if \(f(0)\neq 0\), then \(f\) has a unique reciprocal \(1/f\) such that \(f\cdot 1/f = 1\)

• composition \(f(g(x))\) defined iff \(g(0) = 0\) or \(f\) is a polynomial (cf. \(e^{e^x-1}\) vs. \(e^{e^x}\))

• formal derivative \({\rm D}\): \({\rm D}\sum_n a_nx^n = \sum na_nx^{n-1}\); usual rules for sum, product etc.

• exercise: find all \(f\) such that \({\rm D}f = f\)

<!-- -->

• Rule 1: for a positive integer \(h\), \((a_{n+h})\overset{\text{ogf}}{\longleftrightarrow}(f-a_0-\dots-a_{h-1}x^{h-1})/x^h\)

• Rule 2: if \(P\) is a polynomial, then \(P(x{\rm D})f\overset{\text{ogf}}{\longleftrightarrow}(P(n)a_n)_{n\ge 0}\)

  • example: \((n+1)a_{n+1} = 3a_n+1\) for \(n\ge 0\), \(a_0 = 1\); thus \(f' = 3f + 1/(1-x)\)

  • example: \(\sum_{n\ge 0} {n^2+4n+5\over n!}\); thus \(f=\sum_{n\ge 0} (n^2+4n+5){x^n\over n!} = ((x{\rm D})^2+4x{\rm D}+5)e^x = (x^2+5x+5)e^x\)
    we need \(f(1)=11e\); works because the resulting \(f\) is analytic in a disk
    containing \(1\) in the complex plane (that is, it converges to its Taylor series)

• Rule 3: if \(g\overset{\text{ogf}}{\longleftrightarrow}(b_n)\), then \(fg\overset{\text{ogf}}{\longleftrightarrow}(\sum_{k=0}^n a_kb_{n-k})_{n\ge 0}\) \(\sum_{k=0}^n (-1)^kk = (-1)^n\sum_{k=0}^n k\cdot (-1)^{n-k} = (-1)^n[x^n]{x\over (1-x)^2}\cdot{1\over 1+x} = {(-1)^n\over 4}\left(2n+1-(-1)^n\right)\)

• Rule 4: for a positive integer \(k\), we have \(\displaystyle f^k\overset{\text{ogf}}{\longleftrightarrow}\left(\sum_{n_1+n_2+\dots+n_k=n} a_{n_1}a_{n_2}\dots a_{n_k}\right)_{n\ge 0}\)

  • example: let \(p(n,k)\) be the number of ways \(n\) can be written as an ordered sum of \(k\) nonnegative integers

  • according to R4,
  • \((p(n,k))_{n\ge 0}\overset{\text{ogf}}{\longleftrightarrow}1/(1-x)^k\), so \(p(n,k) = {n+k-1\choose n}\) thanks to [binomial]

• Rule 5: \(\displaystyle {f\over (1-x)}\overset{\text{ogf}}{\longleftrightarrow}\left(\sum_{k=0}^n a_k\right)_{n\ge 0}\)\

  • example: \(\displaystyle (\square_n)_{n\ge 0}\overset{\text{ogf}}{\longleftrightarrow}{1\over 1-x}\cdot (x{\mathrm D})^2 {1\over 1-x} = {x(1+x)\over (1-x)^4}\), so by [binomial], \(\square_n = {n+2\choose 3}+{n+1\choose 3}\)

1. Using Rule 5, prove that \(F_0+F_1+\dots+F_n=F_{n+2}-1\) for \(n\ge 0\) [Wilf 38, example 6].\

2. Solve \(g_n=g_{n-1}+g_{n-2}\) for \(n\ge 2\), \(g_0 = 0\), \(g_{10} = 10\).\

3. Solve \(a_n = \sum_{k=0}^{n-1}a_k\) for \(n > 0\); \(a_0 = 1\). [R16]\

4. Solve \(f_n=2f_{n-1}+f_{n-2}+f_{n-3}+\dots+f_1+1\) for \(n\ge 1\), \(f_0 = 0\) [Knuth 349/(7.41)]\

5. Solve \(g_n = g_{n-1} + 2g_{n-2}+\dots +ng_0\) for \(n> 0\), \(g_0 = 1\). [K7.7]\

6. Solve \(g_n = \sum_{k=1}^{n-1} {g_k + g_{n-k} + k\over 2}\) for \(n\ge 2\), \(g_1 = 1\).

7. Solve \(g_n=g_{n-1}+2g_{n-2}+(-1)^n\) for \(n\ge 2\), \(g_0 = g_1 = 1\). [Knuth 341, example 2]\

8. Solve \(a_{n+2}=3a_{n+1}-2a_n+n+1\) for \(n\ge 0\); \(a_0 = a_1 = 1\). [R24]\

9. Prove that \(\displaystyle \ln {1\over 1-x} = \sum_{n\ge 1} {1\over n} x^n\).

Skipping sequence elements, Catalan numbers

• \((1+x)^r = \sum_{k\ge 0} {r\choose k}x^k\); consider \((1+x)^r(1+x)^s = (1+x)^{r+s}\)

• comparison of coefficients yields \(\sum_{k\ge 0}^n {r\choose k}{s\choose n-k}={r+s\choose n}\) — Vandermonde

• by considering \((1-x)^r(1+x)^r = (1-x^2)^r\), we obtain \(\sum_{k=0}^n {r\choose k}{r\choose n-k}(-1)^k = (-1)^{n/2}{r\choose n/2}[2\mid n]\)

<!-- -->

• why \({1\over 2}(A(x)+A(-x))\overset{\text{ogf}}{\longleftrightarrow}a_0, 0, a_2, 0, a_4, \dots\) works: \({1\over 2}(1^n + (-1)^n) = [2\mid n]\)

• in general, for \(\omega\) being \(r\)-th root of unity, \({1\over r}\sum_{j=0}^{r-1} (\omega^j)^n = {1\over r}\sum_{j=0}^{r-1} e^{2\pi ijn/r} = [r\mid n]\)
  — just a geometric progression, or a consequence of \(t^r-1=(t-1)(t^{r-1}+\dots+t+1)\)

• problem: find \(S_n = \sum_k (-1)^k{n\choose 3k}\)

• if we knew \(f(x) = \sum_k {n\choose 3k}x^{3k}\), we would have \(S_n = f(-1)\)

• for \(A(x) = (1+x)^n\), we have \(f(x) = {1\over 3}\big(A(x) + A(x\omega^1) + A(x\omega^2)\big)\) for \(\omega=e^{2\pi i/3}\) and so \(S_n=f(-1) ={1\over 3}[(1-\omega)^n + (1-\omega^2)^n)] =\)

[= {1\over 3}\left[\left({3-\sqrt3 i\over 2}\right)^n+\left({3+\sqrt3 i\over 2}\right)^n\right] = 2\cdot 3^{\frac{n}{2}-1}\cos\left({\pi n\over 6}\right)]

<!-- -->

• consider the number of possibilities \(c_n\) of how to specify the multiplication order of \(A_0A_1\dots A_n\) by parentheses; let \(C(x)=\sum_{n\ge 0} c_nx^n\)

• divide possibilities by the place of last multiplication; \(c_n = \sum\limits_{k=0}^{n-1} c_kc_{n-1-k}\) for \(n > 0\); \(c_0=1\)

• many ways to deal with the recurrence:

  • shift the recurrence to \(c_{n+1} = \sum_{k=0}^n c_kc_{n-k}\) and use Rules 1 and 3; \({C(x)-1\over x} = C(x)^2\)

  • RHS as a convolution of \(c_n\) with \(c_{n-1}\), i.e. \(C(x)\cdot xC(x)\)

  • RHS as a convolution of \(c_n\) with \(c_n\) shifted by Rule 1, i.e. \(x\cdot C(x)^2\)

  • rewriting through sums and changing the order of summation: \(\sum_{n\ge 1}x^n\sum_{k=0}^{n-1}c_kc_{n-1-k}=\sum_{k=0}^\infty x^kc_k\sum_{n\ge k+1} c_{n-1-k}x^{n-k}= \sum_{k=0}^\infty x^kc_k xC(x)=xC(x)\cdot C(x)\)

• consequently, \(C(x) - 1 = xC(x)^2\) and thus \(C(x) = {1\pm \sqrt{1-4x}\over 2x}=\displaystyle {1\over 2x}\left(1 - \sqrt{1-4x}\right)\)

• we want \(C\) continuous and \(C(0) = 1\), so we choose the minus sign (note that the resulting function below is analytical since \({2n\choose n}/(n+1) < 2^{2n}\); it would be analytical also if we chose the plus sign)

• binomial theorem yields \(\begin{aligned} \sqrt{1-4x} = (1-4x)^{1/2} = \sum_{k\ge 0} {1/2\choose k}(-4x)^k &=& 1+\sum_{k\ge 1}{1\over 2k\cdot (-4)^{k-1}}{2k-2\choose k-1}(-4)^kx^k\\ &=& 1 - \sum_{k\ge 1}{2\over k}{2k-2\choose k-1}x^k \end{aligned}\)

• we used \({1/2\choose k}={1/2\over k}{-1/2\choose k-1} = {1\over 2k(-4)^{k-1}}{2k-2\choose k-1}\) because \({-1/2\choose m}={1\over (-4)^m}{2m\choose m}\)

• therefore, \(C(x)={1\over 2x}\sum_{k\ge 1}{2\over k}{2k-2\choose k-1}x^k = \sum_{n\ge 0}{1\over n+1}{2n\choose n}x^n\)

1. Assume that \(A(x)\overset{\text{ogf}}{\longleftrightarrow}(a_n)\). Express the generating function for \(\sum_{n\ge 0} a_{3n}x^n\) in terms of \(A(x)\).\

2. Compute \(S_n=\sum_{n\ge 0} F_{3n}\cdot 10^{-n}\) (by plugging a suitable value into the generating function for \(F_{3n}\)).\

3. Compute \(\sum_k {n\choose 4k}\).

4. Compute \(\sum_k {6m\choose 3k+1}\).

5. Evaluate \(S_n = \sum_{k=0}^n (-1)^k k^2\).

6. Find ogf for \(H_n = 1 + 1/2 + 1/3 + \dots\).

7. Find the number of ways of cutting a convex \(n\)-gon with labelled vertices into triangles.\

Snake Oil

The Snake Oil method [Wilf 118, chapter 4.3] – external method vs. internal manipulations within a sum.

1. identify the free variable and give the name to the sum, e.g. \(f(n)\)

2. let \(F(x) = \sum f(n)x^n\)

3. interchange the order of summation; solve the inner sum in closed form

4. find coefficients of \(F(x)\)

• Example 0

  • let’s evaluate \(f(n) = \sum_k {n\choose k}\); after Step 2, \(F(x) = \sum_{n\ge 0} x^n \sum_k {n\choose k}\)

  • \[\displaystyle F(x) = \sum_k \sum_n {n\choose k}x^n = \sum_{k\ge 0} {x^k\over (1-x)^{k+1}}={1\over 1-x}\cdot {1\over 1-{x\over 1-x}}={1\over 1-2x}\]

• Example 1 [Wilf 121]

  • let’s evaluate \(f(n) = \sum_{k\ge 0} {k\choose n-k}\); after Step 2, \(F(x) = \sum_n x^n \sum_{k\ge 0} {k\choose n-k}\)

  • \[\displaystyle F(x) = \sum_{k\ge 0} \sum_n {k\choose n-k}x^n = \sum_{k\ge 0}x^k\sum_n {k\choose n-k}x^{n-k} = \sum_{k\ge 0}x^k (1+x)^k = {1\over 1-x-x^2}\]
  • so \(f(n) = F_{n+1}\)

• Example 2 [Wilf 122]

  • let’s evaluate \(f(n) = \sum_{k} {n+k\choose m+2k}{2k\choose k}{(-1)^k\over k+1}\), where \(m\), \(n\) are nonnegative integers \(\begin{aligned} F(x) &=& \sum_{n\ge 0} x^n \sum_{k} {n+k\choose m+2k}{2k\choose k}{(-1)^k\over k+1} \\ &=& \sum_k {2k\choose k}{(-1)^k\over k+1}x^{-k}\sum_{n\ge 0}{n+k\choose m+2k}x^{n+k}\\ &=& \sum_k {2k\choose k}{(-1)^k\over k+1}x^{-k}{x^{m+2k}\over (1-x)^{m+2k+1}}\\ &=& {x^m\over (1-x)^{m+1}}\sum_k {2k\choose k}{1\over k+1}\left({-x\over (1-x)^2}\right)^k \\ &=& {-x^{m-1}\over 2(1-x)^{m-1}}\left(1-\sqrt{1+{4x\over (1-x)^2}}\right) = {x^m\over (1-x)^m} \end{aligned}\)

  • so \(f(n) = {n-1\choose m-1}\)

• Example 6 [Wilf 127]

  • prove that \(\sum_{k} {m\choose k}{n+k\choose m} = \sum_k {m\choose k}{n\choose k}2^k\), where \(m\), \(n\) are nonnegative integers

  • the ogf of the left-hand side is \(L(x) = \sum_{k} {m\choose k} x^{-k}\sum_{n\ge 0}{n+k\choose m}x^{n+k} ={(1+x)^m\over (1-x)^{m+1}}\)

  • we get the same for the right-hand side

1. Prove that \(\sum_k k{n\choose k} = n2^{n-1}\) via the snake oil method.

2. Evaluate \(\displaystyle f(n)=\sum_k k^2{n\choose k}3^k\).\

3. Find a closed form for \(\displaystyle \sum_{k\ge 0} {k\choose n-k}t^k\). [W4.11(a)]\

4. Evaluate \(\displaystyle f(n)=\sum_k {n+k\choose 2k}2^{n-k}\), \(n\ge 0\). [Wilf 125, Example 4]\

5. Evaluate \(\displaystyle f(n)=\sum_{k\le n/2} (-1)^k{n-k\choose k}y^{n-2k}\). [Wilf 122, Example 3]\

6. Evaluate \(\displaystyle f(n)=\sum_{k} {2n+1\choose 2p+2k+1}{p+k\choose k}\). [W4.11(c)]\

7. Try to prove that \(\sum_k {n\choose k}{2n\choose n+k}={3n\choose n}\) via the snake oil method in three different ways: consider the sum \(\sum_k {n\choose k}{m\choose r-k}\) and the free variable being one of \(n\), \(m\), \(r\).

Asymptotic estimates

• Purpose of asymptotics [Knuth 439]

  • sometimes we do not have a closed form or it is hard to compare it to other quantities

  • \(\displaystyle S_n = \sum_{k=0}^n {3n\choose k}\sim 2{3n\choose n}\); \(\displaystyle S_n = {3n\choose n}\left(2-{4\over n} + O\left({1\over n^2}\right)\right)\)

  • how to compare it with \(F_{4n}\)? we need to approximate the binomial coefficient

  • purpose is to find accurate and concise estimates:
    \(H_n\) is \(\sum_{k\ge 1}^n 1/k\)vs.\(O(\log n)\)vs.\(\ln n + \gamma + O(n^{-1})\)

• Hierarchy of log-exp functions [Hardy, see Knuth 442]

  • the class \(\cal L\) of logarithmico-exponential functions: the smallest class that contains constants, identity function \(f(n) = n\), difference of any two functions from \(\cal L\), \(e^f\) for every \(f\in {\cal L}\), \(\ln f\) for every \(f\in {\cal L}\) that is “eventually positive”

  • every such function is identically zero, eventually positive or eventually negative

  • functions in \(\cal L\) form a hierarchy (every two of them are comparable by \(\prec\) or \(\asymp\))

• Notations

  • \(f(n) = O(g(n))\) iff \(\exists c: |f(n)|\le c|g(n)|\) (alternatively, for \(n\ge n_0\) for some \(n_0\))

  • \(f(n) = o(g(n))\) iff \(\lim_{n\to\infty} f(n)/g(n) = 0\)

  • \(f(n) = \Omega(g(n))\) iff \(\exists c: |f(n)|\ge c|g(n)|\) (alternatively, …)

  • \(f(n) = \Theta(g(n))\) iff \(f(n) = O(g(n))\) and \(f(n) = \Omega(g(n))\)

  • basic manipulation: \(O(f)+O(g) = O(|f|+|g|)\), \(O(f)O(g)=O(fg)=fO(g)\) etc.

  • meaning of \(O\) in sums

  • relative vs. absolute error

• Warm-ups

  1. Prove or disprove: \(O(f+g)=f + O(g)\) if \(f\) and \(g\) are positive. [K9.5]

  2. Multiply \(\ln n + \gamma + O(1/n)\) by \(n + O(\sqrt n)\). [K9.6]

  3. Compare \(n^{\ln n}\) with \((\ln n)^n\).

  4. Compare \(n^{\ln\ln\ln n}\) with \((\ln n)!\).

  5. Prove or disprove: \(O(x+y)^2 = O(x^2) + O(y^2)\). [K9.11]

• Common tricks

  • cut off series expansion (works for convergent series, Knuth 451)

  • substitution, e.g. \(\ln(1+2/n^2)\) with precision of \(O(n^{-5})\)

  • factoring (pulling the large part out), e.g. \({1\over n^2+n} = {1\over n^2}{1\over 1+{1\over n}}={1\over n^2}-{1\over n^3}+O(n^{-4})\)

  • division, e.g. \(\displaystyle {H_n\over \ln (n + 1)}= {\ln n + \gamma + O(n^{-1})\over (\ln n)(1+O(n^{-1}))}=1 + {\gamma\over \ln n} + O(n^{-1})\)

  • exp-log, i.e. \(f(x) = e^{\ln f(x)}\)

• Typical situations for approximation

  • Stirling formula: \(\displaystyle n! = \sqrt{2\pi n}\left({n\over e}\right)^n\left(1+{1\over 12n}+{1\over 288n^2}+O(n^{-3})\right)\)

  • harmonic numbers: \(H_n = \ln n + \gamma + {1\over 2n} - {1\over 12n^2} + O(n^{-4})\)

  • rational functions, e.g. \({n\over n+2} = {1\over 1+{2\over n}} = 1-{2\over n}+{4\over n^2}+O(n^{-3})\)

  • exponentials: \(e^{H_n}=ne^\gamma e^{O(1/n)}=ne^\gamma (1+O(1/n))=ne^\gamma + O(1)\)

  • rational function powered to \(n\), e.g. \(\left(1-{1\over n}\right)^n=e^{n\ln \left(1-{1\over n}\right)}= \exp\left(n\left({-1\over n}+O\left(n^{-2}\right)\right)\right) = e^{-1 + O(n^{-1}))} = {1\over e} + O(n^{-1})\)

  • binomial coefficient, e.g. \(2n\choose n\): factorials and Stirling formula $$\begin{aligned} {2n\choose n}={\sqrt{4\pi n}\left({2n\over e}\right)^{2n}(1+O(n^{-1}))\over 2\pi n\left({n\over e}\right)^{2n}(1+O(n^{-1}))^2}= \over \sqrt{\pi n}}(1+O(n^{-1}))

    \end{aligned}$$

• Exercises

  1. Estimate \(\ln(1+1/n)+ \ln(1-1/n)\) with abs. error \(O(n^{-3})\)

  2. Estimate \(\ln(2+1/n)- \ln(3-1/n)\) with abs. error \(O(n^{-2})\)

  3. Estimate \(\lg (n-2)\), abs. error \(O(n^{-2})\)

  4. Evaluate \(H_n^2\) with abs. error \(O(n^{-1})\).

  5. Estimate \(n^3/(2+n+n^2)\) with abs. error \(O(n^{-3})\)

  6. Prove or disprove: [K9.20] (b) \(e^{(1+O(1/n))^2} = e + O(1/n)\) cm (c) \(n! = O\left(((1-1/n)^nn)^n\right)\)

  7. Evaluate \((n+2+O(n^{-1}))^n\) with rel. error \(O(n^{-1})\). [K9.13]

  8. Compare \(H_{F_n}\) with \(F_{\lceil H_n\rceil}^2\) [K9.2]

  9. Estimate \(\sum_{k\ge 0} e^{-k/n}\) with abs. error \(O(n^{-1})\). [K9.7]

  10. Estimate \(H_n^5/\ln (n + 5)\) with abs. error \(O(n^{-2})\).

  11. Estimate \(2n\choose n\) with relative error \(O(n^{-2})\). [A1]

  12. Estimate \(2n+1\choose n\) with relative error \(O(n^{-2})\). [A2]

  13. Compare \((n!)!\) with \(((n-1)!)!\cdot (n-1)!^{n!}\). [K9.2c] (Homework if not enough time is left.)

Estimates of sums and products

• Warm-ups

  1. Let \(f(n) = \sum_{k=1}^n \sqrt k\). Show that \(f(n) = \Theta(n^{3/2})\). Find \(g(n)\) such that \(f(n) = g(n) + O(\sqrt n)\).

  2. Estimate \((n-2)!/(n-1)\) with abs. error \(O(n^{-2})\).

  3. For a constant integer \(k\), estimate \(n^{\underline{k}}/n^k\) with abs. error \(O(n^{-3})\). [A5]\

• Find a good estimate of \(P_n = {(2n-1)!!\over n!}\).

  • obviously \(\displaystyle 1.5^{n-1}\le {1\over 1}\cdot {3\over 2}\cdot {5\over 3}\cdot \dots \cdot {(2n-1)\over n}\le 2^{n-1}\)

  • we split the product into a “small” part (first \(k\) terms, each at least \(3/2\) except the first one) and a “large” part (remaining \(n-k\) terms); then
    \(P_n\ge \left({2k+1\over k+1}\right)^{n-k}\cdot 1.5^{k-1} = Q_n\cdot 1.5^{k-1}\); we estimate \(Q_n\)

  • if we try \(k = \alpha n\), then \(Q_n = 2^{n-\alpha n} \exp \left((n-\alpha n)\ln \left(1-{1\over 2(\alpha n + 1)}\right)\right)=2^{n(1-\alpha)}e^(1+O(n^{-1})),\) so \(P_n \ge (2^{1-\alpha}\cdot 1.5^\alpha)^n \Theta(1)\)

  • if we try \(k = \ln n\), then \(Q_n = \exp\left((n-\ln n)\left[\ln 2 + \ln \left(1-{1\over 2(1+\ln n)}\right)\right]\right);\) if we expand \(\ln\) into Taylor series, the error will be \(1/\ln^k n = \omega(n^{-1})\), so we can get relative error \(O(1)\) at best;
    anyway, if we carry it through, we get \(P_n = \Omega(2^n n^{-c} e^{-0.5n/\ln n})\)

  • If we try \(k = \sqrt n\), then \(\begin{aligned} Q_n &= \exp\left((n-\sqrt n)\left[\ln 2 + \ln \left(1-{1\over 2(1+\sqrt n)}\right)\right]\right)\\ &= 2^{n-\sqrt n}\exp\left((n-\sqrt n)\left[{-1\over 2\sqrt n} + {3\over 8n}-{7\over 24n^{3/2}}+O(n^{-2})\right]\right)\\ &= 2^{n-\sqrt n}\exp\left(-{\sqrt n\over 2} + {7\over 8}-{2\over 3\sqrt n}+O(n^{-1})\right), \end{aligned}\) thus \(P_n \ge 2^n \cdot 0.75^{\sqrt n}\cdot e^{\frac{-\sqrt n}{2}+\frac{7}{8}-\frac{2}{3\sqrt n}} (1+O(n^{-1})) = \Omega\left(2^n c^{\sqrt n}\right)\) for \(c\in (0, 1)\).

  • TODO compare with previous estimate from \(k=\ln n\); which is better?

  • another approach: \(P_n = {(2n)!\over n! 2^n n!} = {2n\choose n}/2^n = {2^n\over \sqrt{\pi n}}(1+O(n^{-1}))\)

• Estimate \(S_n = \sum_{k=1}^n {1\over n^2+k}\) with absolute error (a) \(O(n^{-3})\), (b) \(O(n^{-7})\). [Knuth 458/Problem 4] First approach: \({1\over n^2+k}={1\over n^2(1+k/n^2)}\) etc.; second approach: \(S_n = H_{n^2+n}-H_n\). (DU)

• Sums — gross bound on the tail: \(S_n = \sum_{0\le k\le n} k! = n!\left(1+{1\over n}+{1\over n(n-1)}+ \dots\right)\), all the terms except the first two are at most \(1/n(n+1)\), so \(S_n = n!(1+{1\over n}+n{1\over n(n-1)}) = n!(1+O(n^{-1}))\)

• Sums — make the tail infinite: \(\begin{aligned} n!\sum_{k=0}^n{(-1)^k\over k!} &= n!\left(\sum_{k=0}^\infty{(-1)^k\over k!}-\sum_{k\ge n+1}{(-1)^k\over k!}\right)\\ &= n!\left(e^{-1}-O\left({1\over (n+1)!}\right)\right)= {n!\over e}+O(n^{-1}) \end{aligned}\)

• Estimate \(S_n=\sum_{k=0}^n {3n\choose k}\) with relative error \(O(n^{-2})\). We split the sum into a “small” and a “large” part at \(b\) (which is yet to be determined). \(\begin{aligned} \sum_{k=0}^{n} \binom {3n}k &= \sum_{k=0}^{n} \binom {3n}{n-k} = \sum_{0\leq k<b} \binom {3n}{n-k}+\sum_{b\le k\le n} \binom {3n}{n-k}.\\ \binom{3n}{n-k} &= \binom{3n}{n} \frac{n(n-1)\cdot\ldots \cdot 1}{(2n+1)(2n+2)\ldots(2n+k)} =\\ &= \binom{3n}{n}\cdot\frac{n^k}{(2n)^k}\frac{\prod_{j=0}^{k-1}\left(1-\frac jn\right)}{\prod_{j=1}^k \left(1+\frac j{2n}\right)}=\binom{3n}{n}\cdot\frac{1}{2^k}\cdot\left[1-\frac{3k^2-k}{4n}+O\left(\frac{k^4}{n^2}\right)\right].\\ \sum_{b\le k\le n} \binom {3n}{n-k} &\le n\cdot \binom{3n}{n-b}=\binom{3n}{n}\cdot \frac{1}{2^b} O(n)=\binom{3n}{n}\cdot O\left(n^{-2}\right) \text{if } \sqrt{n} \succ b \geq 3\lg n.\\ \sum_{0\leq k<3\lg n}\frac{1}{2^k} &= 2-\frac{1}{2^{3\lg n}}=2+O(n^{-3}).\\ -\frac{3}{4n}\sum_{0\leq k<3\lg n}\frac{k^2}{2^k} &= \frac{-9}{2n}+O(n^{-3}).\\ +\frac{1}{4n}\sum_{0\leq k<3\lg n}\frac{k}{2^k} &= \frac{1}{2n}+O(n^{-3}).\\ O(n^{-2})\cdot\sum_{0\leq k<3\lg n}\frac{k^4}{2^k} &= O(n^{-2}) \end{aligned}\)

• Estimate \(S_n=\sum_{k=0}^n \binom{4n+1}{k+1}\) with relative error \(O(n^{-2})\). \(\binom{4n+1}{k+1}=\binom{4n}{k+1}+\binom{4n}{k};\) \(S_n=\sum_{k=0}^n \binom{4n+1}{k+1}=\sum_{k=0}^n\binom{4n}{k}+ \sum_{k=0}^n\binom{4n}{k+1}=\sum_{k=0}^n\binom{4n}{k}+\sum_{k=1}^{n+1}\binom{4n}{k};\) \(S_n=2\sum_{k=0}^n\binom{4n}{k}+\binom{4n}{n+1}-\binom{4n}{0}.\) \(Q_n=\sum_{k=0}^n\binom{4n}{k}=\sum_{k=0}^n\binom{4n}{n-k};\) \(\binom{4n}{n-k}=\binom{4n}{n}\cdot\frac{\prod_{j=0}^{k-1}(n-j)}{\prod_{j=1}^{k}(3n+j)}=\binom{4n}{n}\cdot\left(\frac 13\right)^3\cdot\frac{\prod_{j=0}^{k-1}(1-j/n)}{\prod_{j=1}^{k}(1+j/3n)}\) \(Q_n=\sum_{0\leq k\leq 2\log_3 n}\binom{4n}{n-k}+\sum_{2\log_3 n\leq k<n}\binom{4n}{n-k}\) \(\sum_{2\log_3 n\leq k<n}\binom{4n}{n-k}=O\left(n\cdot\binom{4n}{n-\lceil 2\log_3 n\rceil} \right)=O\left(\binom{4n}{n}\cdot\frac 1n \right).\) \(\frac{\prod_{j=0}^{k-1}(1-j/n)}{\prod_{j=1}^{k}(1+j/3n)}=\frac{1-\frac 1n\cdot\sum_{0\leq j<k}j+O\left(\frac{k^4}{n^2}\right)}{1+\frac{1}{3n}\cdot\sum_{1< j\leq k}j+O\left(\frac{k^4}{n^2}\right)} = 1+\frac{2k^2+k}{3n}+O\left(\frac{\log^n}{n^2}\right),\) \(\sum_{0\leq k\leq 2\log_3 n}\binom{4n}{n-k}=\binom{4n}{n}\cdot\sum_{0\leq k\leq 2\log_3 n}\left( \frac 13\right)^k\cdot[1+\frac{2k^2+k}{3n}+O\left(\frac{\log^n}{n^2}\right)]=\) \(=\frac 32\cdot \binom{4n}{n} (1+O(n^{-1})).\) \(\binom{4n}{n+1}= \binom{4n}{n}\cdot\frac{3n}{n+1}=3\cdot\binom{4n}{n}(1+O(n^{-1}));\) \(S_n=6\cdot\binom{4n}{n}(1+O(n^{-1})).\)

• How many bits are needed to represent a binary tree with \(n\) internal nodes?

  • we need just the internal vertices to capture the structure; what is the relation between the number of internal vertices and total number of vertices?

  • imagine labeling the vertices by \(1,2,\dots,n\) in such a way that we get a binary search tree (descendants in the left subtree are smaller, in the right subtree are larger); by summing over possible roots of the tree we get \(t_n = \sum_{i=1}^n t_{i-1} t_{n-i}\); \(t_0 = 1\)

  • this is the same as for Catalan numbers, so \(t_n = {2n\choose n}{1\over n+1}\)

1) Find explicit formulas for the following sequences:

1) \(a_{n+1} = 3a_n+2\) for \(n\ge 0\), \(a_0=0\)

Solution

\(3x/(1-x)(1-3x)\)
\(3^n-1\)

2. \(a_{n+1} = \alpha a_n + \beta\) for \(n\ge 0\), \(a_0=0\)

Solution

\(\beta x/(1-x)(1-\alpha x)\)
\({\alpha^n-1\over \alpha-1}\beta\)

3. \(a_{n+1} = a_n/3 +1\) for \(n\ge 0\), \(a_0=1\)

Solution

\({3/2\over 1-x}-{1/2\over 1-x/3}\)
\({3^{n+1}-1\over 2\cdot 3^n}\)

4. \(a_{n+2} = 2a_{n+1}-a_n\) for \(n\ge 0\), \(a_0=0\), \(a_1=1\)

Solution

\(x/(1-x)^2\)
\(n\)

5. \(a_{n+2} = 3a_{n+1}-2a_n+3\) for \(n>0\), \(a_0=1\), \(a_1=2\)

Solution

\({4\over 1-2x}-{3\over (1-x)^2}\)
\(2^{n+2}-3n-3\)

6. \(a_n = 2a_{n-1}-a_{n-2}+(-1)^n\) for \(n>1\), \(a_0=a_1=1\)

Solution

\({1/2\over (1-x)^2}-{1/4\over 1-x}+{1/4\over 1+x}\)
\({2n+3+(-1)^n\over 4}\)

7. \(a_n = 2a_{n-1}-n\cdot(-1)^n\) for \(n\ge 1\), \(a_0=0\)

Solution

\({x/9-2/9\over (1+x)^2}+{2/9\over 1-2x}\)
\({2^{n+1}-(3n+2)(-1)^n\over 9}\)

8. \(a_n = 3a_{n-1} + {n\choose 2}\) for \(n\ge 1\), \(a_0=2\)

Solution

\(\)
\(\dfrac{1}{8}(19\cdot 3^n-2n(n+2)-3)\)

9. \(a_n = 2a_{n-1}-a_{n-2}-2\) for \(n > 1\), \(a_0=a_{10}=0\)

Solution

\(n(a_1+1-n)\), so with \(a_{10}\), \(a_n=n(10-n)\)

10. \(a_n = 4(a_{n-1}-a_{n-2})+(-1)^n\) for \(n \ge 2\), \(a_0=1\), \(a_1=4\)

Solution

\(\begin{align*} {1+x+x^2\over (1+x)(1-2x)^2} &= {1\over 9}{1\over 1+x} +\left({-5\over 18}\right) {1\over 1-2x} + \left({7\over 6}\right){1\over (1-2x)^2} \end{align*}\)
\({1\over 9}(-1)^n-{5\over 18}\cdot 2^n+{7\over 6}(n+1)\cdot 2^n\)

11. \(a_n = -3a_{n-1}+a_{n-2}+3a_{n-3}\) for \(n\ge 3\), \(a_0=20\), \(a_1=-36\), \(a_2=60\)

Solution

\(\)
\(5(-3)^n+18(-1)^n-3\)

12. \(a_n = -3a_{n-1}+a_{n-2}+3a_{n-3}+128n\) for \(n\ge 3\), \(a_0=0\), \(a_1=0\), \(a_2=0\)

Solution

\(\)
\(8n^2+28n-29-11(-3)^n+40(-1)^n\)

Key Points

Statistics

Overview

Teaching: min
Exercises: min

Questions

Objectives

Basics

Dataset: Any group of values retrieved through a common method/procedure of collection.

Weighted Mean: Weighted mean of a group of numbers with weights given in percentages and scores given in the same range is given by:

Note: Always question about how the weights and categories were collected and why and to what extent it holds importance.

Standard Deviation: A quantity expressing by how much the members of a group differ from the mean value for the group. In more formal terms, it’s the average distance between a datapoint and the mean of the dataset and is given by the equation:

z-score: The z-score is the distance in standard deviation units for any given observation:

Empirical Rule or Three-Sigma Rule: Most of the data points fall within three standard deviations, with 68% would lie within 1 sd, 95% within 2 sd and 99.7% would fall within 3 sd.

Note: It only works for symmetrically distributed data

Percentile Score: The percentile score for any value x in the dataset is given by:

Probabilities

Event: An event is a set of probabilities e.g. When rolling a dice, an event A = (4,2,6) represents the event that an even number will appear.

Sample Space: All possible outcomes for a particular random experiment constitute it’s sample space.

Probability: The number of favorable outcomes divided by the number of all possible outcomes in a given situation.

Types of Probabilities

1. Classical Probability: The flipping of a coin is an example of classical probability because it is known that there are two sides to it and the likelihood of any one of them turning up is 50%. Objective Probabilities are based on calculations and classical probability is a type of Objective Probability.

2. Empirical Probability: The probabilities are the ones that are based on previously known data. For instance, Messi scoring more than ten goals in this season of FIFA is an example of Empirical probability because it is calculated on the basis of Messi’s previous record. This too is an example of Objective Probability since it’s also based on calculations.

3. Subjective Probability: These probabilities are not based on mathematical calculations and people use their opinions, their experiences to assert their viewpoints with some amount of relevant data that makes their own point stronger.

Addition Rule: The addition rule in probabilities is a rule that ensures that an event is not counted twice when calculating probabilities (where P(overlap) is the probability of both E1 and E2 occurring:

Conditional Probability: The probability of occurrence of an event given that some other event has already occurred.

Independent Events: If the probability of two events are completely unrelated. If we can prove that the probability of two occurring together is equal to the product of their individual probabilities, then we can say that the two events are independent.

Random Variable: The result of an experiment that has random outcomes is called a random variable. They can be of two types:

1. Discrete RV: The number of drinks a person will order at Tank is an example of discrete random variable because it has to be a whole number.
2. Continuous RV: The waiting time in line before one can order at a Burger King is an example of continuous variable because there are no fixed values that could be outcomes. The possibilities are infinite and continuous.
3. Binomial RV: When an event has only two possible outcomes, the result is called a Binomial Random Variable.

Probability Density: The curves that represent the distribution of probabilities are called Probability Density curves.

Sampling

There are a few conditions that need to be considered before sampling is done from various sources:

1. Size to Cost ratio: The appropriate size of the sample based upon the cost per data point in the sample
2. Inherent Bias: If any bias was knowingly/unknowingly introduced while creating the sample, it will need to be considered!
3. Quality of Sample

A simple random sample is the gold standard when collecting samples. This means, that any given point during the sample selection process, any individual has the same probability of being chosen as any other individual.

Some alternative sampling methods are:

1. kth point: The first data point is selected and then every kth data point is selected in this method.
2. Opportunity Sampling: The first n values are selected from the total data.
3. Stratified Sampling: The whole sample is broken out into homogenous groups. Then we select few samples from each strata.
4. Cluster Sampling: The whole sample is collected from heterogenous groups with data points with different characteristics. Then we select few samples from each group.

Confidence Intervals

As the name suggests, the confidence intervals present a level of confidence for a given interval.

Hypothesis Testing

The process to be able to test a hypothesis that has been presented.

Visualization Tips

1. Tables: Useful for detailed recording and sharing of actual data
2. Frequency Table: Displays the frequency of each observation in the data set
3. Dot Plots: When you want to convey information that is discrete and individual to each observation.
4. Histograms: When you want to convey frequencies of grouped bins, this could be useful.
5. Pie Charts: Relative Frequency distributions are best represented with Pie Charts. They are also useful for representing distributions among qualitative variables where histograms wouldn’t be a very good measure.

Key Points

Types of Data

Overview

Teaching: min
Exercises: min

Questions

Objectives

Cross-Sectional Data

It consists of variables, that are either quantitative or qualitative for a lot of different observations (usually called ‘cases’) and are taken from some defined population. All the cases are registered at a single point in time, or a reasonably short period of time. The techniques commonly used for this kind of data are ‘t-tests’ analysis of variance or regression depending on the kind and number of variables that are there in the data. It is noteworthy that each observation of a given variable or each set of variables are independent of every other observation in the dataset. The independence of the variables is a critical assumption when modelling cross-sectional data.

It is called cross-sectional data because we are measuring a cross-section of a defined population at a particular point in time or a very short period in time.

Time Series Data

If measurement on variables are taken over or through time. Every variable in a time series dataset is measured at equally spaced time intervals. Usually, the observations are not independent of each other in this case. Time series data can be classified into two types other thank the univariate and multivariate distinction and they are discussed in the time series chapter.

Key Points

Random Variables and Event Probabilities

Overview

Teaching: min
Exercises: min

Questions

Objectives

Random Variables and Event Probabilities

Random variables

Let \(Y\) be the result of a fair coin flip. Not a general coin flip, but a specific instance of flipping a specific coin at a specific time. Defined this way, \(Y\) is what’s known as a random variable, meaning a variable that takes on different values with different probabilities.^[Random variables are conventionally written using upper-case letters to distinguish them from ordinary mathematical variables which are bound to single values and conventionally written using lower-case letters.]

Probabilities are scaled between 0% and 100% as in natural language. If a coin flip is fair, there is a 50% chance the coin lands face up (“heads”) and a 50% chance it lands face down (“tails”). For concreteness and ease of analysis, random variables will be restricted to numerical values. For the specific coin flip in question, the random variable \(Y\) will take on the value 1 if the coin lands heads and the value 0 if it lands tails.

Events and probability

An outcome such as the coin landing heads is called an event in probability theory. For our purposes, events will be defined as conditions on random variables. For example, \(Y = 1\) denotes the event in which our coin flip lands heads. The functional \(\mbox{Pr}[\, \cdot \,]\) defines the probability of an event. For example, for our fair coin toss, the probability of the event of the coin landing heads is written as

[\mbox{Pr}[Y = 1] = 0.5.]

In order for the flip to be fair, we must have \(\mbox{Pr}[Y = 0] = 0.5\), too. The two events \(Y = 1\) and \(Y = 0\) are mutually exclusive in the sense that both of them cannot occur at the same time. In probabilistic notation,

[\mbox{Pr}[Y = 1 \ \mbox{and} \ Y = 0] = 0.]

The events \(Y = 1\) and \(Y = 0\) are also exhaustive, in the sense that at least one of them must occur. In probabilistic notation,

[\mbox{Pr}[Y = 1 \ \mbox{or} \ Y = 0] = 1.]

In these cases, events are conjoined (with “and”) and disjoined (with “or”). These operations apply in general to events, as does negation. As an example of negation,

[\mbox{Pr}[Y \neq 1] = 0.5.]

Sample spaces and possible worlds

Even though the coin flip will have a specific outcome in the real world, we consider alternative ways the world could have been. Thus even if the coin lands heads \((Y = 1)\), we entertain the possibility that it could’ve landed tails \((Y = 0)\). Such counterfactual reasoning is the key to understanding probability theory and applied statistical inference.

An alternative way the world could be, that is, a possible world, will determine the value of every random variable. The collection of all such possible worlds is called the sample space.^[The sample space conventionally written as \(\Omega\), the capitalized form of the last letter in the Greek alphabet.] The sample space may be conceptualized as an urn containing a ball for each possible way the world can be. On each ball is written the value of every random variable.^[Formally, a random variable \(X\) can be represented as a function from the sample space to a real value, i.e., \(X: \Omega \rightarrow \mathbb{R}\).

For each possible world \(\omega \in \Omega\), the variable \(X\) takes on a specific value

\(X(\omega) \in \mathbb{R}\).]

Statistic is the logic of uncertainty.

A sample space is the set of all possible outcomes of an experiment

An event is a subset of sample space

Now consider the event \(Y = 0\), in which our coin flip lands tails. In some worlds, the event occurs (i.e., \(0\) is the value recorded for \(Y\)) and in others it doesn’t. An event picks out the subset of worlds in which it occurs.^[Formally, an event is defined by a subset of the sample space, \(E \subseteq \Omega\).]

Naive definetion of probability

P(A) = #favorable outcome / #possible outcomes

Example: tossing of fair coins

Assumes all outcome equally likely, finite sample space

Sampling table

choose k objects out of n

1. Don’t lose common sense
2. Do check answers, especially by doing simple and extreme cases
3. Label people , objects etc. If have n people, then label them 1,2…n

Example: 10 people, split into them of 6, team of 4 => \(\binom{10}{6}\) 2 teams of 5 => \(\binom{10}{5}\) /2

Problem: pick k times from set of n objects, where order doesn’t matter, with replacement.

Extreme cases: k = 0; k = 1; n = 2

Equiv : how many ways are there to put k indistinguishable particles into n distinguishable boxes?

Axioms of Probability

Non-naive definition

Probability sample consists of S and P, where S is sample space, and P , a function which takes an event \(A\subseteq S\) as input, returns \(P(A) \in [0,1]\) as output.

such that

1. \(P(\phi) = 0, P(S) = 1\)\
2. \(P(U_{n=1}^{\infty}A_n) = \sum_{n=1}^{\infty} P(A_n)\) if \(A1,A2..An\) are disjoint (not overlap)

Birthday Problem

(Exclude Feb 29, assume 365 days equally likely, assume indep. of birth)

k people , find prob. that two have same birthday

If k > 365, prob. is 1

Let k <= 365, \(P(no match) =\frac{ 365 * 364 *... (365 - k + 1) }{365^k}\)

P(match) ~ 50.7%, if k = 23; 97% if k = 50; 99.9999%, if k = 100

\(\binom{k}{2} = \frac{k(k-1)}{2}\) \(\binom{23}{2} = 253\)

Properties of Probability

1. \(P(A^c) = 1 - P(A)​\)\
2. If \(A \subseteq B\) , then \(P(A) \subseteq P(B)\)\
3. \[P(A\cup B) = P(A) + P(B) - P(A\cap B)\]

[P(A\cup B\cup C) = P(A) + P(B) + P(C) - P(A\cap B) - P(A\cap C) - P(B\cap C) + P(A\cap B\cap C)]

Proof:

1. \(1 = P(S) = P(A\cap A^c) = P(A) + P(A^c)\)\
2. \(B = A\cup(B\cap A^c)\) \(P(B) = P(A)+P(B\cap A^c)\)\
3. \[P(A\cup B) = P(A\cap (B\cap A^c)) = P(A) + P(B\cap A^c)\]

General case:

deMontmort’s Problem(1713)

matching problem

n cards labeled 1 to n, flipping cards over one by one, you win if the card that you name is the card that appears.

Let \(A_j\) be the event, ‘‘jth card matches”

\(P(A_j) = 1 / n\) since all position equally likely for card labeled j\

\(P(A_1\cap A_2) = (n-2)! / n! = 1/n(n-1)\)\

…

\(P(A_1\cap … A_k) = (n-k)! / n!\)\

\(P(A_1\cup …A_n) = n*1/n - n(n-1)/2 * 1/n(n-1) + …\)\

\(= 1 - 1/2! + 1/3! - 1/4! … (-1)^n1/n!\) \(\approx 1- 1/e\)

Story proof- proof by interpretation

Ex1 \(\binom{n}{k}\) = \(\binom{n}{n-k}\)

Ex2 n\(\binom{n-1}{k-1}\) = k\(\binom{n}{k}\) pick k people out of n, with one desigenate as president.

Ex3 \(\binom{m+n}{k} = \sum_{j=0}^k \binom{m}{j} \binom{n}{k-j}\) (vandermonde)

Simulating random variables

We are now going to turn our attention to computation, and in particular, simulation, with which we will use to estimate event probabilities.

The primitive unit of simulation is a function that acts like a random number generator. But we only have computers to work with and they are deterministic. At best, we can created so-called pseudorandom number generators. Pseudorandom number generators, if they are well coded, produce deterministic streams of output that appear to be random.^[There is a large literature on pseudorandom number generators and tests for measurable differences from truly random streams.]

For the time being, we will assume we have a primitive pseudorandom number generator uniform_01_rng(), which behaves roughly like it has a 50% chance of returning 1 and a 50% chance of returning 0.^[The name arises because random variables in which every possible outcome is equally likely are said to be uniform.]

Suppose we want to simulate our random variable \(Y\). We can do so by calling uniform_01_rng and noting the answer.

A simple program to generate a realization of a random coin flip, assign it to an integer variable y, and print the result could be coded as follows.^[Computer programs are presented using a consistent pseudocode, which provides a sketch of a program that should be precise enough to be coded in a concrete programming language. R implementations of the pseudocode generate the results and are available in the source code repository for this book.]

import random

y = random.randint(0, 1)
print("y =", y)

The variable y is declared to be an integer and assigned to the result of calling the uniform_01_rng() function.^[The use of a lower-case \(y\) was not accidental. The variable \(y\) represents an integer, which is the type of a realization of a random \(Y\) representing the outcome of a coin flip. In code, variables are written in typewriter font (e.g., y), whereas in text they are written in italics like other mathematical variables (e.g., \(y\)).] The print statement outputs the quoted string y =   followed by the value of the variable y. Executing the program might produce the following output.

y = 1

If we run it a nine more times, it might print

for i in range(9):
    y = random.randint(0,1)
    print("y =", y)

y = 0
y = 1
y = 0
y = 1
y = 1
y = 0
y = 1
y = 1
y = 0

When we say it might print these things, we mean the results will depend on the state of the pseudorandom number generator.

Seeding a simulation

Simulations can be made exactly reproducible by setting what is known as the seed of a pseudorandom number generator. This seed establishes the deterministic sequence of results that the pseudorandom number generator produces. For instance, contrast the program

random.seed(1234)
for n in range(10):
    print(random.randint(0, 1)')
for n in range(10):
    print(random.randint(0, 1))

which produces the output

1 0 1 1 0 1 0 0 0 0  
0 0 0 1 0 1 1 1 0 0 

with the program

random.seed(1234)
for n in range(10):
    print(random.randint(0, 1)')
random.seed(1234)
for n in range(10):
    print(random.randint(0, 1))

which produces

1 0 0 0 0 0 0 1 0 0 
1 0 0 0 0 0 0 1 0 0 

Resetting the seed in the second case causes exactly the same ten pseudorandom numbers to be generated a second time. Every well-written pseudorandom number generator and piece of simulation code should allow the seed to be set manually to ensure reproducibility of results.^[Replicability of results with different seeds is a desirable, but stricter condition.]

Using simulation to estimate event probabilities

We know that \(\mbox{Pr}[Y = 1]\) is 0.5 because it represents the flip of a fair coin. Simulation based methods allow us to estimate event probabilities straightforwardly if we can generate random realizations of the random variables involved in the event definitions.

For example, we know we can generate multiple simulations of flipping the same coin. That is, we’re not simulating the result of flipping the same coin ten different times, but simulating ten different realizations of exactly the same random variable, which represents a single coin flip.

The fundamental method of computing event probabilities will not change as we move through this book. We simply simulate a sequence of values and return the proportion in which the event occurs as our estimate.

For example, let’s simulate 10 values of \(Y\) again and record the proportion of the simulated values that are 1. That is, we count the number of time the event occurs in that the simulated value \(y^{(m)}\) is equal to 1.

M = 100  # Set the value of M
occur = 0

for m in range(1, M+1):
    y = random.uniform(0, 1)
    occur += (y == 1)

estimate = occur / M
print(f"estimated Pr[Y = 1] = {estimate}")

The equality operator is written as ==, as in the condition y[m] == 1 to distinguish it from the assignment statement y[m] = 1, which sets the value of y[m] to 1. The condition expression y[m] == 1 evaluates to 1 if the condition is true and 0 otherwise.

If we let uniform_01_rng(M) be the result of generating M pseudorandom coin flip results, the program can be shortened to

M = 100  # Set the value of M
y = [random.uniform(0, 1) for i in range(M)]
occur = sum([1 for i in y if i == 1])
estimate = occur / M
print(f"estimated Pr[Y = 1] = {estimate}")

A condition such as y == 1 on a sequence returns a sequence of the same length with value 1 in positions where the condition is true. For instance, if

y = (2, 1, 4, 2, 2, 1)

then

y == 2

evaluates to

(1, 0, 0, 1, 1, 0).

Thus sum(y == 1) is the number of positions in the sequence y which have the value 1. Running the program provides the following estimate based on ten simulation draws.

import numpy as np

M = 10
y_sim = np.random.binomial(1, 0.5, M)
for n in range(M):
    print(y_sim[n], end=' ')
print(f"estimated Pr[Y = 1] = {np.sum(y_sim) / M}")

Let’s try that a few more times.

M = 10

for k in range(1, 11):
    y_sim = np.random.binomial(1, 0.5, M)
    for n in range(M):
        print(y_sim[n], end=' ')
    print(f"estimated Pr[Y = 1] = {np.sum(y_sim) / M}")

1 1 0 0 0 1 1 1 0 0 estimated Pr[Y = 1] = 0.5
0 0 1 1 1 1 0 0 1 0 estimated Pr[Y = 1] = 0.5
0 0 0 1 1 0 1 0 1 1 estimated Pr[Y = 1] = 0.5
1 0 0 1 0 0 1 0 1 0 estimated Pr[Y = 1] = 0.4
1 1 0 0 1 0 0 0 0 1 estimated Pr[Y = 1] = 0.4
0 0 1 0 1 1 1 1 1 0 estimated Pr[Y = 1] = 0.6
1 1 0 1 1 1 0 0 0 1 estimated Pr[Y = 1] = 0.6
0 0 0 1 1 1 0 1 0 1 estimated Pr[Y = 1] = 0.5
0 0 0 0 1 0 1 0 0 1 estimated Pr[Y = 1] = 0.3
0 1 0 1 1 1 0 1 1 0 estimated Pr[Y = 1] = 0.6

The estimates are close, but not very exact. What if we use 100 simulations?

M = 100
y_sim = np.random.binomial(1, 0.5, M)
for n in range(M):
    print(y_sim[n], end=' ')
print(f"estimated Pr[Y = 1] = {np.sum(y_sim) / M}")

0 0 0 0 0 1 0 0 0 0 1 0 0 1 1 0 0 0 1 1 0 0 1 1 1 1 0 1 1 1 0 0 0 0 1 0 0 0 1 0 0 0 0 1 1 1 1 0 1 0 1 1 0 1 0 1 0 0 0 1 0 1 1 1 1 1 0 0 0 1 1 0 0 1 0 0 1 1 1 1 1 1 1 1 0 0 1 0 1 1 0 0 1 1 1 1 0 0 0 0 
estimated Pr[Y = 1] = 0.48

That’s closer than most of the estimates based on ten simulation draws. Let’s try that a few more times without bothering to print all 100 simulated values,

for k in range(1, 11):
    y_sim = np.random.binomial(1, 0.5, M)
    print(f"estimated Pr[Y = 1] = {np.sum(y_sim) / M}")

estimated Pr[Y = 1] = 0.52
estimated Pr[Y = 1] = 0.58
estimated Pr[Y = 1] = 0.55
estimated Pr[Y = 1] = 0.37
estimated Pr[Y = 1] = 0.52
estimated Pr[Y = 1] = 0.48
estimated Pr[Y = 1] = 0.53
estimated Pr[Y = 1] = 0.53
estimated Pr[Y = 1] = 0.5
estimated Pr[Y = 1] = 0.53

What happens if we let \(M = 10,000\) simulations?

M = 10000

for k in range(1, 11):
    y_sim = np.random.binomial(1, 0.5, M)
    print(f"estimated Pr[Y = 1] = {np.sum(y_sim) / M}")

estimated Pr[Y = 1] = 0.5074
estimated Pr[Y = 1] = 0.4991
estimated Pr[Y = 1] = 0.5088
estimated Pr[Y = 1] = 0.5015
estimated Pr[Y = 1] = 0.4926
estimated Pr[Y = 1] = 0.4987
estimated Pr[Y = 1] = 0.4959
estimated Pr[Y = 1] = 0.5004
estimated Pr[Y = 1] = 0.4928
estimated Pr[Y = 1] = 0.5032

Now the estimates are very close to the true probability being estimated (i.e., 0.5, because the flip is fair). This raises the questions of how many simulation draws we need in order to be confident our estimates are close to the values being estimated.

Law of large numbers

Visualization in the form of simple plots goes a long way toward understanding concepts in statistics and probability. A traditional way to plot what happens as the number of simulation draws \(M\) increases is to keep a running tally of the estimate as each draw is made and plot the estimated event probability \(\mbox{Pr}[Y = 1]\) for each \(m \in 1:M\).^[See, for example, the quite wonderful little book, Bulmer, M.G., 1965. Principles of Statistics. Oliver and Boyd, Edinburgh.]

To calculate such a running tally of the estimate at each online, we can do this:

M = 100
y = np.zeros(M)
estimate = np.zeros(M)
occur = 0

for m in range(1, M+1):
    y[m-1] = np.random.randint(0, 2)
    occur += (y[m-1] == 1)
    estimate[m-1] = occur / m
    print(f"estimated Pr[Y = 1] after {m} trials = {estimate[m-1]}")

estimated Pr[Y = 1] after 1 trials = 0.0
estimated Pr[Y = 1] after 2 trials = 0.5
estimated Pr[Y = 1] after 3 trials = 0.3333333333333333
estimated Pr[Y = 1] after 4 trials = 0.25
estimated Pr[Y = 1] after 5 trials = 0.2
estimated Pr[Y = 1] after 6 trials = 0.3333333333333333
estimated Pr[Y = 1] after 7 trials = 0.42857142857142855
estimated Pr[Y = 1] after 8 trials = 0.375
estimated Pr[Y = 1] after 9 trials = 0.4444444444444444
estimated Pr[Y = 1] after 10 trials = 0.4
estimated Pr[Y = 1] after 11 trials = 0.36363636363636365
estimated Pr[Y = 1] after 12 trials = 0.3333333333333333
estimated Pr[Y = 1] after 13 trials = 0.3076923076923077
estimated Pr[Y = 1] after 14 trials = 0.2857142857142857
estimated Pr[Y = 1] after 15 trials = 0.3333333333333333
estimated Pr[Y = 1] after 16 trials = 0.375
estimated Pr[Y = 1] after 17 trials = 0.35294117647058826
estimated Pr[Y = 1] after 18 trials = 0.3888888888888889
estimated Pr[Y = 1] after 19 trials = 0.42105263157894735
estimated Pr[Y = 1] after 20 trials = 0.4
estimated Pr[Y = 1] after 21 trials = 0.42857142857142855
estimated Pr[Y = 1] after 22 trials = 0.45454545454545453
estimated Pr[Y = 1] after 23 trials = 0.43478260869565216
estimated Pr[Y = 1] after 24 trials = 0.4166666666666667
estimated Pr[Y = 1] after 25 trials = 0.4
estimated Pr[Y = 1] after 26 trials = 0.4230769230769231
estimated Pr[Y = 1] after 27 trials = 0.4444444444444444
estimated Pr[Y = 1] after 28 trials = 0.42857142857142855
estimated Pr[Y = 1] after 29 trials = 0.41379310344827586
estimated Pr[Y = 1] after 30 trials = 0.4
estimated Pr[Y = 1] after 31 trials = 0.3870967741935484
estimated Pr[Y = 1] after 32 trials = 0.40625
estimated Pr[Y = 1] after 33 trials = 0.3939393939393939
estimated Pr[Y = 1] after 34 trials = 0.4117647058823529
estimated Pr[Y = 1] after 35 trials = 0.4
estimated Pr[Y = 1] after 36 trials = 0.3888888888888889
estimated Pr[Y = 1] after 37 trials = 0.3783783783783784
estimated Pr[Y = 1] after 38 trials = 0.39473684210526316
estimated Pr[Y = 1] after 39 trials = 0.41025641025641024
estimated Pr[Y = 1] after 40 trials = 0.425
estimated Pr[Y = 1] after 41 trials = 0.43902439024390244
estimated Pr[Y = 1] after 42 trials = 0.42857142857142855
estimated Pr[Y = 1] after 43 trials = 0.4186046511627907
estimated Pr[Y = 1] after 44 trials = 0.4090909090909091
estimated Pr[Y = 1] after 45 trials = 0.4
estimated Pr[Y = 1] after 46 trials = 0.391304347826087
estimated Pr[Y = 1] after 47 trials = 0.40425531914893614
estimated Pr[Y = 1] after 48 trials = 0.4166666666666667
estimated Pr[Y = 1] after 49 trials = 0.42857142857142855
estimated Pr[Y = 1] after 50 trials = 0.42
estimated Pr[Y = 1] after 51 trials = 0.43137254901960786
estimated Pr[Y = 1] after 52 trials = 0.4423076923076923
estimated Pr[Y = 1] after 53 trials = 0.4528301886792453
estimated Pr[Y = 1] after 54 trials = 0.46296296296296297
estimated Pr[Y = 1] after 55 trials = 0.4727272727272727
estimated Pr[Y = 1] after 56 trials = 0.48214285714285715
estimated Pr[Y = 1] after 57 trials = 0.47368421052631576
estimated Pr[Y = 1] after 58 trials = 0.4827586206896552
estimated Pr[Y = 1] after 59 trials = 0.4745762711864407
estimated Pr[Y = 1] after 60 trials = 0.4666666666666667
estimated Pr[Y = 1] after 61 trials = 0.47540983606557374
estimated Pr[Y = 1] after 62 trials = 0.46774193548387094
estimated Pr[Y = 1] after 63 trials = 0.47619047619047616
estimated Pr[Y = 1] after 64 trials = 0.46875
estimated Pr[Y = 1] after 65 trials = 0.46153846153846156
estimated Pr[Y = 1] after 66 trials = 0.4696969696969697
estimated Pr[Y = 1] after 67 trials = 0.47761194029850745
estimated Pr[Y = 1] after 68 trials = 0.4852941176470588
estimated Pr[Y = 1] after 69 trials = 0.4782608695652174
estimated Pr[Y = 1] after 70 trials = 0.4714285714285714
estimated Pr[Y = 1] after 71 trials = 0.4788732394366197
estimated Pr[Y = 1] after 72 trials = 0.4861111111111111
estimated Pr[Y = 1] after 73 trials = 0.4931506849315068
estimated Pr[Y = 1] after 74 trials = 0.5
estimated Pr[Y = 1] after 75 trials = 0.49333333333333335
estimated Pr[Y = 1] after 76 trials = 0.4868421052631579
estimated Pr[Y = 1] after 77 trials = 0.4935064935064935
estimated Pr[Y = 1] after 78 trials = 0.5
estimated Pr[Y = 1] after 79 trials = 0.5063291139240507
estimated Pr[Y = 1] after 80 trials = 0.5125
estimated Pr[Y = 1] after 81 trials = 0.5185185185185185
estimated Pr[Y = 1] after 82 trials = 0.5121951219512195
estimated Pr[Y = 1] after 83 trials = 0.5180722891566265
estimated Pr[Y = 1] after 84 trials = 0.5119047619047619
estimated Pr[Y = 1] after 85 trials = 0.5176470588235295
estimated Pr[Y = 1] after 86 trials = 0.5116279069767442
estimated Pr[Y = 1] after 87 trials = 0.5172413793103449
estimated Pr[Y = 1] after 88 trials = 0.5227272727272727
estimated Pr[Y = 1] after 89 trials = 0.5280898876404494
estimated Pr[Y = 1] after 90 trials = 0.5333333333333333
estimated Pr[Y = 1] after 91 trials = 0.5384615384615384
estimated Pr[Y = 1] after 92 trials = 0.5434782608695652
estimated Pr[Y = 1] after 93 trials = 0.5483870967741935
estimated Pr[Y = 1] after 94 trials = 0.5425531914893617
estimated Pr[Y = 1] after 95 trials = 0.5473684210526316
estimated Pr[Y = 1] after 96 trials = 0.5520833333333334
estimated Pr[Y = 1] after 97 trials = 0.5463917525773195
estimated Pr[Y = 1] after 98 trials = 0.5510204081632653
estimated Pr[Y = 1] after 99 trials = 0.5454545454545454
estimated Pr[Y = 1] after 100 trials = 0.54

Recall that the expression (y[m] == 1) evaluates to 1 if the condition holds and 0 otherwise. The result of running the program is that estimate[m] will hold the estimate \(\mbox{Pr}[Y = 1]\) after \(m\) simulations. We can then plot the estimates as a function of the number of draws using a line plot to display the trend.

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import random

np.random.seed(0)

M = 100000
Ms = []
y_sim = [random.randint(0,1) for i in range(M)]
hat_E_Y = []
Ms = []

for i in range(51):
    Ms.append(min(M, 10 ** (i / 10)))
    hat_E_Y.append(np.mean(y_sim[:int(Ms[i])]))

df = pd.DataFrame({'M': Ms, 'hat_E_Y': hat_E_Y})
plot = plt.scatter(df['M'], df['hat_E_Y'])
plt.axhline(y=0.5, color='red')
plt.xscale('log')
plt.xlabel('simulation draws')
plt.ylabel('estimated Pr[Y = 1]')
plt.xlim((1, 100000))
plt.ylim((0, 1))
plt.xticks([1, 50000], ["1", "50,000"])#, "100,000"])
plt.show()

The linear scale of the previous plot obscures the behavior of the estimates. Consider instead a plot with the \(x\)-axis on the logarithmic scale.

Monte Carlo estimate of probability that a coin lands head as a function of the number of simulation draws. The line at 0.5 marks the true probability being estimated. The log-scaled \(x\)-axis makes the early rate of convergence more evident. Plotted on a linear scale, it is clear how quickly the estimates converge to roughly the right value.

import matplotlib.pyplot as plt
import pandas as pd

df = pd.DataFrame({'M': Ms, 'hat_E_Y': hat_E_Y})

plt.plot(df['M'], df['hat_E_Y'], linestyle='-', marker='o')
plt.axhline(y=0.5, color='red')
plt.xscale('log')
plt.xlabel('simulation draws')
plt.ylabel('estimated Pr[Y = 1]')
plt.xlim((1, 50000))
plt.ylim((0, 1))
plt.xticks([1, 50000], ["1", "50,000"])#, "100,000"])
plt.yticks([0, 0.25, 0.5, 0.75, 1.0])
plt.show()

With a log-scaled \(x\)-axis, the values between 1 and 10 are plotted with the same width as the values between \(10\,000\) and \(100\,000\); both take up 20% of the width of the plot. On the linear scale, the values between 1 and 10 take up only \(\frac{10}{100\,000}\), or 0.01%, of the plot, whereas values between \(10\,000\) and \(100\,000\) take up 90% of the plot.

Plotting the progression of multiple simulations demonstrates the trend in errors.

Each of the one hundred grey lines represents the ratio of heads observed in a sequence of coin flips, the size of which indicated on the \(x\)-axis. The line at 0.5 indicates the probability a coin lands heads in a fair coin toss. The convergence of the ratio of heads to 0.5 in all of the sequences is clearly visible.

 
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns

np.random.seed(0)

M_max = 10000
J = 100
I = 47
N = I * J
df2 = pd.DataFrame({'r': [np.nan]*N, 'M': [np.nan]*N, 'hat_E_Y': [np.nan]*N})
pos = 0

for j in range(1, J+1):
    y_sim = np.random.binomial(1, 0.5, size=M_max)
    for i in range(4, 51):
        M = max(100, min(M_max, int((10**(1/10))**i)))
        hat_E_Y = np.mean(y_sim[:M])
        df2.loc[pos, :] = [j, M, hat_E_Y]
        pos += 1

pr_Y_eq_1_plot = sns.lineplot(data=df2, x='M', y='hat_E_Y', hue='r', alpha=0.15, linewidth=2)
pr_Y_eq_1_plot.axhline(y=0.5, color='red', linewidth=2)
pr_Y_eq_1_plot.set(xscale='log', xlabel='simulation draws', ylabel='estimated Pr[Y = 1]',
                    xlim=(100, 10000), ylim=(0.375, 0.625), xticks=[1000, 10000],
                    xticklabels=["1,000", "10,000"])
sns.set_theme(style='ticks')

plt.show()

Continuing where the previous plot left off, each of the one hundred grey lines represents the ratio of heads observed in a sequence of coin flips. The values on the \(x\) axis is one hundred times larger than in the previous plot, and the scale of the \(y\)-axis is one tenth as large. The trend in error reduction appears the same at the larger scale.

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns

np.random.seed(0)

M_max = 1000000
J = 100
I = 47
N = I * J
df2 = pd.DataFrame({'r': [np.nan]*N, 'M': [np.nan]*N, 'hat_E_Y': [np.nan]*N})
pos = 0

for j in range(1, J+1):
    y_sim = np.random.binomial(1, 0.5, size=M_max)
    for i in range(4, 61):
        M = max(100, min(M_max, int((10**(1/10))**i)))
        hat_E_Y = np.mean(y_sim[:M])
        df2.loc[pos, :] = [j, M, hat_E_Y]
        pos += 1

pr_Y_eq_1_plot = sns.lineplot(data=df2, x='M', y='hat_E_Y', hue='r', alpha=0.15)
pr_Y_eq_1_plot.axhline(y=0.5, color='red')
pr_Y_eq_1_plot.set(xscale='log', xlabel='simulation draws', ylabel='estimated Pr[Y = 1]',
                    xlim=(1e4, 1e6), ylim=(0.485, 0.515), xticks=[1e4, 1e5, 1e6],
                    xticklabels=["10,000", "100,000", "1,000,000"])
sns.set_theme(style='ticks')

plt.show()

The law of large numbers^[Which technically comes in a strong and weak form.] says roughly that as the number of simulated values grows, the average will converge to the expected value. In this case, our estimate of \(\mbox{Pr}[Y = 1]\) can be seen to converge to the true value of 0.5 as the number of simulations \(M\) increases. Because the quantities involved are probabilistic, the exact specification is a little more subtle than the \(\epsilon\)-\(\delta\) proofs in calculus.

Simulation notation

We will use parenthesized superscripts to pick out the elements of a sequence of simulations. For example,

[y^{(1)}, y^{(2)}, \ldots, y^{(M)}]

will be used for \(M\) simulations of a single random variable \(Y\).^[Each \(y^{(m)}\) is a possible realization of \(Y\), which is why they are written using lowercase.] It’s important to keep in mind that this is \(M\) simulations of a single random variable, not a single simulation of \(M\) different random variables.

Before we get going, we’ll need to introduce indicator function notation. For example, we write

[\mathrm{I}[y^{(m)} = 1]

\begin{cases} 1 & \mbox{if} \ y^{(m)} = 1 \[4pt] 0 & \mbox{otherwise} \end{cases}]

The indicator function maps a condition, such as \(y^{(m)} = 1\) into the value 1 if the condition is true and 0 if it is false.^[Square bracket notation is used for functions when the argument is itself a function. For example, we write \(\mbox{Pr}[Y > 0]\) because \(Y\) is a random variable, which is modeled as a function. We also write \(\mathrm{I}[x^2 + y^2 = 1]\) because the standard bound variables \(x\) and \(y\) are functions from contexts defining variable values.]

Now we can write out the formula for our estimate of \(\mbox{Pr}[Y = 1]\) after \(M\) draws,

[\begin{array}{rcl} \mbox{Pr}[Y = 1] & \approx & \frac{\displaystyle \mathrm{I}[y^{(1)} = 1] \ + \ \mathrm{I}[y^{(2)} = 1] \ + \ \cdots \ + \ \mathrm{I}[y^{(M)} = 1]} {\displaystyle M} \end{array}]

That is, our estimate is the proportion of the simulated values which take on the value 1. It quickly becomes tedious to write out sequences, so we will use standard summation notation, where we write

[\mathrm{I}!\left[y^{(1)} = 1] + \mathrm{I}[y^{(2)} = 1] + \cdots + \mathrm{I}[y^{(M)} = 1\right] \ =
\sum_{m=1}^M \mathrm{I}[y^{(m)} = 1]]

Thus we can write our simulation-based estimate of the probability that a fair coin flip lands heads as^[In general, the way to estimate an event probability \(\phi(Y)\) where \(\phi\) defines some condition, given simulations \(y^{(1)}, \ldots, y^{(M)}\) of \(Y\), is as \(\mbox{Pr}[\phi(Y)] = \frac{1}{M} \sum_{m = 1}^M \mathrm{I}[\phi(y^{(m)})].\)]

[\mbox{Pr}[Y = 1] \approx \frac{1}{M} \, \sum_{m=1}^M \mathrm{I}[y^{(m)} = 1]]

The form \(\frac{1}{M} \sum_{m=1}^M\) will recur repeatedly in simulation — it just says to average over values indexed by \(m \in 1:M\).^[We are finally in a position to state the strong law of large numbers as the event probability of a limit, \(\mbox{Pr}\!\left\lbrack \lim_{M \rightarrow \infty} \frac{1}{M} \sum_{m = 1}^M \mathrm{I}\!\left\lbrack y^{(m)} = 1 \right\rbrack \ = \ 0.5 \right\rbrack,\) where each \(y^{(m)}\) is a separate fair coin toss. ]

Central limit theorem

The law of large numbers tells us that with more simulations, our estimates become more and more accurate. But they do not tell us how quickly we can expect that convergence to proceed. The central limit theorem provides the convergence rate.

First, we have to be careful about what we’re defining. First, we define the error for an estimate as the difference from the true value,

[\left( \frac{1}{M} \sum_{m=1} \mathrm{I}[y^{(m)} = 1] \right) - 0.5]

The absolute error is just the absolute value^[In general, the absolute value function applied to a real number \(x\) is written as \(|x|\) and defined to be \(x\) if \(x\) is non-negative and \(-x\) if \(x\) is negative.] of this,

[\left| \, \left( \frac{1}{M} \sum_{m=1} \mathrm{I}[y^{(m)} = 1] \right) - 0.5 \, \right|]

The absolute error of the simulation-based probability estimate as a function of the number of simulation draws. One hundred sequences of one million flips are shown.

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from scipy.stats import binom

np.random.seed(1234)

M_max = 1000000
J = 100
I = 47
N = I * J

df2 = pd.DataFrame({'r': [np.nan]*N, 'M': [np.nan]*N, 'err_hat_E_Y': [np.nan]*N})
pos = 0

for j in range(1, J+1):
    y_sim = binom.rvs(1, 0.5, size=M_max)
    for i in range(4, 61):
        M = max(100, min(M_max, int((10**(1/10))**i)))
        err_hat_E_Y = abs(np.mean(y_sim[:M]) - 0.5)
        df2.loc[pos, :] = [j, M, err_hat_E_Y]
        pos += 1

abs_err_plot = sns.lineplot(data=df2, x='M', y='err_hat_E_Y', hue='r', alpha=0.15)
abs_err_plot.axhline(y=0.5, color='red')
abs_err_plot.set(xscale='log', xlabel='simulation draws', ylabel='absolute error',
                  xlim=(10000, 1000000), ylim=(0, 0.015),
                  xticks=[10000, 100000, 1000000], xticklabels=["10,000", "100,000", "1,000,000"])
sns.set_theme(style='ticks')

plt.show()

import numpy as np
from scipy.stats import norm

np.random.seed(1234)

M_max = int(1e6)
J = 300
Ms = 10 ** (np.arange(6, 13, 0.5) / 2)
N = len(Ms)
ys = np.empty((N, J))

for j in range(J):
    z = np.random.binomial(1, 0.5, M_max)
    for n in range(N):
        ys[n, j] = abs(np.mean(z[:int(Ms[n])]) - 0.5)

mean = np.empty(N)
sd = np.empty(N)
quant_68 = np.empty(N)
quant_95 = np.empty(N)

sixty_eight_pct = norm.cdf(1) - norm.cdf(-1)
ninety_five_pct = norm.cdf(2) - norm.cdf(-2)

for n in range(N):
    mean[n] = np.mean(ys[n, :])
    sd[n] = np.std(ys[n, :])
    quant_68[n] = np.quantile(ys[n, :], sixty_eight_pct)
    quant_95[n] = np.quantile(ys[n, :], ninety_five_pct)

Plotting both the number of simulations and the absolute error on the log scale reveals the rate at which the error decreases with more draws.

Absolute error versus number of simulation draws for 100 simulated sequences of \(M = 1,000,000\) draws. The blue line is at the 68 percent quantile and the red line at the 95 percent quantile of these draws. The relationship between the log number of draws and log error is revealed to be linear.

import numpy as np
from scipy.stats import binom, norm
import pandas as pd
import matplotlib.pyplot as plt
import math

np.random.seed(1234)

M_max = int(1e6)
J = 300
Ms = np.power(10, np.arange(6, 13, 0.5) / 2)
N = len(Ms)
ys = np.empty((N, J))
for j in range(J):
    z = np.random.binomial(1, 0.5, size=M_max)
    for n in range(N):
        ys[n, j] = np.abs(np.mean(z[:int(Ms[n])]) - 0.5)

mean = np.empty(N)
sd = np.empty(N)
quant_68 = np.empty(N)
quant_95 = np.empty(N)
sixty_eight_pct = norm.cdf(1) - norm.cdf(-1)
ninety_five_pct = norm.cdf(2) - norm.cdf(-2)
for n in range(N):
    mean[n] = np.mean(ys[n, :])
    sd[n] = np.std(ys[n, :], ddof=1)
    quant_68[n] = np.quantile(ys[n, :], q=sixty_eight_pct)
    quant_95[n] = np.quantile(ys[n, :], q=ninety_five_pct)

fudge = 1e-6
x = 0.5 * np.arange(1, 101)
log_abs_err_plot = plt.plot()
df2 = pd.DataFrame({'M': Ms, 'err_hat_E_Y': mean - 0.5})
for r in range(J):
    plt.plot(Ms, ys[:, r], alpha=0.15, color='gray')
plt.plot(Ms, quant_68, color='blue', linewidth=0.5)
plt.plot(Ms, quant_95, color='red', linewidth=0.5)
plt.xscale('log')
plt.xlim(0.9 * 1e3, 1.1 * 1e6)
plt.xticks([10**3, 10**4, 10**5, 10**6], ['1,000', '10,000', '100,000', '1,000,000'])
plt.ylim(0.8e-6, 2e-5)
plt.yticks([0.00005, 0.0005, 0.005, 0.05], ['0.00001', '0.0001', '0.001', '0.01'])
plt.xlabel('simulation draws')
plt.ylabel('log absolute error + 1e-5')
plt.show()

We can read two points \((x_1, y_1)\) and \((x_2, y_2)\) off of the graph for \(x_1 = 10\,000\) and \(x_2 = 100\,000\) as

print("x[1], y[1] = %7.f, %6.5f\nx[2], y[2] = %7.f, %6.5f" % (
        Ms[3], quant_68[3],
        Ms[5], quant_68[5]))

x[1], y[1] =    5623, 0.00685
x[2], y[2] =   17783, 0.00411

which gives us the following values on the log scale

import math

print(f"log x[1], log y[1] = {math.log(Ms[3]):5.2f}, {math.log(quant_68[3]):4.2f}\nlog x[2], log y[2] = {math.log(Ms[5]):5.2f}, {math.log(quant_68[5]):4.2f}")

log x[1], log y[1] =  8.63, -4.98
log x[2], log y[2] =  9.79, -5.49

Using the log scale values, the estimated slope of the reduction in quantile bounds is

print(f"estimated slope\n(log y[2] - log y[1])\n / (log x[2] - log x[1])  =  {(math.log(quant_68[5]) - math.log(quant_68[3])) / (math.log(Ms[5]) - math.log(Ms[3])) :3.2f}")

estimated slope
(log y[2] - log y[1])
 / (log x[2] - log x[1])  =  -0.44

If we let \(\epsilon_M\) be the value of one of the quantile lines at \(M\) simulation draws, the linear relationship plotted in the figures have the form^[A line through the points \((x_1, y_1)\) and \((x_2, y_2)\) has \(\mbox{slope} = \frac{y_2 - y_1}{x_2 - x_1}.\)]

[\log \epsilon_M = -\frac{1}{2} \, \log M + \mbox{const}.]

When writing “const” in a mathematical expression, the presumption is that it refers to a constant that does not depend on the free variables of interest (here, \(M\), the number of simulation draws). Ignoring constants lets us focus on the order of the dependency. The red line and blue line have the same slope, but different constants.

Seeing how this plays out on the linear scale requires exponentiating both sides of the equation and reducing,

[\begin{array}{rcl} \exp(\log \epsilon_M) & = & \exp( -\frac{1}{2} \, \log M + \mbox{const} ) \[6pt] % \epsilon_M & = & \exp( -\frac{1}{2} \, \log M ) \times \exp(\mbox{const}) % \[6pt] \epsilon_M & = & \exp( \log M )^{-\frac{1}{2}} \times \exp(\mbox{const}) \[6pt] \epsilon_M & = & \exp(\mbox{const}) \times M^{-\frac{1}{2}} \end{array}]

Dropping the constant, this relationship between the expected absolute error \(\epsilon_M\) after \(M\) simulation draws may be succinctly summarized using proportionality notation,^[In general, we write \(f(x) \propto g(x)\) if there is a positive constant \(c\) that does not depend on \(x\) such that \(f(x) = c \times g(x).\) For example, \(3x^2 \propto 9x^2,\) with \(c = \frac{1}{3}\).]

[\displaystyle \epsilon_M \ \propto
\frac{\displaystyle 1}{\displaystyle \sqrt{M}}.]

This is a fundamental result in statistics derived from the central limit theorem. The central limit theorem governs the accuracy of almost all simulation-based estimates. We will return to a proper formulation when we have the scaffolding in place to deal with the pesky constant term.

In practice, what does this mean? It means that if we want to get an extra decimal place of accuracy in our estimates, we need one hundred (100) times as many draws. For example, the plot shows error rates bounded by roughly 0.01 with \(10\,000\) draws, yielding estimates that are very likely to be within \((0.49, 0.51).\) To reduce that likely error bound to 0.001, that is, ensuring estimates are very likely in \((0.0499, 0.501),\) requires 100 times as many draws (i.e., a whopping \(1\,000\,000\) draws).^[For some perspective, \(10\,000\) is the number of at bats in an entire twenty-year career for a baseball player, the number of field goal attempts in an entire career of most basketball players, and the size of a very large disease study or meta-analysis in epidemiology.]

Usually, one hundred draws will provide a good enough estimate of most quantities of interest in applied statistics. The variability of the estimate based on a single draw depends on the variability of the quantity being estimated. One hundred draws reduces the expected estimation bound to one tenth of the variability of a single draw. Reducing that variability to one hundredth of the variability of a single draw would require ten thousand draws. In most applications, the extra estimation accuracy is not worth the extra computation.

Exercise

1. A motorist encounters four consecutive traffic lights, each equally likely to be red or green. Let Z be the number of green lights passed by the motorist before being stopped by a red light. What is the probability distribution of Z?

Solution

Thinking of this problem as 4 buckets each with 2 possibilities paves a clear way to the solution of the problem. There are \(2^4 = 16\) possibilities. There is only one way for them to all be green (\(\xi = 4\)), and again only one way for 3 greens followed by one red (\(\xi = 3\)). Once you have get to \(\xi = 2\), then the last light can be either green or red giving two possibilities then at \(\xi = 1\) you have two lights that can be either red or green giving 4 possibilities. Clearly, then, there are 8 for the case of \(\xi = 0\). \(P(\xi) = \left\{ \begin{array}{rl} \frac{1}{2}, \xi=0 \\ \frac{1}{4}, \xi=1 \\ \frac{1}{8}, \xi=2 \\ \frac{1}{16}, \xi=3 \\ \frac{1}{16}, \xi=4 \end{array} \right. \label{answer4.1}\)

2. Give an example of two distinct random variables with the same distribution function.

Solution

To summarize what we want here: \(\begin{aligned} \xi_1 \neq \xi_2 \\ \Phi_{\xi_1}(x) = \Phi_{\xi_2}(x) \\ \int_{-\infty}^x p_{\xi_1}(x')dx' = \int_{-\infty}^x p_{\xi_2}(x')dx'\end{aligned}\) Since the question does not rule out the possibility that the probability distributions are the same, we can just say that \(\xi_1\) is the outcome of a coin-flip experiment and that \(\xi_2\) is the outcome of the first spin measurement of an EPR experiment.

4. A random variable \(\xi\) has probability density \(p_{\xi}(x) = \frac{a}{x^2 + 1} (-\infty < x <\infty)\). Find

a) The constant a;
b) The distribution function of \(\xi\);
c) The probability \(P{-1 < \xi < 1}\).

Solution

The distribution is clearly not normalized so the first step will be to normalise it. \(\int^{\infty}_{-\infty} \frac{a}{x^2+1}dx = 1 = a \pi\)
\(a = \frac{1}{\pi}\)
Just by definition:
\(\Phi(x) = \int_{x}^{\infty} \frac{1}{\pi(x'^2+1)}dx' = \frac{\arctan x'}{\pi} \bigg|_{-\infty}^x\)
\(= \frac{\arctan x}{\pi}+\frac{1}{2}\) and last but not least:
\(P(\{-1 \leq x \leq 1\}) = \int_{-1}^1 \frac{1}{\pi(x'^2+1)} dx' = \frac{1}{2}\)

###

Solution

Once again we start by normalizing \(1= \int^{\infty}_{0} a x^2 e^{-k x} dx = -\left \frac{e^{-k x} \left(2+2 k x+k^2 x^2\right)}{k^3} \right]^{\infty}_{0}= \frac{2 a}{k^3}\) \(a = \frac{k^3}{2} \label{answer4.5a}\) \(\Phi(x) = \int^{x}_{0} \frac{k^3}{2}x'^2 e^{-k x'}dx' = 1 - \frac{e^{-k x} \left(2+2 k x+k^2 x^2\right)}{2} \label{answer4.5b}\) \(P( \{ 0 \leq x \leq \frac{1}{k} \} = \int_{0}^{\frac{1}{k}} \frac{k^3}{2}x^2 e^{-k x}dx = \frac{2e - 5}{2 e} \label{answer4.4c}\)

###

Solution

\(\Phi(\infty) = 1 \Rightarrow 1 = a + \frac{b \pi}{2}\) \(\Phi(-\infty) = 0 \Rightarrow 0 = a - \frac{b \pi}{2}\) Which you can solve to get \(\begin{aligned} b = \frac{1}{\pi} \\ a = \frac{1}{2}\end{aligned}\) Since \(\Phi\) is the integral of \(p\), we can take the derivative of it and then make sure it’s normalized. \(p(x) = \frac{d \Phi}{dx} = \frac{1}{2 \pi \left(1+\frac{x^2}{4}\right)}\) And the normalization is indeed still correct!

###

Solution

The area of the table is simply \(\pi R^2\) and the are of each of the smaller circles are \(\pi r^2\). The ratio of the sum of the area of the two circles to the total table area will be the chance that one of the circles gets hit: \(p = \frac{2r^2}{R^2}\).

###

Solution

Just like in example 2 in the book, this problem will go by very well if we draw a picture indicating the given criteria:

{width=”\textwidth”}

Where does this come from? \(\begin{aligned} x_1 + x_2 \leq 1 \\ x_2 \leq 1 - x_1\end{aligned}\) which is the straight line. \(\begin{aligned} x_1x_2 \leq \frac{2}{9} \\ x_2 \leq \frac{2}{9 x_1}\end{aligned}\) A little bit of algebra shows that these two lines intersect at \(\frac{1}{3}\) and \(\frac{2}{3}\) so the area underneath the straight line but not above the curved line is \(\begin{aligned} A =\int_0^{\frac{1}{3}}(1-x)dx + \int_{\frac{1}{3}}^{\frac{2}{3}}\frac{2}{9x}dx + \int_{\frac{2}{3}}^{1}(1-x)dx \\ = \frac{5}{18} + \int_{\frac{1}{3}}^{\frac{2}{3}}\frac{2}{9x}dx + \frac{1}{18} \\ = \frac{1}{3} + \int_{\frac{1}{3}}^{\frac{2}{3}}\frac{2}{9x}dx \\ = \frac{1}{3 } + \frac{2 \ln 2}{9} \approx 0.487366\end{aligned}\) And the answer is properly normalized since the initial probability distributions were unity (the box length is only 1).

###

Solution

As we showed in example number 4, \(p_{\eta}\) is the convolution of \(p_{\xi_1}\) and \(p_{\xi_2}\). \(\begin{aligned} p_{\eta}(y) = \int_{-\infty}^{\infty} p_{\xi_1}(y-x)p_{\xi_2}(x)dx\end{aligned}\) I think the integration will look more logical if we stick in Heaviside step-functions. \(\begin{aligned} p_{\eta}(y) = \frac{1}{6} \int_{-\infty}^{\infty}e^{-\frac{y-x}{3}}e^{-\frac{x}{3}}H(y-x)H(x)dx\end{aligned}\) Clearly \(x\) has to be greater than zero and \(y\) must be greater than x, leading to the following limits of integration. \(\begin{aligned} p_{\eta}(y) = \frac{1}{6} \int_{0}^{y}e^{-\frac{y-x}{3}}e^{-\frac{x}{3}}dx \\ = e^{-\frac{y}{3}}\left(1- e^{-\frac{y}{6}} \right)\end{aligned}\) Only when y is greater than zero!

###

Solution

Due to the magic of addition, finding the probability distribution of \(\xi_1 + \xi_2 + \xi_3\) is no different than the probability distribution of \((\xi_1 + \xi_2) + \xi_3\) but we already know what the probability distribtution of the parenthesized quantity is: \(p_{\xi_1 + \xi_2}(y) = \int_{-\infty}^{\infty} p_{\xi_1}(y-x)p_{\xi_2}(x)dx\) Therefore the total combination is \(p_{\xi_1 + \xi_2+ \xi_3}(z) = \int_{-\infty}^{\infty} p_{\xi_1 + \xi_2}(z-y)p_{\xi_3}(y)dx\) \(p_{\xi_1 + \xi_2+ \xi_3}(z) = \int_{-\infty}^{\infty} \left[ \int_{-\infty}^{\infty} p_{\xi_1}(z-y-x)p_{\xi_2}(x)dx \right]p_{\xi_3}(y)dy\) \(p_{\xi_1 + \xi_2+ \xi_3}(z) = \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} p_{\xi_1}(z-y-x)p_{\xi_2}(x)p_{\xi_3}(y)dx dy\) which is just a triple convolution.

###

Solution

\(p_{\xi}(n) = \frac{1}{3^n}\) therefore \(\textbf{E}\xi = \sum_{n=1}^{\infty}\frac{n}{3^n} = \frac{3}{4} \label{answer4.11}\)

###

Since finding a ball in one urn versus the other urn has nothing to do with each other, the events are clearly independent so we can multiply probabilities. The probability of finding a white ball at any given urn is \(p_w = \frac{w}{w+b}\) If you find a white ball on the \(n^{th}\) try then that means you found \(n-1\) black balls before you got to the white ball. \(p_w(n) = \frac{b^{n-1} w}{(w+b)^n}\) \(\textbf{E}n = \sum_{i=1}^{\infty}n p_w(n) = \sum_{i=1}^{\infty}n \frac{b^{n-1} w}{(w+b)^n} = \frac{b+w}{w}\) Which is the total number of balls drawn, subtract one to get the average number of black balls drawn: \(m=\frac{b}{w}\) Now for the variance. to start with we need the average of the square of the random variable. \(\textbf{E}n^2 = \sum_{i=1}^{\infty}n^2 p_w(n) = \sum_{i=1}^{\infty}n^2 \frac{b^{n-1} w}{(w+b)^n} = \frac{(b+w) (2 b+w)}{w^2}\) \(\textbf{D}n = \textbf{E}n^2 - (\textbf{E}n)^2 = \frac{(b+w) (2 b+w)}{w^2} - \left( \frac{b+w}{w} \right)^2 = \frac{b^2+wb}{w^2}\) A note that we don’t need to subtract anything for the variance since shifting a distribution over does not affect its variance: just its average.

###

Solution

\(\textbf{E}\xi = \int_{-\infty}^{\infty} x\frac{1}{2}e^{-|x|} = 0\) because the function is even about \(x=0\). \(\textbf{E}\xi^2 = \int_{-\infty}^{\infty} x^2\frac{1}{2}e^{-|x|} = 2\) Therefore: \(\textbf{D}\xi = \textbf{E}\xi^2 - (\textbf{E}\xi)^2 = 2\)

###

Solution

\(\textbf{E}x = \int_{a-b}^{a+b} \frac{xdx}{2b} = a\) \(\textbf{E}x^2 = \int_{a-b}^{a+b} \frac{x^2dx}{2b} = a^2+\frac{b^2}{3}\) Therefore: \(\textbf{D}x = \frac{b^2}{3}\)

###

Solution

If the distribution function is \(\Phi_{\xi}(x) = a + b \arcsin x, |x| \leq 1\) then it must fulfill the proper boundary conditions as specified by both the problem and the definition of a distribution function. \(\Phi_{\xi}(-1)= 0 = a - b \frac{\pi}{2}\) \(\Phi_{\xi}(1) = 1 = a + b \frac{\pi}{2}\) Some easy algebra gets you \(\begin{aligned} \Phi_{\xi}(x) = \frac{1}{2} + \frac{1}{\pi} \arcsin x\end{aligned}\) Therefore: \(p_{\xi}(x) = \frac{1}{\pi\sqrt{1-x^2}}\) \(\begin{aligned} \textbf{E}x = \int_{-1}^{1} \frac{x dx}{\pi\sqrt{1-x^2}} = 0 \\ \textbf{D}x = \int_{-1}^{1} \frac{x^2 dx}{\pi\sqrt{1-x^2}} = \frac{1}{2}\end{aligned}\)

###

Solution

Each side of the die has the same probability, \(\frac{1}{6}\) \(\textbf{E}x = \sum_{i=1}^6 \frac{i}{6} = \frac{7}{2}\) \(\textbf{E}x^2 = \sum_{i=1}^6 \frac{i^2}{6} = \frac{91}{6}\) Therefore: \(\textbf{D}x = \frac{91}{6} - \left( \frac{7}{2} \right)^2 = \frac{91}{6} - \frac{49}{4} = \frac{35}{12}\)

###

Solution

This problem may seem difficult until you realize that being \(\pm \frac{5}{2}\) away from the mean means you’re at either 1 or 6, meaning that there’s a unity probability of being that far from the mean. Chebyshev’s inequality, however, would have us believe that \(P \{ |x + \textbf{E}x| > \frac{5}{2} \} \leq \frac{4}{25} \frac{91}{6} = \frac{182}{75} = 2.42667\) which is not only far off from the actual answer, it’s also unphysical to have probabilities greater than 1!

###

Solution

We want to consider the probability distribution of \(\xi\) by way of \(\eta\) \(\eta = e^{\frac{a\xi}{2}}\) We know from Chebyshev’s identity:

\(P\{ \eta > \epsilon(\eta) \} \leq \frac{\textbf{E}\eta^2}{\epsilon(\eta)^2}\) Let \(\epsilon\) be the error in \(\xi\) we’re looking for. \(P\{ \xi > \epsilon \} \leq \frac{\textbf{E}(e^{\frac{a\xi}{2}})^2}{(e^{\frac{a\epsilon}{2}})^2}\) \(P\{ \xi > \epsilon \} \leq \frac{\textbf{E}e^{a\xi}}{e^{a\epsilon}}\)

###

Solution

First some initial calculations: \(\begin{aligned} \textbf{E}\xi = \frac{1}{4} \left( -2+-1+1+2 \right) = 0 \\ \textbf{E}\xi^2 = \frac{1}{4} \left( (-2)^2+(-1)^2+1^2+2^2 \right) = \frac{10}{4} = \frac{5}{2} = \textbf{E}\eta \\ \textbf{E}\xi^4 = \frac{1}{4} \left( (-2)^4+(-1)^4+1^4+2^4 \right) = \frac{34}{4} = \frac{17}{2} = \textbf{E}\eta^2\end{aligned}\) Now, we know that \(\begin{aligned} r = \frac{\textbf{E} \left[ (\xi - \textbf{E}\xi)(\eta - \textbf{E}\eta) \right]}{(\textbf{E}\xi^2 - (\textbf{E}\xi)^2)(\textbf{E}\eta^2 - (\textbf{E}\eta)^2)}\end{aligned}\) The denominator \(D\) is the easiest: \(D = \left( \frac{5}{2} - 0 \right)\left( \frac{17}{2} - \frac{5^2}{2^2} \right) = \frac{45}{8}\) Now for the numerator: \(\begin{aligned} \textbf{E} \left[ (\xi - \textbf{E}\xi)(\eta - \textbf{E}\eta) \right] \\ = \frac{1}{16} \left[ \sum_{\xi, \eta} (\xi - \textbf{E}\xi)(\eta - \textbf{E}\eta) \right]\end{aligned}\) If we look at the set we’ll be summing over, we have \({\xi - \textbf{E}\xi} = -2, -1, 1, 2\) and \({\eta - \textbf{E}\eta} = \frac{3}{2}, -\frac{3}{2}, -\frac{3}{2}, \frac{3}{2}\) and clearly, since we have to sum over all possible products, since both distributions are symmetric about 0, the sum will go to zero. \(r=0\)

###

Solution

First some initial calculations: \(\begin{aligned} \textbf{E}x_1 = \int_0^{\frac{\pi}{2}} x_1 \sin x_1 \sin x_2 dx_1dx_2 = 1 \\ \textbf{E}x_2 = 1 \\ \textbf{E}x_1^2 = \int_0^{\frac{\pi}{2}} x_1^2 \sin x_1 \sin x_2 dx_1dx_2 = (\pi - 2) \\ \textbf{E}x_2^2 = (\pi - 2) \\\end{aligned}\) \(\begin{aligned} \textbf{E} \left[ (x_1 - \textbf{E}x_1)(x_2 - \textbf{E}x_2) \right] \\ = \int_0^{\frac{\pi}{2}} \int_0^{\frac{\pi}{2}} \left[(x_1 - 1)(x_2 - 1)\right]\sin x_1 \sin x_2 dx_1dx_2 = 0 \\ \sigma_1 \sigma_2 = \sqrt{(\pi - 2)- 1}\sqrt{(\pi - 2) - 1} = (\pi - 3) \\ r= \frac{0}{(\pi - 3)} = 0\end{aligned}\)

###

Solution

First some initial calculations: \(\begin{aligned} \textbf{E}x_1 = \frac{1}{2}\int_0^{\frac{\pi}{2}} x_1 \sin (x_1 + x_2) dx_1dx_2 = \frac{\pi}{4} \\ \textbf{E}x_2 = \frac{\pi}{4} \\ \textbf{E}x_1^2 = \frac{1}{2}\int_0^{\frac{\pi}{2}} x_1^2 \sin (x_1 + x_2) dx_1dx_2 = -2+\frac{\pi}{2} +\frac{\pi ^2}{8} \\ \textbf{E}x_2^2 = -2+\frac{\pi}{2} +\frac{\pi ^2}{8} \\\end{aligned}\) \(\begin{aligned} \textbf{E} \left[ (x_1 - \textbf{E}x_1)(x_2 - \textbf{E}x_2) \right] \\ = \frac{1}{2}\int_0^{\frac{\pi}{2}} \int_0^{\frac{\pi}{2}} \sin (x_1 + x_2) \left[(x_1 - \frac{\pi}{2})(x_2 - \frac{\pi}{2})\right]dx_1dx_2 = -\frac{1}{16} (\pi -4)^2 \\ \sigma_1\sigma_2 = -2+\pi +\frac{\pi ^2}{2} +\frac{\pi^2}{8} - \frac{\pi^2}{16} = -2+\pi +\frac{\pi ^2}{2} +\frac{\pi^2}{16} \\ r = \frac{-\frac{1}{16} (\pi -4)^2}{-2+\pi +\frac{\pi ^2}{2} +\frac{\pi^2}{16}}\end{aligned}\)

###

Solution

\(\begin{aligned} E\xi = \int_{-\infty}^{\infty} \frac{x}{\pi(1+x^2)}dx = \left \frac{\log \left(x^2+1\right)}{2 \pi } \right]_{-\infty}^{\infty} \rightarrow \infty \\ E\xi^2 = \int_{-\infty}^{\infty} \frac{x^2}{\pi(1+x^2)}dx = \left \frac{x}{\pi }-\frac{\tan ^{-1}(x)}{\pi } \right]_{-\infty}^{\infty} \rightarrow \infty\end{aligned}\) Neither integral actually converges so we cannot define averages or dispersions for this distribution.

Key Points

Multiple Random Variables and Probability Functions

Overview

Teaching: min
Exercises: min

Questions

Objectives

Multiple Random Variables and Probability Functions

Multiple random variables

Random variables do not exist in isolation. We started with a single random variable \(Y\) representing the result of a single, specific coin flip. Suppose we fairly flip the coin three times? Then we can have random variables \(Y_1, Y_2, Y_3\) representing the results of each of the flips. We can assume each flip is independent in that it doesn’t depend on the result of other flips. Each of these variables \(Y_n\) for \(n \in 1:3\) has \(\mbox{Pr}[Y_n = 1] = 0.5\) and \(\mbox{Pr}[Y_n = 0] = 0.5\).

We can combine multiple random variables using arithmetic operations. We have already seen comparison operators in writing the event \(Y = 1\). If \(Y_1, \ldots, Y_{10}\) are random variables representing ten coin flips, then we can define their sum as

[Z = Y_1 + Y_2 + Y_3]

We can simulate values of \(Z\) by simulating values of \(Y_1, Y_2, Y_3\) and adding them.

from random import randint

y1 = randint(0, 1)
y2 = randint(0, 1)
y3 = randint(0, 1)
z = y1 + y2 + y3
print("z =", z)

z = 3

It is easier and less error prone to collapse similar values into arrays and operate on the arrays collectively, or with loops if necessary.

y = []
for n in range(1, 4):
    y.append(randint(0, 1))
z = sum(y)
print("z =", z)

Running this program a few times we get

z = 2
z = 3
z = 1
z = 3
z = 1

We can use simulation to evaluate the probability of an outcome that combines multiple random variables. For example, to evaluate \(\mbox{Pr}[Z = 2]\), we run the simulation many times and count the proportion of results that are two.^[The sum is calculated using notation \(sum(y[m, ])\), which is defined to be \(sum(y[m, ]) = y[m, 1] + ... + y[m, N],\) where N is the number of entries in row m of the variable y.]

M = 10000  # Define the value of M
y = [[randint(0, 1) for n in range(3)] for m in range(M)]  # Generate random integers and store them in a 2D list y
z = [sum(row) for row in y]  # Calculate the sum of each row in y and store them in a list z
Pr_is_two = sum(z == 2) / M  # Calculate the proportion of elements in z that are equal to 2
print("Pr_is_two =", Pr_is_two)

Pr_is_two = 0.3692

As in our other probability estimates, we simulate the variable of interest \(Z\) a total of \(M\) times, yielding \(z^{(1)}, \ldots, z^{(m)}\). Here, that requires simulating \(y_1^{(m)}, y_2^{(m)}, y_3^{(m)}\) and adding them for each \(z^{(m)}\). We then just count the number of times Z is simulated to be equal to 2 and divide by the number of simulations.

Letting \(M = 100,000\), and running five times, we get

import random
M = 10000  # Define the value of M
random.seed(1234)
for i in range(5):
    y = [[randint(0, 1) for n in range(3)] for m in range(M)]  # Generate random integers and store them in a 2D list y
    z = [sum(row) for row in y]  # Calculate the sum of each row in y and store them in a list z
    Pr_is_two = sum([1 for i in z if i == 2]) / M  # Calculate the proportion of elements in z that are equal to 2
    print("Pr_is_two =", Pr_is_two)
 

Pr_is_two = 0.3709
Pr_is_two = 0.3712
Pr_is_two = 0.3707
Pr_is_two = 0.3721
Pr_is_two = 0.3725

Nailing down that final digit is going to require one hundred times as many iterations (i.e, \(M = 10,000,000\) iterations). Let’s see what that looks like.

from numpy.random import binomial
M = 10000000

for k in range(1, 6):
    Z = binomial(3, 0.5, M)
    Pr_Z_is_2 = sum(Z == 2) / M
    print(f"Pr[Z == 2] = {Pr_Z_is_2:.3f}")

Pr[Z == 2] = 0.375
Pr[Z == 2] = 0.375
Pr[Z == 2] = 0.375
Pr[Z == 2] = 0.375
Pr[Z == 2] = 0.375

We can do the same for the other numbers, to get a complete picture of taking the probability of each number of heads in separate coin flips.

M = 10000000

for z in range(0, 4):
    Z = binomial(3, 0.5, M)
    Pr_Z_is_z = s

Pr[Z == 0] = 0.125
Pr[Z == 1] = 0.375
Pr[Z == 2] = 0.375
Pr[Z == 3] = 0.125

What if we flip four coins instead of three?

M = 10000000

for z in range(0, 5):
    Z = binomial(4, 0.5, M)
    Pr_Z_is_z = sum(Z == z) / M
    print(f"Pr[Z == {z}] = {Pr_Z_is_z:.3f}")

Pr[Z == 0] = 0.063
Pr[Z == 1] = 0.250
Pr[Z == 2] = 0.375
Pr[Z == 3] = 0.250
Pr[Z == 4] = 0.063

Discrete random variables

So far, we have only considered random numbers that take a finite number of integer values. A random variable that only takes values in the integers, i.e., values in

[\mathbb{Z} = \ldots -2, -1, 0, 1, 2, \ldots]

is said to be a discrete random variable.^[In general, any countable set of numerical values could be used as values of a discrete random variable. A set of values is countable if each of its members can be assigned a unique counting number in \(\mathbb{N} = 0, 1, 2, \ldots\). The integers \(\mathbb{Z}\) can be mapped to natural numbers \(\mathbb{N}\) by interleaving, \(\begin{array}{rcl}\mathbb{Z} & & \mathbb{N} \\ \hline 0 & \mapsto & 0 \\ -1 & \mapsto & 1 \\ 1 & \mapsto & 2 \\ -2 & \mapsto &3 \\ 2 & \mapsto & 4 \\ & \vdots & \end{array}\)]

Probability mass functions

It’s going to be convenient to have a function that maps each possible outcome in a variable to its probability. In general, this will be possible if and only if the variable is discrete, as defined in the previous section.

For example, if we reconsider \(Z = Y_1 + \cdots Y_4\), the number of heads in four separate coin flips, we can define a function^[We implicitly assume that functions return zero for arguments not listed.]

[\begin{array}{rclll} p_Z(0) & = & 1/16 & & \mathrm{TTTT}
p_Z(1) & = & 4/16 & & \mathrm{HTTT, THTT, TTHT, TTTH}
p_Z(2) & = & 6/16 & & \mathrm{HHTT, HTHT, HTTH, THHT, THTH, TTTH}
p_Z(3) & = & 4/16 & & \mathrm{HHHT, HHTH, HTHH, THHH}
p_Z(4) & = & 1/16 & & \mathrm{HHHH} \end{array}]

There are sixteen possible outcomes of flipping four coins. Because the flips are separate and fair, each possible outcome is equally likely. The sequences corresponding to each count of heads (i.e., value of \(Z\)) are recorded in the rightmost columns. The probabilities are derived by dividing the number of ways a value for \(Z\) can arise by the number of possible outcomes.

This function \(p_Z\) was constructed to map a value \(u\) for \(Z\) to the event probability that \(Z = u\),^[Conventionally, this is written as \(p_Z(z) = \mbox{Pr}[Z = z],\) but that can be confusing with upper case \(Z\) denoting a random variable and lower case \(z\) denoting an ordinary variable.]

[p_Z(u) = \mbox{Pr}[Z = u].]

A function defined as above is said to be the probability mass function of the random variable \(Z\). Every discrete random variable has a unique probability mass function.

Probability mass functions represent probabilities of a discrete set of outcomes. The sum of all such probabilities must be one because at least one of the outcomes must occur.^[More formally, if \(Y\) is a discrete random variable, then \(\sum_u \, p_Y(u) = 1,\) where the summation variable \(u\) ranges over all possible values of \(Y\). We are going to start writing this with the standard overloading of lower and upper case \(Y\) as \(\sum_y \, p_Y(y) = 1.\)]

With large numbers of counts based on simulation, we can more readily apprehend what is going on with a plot. Discrete simulations are typically plotted using bar plots, where the outcomes are arrayed on the \(x\) axis with a vertical bar over each one whose height is proportional to the frequency of that outcome.

Plot of \(M = 100,000\) simulations of the probability mass function of a random variable defined as the number of heads in ten specific coin flips.

import numpy as np
import matplotlib.pyplot as plt

np.random.seed(1234)

M = 100000
u = np.random.binomial(10, 0.5, M)
x = np.arange(0, 11)
y = np.zeros(11)

for n in range(0, 11):
    y[n] = np.sum(u == n)

plt.bar(x, y, color='#F8F8F0', edgecolor='black')
plt.xticks(np.arange(0, 11, 2), labels=['0', '2', '4', '6', '8', '10'])
plt.yticks(np.arange(0, 35000, 10000), labels=['0', '10000', '20000', '30000'])
plt.xlabel('Z')
plt.ylabel('count')
plt.show()

The actual frequencies are not relevant, only the relative sizes. A simple probability estimate from simulation provides a probability for each outcome proportional to its height.^[And proportional to its area because the bars are of equal width.]

This plot can easily be repeated to see what happens as the number of bins grows.

Plot of \(M = 1,000,000\) simulations of a variable \(Z\) representing the number of heads in \(N\) coin flips. Each plot represents a different \(N\). Because the bars are the same width and the \(x\) axes are scaled to the same range in all plots, the total length of all bars laid end to end is the same in each plot; similarly, the total area of the bars in each plot is the same.

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns

np.random.seed(1234)
df = pd.DataFrame()
for N in [5, 10, 15, 20, 25, 30]:
    M = 1000000
    u = np.random.binomial(N, 0.5, size=M)
    x = np.arange(N+1)
    y = np.array([np.sum(u==n) for n in x])
    df = df.append(pd.DataFrame({'N': np.repeat(N, N+1),
                                  'Z': x,
                                  'count': y}))

bar_plot = sns.catplot(data=df, x='Z', y='count', kind='bar',
                       col='N', col_wrap=3, palette='Blues_d')
bar_plot.set_axis_labels('Z', 'Count')
bar_plot.set_xticklabels([str(i) if i%10==0 else '' for i in x])
plt.show()

Dice

Let \(Y\) be a random variable representing a fair throw of a six-sided die. We can describe this variable easily through it’s probability mass function, which is uniform (i.e., assigns each possible outcome the same probability).

[p_Y(y) \ =
\begin{cases} \frac{1}{6} & \mbox{if} \ y \in 1:6 \[4pt] 0 & \mbox{otherwise} \end{cases}]

Games like Monopoly use a pair of six-sided dice and consider the sum of the results. That is, \(Y_1\) and \(Y_2\) are fair six-sided die rolls and \(Z = Y_1 + Y_2\) is the result. Games like Dungeons \& Dragons use a trio of six-sided dice and consider the sum of the results. In that scenario, \(Y_1, Y_2, Y_3\) are the results of fair six-sided die rolls and \(Z = Y_1 + Y_2 + Y_3\). Dungeons and Dragons also uses four six-sided die of which the best 3 are summed to produce a result. Let’s simulate some of these approaches and see what the results look like based on \(M = 100,0000\) simulations.^[The simulations are identical to before, only using 1:6 in the range of uniform variables.] Estimated \(p_Y(y)\) for case of \(Y\) being the sum of three six-sided dice (3d6) or the sum of the highest three of four six-sided dice (3 of 4d6).

import numpy as np
import pandas as pd
from plotnine import *

M = 100000

three_d6 = np.random.randint(1, 7, size=M) + np.random.randint(1, 7, size=M) + np.random.randint(1, 7, size=M)

three_of_4d6 = np.zeros(M)
for m in range(M):
    rolls = np.sort(np.random.randint(1, 7, size=4))
    three_of_4d6[m] = np.sum(rolls[1:])

tot_3d6 = np.zeros(16)
tot_3_of_4d6 = np.zeros(16)
for n in range(3, 19):
    tot_3d6[n - 3] = np.sum(three_d6 == n)
    tot_3_of_4d6[n - 3] = np.sum(three_of_4d6 == n)

att_dice_df = pd.DataFrame({'probability': np.concatenate([tot_3d6 / M, tot_3_of_4d6 / M]),
'roll': np.concatenate([np.arange(3, 19), np.arange(3, 19)]),
'dice': np.concatenate([np.repeat('3d6', 16), np.repeat('3 of 4d6', 16)])})

att_dice_plot = (ggplot(att_dice_df, aes(x='roll', y='probability')) +
geom_bar(stat='identity', colour='black', fill='#F8F8F0') +
facet_wrap('dice', labeller=labeller(dice=label_both)) +
scale_x_continuous(name='y',
breaks=[3, 6, 9, 12, 15, 18],
labels=[3, 6, 9, 12, 15, 18]) +
scale_y_continuous(name="estimated " + "p[Y](y)") +
ggtitle("Distribution of 3d6 vs. 3 of 4d6 rolls") +
theme(panel_spacing=2.0))
att_dice_plot

Dungeons and Dragons also uses 20-sided dice.^[In physical games, an icosahedral die is used. The icosahedron is a regular polyhedron with 20 equilateral triangular faces, giving it the largest number of faces among the five Platonic solids.] The fifth edition of the game introduced the notion of advantage, where two 20-sided dice are rolled and the higher result retained, as well as disadvantage, which retains the lower result of the two dice. Here’s a simulation using \(M = 100,000\). The counts are converted to estimated probabilities on the vertical axis in the usual way by dividing by \(M\).

Estimated \(p_Y(y)\) for case of \(Y\) being a single twenty-sided die (d20), the higher two twenty-sided die rolls (max 2d20), and the lower of two 20-sided die rolls (min 2d20).

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns

M = 100000
d20 = np.random.choice(20, size=M, replace=True)
max_2d20 = np.zeros(M)
min_2d20 = np.zeros(M)

for m in range(M):
    max_2d20[m] = max(np.random.choice(20, size=2, replace=True))
    min_2d20[m] = min(np.random.choice(20, size=2, replace=True))

tot_d20 = np.zeros(20)
tot_max_2d20 = np.zeros(20)
tot_min_2d20 = np.zeros(20)

for n in range(20):
    tot_d20[n] = sum(d20 == (n+1))
    tot_max_2d20[n] = sum(max_2d20 == (n+1))
    tot_min_2d20[n] = sum(min_2d20 == (n+1))

att_d20_df = pd.DataFrame({'probability': np.concatenate([tot_d20/M, tot_max_2d20/M, tot_min_2d20/M]),
                           'roll': np.concatenate([np.arange(1,21), np.arange(1,21), np.arange(1,21)]),
                           'dice': np.repeat(['1d20', 'max 2d20', 'min 2d20'], 20)})

# create a FacetGrid with 3 columns
g = sns.FacetGrid(att_d20_df, col='dice', col_wrap=3, height=4, aspect=1)

# plot each histogram using FacetGrid.map()
g.map(sns.barplot, 'roll', 'probability', color="blue")#F8F8F0")

# set x and y axis labels and titles
g.set_axis_labels('y', 'estimated $p_Y(y)$')
g.fig.subplots_adjust(top=0.9)
g.fig.suptitle('')

# set xticks for each plot
for ax in g.axes.flat:
    ax.set_xticks([0, 4, 9, 14, 19])
    ax.set_xticklabels([1, 5, 10, 15, 20])

plt.show()

The most likely roll is a 20 when taking the best of two rolls and the most likely roll is 1 when taking the worst of two rolls.^[The chance for a 20 when taking the best of two 20-sided die rolls is \(1 - \left(\frac{19}{20}\right)^2 \approx 0.098\); the chance of rolling 1 is \(\left(\frac{1}{20}\right)^2 = 0.0025\). The probabilities are reversed when taking the worst of two 20-sided die rolls.] The min and max plots are mirror images of each other as is to be expected by the consecutive nature of the numbers and the min/max operations.

Spinners and the Bernoulli distribution

Some games, such as All Star Baseball, come with spinners rather than dice. The beauty of spinners is that they can be divided into two areas, one with a 27% chance of occurring in a fair spin and one with a 73% chance of occurring.^[27% is roughly the chance of a hit in an at bat.] Now suppose we have a random variable \(Y\) representing the result of a fair spin. Its probability mass function is

[p_Y(y) \ =
\begin{cases} 0.27 & \mbox{if} \ \ y = 1
0.73 & \mbox{if} \ \ y = 0 \end{cases}]

To simplify our notation, we are going to start defining useful functions that can be used as probability mass functions. Our first example is the so-called Bernoulli^[Named after Jacob Bernoulli (1654–1705), one of several prominent mathematicians in the family.] distribution. We define it as a function with somewhat peculiar notation,

[\mathrm{Bernoulli}(y \mid \theta) \ =
\begin{cases} \theta & \mbox{if} \ \ y = 1
1 - \theta & \mbox{if} \ \ y = 0 \end{cases}]

The vertical bar (\(\mid\)) separates the variate argument \(y\), which we think of as an outcome, from the parameter argument \(\theta \in [0, 1]\), which determines the probability of the outcome. In this case, the variate \(y\) is discrete, and can take on only the values zero and one, so we write \(y \in 0:1\). The parameter \(\theta\), on the other hand, is continuous and can take on any value between zero and one (inclusive of endpoints), so we write \(y \in [0, 1]\).^[Interval notation \([0, 1]\) is used for the set of values \(x\) such that \(0 \leq x \leq 1\). Parentheses are used for exclusive endpoints, so that \((0, 1)\) is taken to be the set of \(x\) such that \(0 < x < 1\).]

This notation allows us to simplify our baseball example. Going back to our example random variable \(Y\) which had a 27% chance of being 1 and a 73% chance of being 0, we can write

[p_Y(y) = \mathrm{Bernoulli}(y \mid 0.27).]

To simplify notation even further, we will say that a random variable \(U\) has a Bernoulli distribution and write

[U \sim \mathrm{Bernoulli}(\theta)]

to indicate that the probability mass function of \(U\) is

[p_U(u) = \mathrm{Bernoulli}(u \mid \theta).]

Cumulative distribution functions

In Dungeons \& Dragons, the players are often concerned with probabilities of rolling higher than a given number (or equivalently, rolling lower than a given number. For example, they may need to roll a 15 to sneak by an orc. Such probabilities are conventionally given in the form of cumulative distribution functions. If \(Y\) is a random variable, its cumulative distribution function \(F_Y\) is defined by

[F_Y(y) = \mbox{Pr}[Y \leq y].]

The event probability on the right is calculated the same way as always in a simulation, by counting the number of simulated values in which the condition holds and dividing by the number of simulations.

We can plot the cumulative distribution function for the straight twenty-sided die roll and the rolls with advantage (best of two rolls) or disadvantage (worst of two rolls). Here’s the result using the same simulations as in the last plot, with \(M = 100,000\).

Cumulative distribution function for three variables corresponding to rolling a single 20-sided die, or rolling two 20-sided dice and taking the best or worst result.

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns

M = 100000
d20 = np.random.choice(20, size=M, replace=True)
max_2d20 = np.zeros(M)
min_2d20 = np.zeros(M)

for m in range(M):
    max_2d20[m] = max(np.random.choice(20, size=2, replace=True))
    min_2d20[m] = min(np.random.choice(20, size=2, replace=True))

tot_d20 = np.zeros(20)
tot_max_2d20 = np.zeros(20)
tot_min_2d20 = np.zeros(20)

for n in range(20):
    tot_d20[n] = np.sum(d20 == (n+1))
    tot_max_2d20[n] = np.sum(max_2d20 == (n+1))
    tot_min_2d20[n] = np.sum(min_2d20 == (n+1))

cum_d20 = np.cumsum(tot_d20)
cum_max_2d20 = np.cumsum(tot_max_2d20)
cum_min_2d20 = np.cumsum(tot_min_2d20)

cum_d20_df = pd.DataFrame({
    'probability': np.concatenate((cum_d20/M, cum_max_2d20/M, cum_min_2d20/M)),
    'roll': np.concatenate((np.arange(1, 21), np.arange(1, 21), np.arange(1, 21))),
    'dice': np.concatenate((np.repeat("1d20", 20), np.repeat("max 2d20", 20), np.repeat("min 2d20", 20)))
})

cum_d20_plot = (
    sns.relplot(data=cum_d20_df, x='roll', y='probability', kind='line', hue='dice',
                facet_kws={'margin_titles': True}, col='dice', col_wrap=3)
    .set(xlabel='y', xticks=[1, 5, 10, 15, 20],
         xticklabels=[1, 5, 10, 15, 20], ylabel=r'estimated $F_Y(y)$')
    .set_titles("{col_name}")
    .fig
)

plt.show()

The plot is rendered as a point plot, though this isn’t quite sensible for discrete distributions—the intermediate values are not real. It’s easy to see that there’s a 50% chance of rolling 5 or lower in a single die throw; with the best of 2 it’s more like a 20% chance and with the worst of 2, it’s more like a 75% chance.

Usually in Dungeons \& Dragons, players care about rolling more than a given number, not less, or they’d have to be subtracting all the time. This is where the complementary cumulative distribution function comes in. For a random variable \(Y\), the complementary cumulative distribution function is

[F^{\complement}_Y(y) \ = \ 1 - F_Y(y) \ = \ \mbox{Pr}[Y > y].]

It’s easier to see with a plot how it relates to the usual cumulative distribution function.

Complementary cumulative distributions for a single 20-sided die, the best of two dice, and the worst of two dice.

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns

M = 100000
d20 = np.random.choice(20, size=M, replace=True)
max_2d20 = np.zeros(M)
min_2d20 = np.zeros(M)

for m in range(M):
    max_2d20[m] = max(np.random.choice(20, size=2, replace=True))
    min_2d20[m] = min(np.random.choice(20, size=2, replace=True))

tot_d20 = np.zeros(20)
tot_max_2d20 = np.zeros(20)
tot_min_2d20 = np.zeros(20)

for n in range(20):
    tot_d20[n] = np.sum(d20 == (n+1))
    tot_max_2d20[n] = np.sum(max_2d20 == (n+1))
    tot_min_2d20[n] = np.sum(min_2d20 == (n+1))

cum_d20 = np.cumsum(tot_d20)
cum_max_2d20 = np.cumsum(tot_max_2d20)
cum_min_2d20 = np.cumsum(tot_min_2d20)

ccum_d20_df = pd.DataFrame({
    'probability': np.concatenate((1 - cum_d20/M, 1 - cum_max_2d20/M, 1 - cum_min_2d20/M)),
    'roll': np.concatenate((np.arange(1, 21), np.arange(1, 21), np.arange(1, 21))),
    'dice': np.concatenate((np.repeat("1d20", 20), np.repeat("max 2d20", 20), np.repeat("min 2d20", 20)))
})

ccum_d20_plot = (
    sns.relplot(data=ccum_d20_df, x='roll', y='probability', kind='line', hue='dice',
                facet_kws={'margin_titles': True}, col='dice', col_wrap=3)
    .set(xlabel='y', xticks=[1, 5, 10, 15, 20],
         xticklabels=[1, 5, 10, 15, 20], ylabel=r'estimated $F_Y(y)$')
    .set_titles("{col_name}")
    .fig
)

Infinite discrete random variables

Consider an experiment in which a coin is tossed until a heads appears. Let the random variable \(U\) be the number of tosses that came up tails before the first head comes up. The legal sequences are H (0 tails), TH (1 tails), TTH (2 tails), and so on. There is no upper limit to how many tails may appear before the first heads.

Here’s some code to create \(M\) simulations of the variable \(U\).^[This code uses a while loop, which repeats as long as its condition evaluates to true (i.e., 1). Here, the condition compares the output of the random number generator directly rather than assigning to an intermediate value. We have also introduced the increment operator +=, which adds the value of the right hand side to the variable on the left hand side.]

import random

M = 10 # replace with the desired value of M
u = [0] * M

for m in range(M):
    while random.randint(0, 1) == 0:
        u[m] += 1

This looks dangerous! The body of a while-loop (here u[m] + 1) is executed iteratively as long as the condition is true.^[If we write while (1 + 1 == 2) we produce what is known as an infinite loop, i.e., one that never terminates.] Shouldn’t we be worried that the random number generator will just continue to throw tails (i.e., 0) so that the program never terminates?^[The answer is “yes,” in general, because programmers are error prone.] In this case, no, because the odds are vanishingly small that \(U\) gets large. For example,

[\begin{array}{rcl} \mbox{Pr}[U < 10] & = & p_U(0) + p_U(1) + \cdots p_U(9) \[6pt] & = & \frac{1}{2} + \cdots + \frac{1}{1024} \[6pt] & \approx & 0.999. \end{array}]

Going further, \(\mbox{Pr}[U < 20] \approx 0.999\,999\), and so on. So there’s not much chance of running very long at all, much less forever.

With the concern of non-termination out of the way, let’s see what we get with \(M = 50\) simulations of \(U\).

import random

def sim_u():
    u = 0
    while random.randint(0, 1) == 0:
        u += 1
    return u

random.seed(1234)
for m1 in range(1, 6):
    for m2 in range(1, 11):
        print(f"{sim_u():4d}", end="")
    print()

0	6	4	0	0	5	2	3	1	0 
0	0	2	3	1	0	3	5	2	0 
1	0	2	1	1	1	0	0	0	0 
0	1	1	0	0	0	2	0	2	0 
0	1	2	1	0	2	2	3	0	0

It’s very hard to discern a pattern here. There are a lot of zero values, but also some large values. For cases like these, we can use a bar plot to plot the values. This time, we’re going to use \(M = 10,000\) to get a better picture of the pattern.

Frequency of outcomes in \(10,000\) simulation draws of \(U\), the number of tails seen before a head in a coin-tossing experiment.

import random
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns

def sim_u():
    u = 0
    while random.randint(0, 1) == 0:
        u += 1
    return u

random.seed(1234)
M = 10000
u = np.array([sim_u() for m in range(M)])
x = np.arange(max(u) + 1)
y = np.array([np.sum(u == n) for n in range(max(u) + 1)])

bar_plot = sns.barplot(x=x, y=y, color="blue")#F8F8F0")
bar_plot.set(xlabel="U", ylabel="count")
bar_plot.set(xticks=np.arange(13), xticklabels=np.arange(13))
bar_plot.set(yticks=[1000, 2000, 3000, 4000, 5000], yticklabels=[1000, 2000, 3000, 4000, 5000])
plt.show()

The \(x\)-axis represents the value of \(U\) and the \(y\)-axis the number of times that value arose in the simulation.^[Despite \(U\) having infinitely many possible values, it will only take on finitely many of them in a finite sample.] Each additional throw of tails appears to cut the probability of occurrence in half result seems to have about half the probability of the previous one. This is what we should expect because each coin toss brings a 50% probability of a tails result. This exponential decay^[Exponential decay means each additional outcome is only a fraction as likely as the previous one.] in the counts with the number of tails thrown is more obvious when plotted on the log scale.

Frequency of outcomes in \(10,000\) simulation draws of \(U\), the number of tails seen before a head in a coin-tossing experiment, this time with the outcome count on the log scale to illustrate the exponentially decreasing probabilities of each successive number of tails.

import random
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns

def sim_u():
    u = 0
    while random.randint(0, 1) == 0:
        u += 1
    return u

random.seed('1234')
M = 10000
u = np.array([sim_u() for m in range(M)])
x = np.arange(max(u) + 1)
y = np.array([np.sum(u == n) for n in range(max(u) + 1)])

log_bar_plot = sns.barplot(x=x, y=y, color="blue")#F8F8F0")
log_bar_plot.set(xlabel="U", ylabel="count")
log_bar_plot.set(xticks=np.arange(13), xticklabels=np.arange(13))
log_bar_plot.set(yticks=[8, 40, 200, 1000, 5000], yticklabels=[8, 40, 200, 1000, 5000])
log_bar_plot.set(yscale="log")
plt.show()

There is a 50% probability that the first toss is heads, yielding a sequence of zero tails, and \(U = 0\). Each successive number of tails is half as likely as the previous, because another tail will have to be thrown, which has a 50% probability.^[In symbols, \(\mbox{Pr}[U = n + 1] \ = \ \frac{1}{2} \mbox{Pr}[U = n].\) ]

Thus the overall probability mass function for \(U\) is^[An elementary result of calculus is that \(\sum_{n = 0}^{\infty} \frac{1}{2^{n + 1}} \ = \ \frac{1}{2} + \frac{1}{4} + \frac{1}{8} + \frac{1}{16} + \cdots \ = \ 1.\)]

[\begin{array}{rclcl} p_U(0) & = & \frac{1}{2} & = & \frac{1}{2}
p_U(1) & = & \frac{1}{2} \times \frac{1}{2} & = & \frac{1}{4}
p_U(2) & = & \frac{1}{2} \times \frac{1}{2} \times \frac{1}{2} & = & \frac{1}{8}
& \vdots &
p_U(u) & = & \underbrace{\frac{1}{2} \times \cdots \times \frac{1}{2}}_{u + 1 \ \mathrm{times}} & = & \left( \frac{1}{2} \right)^{u + 1}
& \vdots & \end{array}]

Even though there are infinitely many possible realizations of the random variable \(U\), simulation may still be used to compute event probabilities, such as \(\mbox{Pr}[U \leq 3]\), by

M = 10000
u = np.array([sim_u() for m in range(M)])
leq3 = np.array([u[m] <= 3 for m in range(M)])

print('est Pr[U <= 3] =', sum(leq3) / M)

Let’s see what we get with \(M = 100,000\),

est Pr[U <= 3] = 0.9394

Writing out the analytic answer involves an infinite sum,

[\mbox{Pr}[U \leq 3] \ =
\sum_{u = 0}^{\infty} p_U(u) \mathrm{I}[u \leq 3].]

We can recognize that all of the terms where \(u > 3\) are zero, so that this reduces to

[\begin{array}{rcl} \mbox{Pr}[U \leq 3] & = & p_U(0) + p_U(1) + p_U(2) + p_U(3) \[8pt] & = & \frac{1}{2}

• \frac{1}{4}
• \frac{1}{8}
• \frac{1}{16} \[8pt] & \approx & 0.938. \end{array}]

Simulation is not that clever. It just blindly simulates values of \(u\), many of which turn out to be larger than three. In the sequence of simulated values, many were larger than three—the histogram summarized a much larger simulation, none of the values of which were that large.

Given that we only ever run a finite number of iterations and thus only ever see a finite number of values, how does simulation get the right answer?^[Spoiler alert! The technical answer, like so much in statistics, is that the central limit theorem kicks in.] Isn’t there some kind of bias from only visiting smallish values? The answer is “no” precisely because the values that are not simulated are so rare. Their total probability mass, when added together is small, so they cannot have much influence on the simulated answer. For simulations to get the right answer,^[Given some qualifications!] they need only visit the typical values seen in simulation, not the extreme values.^[The sparsity problem grows exponentially worse in higher dimensions and uncountably worse with continuous random variables.]

The Chevalier de Méré’s challenge

Antoine Gombaud, the Chevalier de Méré,^[Self appointed!] challenged Blaise Pascal to explain how it was possible that the probability of throwing at least one six in four throws of a single six-sided die is slightly greater than \(\frac{1}{2}\), whereas the probability of throwing two sixes in 24 throws of a pair of six-sided die was slightly less than \(\frac{1}{2}\).^[Smith, D. A., 1929. A Source Book on Mathematics. McGraw-Hill, New York; cited in Bulmer, 1967, p. 26.] We can evaluate these claims by simulation directly.^[Working out the example analytically, there is a \(\frac{35}{36}\) chance of not throwing double six with two six-sided dice, and so \(\left( \frac{35}{36} \right)^{24}\) is the probability of not throwing at least one double six in 24 throws, and so \(1 - \left( \frac{35}{36} \right)^{24} \approx 0.491\) is the probability of throwing at least one double-six in 24 fair throws of a pair of six-sided dice.]