0% found this document useful (0 votes)
2 views111 pages

Course Notes (Complete)

The document is a set of lecture notes for a Statistics Lab course, covering essential topics in probability and statistics, including data description, random variables, and statistical testing. It aims to provide students with a thorough understanding of statistical concepts, supported by examples and Python code for practical application. The content is structured into chapters that detail various statistical methods and their applications in data analysis.

Uploaded by

kagarwal.3115
Copyright
© All Rights Reserved
We take content rights seriously. If you suspect this is your content, claim it here.
Available Formats
Download as PDF, TXT or read online on Scribd
0% found this document useful (0 votes)
2 views111 pages

Course Notes (Complete)

The document is a set of lecture notes for a Statistics Lab course, covering essential topics in probability and statistics, including data description, random variables, and statistical testing. It aims to provide students with a thorough understanding of statistical concepts, supported by examples and Python code for practical application. The content is structured into chapters that detail various statistical methods and their applications in data analysis.

Uploaded by

kagarwal.3115
Copyright
© All Rights Reserved
We take content rights seriously. If you suspect this is your content, claim it here.
Available Formats
Download as PDF, TXT or read online on Scribd

Statistics Lab

Lecture Notes

Verena Wolf1

June 22, 2026

1
[Link]@[Link]
ii
Contents

1 Data Description 1
1.1 Important Terms . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 Data description . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.2.1 Frequency . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.2.2 Measures of Location . . . . . . . . . . . . . . . . . . 3
1.2.3 Measures of Dispersion . . . . . . . . . . . . . . . . . . 5
1.2.4 Measures of Shape . . . . . . . . . . . . . . . . . . . . 6
1.2.5 Standardization of Data . . . . . . . . . . . . . . . . . 7

2 Probabilities and Combinatorics 9


2.1 Probabilities . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
2.1.1 Conditional Probabilities . . . . . . . . . . . . . . . . 11
2.1.2 Bayes’ rule . . . . . . . . . . . . . . . . . . . . . . . . 13
2.1.3 Independent Events . . . . . . . . . . . . . . . . . . . 15
2.2 Combinatorics . . . . . . . . . . . . . . . . . . . . . . . . . . 18

3 Discrete Random Variables 27


3.1 Discrete Random Variables and Probability Distributions . . 28
3.2 Combining and Transforming Random Variables . . . . . . . 30
3.2.1 Joint Probability Distribution and Independence . . . 32
3.3 Expectation and Variance . . . . . . . . . . . . . . . . . . . . 33
3.3.1 Properties of the Expectation and Variance Operator 35
3.3.2 Conditional Expectations . . . . . . . . . . . . . . . . 36
3.4 Higher Order Moments, Covariance and Correlation . . . . . 37
3.4.1 Covariance and Correlation . . . . . . . . . . . . . . . 37
3.4.2 Moment Generating Function . . . . . . . . . . . . . . 39
3.5 Important Discrete Probability Distributions . . . . . . . . . 41

4 Continuous Random Variables 45


4.1 σ-algebras . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
4.2 Continuous Random Variables . . . . . . . . . . . . . . . . . . 48
4.3 Important Continuous Distributions . . . . . . . . . . . . . . 50
4.4 Multivariate Random Variables . . . . . . . . . . . . . . . . . 54

iii
iv CONTENTS

5 Generation of Random Variates 57


5.1 Generating Discrete Random Variates . . . . . . . . . . . . . 58
5.1.1 Interval Method . . . . . . . . . . . . . . . . . . . . . 59
5.2 Inverse Transform Method . . . . . . . . . . . . . . . . . . . . 60
5.2.1 Inverse transform sampling: the discrete case . . . . . 61
5.3 Rejection Sampling . . . . . . . . . . . . . . . . . . . . . . . . 62

6 Laws of Large Numbers 65


6.1 Chebyshev’s inequality . . . . . . . . . . . . . . . . . . . . . . 65
6.2 Weak Law of Large Numbers . . . . . . . . . . . . . . . . . . 66
6.3 Strong Law of Large Numbers . . . . . . . . . . . . . . . . . . 67
6.4 Central Limit Theorem. . . . . . . . . . . . . . . . . . . . . . 67

7 Parameter Estimation 69
7.1 Method of Moments . . . . . . . . . . . . . . . . . . . . . . . 71
7.2 Maximum Likelihood Estimation . . . . . . . . . . . . . . . . 74
7.2.1 Variance of MLE (optional content) . . . . . . . . . . 79
7.3 Bayesian Inference . . . . . . . . . . . . . . . . . . . . . . . . 82
7.3.1 Conjugate families of distributions . . . . . . . . . . . 84
7.3.2 Bayesian point-estimators . . . . . . . . . . . . . . . . 86

8 Statistical Testing 87
8.1 Level α tests: general approach . . . . . . . . . . . . . . . . . 88
8.2 Standard Normal Null Distribution (Z-test) . . . . . . . . . . 90
8.3 T-tests for Unknown σ . . . . . . . . . . . . . . . . . . . . . . 93
8.4 p-Value . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95
8.5 Chi-square Tests . . . . . . . . . . . . . . . . . . . . . . . . . 97
8.5.1 Estimated variance . . . . . . . . . . . . . . . . . . . . 97
8.5.2 Observed Counts . . . . . . . . . . . . . . . . . . . . . 98
8.5.3 Testing independence . . . . . . . . . . . . . . . . . . 101
Preface

The goal of this course is to guide students to a thorough understanding of


probability and statistics, which play a very prominent role for most state-
of-the-art machine learning tools.
Throughout the course, a large number of examples and exercises will help to
get an intuitive understanding of the covered topics. Whenever appropriate,
Python code will be included to explore important concepts using concrete
data sets.
Before using these lecture notes, please carefully read the following hints:

• We use hyperlinks to Wikipedia pages as convenient entry points


for further exploration. Wikipedia is well-suited for getting a quick
overview, but the quality of articles varies, and for serious work you
should always consult the primary literature (textbooks, peer-reviewed
papers) cited there.

• The Python code will be provided as Jupyter Notebooks and we use


Python 3. Please make sure that you have all the necessary packages
installed and up to date.

v
vi CONTENTS
Chapter 1

Data Description

Descriptive statistics refers to a set of tools that is typically used in a first


step of data analysis. By examining a number of statistical quantities and
graphical displays, a preliminary understanding of the data is achieved.

Chapter learning objectives:

• distinguish different types of data (e.g. ordinal vs. nominal)

• describe data by means of different summary statistics for a data set


(e.g. mean, variance, skewness)

• apply and understand standardization of data

1.1 Important Terms


We consider a data set and refer to it as a sample or sample set from a
(typically larger) population, i.e., a population is a set of elements which is
so large or even infinite that we can only analyze subsets of it.
Hence, we draw random subsets from it
which are called samples. Examples of pop-
ulation are all the citizens of a country, all
cancer patients, all sentences of a certain
language, etc. We typically hope that our
sample set is representative for the popula-
tion in the sense that statistical quantities Figure 1.1: A sample from a
of the sample set are similar to those of the population.
population. Often we are interested in cer-
tain attributes or factors of the population elements such as the height of a
person or the frequency of a word.

1
1.2. DATA DESCRIPTION

Quantitative data
Attributes may refer to some quantity of something and have numerical
values such as the attribute height, for instance. Quantitative data may
have a discrete range (e.g. subset of N) or a continuous range (e.g. subset
of R≥0 ). Attributes such as height, weight, length, and time are continuous
data while counts of days, successes, etc are discrete data.

Qualtitative data
Qualitative data are any type of data that are not numerical. Examples
are name, gender, country, social security number, etc. If a qualitative
attribute has an order it is called ordinal . For instance, assume that we
consider a person’s educational experience (with values such as elementary
school graduate, high school graduate, and college graduate). These also
can be ordered and we can assign numbers accordingly such as elementary
school (1), high school (2), and college (3). Even though we can order these
from lowest to highest, the spacing between the values may not be the same
(e.g. the gap between elementary school and high school may be larger
than the gap between high school (2), and college). Another example for an
ordinal attribute is shoe size where there is a clear order but size 42 is not
twice as large as size 21.
If attribute values do not have an implied ordering, the attribute is called
nominal . For instance, the attribute color is a nominal one because we want
to distinguish blue, red, yellow, etc but do not order the colors in a certain
way.

1.2 Data description


1.2.1 Frequency
Assume that we are interested in an attribute that has a finite number of
possible values a1 , a2 , . . . , ak . Assume now that we have a sample set of size
n and that xi ∈ {a1 , a2 , . . . , ak } is the attribute value of the i-th element in
our sample set (e.g. x1 =blue, x2 =red, x3 =red, etc.). Then the number of
xi ’s in the sample set with value aj is called the absolute frequency of aj ,
written
hj := h(aj ).
The relative frequency is then

fj := hj /n.

Hence
k
X k
X
hj = n, fj = 1.
j=1 j=1

2
1.2. DATA DESCRIPTION

Later, we will discuss frequency tables and how to generate them using
Python.

Example 1: Absolute and relative frequency

Assume that we have a sample set of size n = 3 with the attribute ’color’
taking values in A = {red,blue,green} and the samples are x1 = blue, x2 =
red, x3 = red. Then, the absolute frequencies are

hblue = 1
hred = 2
hgreen = 0

and the relative frequencies are

fblue = 1/3
fred = 2/3
fgreen = 0

1.2.2 Measures of Location


To find out where the center of the data is, the following measures of location
are common.
Assume that we consider a numerical attribute (values may be real or integer
numbers) and that xi is the value of the i-th element in our sample of size
n. In other words, we have data x1 , . . . , xn . Then the arithmetic mean or
sample mean is given by
n
1X
x̄ := xi .
n
i=1

Note that the mean is a very natural measure, but it is sensitive to extreme
values. In the case of highly skewed data, the mean is not a good measure
of location of the data.
The (sample) median is a value that splits the data into two parts of equal
size. To compute the median, we sort the data x1 , x2 , . . . , xn in increasing
order. Given the sorted list x(1) , x(2) , . . . , x(n) and define
(
x( n+1 ) : n is uneven,
xmed := 1
2

2 (x( 2 ) + x( n2 +1) ) : n is even,


n

3
1.2. DATA DESCRIPTION

Example 2: Median and Mean


For instance, if the sorted data is 1,2,2,3,4,4,5,5,5,5,6,7 then the median
is 4+5 49
2 = 4.5 while the arithmetic mean is 12 ≈ 4.08.

Compared to the mean, the median is harder to compute (sorting necessary!)


but more resistant to extreme values.
Question:
Is it possible to modify the sequence of integers in the above example so
that the median roughly gives a better measure of location compared to
the mean?

Another measure of location which can also be used if an attribute is nomi-


nal, is the mode of the data given by the attribute value that appears most
often in our data.

Example 3: Mode
The mode of the data set 1,2,2,3,4,4,5,5,5,5,6,7 is 5. We illustrate the
mode together with mean and median in Figure 1.2.

In general, the mode and the median are less sensitive to outliers than the
sample mean. If the last value in our example data set would be 70 (instead
of 7), the mode and the median would be unchanged, while the sample mean
increases to approximately 9.33.

Figure 1.2: Plot of the data including mode, median, and mean.

4
1.2. DATA DESCRIPTION

1.2.3 Measures of Dispersion


Next we consider quantities that measure how much the data spreads (around
its mean). The most common measure it the sample variance
n
1 X
s2 := (xi − x̄)2
n−1
i=1

and the sample standard deviation


v
u n
√ u 1 X
s := s2 = t (xi − x̄)2 .
n−1
i=1

To compute the sample variance, it is often helpful to exploit the identity


n n
1 X 1 X 2 n
s2 = (xi − x̄)2 = (xi ) − x̄2 .
n−1 n−1 n−1
i=1 i=1

Example 4: Sample Variance

The sample variance of the data set 1,2,2,3,4,4,5,5,5,5,6,7 is


 2
1 2  12 49
s2 = 1 + 2 · 22 + 32 + 2 · 42 + 4 · 52 + 62 + 72 − ≈ 3.17
11 11 12

Intuitively, the sample variance approximates the average squared distance


of the observations from the sample mean. By taking its root, we scale back
to the measurement units of the original data.
Also the range of the data is a measure of dispersion. It is defined as the
difference between the largest and the smallest value, i.e., if xmax = maxi xi
and xmin = mini xi then the range is
xmax − xmin .

Example 5: Range
The range of the data set 1,2,2,3,4,4,5,5,5,5,6,7 is 7-1=6. We illustrate
the range together with mean and standard deviation in Figure 1.3.

The Python code for the sample mean, sample variance and range for our
example data set can be found below. We give a naive implementation for
the mean and variance, as well as the built-in functions from numpy. We
always import numpy as np. Note that the variance function from numpy
has an optional parameter ddof, which we choose as 1 at this point. Later,
we will discuss this issue in more detail.

5
1.2. DATA DESCRIPTION

1 # list with data


2 x =[1 ,2 ,2 ,3 ,4 ,4 ,5 ,5 ,5 ,5 ,6 ,7]
3 # mean
4 x b a r = sum ( x ) / l e n ( x )
5 # sample variance
6 s 2 = sum ( [ ( i −x b a r ) ∗∗2 f o r i i n x ] ) / ( l e n ( x ) −1)
7 # range
8 r = max ( x )−min ( x )
9
10 # built - in function from numpy ( np )
11 # mean
12 x b a r = np . mean ( x )
13 # sample variance
14 s 2 = np . v a r ( x , d d o f =1)

Later, we will also discuss bounds such as Chebychev’s Rule that consider
how much of the data will fall past a certain distance (e.g. k standard
deviations) from the mean.

Figure 1.3: Plot of the data including mean, sample standard deviation, and
range.

1.2.4 Measures of Shape

The distribution of the data may be symmetric, skewed to the right or left.
In Figure1.4 we show the left- and right-skewed case (plots taken from [2]).

For the sample skewness, different definitions exist, for reasons that will
become clear when we discuss important properties of estimators such as
unbiasedness. For now, we only consider the most convenient formula to

6
1.2. DATA DESCRIPTION

Figure 1.4: Left- and right-skewed data (here as a distribution plot).

estimate the skewness of the data:


n
m3 1X
g1 = 3/2
, where mk = (xi − x̄)k .
m2 n
i=1

The skewness can be any value, −∞ < g1 < ∞ and its sign indicates the
direction of skewness (negative values correspond to skewness to the left and
positive values to skewness to the right). If g1 ≈ 0, then the distribution is
(nearly) symmetric.

1.2.5 Standardization of Data


Assume you have two or more sets of data from different sources and would
like to compare the data. Then, it is helpful to transform (standardize) each

Figure 1.5: Histogram plot of data before (left) and after standardization
(right).

7
1.2. DATA DESCRIPTION

sample set x1 , . . . , xn as follows:


xi − x̄
zi := ,
s
i.e., consider the differences of a data point xi to the mean x̄ of the corre-
sponding sample set relative to the standard deviation s.
Standardized data sets have mean zero and standard deviation 1 and can
then be compared to other standardized data sets.
In Figure 1.5 we show histogram plots of some data before (left) and after
standardization (right). Note that the data is shifted such that the mean
becomes zero and scaled so that we have a standard deviation of one.

8
Chapter 2

Probabilities and
Combinatorics

“Was ein Punkt, ein rechter Winkel, ein Kreis ist, weiß ich schon
vor der ersten Geometriestunde, ich kann es nur noch nicht präzisieren.
Ebenso weiß ich schon, was Wahrscheinlichkeit ist, ehe ich es
definiert habe.“ (Hans Freudenthal)

This section gives a short primer on probabilities combinatorics. The prob-


ability models that we discuss next are used to describe chance experiments.
They allow to analyze systems of the real world that are subject to uncer-
tainty - even if no data is available. The system may even be a hypothetical
one or no measurements have been made. For instance, if we toss a fair coin
three times, we can determine the probability of getting head only once.
However, if we would like to know whether the coin is fair or not, a number
of sample tosses can give us statistical evidence for or against fairness of the
coin.

Chapter learning objectives:


• get familiar with basic concepts of probability: sample space, events,
probability functions, etc.

• understand conditional probability, independence (including the cor-


responding laws and rules for computing (conditional) probabilities)

• determine numbers of possible outcomes by applying urn models

2.1 Probabilities
Let us consider chance experiments with a countable number of possible
outcomes ω1 , ω2 , . . .. The set Ω = {ω1 , ω2 , . . .} of all outcomes is called the

9
2.1. PROBABILITIES

sample space. Subsets of Ω are called events and by 2Ω we denote the set
of all events.

Example 6: Rolling a die and tossing a coin

If we roll a die, the set of possible outcomes is Ω = {1, 2, . . . , 6}. The


event “number is even” is given by E = {2, 4, 6}.
If we toss a coin and count the number of trials until a head turns up for
the first time, Ω = {1, 2, . . .}. The event “number of trials greater than
10” is given by E = {11, 12, . . .}.

There are typically many different possibilities to define Ω. In the above


example of a sequence of coin tosses, we already chose Ω such that it fits to
the question that we have in mind (what is the probability of more than 10
trials). If, however, we were interested in the probability of having tails in
the first three tosses, we would rather define Ω as

{(ω1 , ω2 , ω3 )|ω1 , ω2 , ω3 ∈ {H, T }}

and ignore later tosses. Obviously, we could also encode H as 1 and T as 0.


We define what a probability is using Kolmogorov’s axioms.
Definition 1: Probability
Assume Ω is a discrete (i.e. finite or countably infinite) and non-empty
sample space. Let P be a function such that P : 2Ω → [0, 1]. The value
P (A) is called the probability of event A if P is such that
1. P (Ω) = 1,

2. for any finite or countably infinite sequence of pairwise disjoint events


A1 , A2 , . . .
it holds that

P (A1 ∪ A2 ∪ . . .) = P (A1 ) + P (A2 ) + . . . .

When we reason about the probability of



certain events, we can use many arguments
from set theory. For instance, if the events C
A, B, and C are as illustrated in Figure 2.1,
B
it holds that A
P (A ∪ B) = P (A) + P (B)
P (A ∪ C) = P (A) + P (C) − P (A ∩ C)
P (Ω \ B) = P (Ω) − P (B)
= 1 − P (B) = P (B̄) Figure 2.1: Using set argu-
ments for the calculation of
P (A \ C) = P (A) − P (A ∩ C). event probabilities.

10
2.1. PROBABILITIES

Here, B̄ denotes the complement Ω \ B of the set B. In addition, it is easy


to show that P (∅) = 0 and P (A) ≤ 1 for all events A.

Example 7: Rolling a die


The probability of the event {1, 2, . . . , 6} is one. Moreover, for each n ∈
{1, 2, . . . , 6} it holds that P ({n}) = 1/6. Thus, P ({2, 4, 6}) = 1/6 + 1/6 +
1/6 = 1/2. Similarly, P ({1, 3, 5}) = 1 − P ({2, 4, 6}) = 1/2.

Example 8: Tossing a coin until heads for the first time

Consider a chance experiment, where a fair coin is tossed until a head


comes up for the first time. We count the number of trials needed. Ω =
{1, 2, . . .}. The probability of the events “exactly n trials” and “more than
three trails” are
P ({n}) = (1/2)n and
P ({4, 5, . . .}) = P (Ω \ {1, 2, 3}) = 1 − (1/2 + 1/4 + 1/8) = 1/8.
P
Note that P (Ω) = ω∈Ω P ({ω}) = 1/2 + 1/4 + 1/8 + . . . = 1.

The triple (Ω, 2Ω , P ) is called a discrete probability space.

2.1.1 Conditional Probabilities


It is possible to restrict our reasoning to certain conditions, e.g. we consider
rolling a die but restrict to those cases where the number is even. Under this
condition, we could ask what the probability of getting number 2 is. Hence,
we already know that the number is even and consider only even outcomes.
Then, clearly with 1/3 we get two pips.
Definition 2: Conditional Probability
Let A, B be events and P (B) > 0. Then

P (A ∩ B)
P (A|B) :=
P (B)

is called probability of A under the condition B. Clearly, this implies that


P (A ∩ B) = P (A|B) · P (B).

Sometimes, we simply write PB (A) instead of P (A|B). It is easy to show


that PB is a probability in the sense of Definition 1.

11
2.1. PROBABILITIES

Example 9: Lung Cancer

We define the events


72
A: person gets lung cancer with P (A) = 200000 = 0.00036,
B: person is a smoker with P (B) = 0.25.
From the people that get lung cancer, 90% are smokers. The experiment
consists in choosing a person at random. For the probability of getting
lung cancer under the condition of being a smoker, we calculate
P (A∩B) P (B|A)P (A) 0.9·0.00036
P (A|B) = P (B) = P (B) = 0.25 = 0.001296.

For the probability of getting lung cancer under the condition of not being
a smoker, we get
P (A∩B̄) P (B̄|A)P (A) (1−0.9)·0.00036
P (A|B̄) = P (B)
= P (B)
= 0.75 = 0.000048.

Thus, the chance of getting lung cancer is around 27 times higher for
smokers compared to non-smokers.

Multiplication rule
Multiplication rule
Let A1 , A2 , . . . , An be events with P (A1 ∩ A2 ∩ . . . ∩ An−1 ) > 0. Then the
following multiplication rule holds:

P (A1 ∩ A2 ∩ . . . ∩ An )
= P (A1 ) · P (A2 |A1 ) · P (A3 |A1 ∩ A2 ) · · · · · P (An |A1 ∩ A2 ∩ . . . ∩ An−1 )

Example 10: Drawing from an urn

Suppose you have an urn with 3 red, 2 blue and 1 green ball. You draw
three times without placing the drawn balls back into the urn. The prob-
ability of drawing blue, green, red (in that order) is
2 1 3
· · .
6 5 4
Note that we will discuss scenarios like this in more detail later.

A1 Ω
Law of total probability A2
Assume that we have a finite or countably infi- B
nite number of events A1 , A2 , . . . that are pairwise
disjoint. Assume further Ω = A1 ∪ A2 ∪ . . . and A3
A4
12
Figure 2.2: The law of
total probability.
2.1. PROBABILITIES

P (Ai ) > 0 for all i. Often, the probability of some


event B is unknown, but the conditional probabil-
ities P (B|Ai ), . . . are known. In this case, P (B)
can be computed using the law of total probability, which states that
X
P (B) = P (B|Ai ) · P (Ai ) .
| {z }
i
=P (B∩Ai )

The corresponding partitioning of Ω is illustrated in Fig. 2.2.

Example 11: Weather Forecast

In a city on 70 from 100 days the weather is good (G) and on 30 from
100 days the weather is bad (Ḡ). The local meteorologist can predict good
weather with 90% accuracy and bad weather with 60% accuracy. According
to the law of total probability the forecast is correct with a probability of

P (M ) = P (M |G) · P (G) + P (M |Ḡ) · P (Ḡ) = 0.9 · 0.7 + 0.6 · 0.3 = 0.81.

2.1.2 Bayes’ rule


The law of total probability is helpful, if we are
interested in P (B). If, however, P (Ai |B) is the
probability of interest, one can use Bayes’ theorem
to compute it. With Bayes’ theorem we can for
any event A with P (B) > 0 express P (A|B) by
P (B|A), that is,
P (B|A) · P (A)
P (A|B) = .
P (B)

Thus, for the above problem of computing P (Ai |B) for pairwise disjoint
events Ai , this gives

P (Ai |B) =
P (B∩Ai ) P (B|Ai )·P (Ai )
P (B) = P (B) = PP (B|A i )·P (Ai )
.
i P (B|Ai )·P (Ai )

Example 12: Noisy Channel


Consider the noisy binary channel illus-
trated on the right. 1−p
If the channel is noise-free (p = 0), zero is 0 0
transmitted from the upper left node to the p
upper right node. Similarly, one is trans- p
mitted from the lower left node to the lower
1 1
right node. 13 1−p
2.1. PROBABILITIES

If the channel is noisy, with probability p > 0 one is transmitted instead


of zero and zero instead of one. Assume that the probability of sending
zero is π0 and the probability of sending one is π1 = 1 − π0 .
We define the events

• A0 : send 0 with P (A0 ) = π0 ,

• A1 : send 1 with P (A1 ) = π1 := 1 − π0 ,

• B0 : receive 0,

• B1 : receive 1,

Then, P (B1 |A0 ) = P (B0 |A1 ) = p and P (B1 |A1 ) = P (B0 |A0 ) = 1 − p.
We calculate
P (B1 ) = P (B1 |A0 ) · P (A0 ) + P (B1 |A1 ) · P (A1 ) = p · π0 + (1 − p) · π1 ,
P (B0 ) = P (B0 |A0 ) · P (A0 ) + P (B0 |A1 ) · P (A1 ) = (1 − p) · π0 + p · π1 .

and the receiver can determine the probability that the transmission was
correct as
P (B1 |A1 )·P (A1 ) (1−p)·π1
P (A1 |B1 ) = P (B1 ) = p·π0 +(1−p)·π1
P (B0 |A0 )·P (A0 ) (1−p)·π0
P (A0 |B0 ) = P (B0 ) = (1−p)·π0 +p·π1

Sometimes, it is advantageous to consider an equation for the ratio (the


odds) P (A|C)/P (B|C), which is straightforwardly derived from Bayes’ the-
orem.
P (A|C) P (C|A) P (A)
=
P (B|C) P (C|B) P (B)

Example 13: Bayesian Odds

There are two urns, one containing 7 red and 3 blue balls, the other
containing 3 red and 7 blue balls. We flip a fair coin to determine, from
which urn we draw 12 balls with replacement. As a result, we get 8 times
a red and 4 times a blue ball. What is the probability that it was the first
urn with predominantly red balls?
Let U1 (U2 ) be the event of selecting the first (second) urn, respectively.
Further, let A be the event of getting 8 times a red and 4 times a blue
ball. Some concepts from combinatorics are necessary to determine the
probabilities P (A|U1 ) and P (A|U2 ). They will be discussed in the next

14
2.1. PROBABILITIES

section. For now, just take it as a fact that

P (A|U1 ) ≈ 0.231
P (A|U2 ) ≈ 0.0078

Moreover, as we choose each urn with 1/2, P (U1 ) = P (U2 ) = 1/2. We


are interested in the ratio between P (U1 |A) and P (U2 |A), because this
tells us, whether U1 or U2 is more likely, given the result A:

P (U1 |A) P (A|U1 ) P (U1 )


= ≈ 29.6.
P (U2 |A) P (A|U2 ) P (U2 )

Hence, given A it is about 29 times more likely that we drew from urn 1.
We get P (U1 |A) by exploiting P (U1 |A)+P (U2 |A) = 1 and the above ratio:
P (U1 |A)
P (U1 |A) = 1 − P (U2 |A) = 1 − 29.642
⇐⇒ P (U1 |A) ≈ 0.967

Thus, with nearly 97%, we drew from urn 1, which is intuitive since we
got many more red balls than blue ones.

2.1.3 Independent Events


Assume that we toss a coin twice. The sample space is {HH, HT, T H, T T }
and P (B) = 2/4 = 1/2 if B is the event that the first toss is H. If A is the
event that the second toss is H, P (A) = 1/2 and P (B ∩ A) = 1/4. If we
now condition on B and compute probabilities that reason only about the
second toss, we should not see any difference.
P (B ∩ A) 1/4
P (A|B) = = = 1/2 = P (A).
P (B) 2/4
The information that the first toss is H has no bearing on the probability
that the second toss is H. Morover, if we condition on Ā, for the same
reason, we get

P (A ∩ B̄) 1/4
P (A|B̄) = = = 1/2 = P (A).
P (B̄) 2/4
Definition 3: Independence
Let 0 < P (B) < 1. The event A is called independent of B if

P (A|B) = P (A|B̄).

15
2.1. PROBABILITIES

An equivalent condition for independence is that P (A|B) = P (A).


Thus, the event A has the same probability no matter whether we condition
on B or on B̄. Our intuition tells us that in this case, event B ”has nothing
to do” with the event A. Another example than the coin tosses would be
that one rolls, say, two dice, one is red and one is blue. Event A is ”red die
gives six” and event B is ”blue die gives even number”. Then, clearly, A
is independent of B and vice versa. However, the following example shows
that independence can also occur between events that are highly related.

Example 14: Independent Events


Consider the case A = Ω. The event Ω is independent of any event B
with 0 < P (B) < 1 since
P (Ω∩B) P (Ω∩B̄)
P (Ω|B) = P (B) =1= P (B)
= P (Ω|B̄).

Assume now that Ω = {1, 2, . . . , 6}, A = {5, 6}, and B = {2, 4, 6}. Then
P (A∩B) P ({6}) 1
P (A|B) = P (B) = P ({2,4,6}) =2· 6 = 1/3

and
P (A∩B̄) P ({5}) 1
P (A|B̄) = P (B)
= P ({1,3,5}) =2· 6 = 1/3,

which means that A = {5, 6} is independent of B = {2, 4, 6}.

Note that one can use the alternative (and equivalent) condition

P (A ∩ B) = P (A) · P (B)

to define independence. In particular, this condition can also be checked if


P (B) = 0.
In summary, if P (B) > 0 we have

(1) P (A|B) = P (A|B̄)


(2) P (A ∩ B) = P (A) · P (B)
(3) P (A|B) = P (A)

From (2) we see that independence is a symmetric property. So, if P (A) > 0,
further equivalent conditions can be derived where the roles of A and B are
reversed. In other words: If the events A and B are independent then
• A and B̄ are independent,

• Ā and B are independent,

• Ā and B̄ are independent.

16
2.1. PROBABILITIES

Remark:
In information theory, the negative log probability I(A) := − log P (A) of
an event A is interpreted as the information content or level of surprise
(the smaller P (A), the larger − log P (A)). For two independent events, the
information content of the combined event simply adds up since
I(A ∩ B) = − log P (A ∩ B)
= − log(P (A) · P (B))
= − log P (A) − log P (B)
= I(A) + I(B)
Later, we will take a closer look at information contents, entropy and other
concepts from information theory.

Pairwise and Mutual Independence


In the case of more then two events, we distinguish pairwise independence
and mutual independence. For events A, B, C, we say that A, B, and C
are pairwise independent if event A is independent of event B, event A
independent of event C, and B independent of event C.

Question:
How can the above definition be extended to n events?

Example 15: Pairwise independent events

We consider a roulette wheel with 36 numbers colored in red (R) or black


(B) according to the following pattern:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18
R R R R R B B B B B R R R R B B B B
36 35 34 33 32 31 30 29 28 27 26 25 24 23 22 21 20 19

Now, let A and B denote the event that a spin of the wheel yields a red
number and an even number, respectively. In addition, let C be the event
that the number is smaller than 19. We check whether A, B, and C
are pairwise independent. We have P (A) = P (B) = 1/2 since half of
the numbers are red and also, half of the numbers are even. Moreover,
P (C) = 1/2 since the numbers from 1 to 18 are the first half of all the
numbers. From the above table, we see that

A ∩ B = {2, 4, 10, 12, 24, 26, 32, 34, 36},


A ∩ C = {1, 2, 3, 4, 5, 10, 11, 12, 13},
B ∩ C = {2, 4, 6, 8, 10, 12, 14, 16, 18},

17
2.2. COMBINATORICS

and thus
P (A ∩ B) = 9/36 = 1/4 = 1/2 · 1/2 = P (A) · P (B),
P (A ∩ C) = 9/36 = 1/4 = 1/2 · 1/2 = P (A) · P (C),
P (B ∩ C) = 9/36 = 1/4 = 1/2 · 1/2 = P (B) · P (C),

which shows pairwise independence.

In the above example, it is the case that

P (A ∩ B ∩ C) = 4/36 = 1/9 6= 1/8 = P (A) · P (B) · P (C).

Hence, a stronger notion of independence is necessary, to describe the rela-


tion between events where the probability of the intersection can be deter-
mined by multiplying the probabilities of all individual events.
Definition 4: Mutual Independence
Let K be a finite set of indices. We say that the events Ak , k ∈ K, are
mutually independent if and only if for all L ⊂ K
Y
P (∩`∈L A` ) = P (A` ).
`∈L

Example 16: Mutual independent events

In the example above with three events, we considered P (A ∩ B), P (A ∩


C), P (B ∩ C) to check for pairwise independence. To determine whether
A, B, and C are mutually independent, we additionally need check if

P (A ∩ B ∩ C) = P (A) · P (B) · P (C),

which is not the case as already stated above.

Question:
Does pairwise independence follow from mutual independence?

2.2 Combinatorics
To determine the probability of certain events, combinatorial considerations
are often useful. We first focus on the elementary events, that is, all sin-
gleton sets {ω} where ω is an outcome. The union of these events forms
Ω and in many chance experiments, all elementary events have the same
probability. This is the case, for instance, in dice games, lottery, and other

18
2.2. COMBINATORICS

games of chance. Assume now that we are interested in a certain event


A = {ω1 , ω2 . . . , ωk } which consists of k outcomes. Then, if all outcomes
in Ω have equal probability, the probability P (A) is given as the ratio of
the number nA of outcomes favorable to A to the number of all possible
outcomes N ,

nA |outcomes favorable to A|
P (A) = = .
N |all possible outcomes|

This probability is also called Laplace probability. If an experiment has two


parts, and there are n1 ways for the first part of the experiment to happen,
and n2 ways for the second part to happen, then there are n1 · n2 ways for
the whole experiment to happen (multiplication principle or rule of prod-
uct). This generalizes in the obvious way to more than two parts. The
examples presented in the sequel are inspired by those mentioned here [3].

Example 17: Laplace Probability


We flip a coin and roll a six-sided die. The coin can land in two ways (heads
or tails). The die can land in six ways (showing 1, 2, 3, 4, 5, or 6). The total
number of possible outcomes in the sample space is therefore 2 × 6 = 12.
If we assume that the outcomes are equally likely, the probability for any
single outcome, such as the coin landing heads and the die showing 4, is
therefore 1/12.
In a different experiment, we flip a coin six times and roll a six-sided die.
The total number of possible outcomes in the sample space is therefore
2 × 2 × 2 × 2 × 2 × 2 × 6 = 26 × 6 = 384. Let A be the event that exactly
one flip lands heads and the die shows 1. The number of outcomes in A is
6, since there are six possibilities for which flip is a head and only one way
for the die to show 1. Hence, if we assume that the outcomes are equally
likely, P (A) = 6/384.

Many chance experiments can be mapped to one of the following urn prob-
lems, where we have an urn containing n distinguishable balls. We draw a
ball from the urn k times. We can do this in two ways: we might replace the
ball drawn each time before drawing the next ball, or we might not replace
the ball (in which case k cannot be bigger than n). We may also consider
the order of balls drawn to matter, or we may consider the draws to be
unordered.

Drawing with Replacement, Ordered Result


Since we replace the balls drawn, each draw can pick any of the n balls.
Since we draw k times, the multiplication principle says that the number of
possible outcomes is n multiplied by itself k times, which is nk .

19
2.2. COMBINATORICS

Example 18: With Replacement, Ordered Result


We draw with replacement two times from an urn containing three balls -
red, green, and blue. There are 32 = 9 possible outcomes: RR, RG, RB,
GR, GG, GB, BR, BG, BB.

Question:
Suppose we have a neural network with k layers, where each layer has n
nodes. How many ways are there to traverse the network if all layers are
fully connected?

Question:
In tic tac toe there are 9 fields, which can either be empty or filled with an
X or an O. How many possible field configurations are there if we ignore
the rules of the game (i.e., we count all assignments of the three symbols to
the nine fields, regardless of whether they could occur in an actual game)?

Drawing without replacement, ordered result


When we don’t replace the balls, the number of possible ball choices goes
down by one after each draw. The multiplication principle then says that
the total number of possible outcomes with k draws from an urn with n
balls is
n!
n(n − 1)(n − 2)...(n − k + 1) = .
(n − k)!
Example 19: Drawing without replacement, ordered result

We draw without replacement two times from an urn containing three balls
- red, green, and blue. There are 3 × 2 = 6 possible ordered outcomes: RG,
RB, GR, GB, BR, BG.

If k = n, we are drawing some permutation of the balls. From the formula

above, the number of possible permutations is n! (remembering that 0! = 1).

Drawing without replacement, unordered result


If we draw k balls without replacement, and don’t look at the order of the
balls drawn, the number of possible results of the experiment decreases,
compared to the result above when we do look at the order. There are k!
ways of ordering k balls (the number of permutations of k items), so the
result with ordering over-counts by this factor. So, dividing the number of
outcomes with ordering by this factor, we get that the number of ways of

20
2.2. COMBINATORICS

drawing k balls without replacement ignoring order from an urn with n balls
is  
n(n − 1)(n − 2) . . . (n − k + 1) n! n
= = .
k! (n − k)!k! k
This is called “n choose k”. Note that
       
n n n n
= and = = 1.
k n−k 0 n

Example 20: Drawing without replacement, unordered result

We draw without replacement two times from an urn containing three balls
- red, green, and blue. There are (3 × 2)/(1 × 2) = 3 possible unordered
outcomes: {R,G}, {R,B}, {G,B}.

21
2.2. COMBINATORICS

Example 21: Drawing without replacement, unordered result


A popular german lottery game is “6 aus 49”, where six balls are drawn
from an urn that contains 49 enumerated balls. In total there are
 
49
= 13 983 816
6

possible outcomes.

Drawing with replacement, unordered result


We can use the result above to find how many possible outcomes there are
when drawing k balls with replacement, when we don’t look at the order the
balls are drawn in. Since we don’t look at the order, all that matters is how
many times each of the n balls in the urn are drawn. We can represent a
count as a sequence of Os, with the number of Os being equal to the count.
We can represent the counts of how many times each of the n balls were
drawn by putting together the sequences of Os representing the counts for
each ball, separating them with Xs. (We choose some order for the n balls;
it doesn’t matter which order.)
For example, if n = 3, with the balls labelled red, green, and blue, and
k = 6, one possible outcome is 2 red, 1 green, and 3 blue. Ordering the
balls as red, green, blue, these counts can be represented by the sequence
OOXOXOOO.
Every set of counts will correspond to a sequence of k + (n − 1) Xs and
Os, in which the number of Xs is exactly n − 1, and every such sequence
will correspond to a set of counts. The correspondence is one-to-one, so we
can count the number of outcomes of the experiment by counting how many
sequences there are of length k + n − 1 with n − 1 of the positions being
occupied by Xs.
The number of ways of putting n − 1 Xs down in a sequence of length
k + n − 1 is the same as the number of ways of choosing n − 1 balls without
replacement from an urn with k + n − 1, ignoring order. We figured that
out above that this is k + n − 1 choose n − 1, i.e.,
 
k+n−1
.
n−1

This is the same as the number of ways of choosing places for the k Os out
of the k + n − 1 positions, which is k + n − 1 choose k.

22
2.2. COMBINATORICS

Example 22: Drawing with replacement, unordered result


We draw with replacement two times from an urn containing three balls
- red, green, and blue. This gives us two Os and two Xs. For instance,
OOXX means that both balls are red, XOXO represents one green and one
blue ball, etc. There are  
2+3−1
=6
2
possible outcomes if we ignore order: {R,R}, {R,G}, {R,B}, {G,G},
{G,B}, {B,B}.

Example 23: Processor jobs


A computer has 6 processors. It is regularly used to run jobs of 4 kinds.
It always runs 6 jobs at a time, so that all the processors will be used, but
there won’t be any processor contention between jobs. The performance of
the computer may depend on what kinds of jobs it is running (e.g., it may
go slowly if two jobs that both access the disk a lot are running). We’re
therefore interested in how many possible job mixes there are, since we may
need to evaluate performance for each job mix.
We can treat this as a problem where we draw k = 6 balls (jobs) with
replacement from an urn with n = 4 balls (kinds of jobs), and we care
only about the numbers of jobs of each kind (there’s no order to jobs). The
answer is therefore    
6+4−1 9
= = 84
4−1 3
possible job mixes.

The above models can be used to calculate Laplace probabilities. Recall


that if we assume equally-likely outcomes, for an event A,

nA |outcomes favorable to A|
P (A) = = .
N |all possible outcomes|

Both, the number of outcomes favorable to A and the number of all possible
outcomes may be computed using one of the formulas above.

Example 24: Birthday problem


A famous probability problem is to find how likely it is that, at a party with
n people, at least two people have the same birthday. Let A be the event
that two or more of the n people have the same birthday. Ac is the event
that all birthdays are distinct. We’ll find P (Ac ), and then get P (A) as
1 − P (Ac ).

23
2.2. COMBINATORICS

We assume equal probability for each day of the year. Hence, P (Ac ) =
#(Ac )/#(S). The number of outcomes in the sample space is #(S) = 365n ,
the same as the number of ways of drawing n balls with replacement from an
urn with 365 balls, paying attention to the order. Using the multiplication
principle, the number of outcomes with no birthdays on the same day is
#(Ac ) = 365 · 364 . . . (365 − n + 1), which is the same as the number of
ways of drawing n balls from an urn with 365 balls without replacement,
paying attention to the order. We use these numbers to compute P (A) =
1 − #(Ac )/#(S). For example, if n = 5, then P (A) = 2.7%.
The following Python code computes the probabilties of the birthday problem
(while preventing overflow problems due to large integers).
1 p l t . p l o t ( [ 1 − np . p r o d ( [ ( 3 6 5 − i ) / 365
2 f o r i in range (n) ] )
3 f o r n in range (100) ] )

Try it out for different n!

Indistinguishable Balls
So far, we assumed that the urn contains n distinguishable balls. But what
if the urn contains r red balls and n − r blue balls and balls of the same
color cannot be distinguished?

Question:
Consider the case of drawing multiple balls (possibly more than one ball
of the same color) with replacement, where balls with the same color can
not be distinguished. Assume that we do not care for the order in which
the balls are drawn. Does the number of possible outcomes differ from the
case with only one ball of each color? Why (not)?
Is your claim also true for the probabilities of the outcomes?

Example 25: With Replacement


Consider an urn with n = 9 balls, where 5 are red and 4 are blue. For k
draws the number N of possible outcomes is 2 · 2 · . . . · 2 = 2k . However, the
probability of, say, two times red and once blue for k = 3 is (5/9)2 · (4/9)
in the ordered case and 3 · (5/9)2 · (4/9) in the unordered case as there are
three possible combinations.

Without replacement things are a little bit more complicated. However, for
small examples it is possible to manually count the possible outcomes:

24
2.2. COMBINATORICS

Example 26: Without Replacement


Consider an urn with n = 9 balls, where 5 are red and 4 are blue.
For k draws the number N of possible outcomes, ignoring the order, is
k 1 2 3 4 5 6 7 8 9
N 2 3 4 5 5 4 3 2 1

If we now consider more than two different colors or care for the order,
things get even more complicated. Luckily, if we are only interested in
the probability of the outcomes (and not in the total number of possible
outcomes), we can use probability trees to easily calculate the probabilities.

Example 27: Without Replacement


Consider again an urn with n = 9 balls, where 5 are red and 4 are blue.
Assume we draw three times from the urn without replacement (k = 3).
The probability of a certain event in the case with considering the order
but without replacement can be obtained from the following probability tree
that lists all possible cases of the three draws.

5 4
9 9

R B

1 1 5 3
2 2 8 8

R B R B

3 4 4 3 4 3 5 2
7 7 7 7 7 7 7 7

R B R B R B R B

To get the probability of a certain event, we traverse the tree along the
path that corresponds to the event (from top to bottom) and multiply the
respective probabilities. For example, for the event RBR we get 59 · 21 · 47 .
If we do not care for the order, we go through all suitable paths and add
the resulting probabilities. For {R, R, B}, i.e., two red and one blue ball in
arbitrary order, we have the paths RRB, RBR and BRR which yield the
total probability 59 · 21 · 47 + 59 · 21 · 47 + 49 · 58 · 47 .
Note that if we consider the case with replacement, every edge that leads
to R has the probability 59 and every edge that leads to B the probability 49 .

25
2.2. COMBINATORICS

26
Chapter 3

Discrete Random Variables


In this chapter, we discuss one of the most important concepts of the course:
random variables. To keep things simple, we first stick to discrete random
variables, which take values in a finite or countably infinite support set.
The relationship to the previously explained idea of having a theoretical
population and a random sample is the following: We will assume that our
population corresponds to a certain probability distribution p and a sample
can then be described as a random variable that follows this distribution.
The distribution p is usually unknown and our goal is to approximate it
(e.g. with some theoretical distribution) and we describe samples drawn
from the population as random variables with distribution p. For example,
let x1 , . . . , x20 be the age of 20 first-year computer science students. This
is the information that we have at hand and the theoretical population
is the set of all first-year computer science students. We approximate it
with, say, a binomial distribution1 . Hence, if X is the random variable that
follows the corresponding binomial distribution, we assume that x1 , . . . , x20
are realizations of X.

Chapter learning objectives:

• understand the mathematical construction underlying discrete random


variables in all its depths

• define appropriate discrete random variables for a variety of problems

• computation of basic properties of discrete random variables: expec-


tation, variance, etc.

• combination and transformation of discrete random variables

• kowledge of the most important discrete probability distributions


1
Later in the course, you will learn how to infer the corresponding parameters of the
binomial distribution.

27
3.1. DISCRETE RANDOM VARIABLES AND PROBABILITY
DISTRIBUTIONS

3.1 Discrete Random Variables and Probability


Distributions
A random variable is used to represent
an outcome of an experiment. Tech-
ω3 X X(Ω)
Ω nically, the fact that this variable is
ω2
X
“random” (that is, it takes values with
ω1 a certain probability), is realized by
using a mapping (as illustrated in the
figure on the left). But why are ran-
dom variables useful if we have sample spaces, events and probabilities?
Often, it is difficult to explicitly define a probability space for a specific
problem. For example, consider the experiment where we throw a die 100
times. Then, we are probably not only interested in the frequencies of the
pips, but also in events such as the sum of the pips being equal to 85. Such
transformations become simple when we work with random variables.
Definition 5: Discrete Random Variable
Let (Ω, 2Ω , P ) be a discrete probability space. A function
X:Ω→R
with
ω 7→ X(ω) ∈ R
is called a discrete (real-valued) random variable on (Ω, 2Ω , P ).

The value x := X(ω) is called realization of the random variable.


Now, that we mapped the outcome of an experiment to R, we need to
assign probabilities to the subsets of X(Ω) = {a ∈ R | ∃ω ∈ Ω : X(ω) = a}.
Since we already have probabilities for events A ⊆ Ω, we define a function
PX : 2X(Ω) → [0, 1] such that, for A ∈ 2X(Ω) ,
PX (A) := P (X −1 (A)) = P ({ω ∈ Ω | X(ω) ∈ A}).
Then (X(Ω), 2X(Ω) , PX ) is a discrete probability space.

Example 28: Rolling two dice


We want to define a random variable for the sum of the two numbers on
the two dice. Let Ω = {(ω1 , ω2 ) | ω1 , ω2 , ∈ {1, . . . , 6}} and X : Ω → R is
such that X(ω1 , ω2 ) = ω1 + ω2 .
Then, the probability of having the number 10 is
1
PX ({10}) = P (X −1 ({10})) = P ({(ω1 , ω2 ) ∈ Ω | X(ω1 , ω2 ) = 10})
= P ({(4, 6)}) + P ({(6, 4)}) + P ({(5, 5)}) = 1/12.

In the sequel, we use the following “shortcuts“ to refer to subset of Ω:

28
3.1. DISCRETE RANDOM VARIABLES AND PROBABILITY
DISTRIBUTIONS

0.2
1

0.15 0.8

0.6
0.1

0.4

0.05
0.2

0
2 3 4 5 6 7 8 9 10 11 12 0
2 3 4 5 6 7 8 9 10 11 12

Figure 3.1: Discrete probability dis- Figure 3.2: Cumulative probability


tribution (rolling two dice). (rolling two dice).

• “X = a” stands for the set {ω ∈ Ω | X(ω) = a}

• “X ≤ a” stands for the set {ω ∈ Ω | X(ω) ≤ a}

• “X < a” . . .

Thus, for instance,


X
P (X ≤ a) = P ({ω ∈ Ω | X(ω) ≤ a}) = P (X = c).
c∈X(Ω),c≤a

The function f : X(Ω) → [0, 1] with f (a) = P (X = a) is called the discrete


probability distribution of X (or probability mass function). It tells us
the probability of each possible event of the form “X = a”. Of similar
importance is the following function related to a random variable X.
Definition 6: Cumulative Probability Distribution
Let X be a discrete (real-valued) random variable. The function F : R →
[0, 1] with
X
F (x) := P (X ≤ x) = P (X = a).
a∈X(Ω),a≤x

is called the cumulative probability distribution of X.

Example 29: Rolling two dice

Assume that X is defined as in Example 28. The discrete probability


distribution and the cumulative probability distribution of X are shown in
Figures 3.1 - 3.2, respectively.

The cumulative probability distribution is monotonically increasing and,


as x → ∞, F (x) approaches 1. The most common discrete probability
distributions will be discussed later.

29
3.2. COMBINING AND TRANSFORMING RANDOM VARIABLES

(! !2 )
ω1 , ω X(Ω)

!1
ω
Z(Ω)
(ω1 , ω2 ) ω1
!1 + ω
ω !2

!2
ω ω1 +ω2

ω2
Y (Ω)

Figure 3.3: Sum of two random variables (rolling two dice).

3.2 Combining and Transforming Random Vari-


ables
It is possible to combine random variables on the same sample space Ω using
operations such as +, −, ·, etc. For instance, we define Z = X + Y as a
random variable on Ω by

Z(ω) = X(ω) + Y (ω) for all ω ∈ Ω.

The probability P (Z = z) = P (X + Y = z) is then well-defined since


X
P (Z = z) = P (ω ∈ Ω | Z(ω) = z) = P ({ω}).
ω∈Ω,X(ω)+Y (ω)=z

Moreover,
X X X
P (Z = z) = P ({ω}) = P (Ω) = 1.
z∈Z(Ω) z∈Z(Ω) ω∈Ω,X(ω)+Y (ω)=z

Example 30: Rolling two dice

Assume that Ω = {1, 2, . . . , 6}2 and X, Y : Ω → R are such that, for


(ω1 , ω2 ) ∈ Ω,
X(ω1 , ω2 ) = ω1 , Y (ω1 , ω2 ) = ω2 .
Then for Z = X + Y we have
1
Z(ω1 , ω2 ) = X(ω1 , ω2 ) + Y (ω1 , ω2 ) = ω1 + ω2

30
3.2. COMBINING AND TRANSFORMING RANDOM VARIABLES

and (see illustration in Fig. 3.3)


X X
P (Z = z) = P ({(ω1 , ω2 )}) = P ({(ω1 , ω2 )}).
(ω1 ,ω2 )∈Ω, (ω1 ,ω2 )∈Ω,
X(ω1 ,ω2 )+Y (ω1 ,ω2 )=z ω1 +ω2 =z

Remark:

In the sequel, we often work with combinations of only two random variables.
However, all these properties and definitions carry over in a straightforward
way to more than two, i.e. a finite number of random variables. Exceptions
or additional concepts for the case of more than two random variables are
explicitly mentioned (e.g. for independence we have pairwise and mutual
independence).

General Transformations of X

In general, we can consider a function g : R → R and tranform a discrete


random variable X : Ω → R with probability mass function fX by defining
Y : Ω → R with Y (ω) := g(X(ω)) for all ω ∈ Ω. Then, if the inverse of g
exists, the probability mass function fY of Y is obtained as follows:

fY (y) = P (Y = y) = P (g(X) = y)
= P ({ω ∈ Ω | g(X(ω)) = y})
= P ({ω ∈ Ω | X(ω) = g −1 (y)})
= P (X = g −1 (y))
= fX (g −1 (y))

Examples for transformations:

X −1
2 · X, , eX
2

Example 31: Diameter and volume of a soccer ball

As manufactured products, soccer balls are subject to small fluctuations


in size. Suppose that in order to be sold a soccer ball must meet a certain
volume criteria. Unfortunately the manufacturer is only able to measure
the diameter of the ball with an accuracy up to 1 mm, which is a discrete
random variable denoted by D. Since diameter and volume are directly
related, we can transform the original random variable D to a new RV

31
3.2. COMBINING AND TRANSFORMING RANDOM VARIABLES

V , which describes the volume via


 3
4 D
V = π .
3 2

Hence, for a concrete realization D(ω), the volume can be directly deter-
mined by the above transformation.

In the sequel, we will often use transformations to shift or scale random


variables. More examples will be discussed in the context of continuous
random variables.

3.2.1 Joint Probability Distribution and Independence


Now, having two random variables on the same probability space, we can
define a joint distribution as follows.
Definition 7: Joint Probability Distribution
Let X, Y be discrete (real-valued) random variables on the same probability
space (Ω, 2Ω , P ). The function PX,Y : X(Ω) × Y (Ω) → [0, 1] with

PX,Y (a, b) := P (X = a ∧ Y = b) = P (X −1 ({a}) ∩ Y −1 ({b}))

is called the joint probability distribution of X and Y . The individual dis-


tributions PX and PY are thenP called the marginal distributions and can be
recovered as P (X = a) = b∈Y (Ω) PX,Y (a, b) and similarly for P (Y = b).

Similar to the definition of independence between events we can define when


two random variables are independent.
Definition 8: Independence of Random Variables
Let X, Y be discrete (real-valued) random variables on the same probability
space (Ω, 2Ω , P ). We call X and Y independent iff for all a ∈ X(Ω), b ∈
Y (Ω)
P (X = a ∧ Y = b) = P (X = a) · P (Y = b).

(Equivalently, if P (Y = b) > 0, this is the same as P (X = a | Y = b) =


P (X = a).)

Note that we use ∧ to consider the intersection of two events, i.e. X =


a ∧ Y = b refers to all outcomes that are mapped to a by X and to b by Y .

Example 32: Rolling two dice

Assume that X and Y are defined as in Example 30. Then X and Y are

32
3.3. EXPECTATION AND VARIANCE

independent, since for any a, b ∈ {1, . . . , 6},

P (X = a ∧ Y = b) = P (X −1 ({a}) ∩ Y −1 ({b}))
= P ({(ω1 , ω2 ) ∈ Ω | ω1 = a}
∩{(ω1 , ω2 ) ∈ Ω | ω2 = b})
= P ({(a, b)}) = 1/36

and
P (X = a) · P (Y = b) = P (X −1 ({a})) · P (Y −1 ({b}))
= P ({(ω1 , ω2 ) ∈ Ω | ω1 = a})
·P ({(ω1 , ω2 ) ∈ Ω | ω2 = b})
= 6/36 · 6/36 = 1/36.

We remark that the above definition can be extended to more than two
random variables by formulating conditions for mutual and pairwise inde-
pendence.
If X1 , . . . , Xn are (real-valued) random variables on the same probability
space, then we call X1 , . . . , Xn (mutually) independent iff for all a1 , a2 , . . . , an ∈
R
P (X1 = a1 , . . . , Xn = an ) = P (X1 = a1 ) · . . . · P (Xn = an ).
We call the random variables X1 , . . . , Xn (pairwise) independent iff all pairs
Xi , Xj , i 6= j are independent in the sense of Definition 8.
In the sequel, we will often consider random variables that are independent
and identically distributed (i.i.d.). This means that the considered random
variables all follow the same probability distribution and are (mutually)
independent.

3.3 Expectation and Variance


Besides the probability distribution of a discrete random variable X, other
values related to X are of interest. The most important ones are the expec-
tation and the variance of X which are defined as
P
• E(X) = x∈X(Ω) x · P (X = x)
P
• V (X) = x∈X(Ω) (x − E(X))2 · P (X = x)

respectively. Note that the sum might not converge, in which case the
expectation/variance does not exist. It is common to write E(X) = ∞ or
V (X) p= ∞ in such a case. The standard deviation σX of X is given by
σX = V (X). Note that V (X) (and σX ) are always non-negative.

33
3.3. EXPECTATION AND VARIANCE

Example 33: Exp. and Var. of Special Random Variables

Let c ∈ R. Assume X(ω) = c for all ω ∈ Ω. Then E(X) = c and


V (X) = 0.

Let A ⊆ Ω and define a random variable IA where



1 if ω ∈ A,
IA (ω) =
0 otherwise.

Then one can compute that E(IA ) = 1 · P (A) + 0 · P (Ā) = P (A) and
V (IA ) = P (A)(1 − P (A)).

Often the expectation of a random variable is denoted by the Greek letter


µ while the variance is often denoted by a σ 2 and σ denotes the standard
deviation. The latter measures the variation using the same ”units” as X
and µ while σ 2 measures in squared units.

One might think that E[X] is some-


thing like ”the best guess for X”
that we can make. However, note
1 2 3 4 5 6 7 8 9 10 11
that even if X is discrete (e.g. X ∈
{1, 2, . . .}) its expectation may not be Figure 3.4: Bimodal discrete proba-
a valid outcome (e.g. 2.6) since it is a bility distribution.
real number. In addition, it might be
a bad guess if the probability distribution of X is, for instance, bimodal as
illustrated in Fig. 3.4.

Example 34:

We consider first a standard fair die with 6 sides. For such a die the
P
6
expectation is given by E(X) = k · 16 = 3.5 and the variance is V (X) =
k=1
P
6
2 1
(k − 3.5) · 6 = 2.9. Now let us consider the second fair die where
k=1
number on sides 3 and 4 are substituted by 1 and 6 (so the die has sides
1, 2, 1, 6, 5, 6). For this die the expectation is 3.5 as before but the variance
is much larger, V (Y ) = 4.9. Distribution of both dice are shown on
the figure below (blue colors) together with the corresponding standard
deviations (green colors give expectation +/− standard deviation).

34
3.3. EXPECTATION AND VARIANCE

E(X) E(Y )

1 2 3 4 5 6 1 2 3 4 5 6

3.3.1 Properties of the Expectation and Variance Operator


We are now able to state some properties of the expectation and variance.

Properties of the expectation:


Let X, Y be discrete random variables on the same probability space and
assume that E(X) and E(Y ) exist. Then, for a, b ∈ R
P
1. E(X) = ω∈Ω X(ω)P ({ω}),

2. E(a · X + b) = a · E(X) + b,

3. E(X + Y ) = E(X) + E(Y ).

4. If X and Y are independent, then E(X · Y ) = E(X) · E(Y ).

Proof. (1.)
X
E(X) = x · P (X = x) (definition of E(X))
x∈X(Ω)
X
= x · P ({ω | X(ω) = x}) (definition of P (X = x))
x∈X(Ω)
X X
= x· P ({ω}) (2nd axiom of Def. 1)
x∈X(Ω) ω∈Ω∧X(ω)=x
X
= X(ω)P ({ω}) (X is a function)
ω∈Ω

The remaining proofs are left as an exercise.

Remark:
We saw that expectation is linear. However, in general expectation does not
factor nicely E(XY ) 6= E(X)E(Y )!

Properties of the variance:


We list the properties of the variance operator without proof: Let X, Y be
discrete random variables on the same probability space and assume that
E(X), E(Y ), V (X), and V (Y ) exist. Then

35
3.3. EXPECTATION AND VARIANCE

1. V (a · X + b) = a2 · V (X) for a, b ∈ R,

2. V (X) = E(X 2 ) − E(X)2

3. V (X + Y ) = V (X) + V (Y ) + 2 · (E(X · Y ) − E(X) · E(Y )).

4. If X and Y are independent, then V (X + Y ) = V (X) + V (Y ).

For more than two independent random variables, we generalize the above
properties of the expectation and variance, we have: If X1 , . . . , Xn are (mu-
tually) independent, then

E(X1 · . . . · Xn ) = E(X1 ) · . . . · E(Xn )

and
V (X1 + . . . + Xn ) = V (X1 ) + . . . + V (Xn ).

3.3.2 Conditional Expectations


In Section 2 we have seen that conditional probabilities are probabilities,
i.e., they fulfil the axioms of probability. Hence, we can consider the cor-
responding ”conditional” probability space and compute also expectations
using conditional probabilities. Suppose we have two discrete random vari-
ables X and Y . Assume that the event X = x has positive probability,
i.e. P (X = x) > 0. Then, then expectation of Y conditioned on X = x is
defined as X
E[Y |X = x] = y · P (Y = y | X = x).
y

In the same way, we could compute the mean of X given Y = y (just switch
the roles of X and Y ).

Example 35: Conditional Expectation

Consider the rolling of two 6-sided fair dice D1 and D2 and two random
variables X and Y where X = value of D1 + D2 and Y = value of D2 .
We compute
P
E[X | Y = 6] = x x · P (X = x | Y = 6)
1
= 6 (7 + 8 + . . . + 12) = 57/6 = 9.5.

Note that this intuitively makes sense, as it is equal to E[value of D1 ] +


E[Y | Y = 6] = 3.5 + 6.

Note that the properties of the expectation carry over to conditional expec-
tations as we are only considering a different probability function (namely,
the conditional probability, which is also a probability in the sense of Def. 1).

36
3.4. HIGHER ORDER MOMENTS, COVARIANCE AND
CORRELATION

We remark that it is possible to construct new random variables by leaving


the value in the condition unspecified, i.e. considering Z = E[X | Y ]. The
conditional expectation as a random variable is very useful and popular
construct, but beyond the scope of this course.

3.4 Higher Order Moments, Covariance and Cor-


relation
In the previous section we have used different operators (e.g. +) to combine
several random variables and considered the expectation of combinations of
random variables. Thus, we already used the following general formula for
the expectation of a function g
X
E(g(X)) = g(x) · P (X = x).
x

(Again the sum might not converge, in which case E(g(X)) does not exist.)
Note that this formula can also be applied if g is not a one-to-one function,
e.g. if g(x) = x2 . In the case g(x) = xi , we call E(X i ) the i-th moment of
X. Note that the moments of a distribution are related to its skewness and
its kurtosis. The former characterizes the degree of asymmetry while the
latter characterizes the flatness or peakedness of the distribution.

3.4.1 Covariance and Correlation


While expectation, variance and standard deviation describe properties of a
single random variable, the covariance and the correlation measure the level
of dependence between variables:

• The covariance of two random variables X and Y is defined by


COV (X, Y ) = E([X − E(X)][Y − E(Y )]).

• Assuming V AR(X), V AR(Y ) > 0, the Pearson correlation coefficient


of X and Y is a value between -1 and 1 and defined by
COV (X, Y )
COR(X, Y ) = p p .
V (X) V (Y )

The covariance is the expected product of the deviations of X and Y from


their respective means. The correlation can be seen as a scaled version of
the covariance having the same sign. Note that both are symmetric, i.e.
COV (X, Y ) = COV (Y, X)
and
COR(X, Y ) = COR(Y, X).

37
3.4. HIGHER ORDER MOMENTS, COVARIANCE AND
CORRELATION

If X and Y are positively correlated, COR(X, Y ) >


0, large values of X (positive deviations from
Y
E(X)) generally correspond to large values of
Y . Similarly, small values of X imply small
values of Y . If X and Y are negatively corre- X
lated, COR(X, Y ) < 0, large values of X gener-
ally correspond to small values of Y (example is Figure 3.5: Positively
given in Fig. 3.5). If X and Y are uncorrelated, correlated X and Y .
COR(X, Y ) = 0, we know that there is no linear
dependency between X and Y . This does, however, not imply that they are
independent since other forms of dependence are possible as the following
example shows.

Example 36: Uncorrelated Variables

Assume that X is a discrete random variable such that P (X = −1) =


P (X = 0) = P (X = 1) = 31 and Y = X 2 . Then E(X) = 0 and E(Y ) =
2/3. Thus,
P
COV (X, Y ) = x,y P (X = x ∧ Y = y)(x − E(X))(y − E(Y ))
= 13 (−1 − 0)(1 − 23 ) + 13 (1 − 0)(1 − 23 )
= 0

Thus, X and Y are uncorrelated even though Y is a function of X (the


strongest form of dependence). Clearly, X and Y are not independent
since, for instance,
1 1 2 2
P (X = 1 ∧ Y = 1) = 6= P (X = 1)P (Y = 1) = · = .
3 3 3 9

Next, we list some important properties of the covariance. Let X and Y be


two discrete random variables.

• COV (X, Y ) = COV (Y, X)

• COV (X, Y ) = E(X · Y ) − E(X) · E(Y )

• COV (X, X) = V AR(X)

• If X and Y are independent, then COV (X, Y ) = 0.

Finally, we note that since the correlation is scaled and


has the property −1 ≤ COR(X, Y ) ≤ 1, we have a max- Y
imal positive (negative) correlation at COR(X, Y ) = 1
(COR(X, Y ) = −1). In this case all possible realizations
for (X, Y ) lie on a straight line. Thus, given X = x the X

38
3.4. HIGHER ORDER MOMENTS, COVARIANCE AND
CORRELATION

value Y = y is fixed by the line. The figure on the right


illustrates this case.
Note that other correlation tests exist such as the Spearman rank correlation,
which tests whether one variable is a monotone function of the other.

3.4.2 Moment Generating Function

When we discussed transformations of random variables, one example was


the transformation g(X) = eX . This is a very popular transformation as it
can be turned into the powerful tool of moment generating functions, which
are useful for deriving equations for the moments of many popular random
variables.
The moment generating function for a random variable X is defined by

 
tX (tX)2 (tX)3
MX (t) = E[etX ] = E 1 + + + + ... , t ∈ R.
1! 2! 3!

Here, the variable t is basically a placeholder that determines the degree k


(from (tX)k ). For the moment generating function, the coefficient of tk is
the kth moment E[X k ]:

" ∞
# ∞
X X k tk X tk
tX
MX (t) = E[e ]=E = E[X k ]
k! k!
k=0 k=0

If MX (t) exists on an open interval around t = 0, then the kth moment


equals the kth derivative of the moment generating function, evaluated at
t = 0:

dk MX (t)
E[X k ] =
dtk t=0

Another usefull equation for computing moment generating functions is

X
MX (t) = P (X = x) · etx .
x∈X(Ω)

39
3.4. HIGHER ORDER MOMENTS, COVARIANCE AND
CORRELATION

Proof.

MX (t) = E[etX ] (definition of MX (t))


"∞ #
X X k tk
=E (definition of etX )
k!
k=0

X tk
= E[X k ] (linearity of E)
k!
k=0
 

X X
 tk
= xk · P (X = x) (definition of E[g(X)])
k!
k=0 x∈X(Ω)

!
X X tk
= xk · P (X = x) (reorder sums)
k!
x∈X(Ω) k=0
X ∞
X (xt)k
= P (X = x) (rearrange)
k!
x∈X(Ω) k=0
X P∞ xk
= P (X = x) · etx (ex = k=0 k! )
x∈X(Ω)

Example 37: MX (t) of the Bernoulli Distribution

A random variable X is Bernoulli-distributed with parameter p ∈ [0, 1] if


P (X = 1) = p and P (X = 0) = 1−p. Thus, X is the indicator variable of
a single trial that has outcome “success” (encoded as 1) with probability
p and “failure” (encoded as 0) with probability 1 − p. We consider its
moment generating function:

MX (t) = E[etX ]
= P (X = 0) · et·0 + P (X = 1) · et·1
= (1 − p)e0 + pet
= 1 − p + pet

Next, we consider

dMX (t) d2 MX (t)


E[X] = = pe0 = p and E[X 2 ] = = pe0 = p.
dt t=0 d2 t t=0

Thus, E[X] = p and V [X] = E[X 2 ] − E[X]2 = p − p2 = p(1 − p).

An important property of the moment generating function is that if two

40
3.5. IMPORTANT DISCRETE PROBABILITY DISTRIBUTIONS

distributions have the same moment generating function, then they have
the same distribution, i.e. MX = MY implies PX = PY .

3.5 Important Discrete Probability Distributions


Let us now consider some important discrete probability distributions.

• Bernoulli distribution: Consider again an experiment where there


are only two possible outcomes (e.g. flipping a coin). This is called
a Bernoulli trail. Let us use 0 (failure) and 1 (success) for the two
possible outcomes. The probability of 1 is P (X = 1) = p for some
fixed p ∈ [0, 1] and P (X = 0) = 1 − p. We say that X is Bernoulli
distributed with parameter p, shorthand

X ∼ Bernoulli(p).

As we have seen above, the mean of the Bernoulli distribution is


E(X) = 0(1 − p) + 1(p) = p.

• Binomial distribution: We consider a sequence of n independent


Bernoulli trials and count the number X of successes. Then X follows
a binomial distribution, i.e. for k ∈ {0, 1, . . . , n}
 
n k
P (X = k) = p (1 − p)n−k
k

where the first factor counts all possible orderings for k successes
among n trials, the second one gives the probability that in the in-
dependent trials we have k successes and the last factor describes the
probability of (the remaining) n−k failures. Note that one can express
X as the sum of n independent Bernoulli variables X1 , . . . , Xn :

X = X1 + . . . + Xn .

Thus, E(X) = E(X1 ) + . . . + E(Xn ) = p + . . . + p = n · p. Please use


the Wikipedia link to view the plot of the distribution and learn more
about the Binomial distribution.

• Geometric Distribution: Consider again an experiment where the


probability of an event A is P (A) = p > 0. We repeat this experiment
independently until A occurs for the first time. Let X be the random
variable that describes the number of trials until A occurs for the first
time, so X takes values in {1, 2, . . .}. Then X is called geometrically
distributed. We have P (X = i) = (1 − p)i−1 · p and E(X) = 1/p. The
variance of X is V (X) = (1 − p)/p2 .

41
3.5. IMPORTANT DISCRETE PROBABILITY DISTRIBUTIONS

Note: An alternative convention defines the geometric distribution as


the number of failures before the first success, with support {0, 1, . . .}
and E(X) = (1 − p)/p. Both conventions appear in the literature, so
always check which one is used.
Please use the Wikipedia link to view the plot of the distribution and
learn more about the geometric distribution.
• Poisson Distribution Consider a call center where on average µ = 6
calls per minute arrive. Let X be the random variable that represents
k
the number of calls in the next minute and assume P (X = k) = µk! e−µ .
Then X is called Poisson distributed and has expectation

X ∞
X µk−1
µk −µ
E(X) = k· e = µe−µ = µe−µ eµ = µ
k! (k − 1)!
k=0 k=1

and variance V (X) = µ (without proof). Please use the Wikipedia link
to view the plot of the distribution and learn more about the Poisson
distribution.

Note that the Poisson distribution is the limit of the binomial distribution,
when n is large and p is small. This can be shown as follows:
Let µ = np. Starting from a binomial distribution we then get
 
n k
P (X = k) = p (1 − p)n−k
k
n · (n − 1) · . . . · (n − k + 1)  µ k  µ n−k
= 1−
k! n n
µ n
n · (n − 1) · . . . · (n − k + 1) µk 1 − n
= 
nk k! 1 − µ k
n
µk −µ
≈ e ,
k!
where we used the following limits for fixed k as n → ∞ (with µ = np held
k n
constant): n·(n−1)·...·(n−k+1)
nk
→ 1, 1 − nµ → 1 and 1 − nµ → e−µ .

Example 38: Binomial vs Poisson distribution

Suppose a chip is defective in 10% of the time. You have 10 chips and
the number of defective chips is described by the RVs Y or X. We would
like to know the probability that no more than 1 chips are defective. For
the binomial distribution Y ∼ binomial(10, 0.1) we get
P (Y ≤ 1) = P (Y = 0) + P (Y = 1)
   
10 10
= 0.10 0.910 + 0.11 0.99 ≈ 0.7361.
0 1

42
3.5. IMPORTANT DISCRETE PROBABILITY DISTRIBUTIONS

For the Poisson distribution X ∼ P oisson(µ = 10 · 0.1 = 1) on the other


hand, we get

P (X ≤ 1) = P (X = 0) + P (X = 1)
10 11
= e−1 + e−1 ≈ 0.7358.
0! 1!
Note that even for a moderate size of n = 10 the results are already very
similar.

43
3.5. IMPORTANT DISCRETE PROBABILITY DISTRIBUTIONS

44
Chapter 4

Continuous Random
Variables

So far, we considered only random variables that take discrete values such
as 0, 1, 2, . . .. However, many chance experiments can only be described by
means of continuous random variables. For instance, consider a dartboard
with one foot in radius. The experiment consists of throwing a dart at the
board and the outcome is the point at which it hits. Then, any point in

S = {(x, y) ∈ R2 : x2 + y 2 ≤ 1}

is a possible outcome. Hence, the sample space is uncountably large. In this


case, the defintion of random variables (in fact even that of a probability)
gets more complicated.

Chapter learning objectives:


• understand the mathematical construction underlying continuous ran-
dom variables in all its depths
• define appropriate continuous random variables for a variety of prob-
lems
• computation of basic properties of discrete random variables: expec-
tation, variance, etc.
• kowledge of the most important continuous probability distributions
• extension of all these concepts to the multivariate case

4.1 σ-algebras
Consider a chance experiment where the sample space Ω contains uncount-
ably many elements. In this case, assigning probabilities to all elements in

45
4.1. σ-ALGEBRAS

2Ω poses problems. Assume, for instance, we randomly choose a real num-


ber in [0, 1]. If all numbers are equally likely to occur, we have to assign
probability zero to each since their “sum” must be one (note that the sum
over uncountably many nonzero values is undefined). Instead of giving a
new definition of probabilities, we restrict ourselves to a set of events, called
σ-algebra, for which we can define probabilities as in the discrete case.
If we define a random variable as a function X : Ω → R, we want to reason
about the probability of events such as X = x for any x ∈ R or a < X ≤ b
for any interval (a, b]. From the properties of probabilities (see Def. 1), we
then know the probability of the disjoint union of countably many of such
sets as well as the probability of the complement of such a set1 .
Definition 9: σ-algebra
A set F ⊆ 2Ω is called a σ-algebra if
1. Ω ∈ F,

2. A ∈ F implies Ā ∈ F.

3. If A1 , A2 , . . . ∈ F is a sequence of sets then

A1 ∪ A2 ∪ . . . ∈ F.

The most important example of a σ-algebra, which is needed in the sequel,


is the σ-algebra that is generated by a set E ⊆ 2Ω . We define the smallest
σ-algebra that contains E by

σ(E) := ∩{F ⊃ E : F is a σ-algebra}.

Example 39: Generated σ-algebra

For simplicity, we consider the finite set Ω = {1, 2, . . . , 6}. Let E =


{{2, 6}, {5, 6}}. Then

σ(E) = {{2, 6}, {5, 6}, {1, 3, 4, 5}, {1, 2, 3, 4}, {1, 2, 3, 4, 5}, {6}, . . .}.

If Ω is finite, the idea is to construct a sequence of sets in an iterative


fashion, where we start with E and obtain the next set from the previous
one by joining and complementing elements of the current set. If no
new element can be constructed by union or complement operations, the
current set equals σ(E).

Net we assume Ω = R and

E = {(a, b] : a, b ∈ R, a ≤ b}.
1
From the two conditions in Def. 1, one can easily derive that P (Ā) = 1 − P (A).

46
4.1. σ-ALGEBRAS

The σ-algebra B := σ(E) is called the Borel algebra on the reals. It contains
all subsets, called Borel sets, of 2R that can be obtained from E by countable
union and complement operations. Note that also intervals of the form
(a, b) or [a, b] are Borel sets. A similar construction is possible for Ω = Rn .
Intuitively, the Borel sets are those sets, for which we can assign a “volume”,
“area” or “size”. Subsets of Rn that are not Borel sets are only of theoretical
interest since for practical applications they are not of importance.
We are now able to give a more general definition of a probability space (i.e.
also for sample sets with uncountably many elements).

Definition 10: Probability Space


Let Ω be a nonempty set and let F ⊆ 2Ω be a σ-algebra. A probability space
is a triple (Ω, F, P ) where the probability measure P : F → [0, 1] is such
that

• P (Ω) = 1 and,

• if A1 , A2 , . . . ∈ F is a sequence of pairwise disjoint sets, then

P (A1 ∪ A2 ∪ . . .) = P (A1 ) + P (A2 ) + . . . .

Example 40: Discrete Probability Space

Any discrete probability space (Ω, 2Ω , P ) fulfills Def. 10, since 2Ω is a


σ-algebra and P is as in Def. 1.

Example 41: One-dimensional Interval

Let Ω = R, F = B and a < b, a, b ∈ R. Consider a probability measure


P : F → [0, 1] such that
( min(y,b)−max(x,a)
b−a if max(x, a) < min(y, b)
P ({ω ∈ Ω | x < ω ≤ y}) =
0 otherwise.

Similar to how we extended the set E = {(a, b] : a, b ∈ R, a ≤ b} to a


σ-algebra, we can show that if P is a probability measure and defined as
above for all half-open intervals, its value for the remaining sets in B is
uniquely determined.

47
4.2. CONTINUOUS RANDOM VARIABLES

4.2 Continuous Random Variables


ω3
Definition 11: Real-valued Ran- Ω X X(Ω)
ω2
dom Variable
X A
Let (Ω, F, P ) be a probability space. ω1
A real-valued random variable on
(Ω, F, P ) is a function X : Ω → R
such that for all A ∈ B Figure 4.1: Is the inverse image of
A w.r.t. X an element of F?.
X −1 (A) = {ω ∈ Ω | X(ω) ∈ A} ∈ F.

The above definition ensures that if we want to know the probability that X
is in some Borel set A, we can consider the inverse image of A with respect
to X, for which we know its probability.
Clearly, we can define a probability measure PX : B → [0, 1] by setting
PX (A) := P ({ω | X(ω) ∈ A}) = P (X −1 (A)) and use similar notations as in
the discrete case (e.g. P (a < X ≤ b)).
Definition 12: Cumulative Probability Distribution
Let X be a real-valued random variable on (Ω, F, P ). The function F :
R → [0, 1] with x 7→ F (x) := P (X ≤ x) is called the cumulative probability
distribution (CDF) of X.

Example 42: One-dimensional Interval

Assume that X is a randomly chosen point in the interval [a, b] and


(Ω, F, P ) is as in Ex. 41. Then
 y−a
 if y ∈ [a, b],
 b−a
F (y) = P (X ≤ y) = 1 if y > b,


0 otherwise.

We call a random variable X on (Ω, F, P ) discrete if X(Ω) is a discrete


set (finite or countably infinite). We call X continuous if X(Ω) contains
uncountably many elements and there exists a non-negative and integrable
1
function f : R → R≥0 , called density, with
Z x
F (x) = P (X ≤ x) = f (y) dy.
−∞

Clearly, Z ∞
f (y) dy = 1
−∞

48
4.2. CONTINUOUS RANDOM VARIABLES

f (y) F (x)

1
b−a

y x
a b a b

Figure 4.2: Density and cumulative probability distribution of a random


variable X.

since F (∞) = 1, but note that P (X = y) 6= f (y).

Example 43: One-dimensional Interval

Assume that X is as in Ex. 42. Then X is a continuous random variable


since the (constant) function f with
( 1
b−a if y ∈ [a, b],
f (y) =
0 otherwise.

is the density of X. We can verify this by calculating


Z x Z x  x
1 y x−a
F (x) = f (y) dy = dy = = = P (X ≤ x)
−∞ a b−a b−a a b−a

for x ∈ [a, b], F (x) = 0 for x < a, and F (x) = 1 for x > b. Figure 4.2
shows a plot of the functions f and F .

The distribution in the above example is called the (continuous) uniform


distribution.
We summarize the most important properties of the CDF F (x) and the
density f (x) in the following table:

F (x) f (x)
domain R R
codomain [0, 1] R≥0
monotonicity increasing not necessarily
Rx
relation = −∞ 1 f (y)dy = dFdx(x)

The expectation and variance of a continuous random variable X with den-

49
4.3. IMPORTANT CONTINUOUS DISTRIBUTIONS

sity f are defined as


Z ∞ Z ∞
E(X) = x · f (x) dx and V (X) = (x − E(X))2 · f (x) dx.
−∞ −∞

Note that these integrals may not exist, in which case the expectation/vari-
ance is undefined.

What from the discrete case carries directly over to the continuous
case?
Many results for discrete random variables carry over to the continuous
setting. For instance, properties of the expectation and variance (e.g. E(X +
Y ) = E(X) + E(Y )) or results concerning the combination/transformation
of random variables, joint distributions, stochastic independence, etc. Also,
we define conditional expectation, covariance, correlation, and higher-order
moments for continuous random variables in an equivalent way. All results
of the previous sections carry over to the continuous case (as long as the
corresponding integrals exist). We do not repeat these definitions here as
the only things that change compared to the discrete case are that we replace
the sum by an integral and the probability of a value x by the density f (x).

4.3 Important Continuous Distributions


Besides the uniform distribution that we have seen above, there are some
other important continuous distributions:

Uniform Distribution
We already described the uniform distribution in the examples above. In
summary we get the following properties if X is uniformly distributed on
the interval (a, b) (denoted by X ∼ U (a, b)):

• Its density is constant on the interval (a, b), that is, f (x) = 1/(b − a)
for x ∈ (a, b) and f (x) = 0 otherwise.

• If we want to know the probability that X falls into an subinterval of


(a, b) of length d we get

P (X ∈ (c, c + d)) = d/(b − a), (c, c + d) ⊆ (a, b)

which means that it is proportional to the length d (and independent


of c, cf. Figure 4.2).

• The expectation is given by the midpoint of the interval E(X) =


(a + b)/2.

50
4.3. IMPORTANT CONTINUOUS DISTRIBUTIONS

Exponential Distribution
Let λ > 0. We say that a continuous random variable X is exponentially
distributed with parameter λ (denoted by X ∼ Exp(λ)) if the density of X
is such that, for t ∈ R,
(
λ · e−λt if t ≥ 0,
f (t) =
0 otherwise.
The cumulative probability distribution of X is then given by
( Rx Rx −λt dt = 1 − e−λx
−∞ f (t) dt = 0 λ · e if x ≥ 0,
F (x) =
0 otherwise.
The exponential distribution is often used to describe waiting times or in-
terarrival times since in many real-world systems where a sequence of rare
events plays an important role, the assumption that the time between these
events is exponentially distributed is used. For instance, the time of ra-
dioactive decay is exponentially distributed. This is related to the fact that
these events are assumed to occur spontaneously and that the exponential
distribution has the so-called memoryless property:
Assume that the random variable X describes a waiting time and we already
know that X ≥ t and ask for the probability that X ≥ t + h (we have to
wait for additional h time units). Then, if X is exponentially distributed we
can verify that
P (X ≥ t + h | X ≥ t) = P (X ≥ h), t, h > 0.
The exponential distribution is the only continuous distribution that is mem-
oryless. In fact, it is possible to derive from the memoryless property that
the distribution of X must be exponential. Similarly, one can show that the
geometric distribution is the only discrete distribution that is memoryless.

Example 44: Bus waiting paradox


Assume that Bob arrives at a bus stop where buses arrive on average
every 10 min. The interarrival times X of the buses are exponentially
distributed. Bob asks one of the people waiting there how long he is already
waiting. Whatever the answer of that person is (e.g. t = 20 min or 1 min),
it does not change the probability that Bob has to wait for h minutes for
the next bus.

Normal distribution
The normal distribution is one of the most important continuous distri-
butions since it naturally arises when we sum up many random variables.

51
4.3. IMPORTANT CONTINUOUS DISTRIBUTIONS

Therefore, measurement errors and fluctuations are often (approximately)


normally distributed. Also measurements of the length or size of certain
objects give often normally distributed results.
A normally distributed random variable X with mean µ and variance σ 2
(denoted by X ∼ N (µ, σ 2 )) has the density

1 −(x − µ)2
f (x) = √ exp ,
σ 2π 2σ 2
The density of the standardized normal distribution with mean µ = 0 and
variance σ 2 = 1 is then defined as
1 −x2
φ(x) = √ e 2

and its cumulative distribution function is
Z x
1 −z 2
Φ(x) = √ e 2 dz.
−∞ 2π

Assume we are interested in the probability that X ∼ N (µ, σ 2 ) falls into a


certain interval (a, b). Since the solution of the above integral is difficult to
compute, we can standardise X and use the fact that we have a table with
the values for Φ.
Definition 13: Standardized Random Variables.
For any random variable X with finite E(X) = µ and finite V (X) = σ 2 , we
define its standardised version
X −µ
X∗ = .
σ
Then E(X ∗ ) = 0 and V (X ∗ ) = 1 since

 
∗ X −µ E(X) − µ
E(X ) = E = =0
σ σ
and  
∗ X −µ V (X)
V (X ) = V = = 1.
σ σ2
For X ∼ N (µ, σ 2 ) this means that we can find the cumulative distribution
of X ∗ in a standard normal table and use it to compute that of X. Using
X ∗ = X−µ
σ ⇐⇒ X = X ∗ σ + µ we get
x−µ
P (X ≤ x) = P (X ∗ σ + µ ≤ x) = P (X ∗ ≤ ).
σ

Example 45: Computing Normal Probabilities

52
4.3. IMPORTANT CONTINUOUS DISTRIBUTIONS

The average salary of a full-time employee in Germany is µ = 45240


Euros per year (in 2017, brutto). Assuming a normal distribution and a
standard deviation of σ = 10000 Euros, we can compute the proportion
of employees that earn between, say 30 000 and 50 000 Euros per year as
follows:
 
P (30000 < X < 50000) = P 30000−µ σ < X−µ
σ < 50000−µ
σ

= P (−1.524 < X ∗ < 0.476)

= Φ(0.476) − Φ(−1.524)

≈ 0.68439 − 0.06426 = 0.6203.

Next, we can compute the income limit of the poorest 3% of the employees,
i.e. find x such that P (X < x) = 0.03:

P (X < x) = P X < x−µ σ = Φ( x−µ
σ ) = 0.03

Thus,

x = µ + σΦ−1 (0.03) = 45240 + 10000 · (−1.88) = 26440

where we used that Φ(−1.88) = 0.03.

Finally, we note that if a random variable X is normally distributed with


expectation µ̂ and variance σ̂ 2 , written X ∼ N (µ̂, σ̂ 2 ), then, for a, b ∈ R,
aX + b has distribution N (aµ̂ + b, (aσ̂)2 ).

Gamma distribution
The Gamma distribution generalizes a number of other distributions (e.g.
exponential) and has a rather flexible shape compared to other continuous
distributions. Therefore it can be used to fit very different types of data.
A random variable X is Gamma distributed with parameters α and β (writ-
ten X ∼ Gamma(α, β)) if its density is given by
 
1 1 α−1 −x
f (x) = x exp , x > 0, α > 0, β > 0
Γ(α) β α β

and f (x) = 0 whenever x ≤ 0. Thus, it assigns positive probability only to


positive values of X and also the parameters α (for the shape) and β (for the
scaĺe) must be positive. The gamma function Γ(α) is used as a normalization
factor of the density. Often α is a positive integer and then Γ(α) = (α − 1)!
which means that the density (and its integral) is easy to evaluate (for plot of

53
4.4. MULTIVARIATE RANDOM VARIABLES

the Gamma densities please visit [Link]


commons/e/e6/Gamma distribution [Link]).
How do the parameters α and β influence the mean and the variance of
the distribution? It can be shown that the mean of X ∼ Gamma(α, β) is
µ = αβ and the variance is σ 2 = αβ 2 . Further interesting properties of the
Gamma distribution are:
• For large α and β, the distribution is very similar to a normal distri-
bution (with mean µ = αβ and variance σ 2 = αβ 2 )
• If α = 1 and β = 1/λ, then X ∼ Exp(λ).
• If α = ν/2 and β = 2 then the distribution is a chi-squared distribution
with ν degrees of freedom (where ν is a positive integer). In that case,
we have the distribution of a sum of the squares of ν independent
standard normal random variables. The chi-squared distribution is
often used for hypothesis testing (goodness of fit tests).

4.4 Multivariate Random Variables


So far we only considered one dimensional RVs, i.e. the outcome of the
mapping X leads to a one dimensional result. Imagine now, that during a
measurement you obtain multiple quantities at once. For example during a
health screening, you obtain the age, height, weight, . . . for each participant.
In order to properly describe the outcome, we need the following definition:
Definition 14: Multivariate Random Variable
Let Ω be a sample space and

X1 , X2 , . . . , Xn : Ω → R

n one dimensional RVs. The vector X = (X1 , X2 , . . . , Xn ) is then called


an n-dimensional random variable or an n-dimensional random vector.
A vector of realizations of the corresponding one dimensional RVs x =
(x1 , x2 , . . . , xn ) is called realization of X.

If the sample space of X = (X1 , . . . , Xn ) is countable, we call X a discrete


n-dimensional RV. Since each entry Xi in the vector of our multivariate
RV X is a one dimensional RV, we can adopt Definition 7 for the joint
distribution of the vector entries. The generalization to more than two
entries is straightforward.
In a similar fashion to the one dimensional case, we can also define the
cumulative probability distribution for multivariate RVs.
Definition 15: Cumulative Probability Distribution
The cumulative distribution function F : Rn → [0, 1] of a RV X = (X1 , . . . , Xn )
is defined as
F (x) = P (X1 ≤ x1 , . . . , Xn ≤ Xn ).

54
4.4. MULTIVARIATE RANDOM VARIABLES

In an analogous way we lift properties such as expectation and covariance

to the multi-dimensional case.


To ease the notation, we consider the case n = 2 in the following and call
the two one dimensional RVs X and Y , i.e., our multivariate RV is (X, Y ).
Suppose now, that there are N possible realizations x1 , . . . , xN of X and K
possible realizations y1 , . . . , yK of Y . The possible values i ∈ {1, . . . , N } and
j ∈ {1, . . . , K} then form the index sets I and J. Note that N and K might
be (countably) infinite and that N and K might in general be different.
Besides the probability of a certain realization P (X = xi , Y = yj ) obtained
from the joint distribution, other interesting quantities are P (X = xi ) and
P (Y = yj ), i.e, we keep one of the variables fixed, while the other still
changes. This leads to the following definition.
Definition 16: Marginal Distribution
Given a RV (X, Y ) and its joint distribution P (X = x, Y = y), the marginal
distribution of X is defined as
X
P (X = xi ) = P (X = xi , Y = yj )
j∈J

for all i ∈ I. In a similar fashion the marginal distribution of Y is given by


X
P (Y = yj ) = P (X = xi , Y = yj )
i∈I

for all j ∈ J. For more than two dimensions, we sum over all variables but
one. For two dimensional RV a well-arranged representation for the joint

and marginal distribution is a so-called contingency table.

Example 46: Contingency Table

Consider the RV (X, Y ) where X has two possible realizations x1 and x2


and Y has the three possible realizations y1 , y2 , y3 .
Y X x1 x2
y1 0.1 0.25 0.35
y2 0.1 0.4 0.5
y3 0.05 0.1 0.15
0.25 0.75
The values form the joint distribution are in the inside of the table, for
example P (X = x2 , Y = y2 ) = 0.4. The sum of each column gives the
marginal distribution for X, the sum of each row for Y . For example
P (X = x1 ) = 0.1 + 0.1 + 0.05 = 0.25 and P (Y = y3 ) = 0.05 + 0.1 = 0.15.
Note that the sum of all entries within the table, as well as the sum of
the respective marginal probabilities is 1.

55
4.4. MULTIVARIATE RANDOM VARIABLES

In the continuous case (i.e. Ω is not discrete), we have a continuous n-


dimensional RV with an n-dimensional joint density f : Rn → R≥0 and,
compared to the discrete case, we replace summations by integrals in the
definitions given above.

Example 47: Bivariate continuous distribution

Consider the function


(
c · exp(−(x + 2y)), if 0 ≤ x, y ≤ ∞,
f (x, y) =
0, else.

How do we have to choose c such that f becomes a density? We determine


c such that the integral becomes 1, i.e.
Z Z Z ∞Z ∞
1= f (x, y)dxdy = c · exp(−(x + 2y))dxdy
0 0
Z ∞
c
=c e−2y dy = ⇒ c = 2.
0 2

We can now compute the marginal distributions


Z ∞ Z ∞
−x
fX (x) = 2 exp(−(x + 2y))dy = e 2e−2y dy = e−x ,
Z0 ∞ 0
Z ∞
−2y
fY (y) = 2 exp(−(x + 2y))dx = 2e e−x dx = 2e−2y ,
0 0

and the expectations


Z ∞Z ∞ Z ∞
E(X) = 2x exp(−(x + 2y))dxdy = xe−x dx = 1
Z0 ∞ Z0 ∞ Z 0∞
E(Y ) = 2y exp(−(x + 2y))dxdy = 2ye−2y dy = 0.5
0 0 0
Z ∞Z ∞
E(XY ) = 2xy exp(−(x + 2y))dxdy
Z0 ∞ 0 Z ∞
−x
= xe dx 2ye−2y dy = E(X) · E(Y ) = 0.5.
0 0

From this it is obvious that

Cov(X, Y ) = 0

as expected, since X and Y are independent, because f (x, y) = fX (x) ·


fY (y).

56
Chapter 5

Generation of Random
Variates

The generation of independent


samples, i.e. random variates,
that follow a specific distribu-
tion is an important part of the
popular sampling methods in statis-
tics. All common programming
languages provide methods to
generate pseudo random num-
bers, i.e. numbers that look as
if they random, but are gener-
ated using a deterministic algo-
rithm, a so-called pseudo ran-
dom number generator. On the
right, you see the results of the
Mersenne Twister, a popular pseudo
random number generator. It
generates numbers that are uni-
formly distributed on the area (0, 1) × (0, 1). The screenshot is taken from
the WOLFRAM Demonstrations Project.
In the sequel, we will not discuss algorithms to generate pseudo random
numbers but focus on the transformation of uniformly distributed numbers
such that the results are variates that follow a certain distribution.

Chapter learning objectives

• how to exploit relationships between random variables for the genera-


tion of random variates

• understand the concept of the inverse transform method for the con-

57
5.1. GENERATING DISCRETE RANDOM VARIATES

tinuous and the discrete case


• understand the concept of rejection-based approaches
• implemention of concrete random variate sampling algorithms

5.1 Generating Discrete Random Variates


Assume we want to generate random numbers that are realizations of a
discrete random variable X. For discrete distributions there are two possi-
bilities.

a) We either use a specific relationship between X and U ∼ U (0, 1) or


b) We split the interval (0, 1) into disjoint intervals and do a case distinc-
tion as explained below.

Example 48: Generating Bernoulli distributed numbers

Let U ∼ U (0, 1) and p ∈ (0, 1). We define



1 if U < p,
X=
0 if U ≥ p.

Then P (X = 1) = P (U < p) = p and P (X = 0) = 1−p. Thus, X follows


the Bernoulli distribution with probability p.
Algorithmically, we first generate U , set X accordingly and return X.
A Python code for this is (after setting the value for p):

1 i m p o r t numpy a s np
2 U = np . random . r a n d ( )
3 X = i n t (U < p )

Example 49: Generating binomially distributed numbers

Let U1 , . . . , Un be independent and U (0, 1)-distributed. Moreover, let p ∈


(0, 1). We define 
1 if Ui < p,
Xi =
0 if Ui ≥ p,
Pn
and Y = i=1 Xi . Then X1 , X2 , . . . , Xn are Bernoulli distributed (with
parameter p) and thus their sum Y is binomially distributed with param-
eters n and p (see Section 3.5).
Algorithmically, P we generate U1 , . . . , Un , set X1 , X2 , . . . , Xn accordingly
and return Y = ni=1 Xi .
A Python code for this is (after setting the values for n and p):

58
5.1. GENERATING DISCRETE RANDOM VARIATES

1 i m p o r t numpy a s np
2 U = np . random . r a n d ( n )
3 X = np . sum (U < p )

Example 50: Generating geometrically distributed numbers

Let U1 , U2 , . . . be independent and U (0, 1)-distributed. Moreover, let p ∈


(0, 1). We define
X = min{Ui < p}.
i

Then, clearly, P (X = 1) = p and it is easy to see that P (X = k) =


(1 − p)k−1 p.
Algorithmically, we set X = 1 and perform a while-loop in which we
generate U and check whether U > p (loop condition). If this is the case
we set X = X + 1. After the loop we return X. A Python code for this
is (after setting the value for p):

1 i m p o r t numpy a s np
2 X = 1
3 w h i l e ( np . random . r a n d ( ) > p ) :
4 X += 1
5 X

5.1.1 Interval Method


The splitting of the interval (0, 1) can be done for any discrete random
variable X. Assume that X takes the values x0 , x1 , . . . and P (X = x0 ) = p0 ,
P (X = x1 ) = p1 , . . . . Then we divide the interval (0, 1) into the subintervals

A0 = (0, p0 ),
A1 = [p0 , p0 + p1 ),
A2 = [p0 + p1 , p0 + p1 + p2 ).
...

Then we let X = xi if U falls into interval Ai . Note that since Ai has length
pi the probability that U falls into Ai is exactly pi . Hence X has the desired
distribution.
Algorithmically, the naive approach to find the index i with

p0 + . . . + pi−1 ≤ U < p0 + . . . + pi

is performed as follows:
1. initialize c = p0 for the cumulative value,

2. iteratively set c = c + pk in the k-th execution of a while loop that is


only executed if U ≥ c.

59
5.2. INVERSE TRANSFORM METHOD

If U < c holds after i iterations, then X = xi .


A Python code for this is (when pi is stored in p(i) as in Python array indices
start at 0):
1 i m p o r t numpy a s np
2 c = p [ 0 ] # set to p_0
3 i = 0
4 U = np . random . r a n d ( )
5 w h i l e (U >= c ) :
6 i += 1
7 c += p [ i ] # set to p_i and add to c
8 X = i

There is a simple way to improve the efficiency of the above approach. If


we sort the values x0 , x1 , . . . according to their probability p0 , p1 , . . . we can
ensure that large intervals (high probabilities) are considered first. Hence it
is more likely that the number of iterations is small since U is more likely
to fall in one of the first intervals.
In the case that X can take infinitely many values with positive probability,
the intervals can be constructed on the fly. However, care has to be taken in
how the intervals are traversed in order to ensure convergence of the while
loop.
We remark that Python also offers functions for generating random variates
for many well-known continuous and discrete distributions (e.g. package
[Link]).

5.2 Inverse Transform Method


Several approaches exist for the gen-
eration of independent samples that
follow certain probabilities distribu-
tions.
Here, we present the most common
one: inverse transform sampling, also
called the inverse transform method.
It can be used for any distribution
(discrete or continuous) whenever the
inverse F −1 of the cumulative distri-
bution function (CDF) F exists.
The main idea of the method is based
on the following intuition: If we plot the CDF FX of a RV X and randomly
choose a point y on the [0, 1]-y-axis (all points are equally likely), then the
corresponding value FX−1 (y) = x on the x-axis will have distribution FX . In
the above plot (source: Wikipedia), this is illustrated for the exponential
distribution with parameter λ.

60
5.2. INVERSE TRANSFORM METHOD

Example 51: Generating exponentially distributed numbers

We have the CDF F (x) = 1 − e−λx for X ∼ Exp(λ). Thus, y = F (x) is


a number between 0 and 1. Solving for x yields

y := F (x) = 1 − e−λx
⇐⇒ 1 − y = e−λx
⇐⇒ ln(1 − y) = −λx
⇐⇒ x = − λ1 ln(1 − y) = F −1 (y).

Thus, if U is uniformly distributed on (0, 1), then the random variable


X = F −1 (U ) = − λ1 ln(1 − U ) is exponentially distributed with parameter
λ. Moreover, since 1−U has the same distribution as U , the random vari-
able X = − λ1 ln(U ) must also be exponentially distributed with parameter
λ.

Next, we formulate the inverse transform sampling algorithm for any given
CDF FX with inverse FX−1 :

1. Generate U ∼ U (0, 1).

2. Return X = F −1 (U ).

Next, we prove that X = F −1 (U ) has the cumulative probability distribution


F as desired: First observe that

P (X ≤ x) = P (F −1 (U ) ≤ x) = P (U ≤ F (x))

where we used that F is monotone in the last step. Finally, we note that
because U is uniformly distributed, P (U ≤ y) = y for any y ∈ [0, 1] and
in particular for y = F (x). Hence, P (U ≤ F (x)) = F (x) and thus P (X ≤
x) = F (x) holds.

5.2.1 Inverse transform sampling: the discrete case


In Section 5.1.1 we developed an interval method to sample discrete distri-
butions. How does this relate to the idea of inverse transform sampling?
For the cumulative distribution function FX of a discrete random variable
X we define FX−1 (y) = inf{x | FX (x) ≥ y}. But then, inverse transform
sampling means that we look for the interval into which U ∼ U (0, 1) falls on
the y-axis (which corresponds to the while-loop in Section 5.1.1) and choose
the smallest x such that FX (x) ≥ U . Hence, inverse transform sampling for
discrete distributions is identical to the interval method discussed before.

61
5.3. REJECTION SAMPLING

5.3 Rejection Sampling


Random variates can also be generated using acceptance-rejection sampling.
The idea is similar to that of Monte-Carlo methods for computing compli-
cated integrals and is best explained starting with the example of an expo-
nential distribution with parameter λ = 1.

Example 52: Generating exponentially distributed numbers

Assume we want to generate X ∼ fcx) e-


×

Ät
"" ""
=

exp(1). The corresponding density is


f (x) = e−x . Now, consider a rect-
angle, which includes most of f (if
we choose, say, [0, 8] × [0, 1] as illus-
trated on the right). We now repeat-
edly draw uniformly distributed num-
bers U1 ∼ U (0, 1) and U2 ∼ U (0, 8).
0

The numbers U1 and U2 give use points (U2 , U1 ) that are uniformly dis-
tributed within the 2-dimensional rectangle. Whenever the pair (x, y) =
(U2 , U1 ) is below the density f (black points), i.e. y ≤ f (x), we keep x
as a sample and otherwise, we reject x (illustrated as grey pairs (x, y)).
Observe that if x is small, f (x) is large and it is very likely that we keep x
as the point (x, y) falls below the red line. Indeed, it is easy to show that
the samples that we keep have the desired distribution f (x) = e−x (up to
the approximation we made by truncating the x-axis at 8).

For the general method, we first define what the body Bh of a nonnegative
integrable function h on Rd is.
Bh = {(x, y) : x ∈ Rd , 0 ≤ y ≤ h(x)}.
Now, observe the following: if we sample pairs X - h

(X, Y ) that are uniformly distributed on Bh , then µ


X has a density that is proportional to h. This is
simply because in those regions where h is large,
we will have more samples (X, Y ) and only few
where h is small (see illustration on the right).
×

Moreover, if we generate U ∼ U (0, 1) and X, in-


dependently of U such that its density is propor-
tional to h, then (X, U h(X)) is uniformly distributed on Bh . Intuitively,
this is because in regions where h is large, we have many realizations for X
and we use U to choose one of them (uniformly distributed).
Now, consider the probability density f of the random variable for which

62
5.3. REJECTION SAMPLING

we want to sample variates. We choose another function g that majorizes


f , i.e. that is at least as large as f at all points x ∈ R, g(x) ≥ f (x) for
all x. Let c > 0 be the scaling constant such that 1c g(x) is a density (you
should know how to determine c!). We then generate U ∼ U (0, 1) and X
such that it has density 1c g(x). If U g(X) ≤ f (X) then we keep X. Other-
wise, we reject it. Using pseudocode, we summarize the algorithm as follows:

i) Generate U ∼ U (0, 1).

ii) Generate X, independent of U , distributed according to density 1c g(x).

iii) If U g(X) ≤ f (X) then return X (’accept’); otherwise go back to i)


(’reject’);

Using the arguments before, we note that we sample uniformly points (X, U g(X))
on Bg , but partition Bg into Bf and its complement Bg \ Bf . We accept
whenever we hit the area Bf . Hence, the x-components of our samples are
distributed according to f . Moreover, the acceptance probability is

P (acceptance | X = x) = P (U g(X) ≤ f (X) | X = x)


 
(X)
= P U ≤ fg(X) |X=x
f (x)
= g(x)

Unconditioning and recalling that X has density 1c g(x) yields


Z +∞ Z +∞
f (x) 1 1 1
P (acceptance) = g(x)dx = f (x)dx = .
−∞ g(x) c c −∞ c

Hence, in order to achieve a high probability for acceptance, we must choose


g such that only a small scaling constant c is needed, i.e. g must majorize
f but only slightly such that only little scaling to a density is needed.

Example 53: Generating normally distributed numbers

Assume we want to generate X ∼


N (µ, σ). Since by shifting and scal-
ing, we can express X as X = σZ +
µ, where Z ∼ N (0, 1), it is enough to
concentrate on the standard normal
distribution.

63
5.3. REJECTION SAMPLING

Another simplification is that it is enough, if we can sample the absolute


value |Z| because the density is symmetric. We can simply sample the
sign by choosing + and − each with probability 1/2. (Sample U ∼ U (0, 1)
and choose +Z if U < 1/2 and −Z if U ≥ 1/2.) The density of the
non-negative RV |Z| is
2 2
f (x) = √ e−x /2 , x ≥ 0.

Next, we need a function g that majorizes f and is proportional to a
density, from which we can sample easily. We choose the density r(x) =
e−x because we can use the inverse transform method to sample X ∼
exp(1) and it has the same exponential form as f . Now, we need to scale
up
p r to a x−x
function g that majorizes f . Observe that the ratio f (x)/r(x) =
2 /2
2/π · e has its maximal value at x = 1 (you find this whenp you
set the derivative of f (x)/r(x) to zero!). Then, we determine c = 2e/π
by setting f (1) = c · r(1) because at x = 1 we want that f (x) = g(x) (the
curves coincide at this point). Checking
2 2 p
f (x) = √ e−x /2 ≤ 2e/π · e−x = c · r(x), for all x ≥ 0

we find a function g(x) = c · r(x) that majorizes f and apply the following
rejection algorithm to generate Z ∼ N (0, 1):

i) Generate U1 ∼ U (0, 1).

ii) Generate X exponentially distributed with λ = 1 by first generating


U2 ∼ U (0, 1) and then setting X = − ln(U2 ).

iii) If U1 g(X) ≤ f (X) then return |Z| := X (’accept’); otherwise go


back to i) (’reject’);

iv) Generate U3 ∼ U (0, 1). Set Z := |Z| if U3 < 1/2 and Z := −|Z|
otherwise.
f (X)
Note that the algorithm can be simplified when noting that g(X) =
2
e−(X−1) /2 and thus
f (X) 2 /2
U g(X) ≤ f (X) ⇐⇒ U≤ g(X) = e−(X−1)
⇐⇒ − ln(U ) ≥ (X − 1)2 /2.

Since − ln(U ) is exponentially distributed with rate 1, we simply need two


of those generated in step ii), say X1 and X2 and only have to check
whether X2 ≥ (X1 − 1)2 /2 in step iii) and set |Z| := X1 if so.

64
Chapter 6

Laws of Large Numbers

In this chapter, we will discuss some fundamental results of probability the-


ory: the laws of large numbers. They formalize and generalize what is
intuitive for most of us: if we consider many independent replications of the
same chance experiment, we can, based on the outcomes, approximate the
characteristics of the average outcome. E.g. consider flipping a fair coin
very often. Then, we expect that approximately half of the time we see
heads. More formally, consider n realizations x1 = X(ω1 ), x2 = X(ω2 ), . . . ,
xn = X(ωn ) of a random variable X with finite expectation µ. They could,
for instance, be generated by repeating an experiment n times under equal
conditions. Intuitively, for large n the sample mean approximates µ, i.e.,
n
1X
x̄ = xi ≈ µ.
n
i=1

In this chapter, we discuss important aspects related to this approximation.

Chapter learning objectives

• understand the consequences of the laws of large numbers and the


central limit theorem

• apply Chebyshev’s inequality

• apply the central limit theorem

6.1 Chebyshev’s inequality


In order to determine the quality of the above approximation, we recall
Chebyshev’s Inequality.

65
6.2. WEAK LAW OF LARGE NUMBERS

µ−a µ µ+a
Let X be a random variable with fi-
nite expectation E(X) and finite vari-
ance V (X). Assume that besides E(X) {ω ∈ Ω | |X(ω) − µ| < a}
and V (X), nothing is known about X
(e.g. the cumulative probability distribu-
tion). Our aim is to reason about the
deviation of X from its expectation. Ac-
cording to Chebyshev’s Inequality (see
figure on the right for an illustration), for any a > 0

V (X)
P (|X − E(X)| ≥ a) ≤ . (6.1)
a2
Note that the proof of the theorem is quite straight forward.

6.2 Weak Law of Large Numbers


In order to apply the inequality above, we assume that x1 , x2 , . . . , xn are
realizations of independent random variables X1 , X2 , . . . , Xn , all having the
same distribution as X. Define
n
1X
Zn = Xi .
n
i=1
P
Then E(Zn ) = n1 ni=1 E(Xi ) = µ and, if X has finite variance V (X) = σ 2 ,
n
! n
1X 1 X n · σ2 σ2
V (Zn ) = V Xi = 2 V (Xi ) = = .
n n n2 n
i=1 i=1

(See also page 36 for the properties of the variance operator.) Thus, Eq. (6.1)
gives us, for any  > 0,
σ2
P (|Zn − µ| ≥ ) ≤ .
n · 2
σ2
If  is given, we can make n·2
arbitrarily small by increasing n. Thus, for
any  > 0,
lim P (|Zn − µ| ≥ ) = 0 (or, equivalently, lim P (|Zn − µ| < ) = 1),
n→∞ n→∞

which is known as the weak law of large numbers 1 . Here, the sequence 1
{Zn }n≥1 of random variables converges “weakly” since only the correspond-
ing probabilities converge.

Example 54: Bernoulli Trials


1
Note that, as opposed to the strong law of large numbers, for the weak law of large
numbers pairwise independence of X1 , . . . , Xn is already sufficient and (mutual) indepen-
dence is not necessary.

66
6.3. STRONG LAW OF LARGE NUMBERS

Consider a sequence of n Bernoulli trials and for 1 ≤ j ≤ n, let Xj be


the random variable that is 1 if the j-th trial is success and 0 otherwise.
Moreover, let P (Xj = 1) = p for all j. Then E(X1 ) = E(X2 ) = . . . =
E(Xn ) = p. According to the weak law of large numbers, for any  > 0,
 
Xn
1
P Xj − p <  → 1
n
j=1

as n → ∞.

6.3 Strong Law of Large Numbers


An even stronger result than the weak law of large numbers provides the
strong law of large numbers, which states that
 
P lim |Zn − µ| = 0 = 1.
n→∞

Here, we measure the probability of all outcomes ω ∈ Ω with limn→∞ |Zn (ω)−
µ| = 0 and find that this event occurs with probability one. The strong law
of large numbers implies the weak law of large numbers.
We have seen two laws that match our initial intuition that
n
1X
x̄ = xi ≈ µ.
n
i=1

In the following we will see that even more can be derived about sums of
random variables.

6.4 Central Limit Theorem.


P
Recall that Zn = n1 ni=1 Xi , where the Xi are independent and identically
distributed random variables with finite expectation µ and finite variance
σ 2 > 0. According to the central limit theorem, the distribution of

Zn − E(Zn ) Zn − µ
Zn∗ = p = √ .
V (Zn ) σ/ n

converges to the normal distribution with expectation 0 and variance 1 (stan-


dard normal distribution). Thus, for x, y ∈ R, x < y,
Z y
1 2
lim P (x < Zn∗ < y) = √ e−0.5t dt
n→∞ 2π x

67
6.4. CENTRAL LIMIT THEOREM.

n=4 n=10 n=50


1 1 1
distr ZnStar distr ZnStar distr ZnStar
0.8 normal distr 0.8 normal distr 0.8 normal distr

0.6 0.6 0.6

0.4 0.4 0.4

0.2 0.2 0.2

0 0 0
−4 −2 0 2 4 −4 −2 0 2 4 −4 −2 0 2 4

Figure 6.1: The distribution of Zn∗ approaches the standard normal distri-
bution as n increases.

and Zn is approximately normally distributed with expectation µ and vari-



ance σ 2 /n, since Zn = σ/ n · Zn∗ + µ and Zn∗ ∼ N (0, 1).

Example 55: Approximation of the binomial distribution

We have already seen that for a sequence of n independent Bernoulli trials


the number Sn of successes follows a binomial distribution. We can write
Sn as
Sn = X1 + . . . + Xn
where X1 , . . . , Xn are independent Bernoulli random variables with pa-
rameter p. If n is large, we can approximate the binomial distribution
of Sn by a normal distribution. The central limit theorem tells us that
Zn = n1 Sn is approximately normally distributed with mean p (of the
Bernoulli random variables) and variance p(1−p)n . Therefore Sn must ap-
proximately follow the distribution N (np, np(1−p)). Fig. 6.1 shows a plot
of the cumulative probability distribution of Zn∗ . It can be seen that the
distribution (in yellow) approaches N (0, 1) (red curve) as n increases.

68
Chapter 7

Parameter Estimation

We consider an unknown distribution of which we have a collection of sam-


ples and we know the family of the distribution (e.g. Poisson) but some of
the parameters are unknown and have to be estimated. This problem is an
instance of statistical inference, which differs from the methods of descrip-
tive statistics (see Section 1) in that we do not compute numbers to describe
our data but we infer (i.e., estimate) the parameters of the process (here this
is our theoretical distribution) that is generating our data.
Values that we estimate based on the data that we have are called statistical
estimators. In this chapter, we will consider estimators for the mean of
the distribution and for its standard deviation. They are random variables
computed from a collection of data and are subject to a sampling error.
We will also construct random intervals that contain the true value of the
parameter with high probability.

Chapter learning objectives


• understand the concept of an estimator and important properties of
estimators
• understand and apply popular methods for estimating parameters:
method of moments, maximum likelihood estimation, Bayesian infer-
ence
• determine standard deviations of the corresponding estimators
• know the advantages and disadvantages of the different estimation
methods

Example 56: Photon Count


Consider a telescope observing the photon flux of a certain star during
certain time intervals of fixed length. We assume that the distribution of

69
the number of observed photons in the interval [t, t+h] is constant in t, i.e.
we have the same distribution for all intervals of length h. Assume further
that n independent measurements are performed for certain intervals of
length h, i.e. we have n numbers x1 , . . . , xn for the photon flux within
intervals of length h. Clearly, we will in reality never have exact and
truly independent measurements but let us assume for now that

a) the measurements have been taken over n days for, say, h = 1 hour
every day,

b) the telescope gives a very accurate measurement.

From a) we get approximate independence and from b) we derive that


measurement errors may be neglected. Random photon counts can well be
described by a Poisson distribution with some parameter λ.
We know that the mean of the Poisson distribution is E(X) = λ and the
variance is V (X) = λ as well. Shall we now use the sample mean
n
1X
xi = x̄
n
i=1

to fit E(X) = λ or shall we use the sample variance


n
2 1 X
s = (xi − x̄)2 ?
n−1
i=1

as defined in Section 1.2.3 to fit V (X) = λ? Which of the two will give
us a better estimate for λ? What does ”better” mean in this context? We
give answers to these questions in the next sections.

Recall that we use a lower case ’x’ for concrete realizations (real values!)
of a random variable X. When we discuss independent observations that
follow the distribution of X, then we can take two different views:

View 1: We have random samples of X, i.e. independent realizations


x1 , . . . , xn ∈ R. If we compte the corresponding mean x̄, then x̄ is a real
value.

View 2: We look at the observation process in a more abstract view


and describe the n independent observation samples as random variables
X1 , . . . , Xn with the following properties:

• they are independent because we ensured independence during the


observation process. E.g. assume you consider data for the height of a

70
7.1. METHOD OF MOMENTS

set of individuals. If these individuals are relatives, your independence


assumption is no longer true as this might bias the data towards a
higher or lower height than you have on average in the total population
that you consider.

• the follow the same distribution, namely the distribution of the random
variable X which represents the random variable that describes a single
data point (e.g. height of a person). Hence, the all have, for instance,
the same expectation:

E[X] = E[X1 ] = . . . = E[Xn ].

• any combination or transformation of the Xi is again a random vari-


able. E.g. their mean
n
1X
X̄ = Xi
n
i=1

is a random variable! This is intuitive, if you imagine that you repeat


your observation process several times. Your data might change and
the mean of your data as well. In the sequel, the properties (such as its
mean and variance) of such an estimator will be an important aspect.

7.1 Method of Moments


Given data points x1 , . . . , xn of n independent measurements (view 1!), the
idea of the method of moments is to adjust the moments of the distribution
such that they fit the sample moments of the data. In general, the k-th
sample moment is defined as
n
1X
mk = (xi )k .
n
i=1

Note that we often write x̄ for m1 . For k > 1 it is often a good alternative
to consider the k-th sample central moments
n
1X
m̄k = (xi − x̄)k
n
i=1

since non-central moments may become very large. Note that for the sample
variance two definitions exist, the one for m̄2 and the one for s2 . The
difference is that for m̄2 we divide the sum by n while for s2 we divide it by
n−1. The reason is that if we view the sample variance as an estimator of the
true variance of a random variable X based on the data points X1 , . . . , Xn

71
7.1. METHOD OF MOMENTS

(view 2!) where the Xi are i.i.d. (same distribution as X with mean µ and
variance σ 2 ), then the estimator
n
2 1 X 1X
S = (Xi − X̄)2 , where X̄ = Xi ,
n−1 n
i=1 i

is unbiased (E(S 2 ) = σ 2 ) while


n
2 1X
S̃ = (Xi − X̄)2
n
i=1

is not (E(S̃ 2 ) 6= σ 2 ). We will come back to this issue when we discuss the
pros and cons of the method of moments.

Now we assume that we consider a theoretical distribution with exactly K


unknown parameters θ1 , . . . , θK . Then the theoretical moments µ1 , µ2 , . . .
of the distribution (often called population moments) are functions of these
parameters, i.e. E(X k ) = µk (θ1 , . . . , θK ). For example, the exponential
distribution has a parameter θ1 = λ > 0 that determines the mean E[X] =
µ1 (θ1 ) = λ−1 .
We may also simply write µ instead of µ1 (θ1 , . . . , θK ). Clearly, the central
moments E((X − µ)k ) = µ̄k (θ1 , . . . , θK ), k = 2, 3, . . . are also a function of
θ1 , . . . , θ K .
To determine these parameters we construct an equation system

µ1 (θ̂1 , . . . , θ̂K ) = m1
µ̄2 (θ̂1 , . . . , θ̂K ) = m̄2
...
µ̄K (θ̂1 , . . . , θ̂K ) = m̄K

and solve for θ̂1 , . . . , θ̂K . Note that it is also possible to equate the non-
central moments instead. Also note that in general, it may be necessary to
add equations for higher moments if there is no unique solution for θ̂1 , . . . , θ̂K
with K equations (for instance, because some equations are linearly depen-
dent).

Example 57: Fitting the mean of a Poisson distribution

Revisiting Example 56 we have the equation µ = m1 . We know that


the parameter λ is equal to the theoretical mean µ. Hence, we choose
λ̂ = x̄ = m1 as an estimator for the parameter λ.

72
7.1. METHOD OF MOMENTS

Example 58: Fitting the parameters of a normal distribution

Assume that our data points x1 , . . . , xn are IID realizations of a random


variable X that follows a normal distribution. The normal distribution
has parameters θ1 (equal to the theoretical mean µ) and θ2 (equal to the

standard deviation σ). Also, µ̄2 = θ2 since µ̄2 = σ 2 .

We use the two equations µ = m1 and µ̄2 = m̄2 . Thus, we directly get

as a solution that θ1 = m1 = x̄ (from the first equation) and θ2 = m̄2
(from the second equation).
As a variant, assume that it is known that X has mean µ = 0. Then, for
the single remaining parameter σ we do not set K = 1 and consider only
the first equation µ1 (µ, σ) = m1 since the first moment µ1 (µ, σ) is equal to
the first parameter µ and thus does not give any constraints on σ. Instead
we consider the second equation µ̄2 (µ, σ) = m̄2 since µ̄2 (µ, σ) = σ 2 and

thus σ = m̄2 .

The method of moments has a number of disadvantages compared to the


estimation methods considered in the sequel. The estimators may be biased
and thus may have a systematic estimation error. An example is the bias of
the sample variance

Example 59: Bias of the sample variance

As already discussed above, according to the method of moments, we would


estimate the variance of the true distribution by
n
1X
s̃2 = (xi − x̄)2
n
i=1

instead of using the unbiased estimator


n
2 1 X
s = (xi − x̄)2 .
n−1
i=1

Intuitively, s̃2 underestimates the variance since in its definition we use


an estimator for the mean instead of the true mean and this estimator
must lie in the center of the samples (while the true mean might not) and
therefore the deviation from it is smaller than the deviations from the true
mean. Consider for instance a normal distribution with mean 0 and vari-
ance 100. For a small sample size of, say, n = 3 one might get the data
points x1 = −1.9412, x2 = −21.3836, x3 = −8.3959 using the following
Python command x = [Link](0,100,3). The sample mean
x̄ = −10.5736 is far away from the true mean. Obviously, the deviations

73
7.2. MAXIMUM LIKELIHOOD ESTIMATION

from the true mean 0 are much larger than the deviations from x̄ and thus
s̃2 is only [Link](x) = 65.3717 while [Link](x,ddof=1) = S 2 = 98.0576
is closer to the true standard deviation (note that the Python command
[Link](x) uses the biased estimator!).
If we could consider the deviations from the (usually unknown) true mean
0 the estimator with the factor n1 becomes unbiased:
1 Pn  1 Pn
E n i=1 (Xi − E(X))2 ) = n i=1 E((Xi − E(X))2 )
1
= nn · V (Xi ) = V (X)

The method of moments has the drawback that it does not take into account
information about more moments than needed to get a unique solution for
the equation system, i.e. if we have a single parameter, we fully rely on only
the mean or only the variance to estimate it. For instance, in Example 56
we only use x̄ and do not additionally take into account the variation of the
data, i.e. s2 or s̃2 .
The generalized method of moments does take into account this information
by considering cost functions instead of equating the population moments
of the distribution and the sample moments.

7.2 Maximum Likelihood Estimation


Again, we assume that n independent observations X1 , . . . , Xn are given as
well as the family of the distribution. The maximum likelihood method is
the most common way of estimating parameters of a distribution. Recall
that an estimator is a function of the data. Here, we consider maximum
likelihood estimators.
Let θ be the unknown parameter and let Pθ (·) be the corresponding discrete
probability mass function, i.e. the Xi are all independent and identically
distributed with P (Xi = x) = Pθ (x), x ∈ R. Therefore, the probability of
observing X1 = x1 , . . . , Xn = xn , assuming that θ is the true value, is given
by
Lθ (x1 , . . . , xn ) = Pθ (x1 )Pθ (x2 )...Pθ (xn ).

Lθ (x1 , . . . , xn ) is also called the likelihood of the data and we say that θ̂ is
a maximum likelihood estimator (MLE) of θ if for all possible θ

Lθ (x1 , . . . , xn ) ≤ Lθ̂ (x1 , . . . , xn ).

Note that the MLE of θ is, in general, not unique. Intuitively, the MLE is
the parameter value that “best explains” the observed data.

74
7.2. MAXIMUM LIKELIHOOD ESTIMATION

To maximize the likelihood, we can consider its derivatives w.r.t. θ and find

parameter values where ∂θ Lθ (x1 , . . . , xn ) equals 0 (in some cases such points
do not exist and we find θ̂ at the boundary of the set of possible values for
θ).
Often the maximum of the log-likelihood

ln Lθ (x1 , . . . , xn ) = ln Pθ (x1 ) + ln Pθ (x2 ) + ... + ln Pθ (xn ),

which is equal to that of the likelihood (as the logarithm is strictly mono-
tonically increasing), is easier to compute. More concretely, θ̂ maximizes Lθ
if and only if θ̂ maximizes ln Lθ .

Example 60: MLE of the geometric distribution

Assume we have collected data that represents inter-arrival times of cus-


tomers (in minutes). Since the inter-arrival times are approximately ge-
ometrically distributed (we hypothesize this after looking at certain sum-
mary statistics such as the range of the data, the mean, the skewness,
the variance or a histogram plot of the data), we wish to fit the data to
a geometric distribution with parameter θ, i.e. we want to find a θ that
best explains the data. The likelihood of the data x1 , . . . , xn is in this case
given by

Lθ (x1 , . . . , xn ) = Pθ (x1 )Pθ (x2 )...Pθ (xn )

= θ(1 − θ)x1 · θ(1 − θ)x2 · . . . · θ(1 − θ)xn


n
P
xi
= θn (1 − θ)i=1

since the probability to observe xi is Pθ (xi ) = θ · (1 − θ)xi . Note, that


we use the alternative definition of the geometric distribution here, where
we only count the number of unsuccessful trials until the first success.
Since Lθ (x1 , . . . , xn ) is a function of θ, we shortly write L(θ) and omit
the dependence on the data. Considering the log-likelihood yields
n
X
ln L(θ) = n · ln θ + xi ln(1 − θ)
i−1

We maximize ln L(θ) by setting the derivative to zero:

d ln L(θ) n (−1) X n 1 X
= + xi = − xi
dθ θ 1−θ θ 1−θ
i i

! n 1 X
0= − xi
θ 1−θ
i

75
7.2. MAXIMUM LIKELIHOOD ESTIMATION

X n(1 − θ) 1 1
⇔ xi = ⇔ x̄ = − 1 ⇔ θ̂ =
θ θ x̄ + 1
i

1 P
n
where x̄ = n xi . Next we find that θ̂ is really a maximizer:
i=1

d2 ln L(θ) n 1 X
= − − xi < 0
dθ2 θ2 (1 − θ)2
i

1
Hence, the MLE of θ is θ̂ = x̄+1 . If we use the alternative definition of the
geometric distribution here, where we count the number of unsuccessful
trials and the first success, how would the MLE change? Hint: this means
that we transform X to X + 1.

Note that in the above example, since the mean of the geometric distribu-
tion is (1 − θ)/θ we matched the mean x̄ of the data and the mean of the
distribution. Note that this is not always the case for a maximum likelihood
estimator, i.e. the estimator of the method of moments is not always the
same as the MLE.
Assume now that we have hypothesized a continuous distribution for our
data. In that case, the likelihood function is

Lθ (x1 , . . . , xn ) = fθ (x1 ) · fθ (x2 ) · ... · fθ (xn )

where fθ is the density of the distribution that depends on the parameter


θ. Again, since Lθ (x1 , . . . , xn ) is a function of θ, we shortly write L(θ) and
omit the dependence on the data. An MLE θ̂ of θ is defined as a value that
maximizes L(θ) (over all permissible values of θ).
The reason why the likelihood is a density in the continuous case is that the
probability of a single outcome is zero in the continuous case. Hence one
could consider
Z xi +h
P (xi − h < X < xi + h) = fθ (x)dx
xi −h

for the outcome xi and some very small h. But this quantity is approximately
equal to 2hfθ (xi ) and thus proportional to the factors of the likelihood de-
fined above. Hence, this approach would result in the same MLE.

Example 61: MLE of the exponential distribution


Assume that data were collected on the inter-arrival times for cars in a
drive-up banking facility. Since histograms of the data show exponential
curve, we hypothesize an exponential distribution. Note that some distri-
butions generalize the exponential distribution but have more parameters

76
7.2. MAXIMUM LIKELIHOOD ESTIMATION

(e.g. gamma distribution). They would provide a fit which is at least as


good as the fit of the exponential distribution.
For the exponential distribution the parameter is θ = λ, λ > 0 and the
density is fλ (x) = λe−λx . Hence, the likelihood of the samples x1 , . . . , xn
is
Yn P
−λ xi
−λxi n
L(λ) = λe =λ ·e i
i=1

and the log-likelihood is


X
ln L(λ) = n · ln λ + (−λ xi )
i

and its derivative w.r.t. λ is


d n X
ln L(λ) = − xi
dλ λ
i

Setting the derivative to zero yields

! n X 1 1X
0= − xi ⇔ = xi = x̄
λ λ n
i i

and thus λ = 1/X̄. We find that λ is a maximizer:

d2 n
2
ln L(λ) = − 2 < 0
dλ λ

and finally get λ̂ = 1/x̄.

We remark that for some distributions, the log-likelihood function may not
be useful and also, finding a maximum by setting the derivative to zero is
not always possible. In general, if an analytic solution is not possible, global
optimization methods have to be applied in order to determine, which of
several local optima of the likelihood has the highest value.
We list some important properties of Maximum Likelihood Estimators (some
of the properties require mild “regularity” assumptions, e.g. likelihood func-
tion must be differentiable, support of distribution does not depend on θ.)

• Uniqueness: For most common distributions, the MLE is unique; that


is, L(θ̂) is strictly greater than L(θ) for any other value of θ.

• Asymptotic unbiasedness: θ̂ is a function of the samples x1 , ..., xn and


thus a random variable if Xi = xi is not given, i ∈ {1, 2, . . . , n}.
Assume now that X1 , ..., Xn have distribution Pθ and θ̂n is the MLE

77
7.2. MAXIMUM LIKELIHOOD ESTIMATION

of θ based on the observations X1 , ..., Xn . Then

lim E[θ̂n ] = θ.
n→∞

Therefore, we say that θ̂n is asymptotically unbiased.

• Invariance: If θ̂ is an MLE of θ and if g is a one-to-one function1 , then


g(θ̂) is an MLE of g(θ). To see that this holds, we note that


z }| {
L(θ) = L(g −1 (g(θ))

d = g(θ̂) (where g(θ)


are both maximized by θ̂, so g(θ) d is the MLE of
g(θ) interpreted as the unknown parameter).

Example 62: Invariance

Assume that x1 , x2 , ..., xn are samples from a Bernoulli distribution with


parameter θ = p. Then the likelihood is
n
Y P P
xi n− xi
[θxi + (1 − θ)χ=0 (xi )] = θ i · (1 − θ) i

i=1

1 if x = 0,
where χ=0 (x) =
0 otherwise.

Hence the log-likelihood is


X X
ln L(θ) = xi ln(θ) + (n − xi ) ln(1 − θ)
i i

and setting its derivative to zero yields


∂ X X !
L(θ) = xi /θ − (n − xi )/(1 − θ) = 0
∂θ
i i
P P P P
⇒ (1 − θ̂) xi = θ̂(n − xi ) ⇔ xi = n · θ̂ ⇔ θ̂ = xi /n
i
P i i i
Next we find that θ̂ = xi /n =: x̄ is a maximizer:
i

∂2 X X
2
L(θ) = − xi /θ2 − (n − xi )/(1 − θ)2 < 0 for any θ ∈ [0, 1]
∂θ
i i

1
A function is one-to-one if every element of the range of the function corresponds to
exactly one element of the domain.

78
7.2. MAXIMUM LIKELIHOOD ESTIMATION

Thus, θ̂ = x̄.
Next, we assume that the variance g(θ) = θ(1 − θ) is our parameter (note
that in the permissible range of θ, g is one-to-one). Define L̃ such that
L̃ (g(θ)) = L(θ), i.e., we consider the same values for the likelihood but
the likelihood function is different as it is a function in g(θ). Setting the
derivative to zero yields
∂ ∂ ∂θ !
L̃(g) = L̃(g) · = 0.
∂g ∂θ |{z} ∂g
=L(θ)

∂ ∂
Since ∂θ L (g(θ)) |θ=θ̂ = 0 it holds that ∂g L̃(g) = 0 if we choose ĝ = g(θ̂) =
x̄(1 − x̄).

• Asymptotically normally distributed: For n → ∞, θ̂n converges in


distribution to N (θ, δ(θ)) (normal distribution with mean θ and vari-
ance δ(θ)), where θ̂n is our estimation based on n observations, θ is
“true” parameter of the distribution of the observations. Moreover,
  2 −1
d ln L(θ)
δ(θ) = − E
dθ2

One can even show that for any other estimator θ̃, which converges
in distribution to N (θ, σ 2 ), we have δ(θ) ≤ σ 2 (variance is greater or
equal). Therefore, MLEs are called best asymptotically normal. For
large n, we can estimate the probability that θ̂ deviates from the true
value by more than  (see confidence intervals; later).

• strongly consistent: MLEs are strongly consistent, i.e,

P ( lim θ̂n = θ) = 1.
n→∞

In the case of more than one parameter the MLE is found in a very similar
way, i.e. we find the vector θ̂ = (θ̂1 , ..., θ̂m ) that maximizes L(θ) or ln L(θ).

For instance, for two parameters α and β we try to solve ∂α ln L(θ) = 0 and

∂β ln L(θ) = 0 simultaneously for α and β where θ = (α, β). (See exercises
for an example.)

7.2.1 Variance of MLE


Before we consider the variance of the maximum likelihood estimator, we
will discuss the variance of estimators in general.

79
7.2. MAXIMUM LIKELIHOOD ESTIMATION

It is very important to analyze how good the estimated value given by some
estimator θ̂ is. If V (θ̂) is large, then our estimation is of bad quality and
might lead to wrong conclusions about the real system. Typically, when
results ofpparameter estimations are reported, the (estimated) standard de-
viations V (θ) are given as well.

Example 63: Poisson distribution

In the previous sections we found that both the method of moments and
the maximum likelihood approach gives the estimator

θ̂ = X̄

for the Poisson distribution, i.e. the best value for the unknown mean of
the Poisson distribution is the sample mean. In Section ?? we already
found that the variance of X̄ is σ 2 /n where σ 2 is the variance of the
distribution of the Xi and thus for the Poisson distribution we get

V (θ̂) = V (X̄) = θ∗ /n.

Here, θ∗ is the ’true’ (unknown) value that we estimate with θ̂ = X̄.


Thus, we can estimate V (θ̂) as X̄/n.
Alternatively, we can estimate σ 2 by computing the sample variance S 2
(as explained in Example 59) and estimate V (θ̂) as S 2 /n.

Often, analytic formulas for V (θ̂) cannot be derived. In this case, one can
either use the bootstrap method or estimate V (θ̂) based on the properties
of the estimator.

Maximum Likelihood Method. To estimate the variance of a maxi-


mum likelihood estimator, we exploit the fact that MLEs are asymptotically
normally distributed. Thus, for large sample sizes and a single parameter
θ ∈ R, we have approximately a variance of
  2 −1
∗ d ln L(θ∗ , X1 , . . . , Xn )
δ(θ ) = − E .
dθ2

In the above expression we wrote L(θ, X1 , . . . , Xn ) for the likelihood to em-


phasize that L is a function of θ and X1 , . . . , Xn . Moreover, here θ∗ is
the true value of θ. For an estimator based on the data we simply replace
X1 , . . . , Xn by their corresponding realizations x1 , . . . , xn and the (unknown)
true value of θ by the estimate θ̂, i.e.
!−1
∂ 2 ln L(θ̂, x1 , . . . , xn )
− .
∂θ2

80
7.2. MAXIMUM LIKELIHOOD ESTIMATION

Intuitively, the second derivative tells us the curvature of the likelihood and
if it is ’flat’ at θ̂ then our estimated value might not be very accurate and
the variance (negative inverse) is large. Then either we do not have enough
samples or the parameter is difficult to identify (see also the example after
next below). We are not very confident about our estimated value since
perturbing θ̂ slightly also gives us a similar likelihood.
For several parameters, the same approach is used, i.e. the negative diagonal
entries of the inverse of the Hessian matrix give estimates for the variances
of the parameters.

Example 64: Parameters of the normal distribution


Assume that we have n i.i.d. samples X1 , . . . , Xn that follow a normal
distribution with (unkown) mean µ and (unkown) variance σ 2 . Hence,
θ = (µ, σ 2 ).
The likelihood of the data is
Yn  
1 1
L(X1 , . . . , Xn ; θ) = √ exp − 2 (Xi − µ)2
2πσ 2 2σ
i=1

. Thus, log-likelihood of the data is


n
n n 1 X
ln L(X1 , . . . , Xn ; θ) = − ln(2π) − ln(σ 2 ) − 2 (Xi − µ)2 .
2 2 2σ
i=1

Next we compute the derivatives w.r.t. µ and σ2:


∂ 1 Pn
∂µ ln L(X1 , . . . , Xn ; θ) = σ 2 i=1 (Xi − µ)
Pn

∂σ 2
ln L(X1 , . . . , Xn ; θ) = − n2 (σ 2 )−1 + 21 (σ 2 )−2 i=1 (Xi − µ)2 .
Setting the two derivatives to zero yields from the first equation
n
X n
X n
1X
(Xi − µ) = 0 =⇒ Xi − nµ = 0 =⇒ µ̂ = Xi = X̄.
n
i=1 i=1 i=1

Inserting this into the second equation and multiplying by 2σ 4 gives


n
X n
2 2 2 1X
−nσ + (Xi − µ̂) = 0 =⇒ σ̂ = (Xi − X̄)2 .
n
i=1 i=1

Thus, the MLE of σ is a biased estimator.


To prove that µ̂ and σ̂ 2 yield a maximum we have to consider the Hessian
matrix.
 ∂2 ∂2

∂µ∂µ ln L ∂µ∂σ 2
ln L
H(µ, σ) =  2 
∂ ∂2
∂σ 2 ∂µ
ln L ∂σ2 ∂σ2 ln L

81
7.3. BAYESIAN INFERENCE

 Pn 
− σn2 −(σ 2 )−2 i=1 (Xi − µ)
= P Pn

−(σ 2 )−2 ni=1 (Xi − µ) n 2 −2
2 (σ ) − (σ 2 )−3 i=1 (Xi − µ)
2

∂2
We have a local maximum at θ̂ = (µ̂, σ̂ 2 ) if ∂µ∂µ ln L < 0 and if the
determinant
∂2 ∂2 ∂2 ∂2
det(H(µ̂, σ̂ 2 )) = ln L · ln L − ln L · ln L
∂µ∂µ ∂σ 2 ∂σ 2 ∂µ∂σ 2 ∂σ 2 ∂µ

is positive at θ̂. Next we insert concrete numbers to simplify the computa-


tion. We generate 50 random numbers that are normally distributed with
mean zero and variance one. We get the estimates

µ̂ = −0.1229 and σ̂ 2 = 0.9903

and the Hessian is (approximately)


" #
−50.4906 −0.0000
.
−0.0000 −25.4930

(Note that its determinant is indeed positive.) The inverse of the Hessian
is " #
−0.0198 0.0000
.
−0.0000 −0.0392
Thus, the estimated variances are 0.0198 for µ̂ and 0.0392 for σ̂ 2 , yielding
standard errors of 1 = 0.1407 and 2 = 0.1981. Note that the true mean
and variance are elements of the intervals [µ̂−1 , µ̂+1 ] and [σ̂ 2 −2 , σ̂ 2 +
2 ].

If the standard errors are extremely high, this shows that the log- likelihood
function is ’very flat’ around the minimum and thus the estimated values
may not be close to the true values of the parameters. We then have the
problem of parameters that are not identifiable. However, sometimes at
least ratios of parameters are estimated very accurately.

7.3 Bayesian Inference


In the previous sections, we followed the so called frequentist approach when
estimating population parameters. We worked with distributions, expecta-
tions, and variances of the given random samples (population data) as well
as of estimators and other statistics. The likelihood Lθ (x1 , . . . , xn ) of the
data given the parameter θ is a central function in the frequentist approach.

82
7.3. BAYESIAN INFERENCE

In this section, we will follow the Bayesian approach, which differs from the
frequentist approach in that θ is treated as a random variable2 and has a
certain probability distribution. Hence, not only the data is a source of un-
certainty but also θ. Its distribution π(θ), called prior distribution, reflects
our ideas, beliefs, and past experiences about θ before we make use of the
data x1 , . . . , xn , i.e., before we perform inference based on the data.

For example, our prior belief could be that it is likely


that the true value of θ ∈ [0, 1] lies in a certain in-
terval [a, b] ⊂ [0, 1]. However, we cannot completely
()

exclude the cases θ < a and θ > b. Hence, we could


assume a prior distribution π(θ) as depicted in the
figure on the left for [a, b] = [0.3, 0.7].
0.3 0.7
The information given by the population data may
lead to a change of our prior beliefs, i.e., after the Bayesian inference we
obtain a posterior distribution which tells us how likely different values of θ
are when we know the data. The main advantage of this approach is that
it gives also meaningful results in the case of only few samples since in this
setting we do not rely on asymptotic estimator properties (which require
a large number of samples). On the other hand, the posterior distribution
does not give a fixed value or confidence interval for θ but only a distribu-
tion. In the above example, we could, for instance, compute the posterior
probability of the interval [0.3, 0.7] and compare it to the probability of the
prior to see whether the data lead to an increase of the probability that
θ ∈ [0.3, 0.7].
Hence, to perform Bayesian inference, we need besides the observed data
a prior distribution π(θ) for θ. Assume that the data has likelihood
Lθ (x1 , . . . , xn ) = L(x1 , . . . , xn | θ) (it is helpful to make the condition on
θ explicit at this point by writing ’| θ’). Note that this condition means
that our random parameter – let us call it Θ – takes the concrete value θ.
Since π(θ) is the probability that Θ = θ, we get the posterior distribution
π(θ | x1 , . . . , xn ) as
L(x1 , . . . , xn | θ)π(θ)
π(θ | x1 , . . . , xn ) = (7.1)
f (x1 , . . . , xn )
according to Bayes’ Theorem. The denominator f (x1 , . . . , xn ) is the
marginal probability of the data, i.e., independent of a concrete choice of
θ = Θ. It can be computed using the law of total probability.
X
f (x1 , . . . , xn ) = L(x1 , . . . , xn | θ)π(θ) (7.2)
θ
2
In the previous sections we used upper case letters to denote random variables. In
this section, we use θ although it is a random variable as this is the standard notation in
the literature.

83
7.3. BAYESIAN INFERENCE

If the prior is a continuous distribution, we compute the marginal as


Z
f (x1 , . . . , xn ) = L(x1 , . . . , xn | θ)π(θ) dθ. (7.3)
θ

Example 65: Quality inspection3

A manufacturer claims that the shipment contains only 5% of defective


items, but the inspector feels that in fact it is 10%. We have to decide
whether to accept or to reject the shipment based on θ, the proportion of
defective parts. Before we see the real data, let?s assign a 50-50 chance
to both suggested values of θ, i.e.,

π(0.05) = π(0.10) = 0.5.

A random sample of 20 parts has 3 defective ones.


We calculate the posterior distribution of θ as follows: First, we find that
the probability of X defective parts in a sample of n parts is binomially
distributed with parameters n and θ. Hence,
 
n x
L(x | θ) = θ (1 − θ)n−x
x

and therefore

L(3 | 0.05) ≈ 0.0596 and L(3 | 0.10) ≈ 0.1901.

For the marginal distribution we compute

f (x) = 0.5 · L(x | Θ = 0.05) + 0.5 · L(x | Θ = 0.10) ≈ 0.12485.

Finally, the posterior probabilities are


L(X=3|0.05)π(0.05)
π(0.05 | X = 3) = f (X=3) ≈ 0.2387,
L(X=3|0.10)π(0.10)
π(0.10 | X = 3) = f (X=3) ≈ 0.7613,

which indicate a threefold higher chance that the proportion of defective


parts is 10% compared to the 5% claimed by the manufacturer.

7.3.1 Conjugate families of distributions


It is often an important advantage of Bayesian approaches if we can deter-
mine the posterior distribution analytically, since otherwise we may not be
3
Taken from [1].

84
7.3. BAYESIAN INFERENCE

able to efficiently compute the posterior distribution. We say that a family


of prior distributions is conjugate to a model if the posterior distribution
belongs to the same family as the prior.

Example 66: Gamma prior is conjugate to a Poisson model

We consider realizations x1 , . . . , xn of a Poisson distribution with param-


eter θ as our model. Hence,
n
Y n
Y e−θ θxi
f (x1 , . . . , xn | θ) = f (xi | θ) = .
xi !
i=1 i=1

Next, we assume for θ a Gamma(α, λ) distribution as a prior, i.e.,

λα α−1 −λθ
π(θ) = θ e ∼ θα−1 e−λθ ,
Γ(α)

where ∼ means “proportional to”. In the sequel, we will determine the


posterior only up to a constant coefficientQ to simplify the derivation.
Therefore, we will also drop the factor ni=1 x1i ! of the density f of the
data, which is constant in θ, i.e., for ~x = (x1 , . . . , xn )
n
Y n
Y
e−θ θxi P
f (~x | θ) = ∼ e−θ θxi = e−nθ θ xi
.
xi !
i=1 i=1

Now, the posterior is given by

π(θ | ~x) ∼ f (~x | θ)π(θ)


P  
∼ e−nθ θ xi θα−1 e−λθ
P
= θα+ xi −1 e−(λ+n)θ .

P
Thus, the posterior is a Gamma(α + xi , λ + n) distribution. Note that
the missing constant factor of the distribution can be uniquely determined
as π(θ | ~x) is a proper density.
It is interesting to observe how the mean and the variance of the Gamma
distribution is adjusted, when we take the data into account: For the prior
we Phave E[Θ] = αλ and V P [Θ] = λα2 , while the posterior has E[Θ | ~x] =
α+ xi α+ xi
λ+n and V [Θ | ~x] = (λ+n)2
.

A list of further conjugate prior distributions can be found on Wikipedia.

85
7.3. BAYESIAN INFERENCE

7.3.2 Bayesian point-estimators


From the posterior distribution different estimators can be determined. The
most common estimator is the posterior mean, which gives the average value
of θ conditioned on the data, i.e.,
X
θ̂M = E[Θ | ~x] = θ π(θ | ~x)
θ

if the posterior/prior is a discrete distribution. In the continuous case, the


sum is replaced by an integral over θ.
For this estimator, the variance based on the posterior distribution is com-
puted as
V [Θ | ~x] = E[(θ̂M − Θ)2 | ~x].
Other popular point estimators are the posterior median and the maximum
a posteriori (MAP),
θ̂MAP = arg maxθ π(θ | ~x)
which is given by the maximum of the posterior distribution. Note that
the MAP estimator is not a good choice if the maximum of the posterior
distribution is not representative for the distribution, e.g. if the posterior
distribution is bimodal and θ̂MAP is determined by a narrow peak, while
other very different values of θ are also quite likely.

Example 67: Quality inspection (continued)


Recall the posterior distribution of Example 65,

π(0.05 | x) ≈ 0.2387,
π(0.10 | x) ≈ 0.7613.

The posterior mean is


X
θ̂M = E[Θ | x] = θ π(θ | x) ≈ 0.05 · 0.2387 + 0.1 · 0.7613 ≈ 0.0880,
θ

which is much closer to the inspector’s estimate. The corresponding vari-


ance is

V [θ | x] = E[(θ̂M − Θ)2 | ~x]


≈ (0.088 − 0.05)2 · 0.2387 + (0.088 − 0.1)2 · 0.7613
≈ 0.0004,

which means that the standard deviation is about 0.02. Obviously, the
MAP estimator would be θMAP = 0.1 as the posterior probability at 0.1
is higher than that at 0.05.

86
Chapter 8

Statistical Testing

Carrying out statistical test that are related to the observations of the real
system and/or the model is very important when we make claims or state-
ments about the system. A popular class of tests are hypothesis tests that
are used to verify statistical hypotheses.
In a statistical hypothesis test, we usually first formulate a (null) hypothesis
H0 and an alternative hypothesis HA . The two statements H0 and HA must
be mutually exclusive and the test can either accept or reject H0 (in favor
of HA ). The null hypothesis is either
• an equality
• the absence of an effect or some relation
Note that this leads to null hypotheses that are in many cases not the
same as the statement that we want to verify. The reason for the above
constraint is that intuitively, we need an equality (or absence of an effect or
some relation) to fix the distribution that we consider during the test.
We begin with a motivating example.

Example 68: Defective Products

Assume that a manufacturer claims that at most 3% of his products are


defective. We want to verify this statement to decide whether we accept
the shipment of the products or not. We define

H0 : fraction of defective products is equal to 3%


HA : fraction of defective products is greater than 3%

where for HA we used the right-tail alternative. If we have evidence that


H0 is rejected in favor of HA then we reject the shipment.
Note that HA : ’fraction of defective products is smaller than 3%’ is not
useful since then we will always accept the shipment - no matter if we
have evidence for H0 or for HA .

87
266 Probability and Statistics for Computer Scientists

where µ1 is the average number of concurrent users last year, and µ2 is the average number
of concurrent users this year. Depending on the situation, we may replace the two-sided
(1)
alternative HA : µ2 − µ1 ̸= 2000 with a one-sided alternative HA : µ2 − µ1 < 2000 or
(2)
8.1. LEVEL α TESTS: (1)
GENERAL APPROACH
HA : µ2 − µ1 > 2000. The test of H0 against HA evaluates the amount of evidence that
(2)
the mean number of concurrent users changed by fewer than 2000. Testing against HA , we
see if there is sufficient evidence to claim that this number increased by more than 2000. ♦
Example 69: Concurrent Users
Assume that we want to verify the statement that the average number of
concurrent
Example 9.24. To verify users ofofdefective
if the proportion an online PCisgaming
products platform
increased by 2000 this
at most 3%, we test
[Link] define
0 : p = 0.03 vs HA : p > 0.03,

where p is the proportion of defects in the whole shipment. H0 : µ2 -µ1 =2000


Why do we choose the right-tail alternative HA : p > 0.03?HThat -µ1 6= 2000
A : µis2because we reject the
shipment only if significant evidence supporting this alternative is collected. If the data
where µ1 is the average number of concurrent users♦ of
suggest that p < 0.03, the shipment will still be accepted.
the last year and µ2
the average number of this year. Here, HA is called a two-sided alternative
since it covers both cases µ2 -µ1 > 2000 and µ2 -µ1 < 2000.
9.4.2 Type I and Type II errors: level of significance

Fromwethe
When testing hypotheses, twothat
realize examples
all we seeabove we see
is a random that Therefore,
sample. there can be two-sided alternatives,
with
all the best statisticsone-sided,
skills, our decision to accept
left-tail or to reject H0(H
alternatives may still
is be
µ wrong.
< µ That
A 0 and one-sided, right-tail
)
would be a sampling error (Section 8.1).
alternatives (HA is µ > µ0 ), where H0 is µ = µ0 .
Four situations are possible,
The outcome of our test depends on a
finite random sample and thus we may
Result of the test
always take a wrong decision. The
Reject H0 Accept H0 four situations depicted on the left are
possible. Our goal is to keep each of
H0 is true Type I error correct
the two errors small. Thus, a good
H0 is false correct Type II error test only results in a wrong decision if
the sample is not very representative
(i.e.
In two of the four cases, the test results in a correct decision. Either we extreme).
accepted a true Often the type I error
hypothesis, or we rejected a false hypothesis. The other two situations are sampling errors.
is seen as more dangerous since it corresponds to ’convicting an innocent
defendant’ or ’sending a healthy patient to a surgery’. Therefore, we fix the
DEFINITION 9.7
probability α of a type I error which is also called the significance level of
A type I error occurs when we reject the true null hypothesis.
the test.
A type II error occurs when we accept the false null hypothesis.
α = P {reject H0 | H0 is true}
The probability of rejecting a false hypothesis (avoid a type II error) is the
Each error occurs with a certain
power of probability
the test that
andweahope to keep small.
function of theA good test resultsθ
parameter
about which we make
in an erroneous decision only if the observed data are somewhat extreme.
our hypothesis:
p(θ) = P {reject H0 | θ; HA is true}
Typically, α is chosen very small, e.g. α ∈ {0.01, 0.05, 0.10} such that the
type I error is kept small and we only reject H0 with a lot of confidence.

8.1 Level α tests: general approach


To test H0 against HA we perform the following steps:
1. Compute a test statistic T which is a function of the sample (or an
estimator) and thus a random number. The distribution of T , given
H0 is true, is known.

88
8.1. LEVEL α TESTS: GENERAL APPROACH

2. Consider the null distribution F0 of T , given H0 is true, and find the


portion that corresponds to α, i.e. the part of the area below the
density curve that gives α and is thus the region where we reject H0 .
268 The remaining
Probability partfor
and Statistics (ofComputer
area 1-α) is the acceptance region.
Scientists

In the illustration on the left,


f0 (T ) ✻ T is expected to be large if HA
Rejection is true. Therefore, the rejec-
region, tion region is at the right tail
probability α
Acceptance when H0 is true
of the distribution. In a two-
region, sided test, we consider por-
probability (1 − α) tions of α/2 at both tails of the
when H0 is true ✠ ✲ distribution.
T
We always have that
FIGURE 9.6: Acceptance and rejection regions.

P {T ∈ acceptance region | H0 is true} = 1 − α


These regions are selected in such a way that the values of test andstatistic T in the rejection
region provide a stronger support of HA than the values T ∈ A. For example, suppose that
P {T
T is expected to be large if HA is true. Then ∈ the
rejection
rejection region | H0 is true}
regioncorresponds = α.
to the right
tail of the null distribution F0 (Figure 9.6).
3. look
As another example, Accept H0 9.3
at Figure if Tonbelongs
p. 249. Ifto
thethe
nullacceptance
distribution ofregion and reject it otherwise.
T is Standard
Normal, then the areaIt
between (−z α/2 ) and z α/2 equals exactly (1 − α). The interval
is important to mention that if we accept H0 we cannot say that
’with probability 1−
A = (−zα/2 α the
, zα/2 ) hypothesis H0 is true’. The reason is that
H0 is not random, i.e. it either holds with probability one or it does
can serve as a level α acceptance region for a two-sided test of H0 : θ = θ0 vs HA : θ ̸= θ0 .
not hold
The remaining part consists of twowith probability
symmetric tails, one. The only correct interpretation is that
if we reject H0 then the data provides sufficient evidence against H0
R = Ā = (−∞, −zα/2 ] ∪ [zα/2 , +∞);
and in favor of HA . Either H0 is really not true or our data is not
this is the rejection region.
representative - which happens with probability α. If we accept H0
Areas under the density then the
curve aredata does not
probabilities, and provide
we concludesufficient
that evidence to reject H0 . In the
absence of sufficient
P {T ∈ acceptance region evidence,
| H0 } = 1 by
− α default, we accept H0 .

and
Let us now P {Tin∈detail reason
rejection region about
| H0 } = finding
α. the rejection region with area α
since there are many regions with area α below the density curve. The best
choice is an area that ensures that the type II error is small. Thus, we
Step 3: Result and its interpretation
choose the rejection region such that it is likely that T falls into this region
if HA isH0true.
Accept the hypothesis if the This will maximize
test statistic T belongs tothe
the power
acceptanceof region.
the test, i.e. the probability
Reject
H0 in favor of of
the rejecting
alternative H HA0ifgiven HAtoisthetrue.
T belongs Often,
rejection [Link] results in the a choice of the
Our acceptancerejection region
and rejection regionsasguarantee
illustrated insignificance
that the Figure 8.1 levelfor a normally
of our test is distributed test
statistic. Usually, the test statistic is defined such that
Significance level = P { Type I error }
= P { Reject | H0 }
• the right-tail alternative forces T to be large,
= P {T ∈ R | H0 }
= α. (9.13)
• the left-tail alternative forces T to be small,
Therefore, indeed, we have a level α test!
• the two-sided alternative forces T to be either large or small.

89
8.2. 270
STANDARD NORMAL NULL DISTRIBUTION (Z-TEST)
Probability and Statistics for Computer Scientists

Reject
if T is here
Accept Accept
if T is here if T is here
✢ ✲ ❫ ✲
0 zα T −zα 0 T
(a) Right-tail Z-test (b) Left-tail Z-test

Reject Reject
if T is here Accept if T is here
if T is here

❯ ☛ ✲
−zα/2 0 zα/2 T

(c) Two-sided Z-test

FIGURE 9.7: Acceptance and rejection regions for a Z-test with (a) a one-sided right-tail
Figure 8.1: Acceptance
alternative; (b) a one-sidedand rejection
left-tail regions
alternative; for a normally
(c) a two-sided alternative. distributed test
statistic. a) one-sided right-tail alternative; b) one-sided left-tail alternative;
c) two-sided alternative.
9.4.5 Standard Normal null distribution (Z-test)

An important case, in terms of a large number of applications, is when the null distribution
of the test statistic is Standard Normal.

8.2The test
Standard Normal
in this case is called Null
a Z-test, and the testDistribution (Z-test)
statistic is usually denoted by Z.
(a) A level α test with a right-tail alternative should
For a large number of applications,
! the null distribution of T (the distribu-
reject H0 if Z ≥ zα
tion of T given H0 is true) is standard
accept H0 if normal.
Z < zα Then the test is (9.14)
called a
Z-test. Usually, one the following cases applies:
The rejection region in this case consists of large values of Z only,

• we consider sample R
means of normally
= [zα , +∞), distributed
A = (−∞, zα ) data,

•(see
weFigure 9.7a).
consider sample means of arbitrarily distributed data where the
Under
number of samplesZ is
the null hypothesis, belongs to A and we reject the null hypothesis with probability
large,
P {T ≥ zα | H0 } = 1 − Φ(zα ) = α,
• we consider sample proportions of arbitrarily distributed data where
making the probability
the number of false rejection
of samples (type I error) equal α .
is large,
For example, we use this acceptance region to test the population mean,
• we consider differences of sample means or sample proportions where
H0 : µ = µ0 vs HA : µ > µ0 .
the number of samples is large.

In all of these cases a Z-test can be used.


Let zα be the number such that P (Z > zα ) = Φ(zα ) = α if Z is a standard
normally distributed random variable. Then we reject H0 if

• Z ≥ zα for a test with right-tail alternative,

• Z ≤ −zα for a test with left-tail alternative,

90
8.2. STANDARD NORMAL NULL DISTRIBUTION (Z-TEST)

• |Z| ≥ zα/2 for a test with two-sided alternative.

Note that each time we have

P {reject H0 | H0 is true} = 1 − Φ(zα ) = α.

Example 70: A test for the mean

Assume that we have used X̄ to estimate the unknown mean µ0 of a


normal distribution based on n = 100 independent samples. Since the
samples are normally distributed, X̄ is normally distributed too and we
know that E[X̄] = µ0 (see previous chapter) and that V [X̄] = σ 2 /n.
Assume further that σ = 800 is known and that the data is such that
X̄ = 5200. We would like to verify (with a significance of 5%, α = 0.05)
whether the true mean is greater than 5000, i.e. H0 : µ0 = 5000 and HA :
µ0 > 5000 (right-tail alternative).

1. We first compute the test statistic

X̄ − µ0 5200 − 5000
Z= √ = √ = 2.5.
σ/ n 800/ 100

2. The critical value is zα = 1.645 (from the table of the standard


normal distribution). Thus, we should reject H0 if Z ≥ 1.645 and
accept it otherwise.

3. Since Z falls into the rejection region, we have enough evidence to


reject H0 and support the alternative hypothesis HA : µ0 > 5000.

Note that if, for instance, X̄ = 5100, we would get Z = 1.25 and not have
enough evidence to reject H0 and believe in the alternative hypothesis that
µ0 > 5000.

Example 71: Two-Sample Z-test of proportions

A quality inspector finds 10 defective parts in a sample of n = 500 parts


received from manufacturer A. Out of m = 400 parts from manufacturer
B, she finds 12 defective ones. A computer-making company uses these
parts in their computers and claims that the quality of parts produced by
A and B is the same. At the 5% level of significance, do we have enough
evidence to disprove this claim?

We test H0 : pA = pB , or H0 : pA − pB = 0, against HA : pA 6= pB
where pA (pB ) is the portion of defective parts from manufacturer A (B),

91
8.2. STANDARD NORMAL NULL DISTRIBUTION (Z-TEST)

respectively.
This is a two-sided test because no direction of the alternative has been
indicated.

1. We first compute the values of the estimated portions


10 12
p̂A = = 0.02 and p̂B = = 0.03.
500 400
These estimators are asymptotically normally distributed (sum of
n and m Bernoulli variables divided by n and m, respectively)
where the means are the true portions pA and pB and the vari-
ances are pA (1 − pA )/n and pB (1 − pB )/m. Obviously, p̂A − p̂B is
also asymptotically normally distributed with mean zero and vari-
ance pA (1 − pA )/n + pB (1 − pB )/m. Thus, inserting the estimated
values we standardize and get
p̂A − p̂B
Z=p = −0.945.
pA (1 − pA )/n + pB (1 − pB )/m

2. The critical value is zα/2 = 1.96 (from the table of the standard
normal distribution). Since this is a two-sided test we should reject
H0 if |Z| ≥ 1.96 and accept it otherwise.

3. Since Z falls into the acceptance region, we do not have enough


evidence to reject H0 . Although the sample proportions of defective
parts are unequal, the difference between them appears too small to
claim that population proportions are different.

Alternatively, we can consider a single estimator p̂ for the overall propor-


tion of defective products since we assume that H0 (pA = pB ) holds. But
with p := pA = pB this implies that

V [XA ] = V [XB ] = p(1 − p)

if XA and XB are Bernoulli distributed with parameter p. Hence, if we


estimate the common portion of defective parts as
number of defective parts np̂A + mp̂B
p̂ = = = 0.0244
total number of parts n+m
then we can use it to replace the unknown probability p in the variance
estimators V [p̂A ] = p(1−p)
n and V [p̂B ] = p(1−p)
m . Hence,

p̂(1 − p̂) p̂(1 − p̂)


V [p̂A − p̂B ] = +
n m

92
8.3. T-TESTS FOR UNKNOWN σ
Statistical Inference I 273

Null Parameter,
If H0 is true: Test statistic
hypothesis estimator
θ̂ − θ0
H0 θ, θ̂ E(θ̂) Var(θ̂) Z=!
Var(θ̂)

One-sample Z-tests for means and proportions, based on a sample of size n

σ2 X̄ − µ0
µ = µ0 µ, X̄ µ0 √
n σ/ n

p0 (1 − p0 ) p̂ − p0
p = p0 p, p̂ p0 !
n p0 (1−p0 )
n

Two-sample Z-tests comparing means and proportions of two populations,


based on independent samples of size n and m

µX − µY , 2
σX σ2 X̄ − Ȳ − D
µX −µY = D D + Y !
X̄ − Ȳ n m 2
σX 2
σY
n + m

p1 − p2 , p1 (1 − p1 ) p2 (1 − p2 ) p̂1 − p̂2 − D
p1 −p2 = D D + !
p̂1 − p̂2 n m p̂1 (1−p̂1 ) p̂2 (1−p̂2 )
n + m

$ p̂1 − p̂2
" # " #
p1 − p2 , 1 1 1 1
p(1 − p) + , p̂(1 − p̂) +
p1 = p2 0 n m n m
p̂1 − p̂2 where p = p1 = p2
np̂1 + mp̂2
where p̂ = n+m

TABLE 9.1: Summary of Z-tests.


Figure 8.2: A summary of Z-tests (from [1])

Example 9.26 (Two-sample Z-test of proportions). A quality inspector finds 10


defective parts in a sample of 500 parts received from manufacturer A. Out of 400 parts
from manufacturer B, she finds 12 defective ones. A computer-making company uses these
and
partsget
in their computers and claims that the quality of parts produced by A and B is the
do−
same. At the 5% level of significance,p̂A wep̂have
B enough evidence to disprove this claim?
Z=q = −0.966.
p̂(1−p̂)
+ p̂(1−p̂)
Solution. We test H0 : pA = pB , ornH0 : pA −mpB = 0, against HA : pA ̸= pB . This is a
two-sided test because no direction of the alternative has been indicated. We only need to
Again, Z falls into the acceptance region and we do not have enough
verify whether or not the proportions of defective parts are equal for manufacturers A and
evidence
B. to reject H0 .
Step 1: Test statistic. We are given: p̂A = 10/500 = 0.02 from a sample of size n = 500;
p̂B = 12/400 = 0.03 from a sample of size m = 400. The tested value is D = 0.

8.3 T-tests for Unknown σ

In the previous section we used an estimator for the unknown true variance
σ 2 . In the special case that our data X1 , . . . , Xn is normally distributed

93
8.3. T-TESTS FOR UNKNOWN σ

with mean µ and variance σ 2 and we estimate the unknown mean with
n
1X
X̄n = Xi
n
i=1

we can make use of a result from statistics that tells us the following:
We know that E[X̄n ] = µ and thus
X̄n − µ
Z= p
σ 2 /n
must follow a standard normal distribution (we standardized it by subtract-
ing the mean and dividing by the standard deviation!). However, since σ 2
is unknown we estimate it using
n
1 X
Sn2 = (Xi − X̄n )2 .
n−1
i=1

It can be shown that


X̄n − µ
T =p
Sn2 /n
follows a Student’s t-distribution with n − 1 degrees of freedom. Thus, once
µ is fixed (because of our hypothesis H0 ) and the data X1 , . . . , Xn is given we
can compute T and check whether it falls into the acceptance or rejection
region by looking at the values tα , −tα or tα/2 depending on the type of
alternative. Note that tα is the value such that
P (T > tα ) = α
and it can be found in the table of the Student’s t-distribution (in the same
way we find tα/2 ). Note that the degrees of freedom n − 1 is a parameter
of the Student’s t-distribution since the distribution of T changes when the
276samples size n changes.
Probability and Statistics for Computer Scientists

Hypothesis Test statistic Degrees of


Conditions
H0 t freedom

Sample size n; X̄ − µ0
µ = µ0 t= √ n−1
unknown σ s/ n

Sample sizes n, m;
unknown but equal X̄ − Ȳ − D
µX − µY = D t= ! n+m−2
standard deviations,
sp n1 + m 1
σX = σY

Sample sizes n, m;
Figure 8.3: Aunequal
summary of T-tests Satterthwaite
unknown, X̄ − Ȳ (from
− D [1]).
µX − µY = D t= ! approximation,
standard deviations, s2X s2Y
+ formula (9.12)
σX ̸= σY n m

94
TABLE 9.2: Summary of T-tests.

at a significance level α = 0.01. From Example 9.19, we have sample statistics n = 18,
X̄ = 0.29 and s = 0.074. Compute the T-statistic,
¯
8.4. P-VALUE

8.4 p-Value Statistical Inference I 281

So far, we were testing hypotheses by (a) Small α


means of acceptance and rejection re- Accept H 0

gions where we need to know the sig-


nificance level α in order to conduct Accept
H 0

a test. However, there is no system- Zobs ✲


0 z α
atic way of choosing α, the probabil-
ity of making type I error. Of course,
when it seems too dangerous to reject (b) Large α,
same Zobs
a true H0 , we choose a low significance Reject H 0

level. But how low? Should we choose Accept


H 0
α = 0.01? Perhaps, 0.001? Or even Zobs ✲
0 z α
0.0001? If we choose α small enough
we can always make sure that H0 isFIGURE accepted
Under
9.10: (a) Under a low level of significance α, we accept the null hypothesis. (b)
a high level ofas illustrated
significance, we reject it. on the right.
280 Probability and Statistics for Computer Scientists
Also, if our observed test
P-value
statistic Zobs belongs to a re-
Using a P-value approach, we try not to rely on the level of significance. In fact, let us try
to test a hypothesis jection
using all levelsregion but it is ”too
of significance!
Zobs
Accept
close to call” as illustrated on
Considering all levels of significance (between 0 and 1 because α is a probability of Type I
error), we notice:
H0 Case 1. If a level ofthe left,is very
significance thenlow, wehow donullwe
accept the report
hypothesis (see Figure 9.10a).
A low value of
the result?
α = P { reject the null hypothesis when it is true }
makes it very unlikely to reject the hypothesis because it yields a very small rejection region.
❘ ✲
The right-tail area Formally,
above the rejectionwe regionshould
equals α. reject the
−zα/2 0 zα/2
null hypothesis, but practi-
Case 2. On the other extreme end, a high significance level α makes it likely to reject the null
hypothesis and corresponds to a large rejection region. A sufficiently large α will produce
FIGURE 9.9: This test is “too close to call”: formally we reject the null
suchhypothesis although
the Z-statistic is almost at the boundary. Figure 9.10b).
cally,
a large rejection
we realize that a slightly
region that will cover our test statistic, forcing us to reject H (see
0

different significance level α could have expanded


Conclusion: the acceptance
there exists a boundary region
value between α-to-accept just
(case 1) and α-to-reject (case
2). This number is a P-value (Figure 9.11).
Example enough to cover Zobs and force us to accept H0 . Is there a statistical measure
9.36 (Unauthorized use of a computer account, continued). A 99%
confidence interval for the mean time between keystrokes is
to quantify how far away we are from a ”too close to call”?
[0.24; 0.34]
The idea is to try to test a hypothesis using all levels of significance. Then
(Example 9.19 on p. 260). Example 9.28 on p. 275 tests whether the mean time is 0.2
seconds,we
whichhave two
would be cases:
consistent with the speed of the account owner. The interval does
not contain 0.2. Therefore, at a 1% level of significance, we have significant evidence that
1) Very small values of α make it very unlikely♦ to reject the hypothesis
the account was used by a different person.
because they yield very small rejection regions. 2) High significance levels α
will make it likely to reject H0 and corresponds to a large rejection region.
9.4.10 We will be forced to reject H0 . The P-value is the boundary value between
P-value
the accept case 1) and reject case 2). Thus,
How do we choose α?

the p-value is the lowest significance level α that forces rejection of H0 and
So far, we were testing hypotheses by means of acceptance and rejection regions. In the last
section, we learned how to use confidence intervals for two-sided tests. Either way, we need
also the highest significance level α that forces acceptance of H0 .
to know the significance level α in order to conduct a test. Results of our test depend on it.
How do we choose α, the probability of making type I sampling error, rejecting the true
hypothesis? Of course, when it seems too dangerous to reject true H0 , we choose a low
Usually α ∈ [0.01, 0.1] (although there are exceptions). Then, a P-value
significance level. How low? Should we choose α = 0.01? Perhaps, 0.001? Or even 0.0001?
greater than 0.1 exceeds all natural significance levels, and the null hypoth-
Also, if our observed test statistic Z = Zobs belongs to a rejection region but it is “too
close to call” (see, for example, Figure 9.9), then how do we report the result? Formally,
esis should be accepted. Conversely, if a P-value is less than 0.01, then it is
we should reject the null hypothesis, but practically, we realize that a slightly different
smaller than all natural significance levels, and the null hypothesis should
significance level α could have expanded the acceptance region just enough to cover Zobs
and force us to accept H0 .
Supposebethatrejected. Only
the result of our test is if the important.
crucially P-valueForhappensexample, the to fall
choice of a between 0.01 and 0.1, we
businessreally have to think about the level of significance. This is the ”too close to
strategy for the next ten years depends on it. In this case, can we rely so heavily
on the choice of α? And if we rejected the true hypothesis just because we chose α = 0.05
instead of α = 0.04, then how do we explain to the chief executive officer that the situation
was marginal? What is the statistical term for “too close to call”?
95
Important:

Pr (observation | hypothesis) ≠ Pr (hypothesis | observation)

The probability of observing a result given that some hypothesis


is true is not equivalent to the probability that a hypothesis is true
given that some result has been observed.

Using the p-value as a “score” is committing an egregious logical error:


8.4. P-VALUE
the transposed conditional fallacy.

More likely observation

Statistical Inference I 283


P-value
Probability density
Computing P-values

Here is how a P-value can be computed from data.


Very un-likely Very un-likely
Let us look at Figure 9.10 again. Start from Figure 9.10a, gradually increase α, and keepobservations
observations
your eye at the vertical bar separating the acceptance and rejection region. It will move to
the left until it hits the observed test statistic Zobs . At this point, our decision changes,
and we switch from case 1 (Figure 9.10a) to case 2 (Figure 9.10b). Observed
Increasing α further, we
pass the Z-statistic and start accepting the null hypothesis. data point
What happens at the border of α-to-accept and α-to-reject? Definition 9.9 says that this
borderline α is the P-value,
P = α. Set of possible results

Also, at this border our observed Z-statistic coincides with the critical value zα ,
Figure 8.4: Interpretation of the
A p-value (shaded p-value
green area) isas
thethe probability
probability that we observe
of an observed
Zobs = zα , result assuming that the null hypothesis is true.
a test statistic (or
T more
thatextreme)
is at least as extreme as Tobs given that H0 is true.
and thus,
P = α = P {Z ≥ zα } = P {Z ≥ Zobs } .
In this formula, Z is any Standard Normal random variable, and Zobs is our observed test
statistic, which is a concrete number, computed from data. First, we compute Zobs , then
use Table A4call”. A good decision is to collect more data until a more
to calculate definitive answer
can be obtained. P {Z ≥ Zobs } = 1 − Φ(Zobs ).

We
P-values for compute
the left-tail and forpthe
bytwo-sided
fixingalternatives
Zobs =arezcomputed
α and similarly,
selecting p =
as given α such that for a
in Table 9.3.
one-sided right-tail alternative we have
This table applies to all the Z-tests in this chapter. It can be directly extended to the case
of unknown standard deviations and T-tests (Table 9.4).
p = α = P (Z ≥ zα ) = P (Z ≥ Zobs ) = 1 − Φ(Zobs )
Understanding P-values
where Z is standard normally distributed and Zobs is the test statistic. The
Looking atcomputation
Tables 9.3 and 9.4,ofwepsee
is that
similar
P-valuefor theprobability
is the one-sided left-tail
the two-sided case
of observing and
a test
statistic at least as extreme as Zobs or tobs . Being “extreme” is determined by the alterna-
as well as for T-tests.
tive. For a right-tail alternative, large numbers are extreme; for a left-tail alternative, small
We summarize the com-
Hypothesis Alternative putation of the p-value for
P-value Computation
H0 HA Z-test on the left where
right-tail we distinguish the three
P {Z ≥ Zobs } 1 − Φ(Zobs )
θ > θ0 different cases for the al-
θ = θ0
left-tail
P {Z ≤ Zobs } Φ(Zobs ) ternative hypothesis HA .
θ < θ0
two-sided From the definition of the
P {|Z| ≥ |Zobs |} 2(1 − Φ(|Zobs |))
θ ̸= θ0
p-value, it is also clear
that
TABLE 9.3: P-values for Z-tests.

p = P (observing test statistic T that is at least as extreme as Tobs | H0 ).

In Figure 8.4 we illustrate this interpretation (the green shaded area is the
p-value here and Tobs is the observed data point). Thus, it is wrong to say
that the p-value tells us something about the probability that H0 is true
(given the observation)! A high p-value tells us that the observed or even
more extreme values of Zobs is not so unlikely (given H0 ), and therefore, we

96
8.5. CHI-SQUARE TESTS

see no contradiction with H0 and do not reject it. Conversely, a low p-value
signals that such an extreme test statistic is unlikely if H0 is true. Since we
really observed it, our data are not consistent with H0 and we reject it.

Example 72: Two-Sample Z-Test of Proportions (revisited)


In Example 71 we computed a test statistic of Zobs = −0.945 for a two-
sided test which compares the quality of parts produced by two different
manufactures. We compute p as

p = P (|Z| ≥ | − 0.945|) = 2(1 − Φ(0.945)) = 0.3472.

This p-value is quite high and indicates that the null hypothesis should not
be rejected. If H0 is true then the chance of observing a value for Z that
is as extreme or more extreme than Zobs is 34%. This is no contradiction
with the assumption that H0 is true.

8.5 Chi-square Tests


Chi-square tests are another popular class of statistical tests that can be used
to support statements about estimated variances, about estimated counts,
about a family of distributions and about independence.

8.5.1 Estimated variance


Recall that an unbiased estimator for the population variance σ 2 is given by
the sample variance
n
1 X
S2 = (Xi − X̄)2 .
n−1
i=1
We would like to know the distribution of S 2 to perform test for the hy-
pothesis H0 : σ 2 = σ02 . Since the summands (Xi − X̄)2 are not independent
we cannot directly use the central limit theorem. In fact, for moderate or
small n the distribution is not Normal at all (for large n a generalization of
the central limit theorem gives an approximate normal distribution). If the
Xi are independent and normally distributed observations with V (Xi ) = σ 2
then the distribution of
n  2
(n − 1)S 2 X Xi − X̄
=
σ2 σ
i=1
is Chi-square with n − 1 degrees of freedom. The Chi-square distribution is
a special case of the Gamma distribution and its density is given by
1
f (x) = xν/2−1 e−x/2 , x > 0,
Γ(ν/2)2ν/2

97
8.5. CHI-SQUARE TESTS

and f (x) = 0 whenever x ≤ 0. Here ν > 0 is a parameter that equals the


degrees of freedom. Note that the mean is equal to ν and the variance is 2ν.
Thus, for a level-α test with H0 : σ 2 = σ02 we compute

(n − 1)S 2
χ2obs =
σ02

and compare the value χ2obs with the critical values of the Chi-square distri-
bution. More concretely,

• for a right-tail alternative HA : σ 2 > σ02 , we reject H0 if χ2obs ≥ χ2α ,

• for a left-tail alternative HA : σ 2 < σ02 , we reject H0 if χ2obs ≤ χ21−α ,

• for the two-sided alternative HA : σ 2 6= σ02 , we reject H0 if either


χ2obs ≥ χ2α/2 or χ2obs ≤ χ21−α/2

where χ2α is such that α = P (χ2 > χ2α ). We summarize the chi-square test
for the population variance below and omit an example here since this test
Statistical Inference I 291
works exactly as the Z and T tests.
Null Alternative Test Rejection
P-value
Hypothesis Hypothesis statistic region
! "
σ 2 > σ02 χ2obs > χ2α P χ2 ≥ χ2obs
(n − 1)s2 ! "
σ 2 = σ02 σ 2 < σ02 χ2obs < χ2α P χ2 ≤ χ2obs
σ02 # ! "
χ2obs ≥ χ2α/2 or 2 min P χ2 ≥ χ2obs ,
σ 2 ̸= σ02 ! "$
χ2obs ≤ χ21−α/2 P χ2 ≤ χ2obs

TABLE 9.5: χ2 -tests for the population variance


8.5.2 Observed Counts
this case,
The chi-square test can #also
α $ be used
% to compare
& observed and expected
P =2 = 2P χ2 ≤ χ2obs = 2F (χ2obs ). (9.27)
2
counts via the chi-square statistic
So, the P-value is either given by (9.26) or by (9.27), depending on which one of them is
2
N
smaller and which boundary hits χ Xobs first. We can write this2in one equation as
2 (Obs(k) − Exp(k))
' % χ = & % &( % & (8.1)
Exp(k)
P = 2 min P χ2 ≥ χ2obs , P χ2 ≤ χ2obs = 2 min F (χ2obs ), 1 − F (χ2obs ) ,
k=1
where F is the cdf of χ2 distribution with ν = n − 1 degrees of freedom.
where N is the number of categories or groups of data defined depending on
Testing
our procedures
testing that we
problem. have just
Obs(k) derived
is the are summarized
observed numberinofTable 9.5. The
samples same
that tests
belong
can also be used for the standard deviation because testing σ 2 = σ02 is equivalent to testing
toσ category
= σ0 . k, and Exp(k) = E[Obs(k) | H0 is true] is the expected number
of samples that fall into category k if the null hypothesis H0 is true. This
isExample
always a9.41.
one-sided, 2 indicate
Refer to right-tail
Example 9.40 test,
on [Link]
288. The only low values
90% confidence of χconstructed
interval
therethe
that contains the suggested
observed countsvalue
are ofclose
σ = 2.2. Then, by
to what weduality
expectbetween
themconfidence intervals
to be under the
and tests, there should be no evidence against this value of 2σ. Measure the amount of
null hypotheses. On the contrary, large values of χ occur if the observed
evidence against it by computing the suitable P-value.
counts are far from the expected counts, which indicates inconsistency of the
data and H0 . Therefore, the test can be seen as a “frequency comparison”
Solution. The null hypothesis H0 : σ = 2.2 is tested against HA : σ 2 ̸= 2.2. This is a two-
sided test because we only need to know whether the standard deviation equals σ0 = 2.2 or
not. 98
Compute the test statistic from the data in Example 9.40,

(n − 1)s2 (5)(6.232)
χ2obs = = = 6.438.
σ02 2.22

Using Table A6 with ν = n − 1 = 5 degrees of freedom, we see that χ20.80 < χ2obs < χ20.20 .
Therefore, % & % &
8.5. CHI-SQUARE TESTS

(compare the data histogram with a fitted density or mass function). A


level α rejection region for this chi-square test is [χ2α , +∞) and the p-value
is computed as
p = P (χ2 ≥ χ2obs ).
It can be shown that for large n
Obs(k) − Exp(k)
p
Exp(k)
follow an standard normal distribution and thus the null distribution of χ2
converges to the Chi-square distribution with (N − 1) degrees of freedom, as
the sample sizePincreases to infinity. We subtract one degree from N because
we know that k Obs(k) = n and once any N − 1 are known, the remaining
one is uniquely determined.
As a rule of thumb, we require an expected count of at least 5 in each
category, i.e. Exp(k) ≥ 5 for all k and at least 3 bins (N ≥ 3). Moreover, it
is often beneficial if
Exp(1) ≈ . . . ≈ Exp(N ).
After defining the bins and computing χ2 we can use the Chi-square distribu-
tion to construct rejection regions and compute p-values to decide whether
H0 can be rejected or not.

Example 73: Testing a distribution


We want to test whether the observed data belong to a particular dis-
tribution. For example, we may want to test whether a sample comes
from the Normal distribution, whether interarrival times are Exponential
and counts are Poisson, whether a random number generator returns high
quality Standard Uniform values, or whether a die is unbiased.
In general for n samples X1 , . . . , Xn with distribution F we test H0 : F =
F0 vs HA : F 6= F0 , where F0 is some given distribution.
To conduct the test, we take all possible values of X under F0 (the sup-
port of F0 ) and split them into N bins B1 , . . . , BN such that the expected
number of samples, that fall into a bin, is at least 5. It could be that the
left-most interval is B1 = (−∞, a1 ] or the right-most interval is [aN −1 , ∞)
or both. The observed count for the k-th bin Bk is
Obs(k) = |{i = 1, . . . , n : Xi ∈ Bk }|,
where for a set A we use |A| to denote the number of elements of A. Hence
PN
Obs(k) = n. If H0 is true and all Xi have the distribution F0 , then
k=1
Obs(k), the number of ”successes” in n trials, has Binomial distribution
with parameters n and
pk = F0 (Bk ) = P (Xi ∈ Bk | H0 ).

99
8.5. CHI-SQUARE TESTS

Then, the corresponding expected count is the expected value of this Bi-
nomial distribution, Exp(k) = npk . We compute χ2 as defined in Eq. 8.1
and conduct the test.
Let us now become more concrete: Suppose that after losing a large
amount of money, an unlucky gambler questions whether the game was
fair and the die was really unbiased. The last 90 tosses of this die gave
the following results

Score 1 2 3 4 5 6
Frequency 20 15 12 17 9 17

We test H0 : F = F0 vs HA : F 6= F0 , where F is the distribution of


the number of dots on the die, and F0 is such that P (X = i) = 1/6 for
i ∈ {1, 2, . . . , 6}. We choose Bi = {i} and then the observed counts are
given by the table above and the expected counts are n/6 = 90/6 = 15 for
all bins (all more than 5). We compute
P (Obs(k)−Exp(k))
2

χ2 = N k=1 Exp(k)
(20−15)2 (15−15)2 (12−15)2 (17−15)2 (9−15)2 (17−15)2
= 15 + 15 + 15 + 15 + 15 + 15 = 5.2
and find (from the table of the chi-square distribution for N − 1 = 5
degrees of freedom) that the p-value

p = P (χ2 ≥ 5.2) ∈ [0.2, 0.8].

Hence there is no significant evidence to reject H0 , and therefore, no


evidence that the die was biased.

Assume that we have the same situation as in the previous example, i.e., we
suppose that X1 , ..., Xn have distribution F0 but we also used X1 , ..., Xn to
fit the parameters of F0 (e.g. we computed the MLE θ̂). Assume F0 has m
parameters (m ≥ 1).

One can show that in this case, if H0 is true, then χ2 converges to a chi-
square distribution with k − m − 1 degrees of freedom. Thus, we have to
subtract the number of estimated parameters from the degrees of freedom
and defined the rejection region accordingly. This type of test is called
goodness of fit test.

In general, it is a major drawback of the test that there is no clear pre-


scription for the interval selection and one can come to different conclusions
depending on how the intervals are chosen. However, it can be applied to
any hypothesized distribution and is often used to test whether artificially
generated pseudo-random numbers are U (0, 1)-distributed and independent.

100
8.5. CHI-SQUARE TESTS

8.5.3 Testing independence


In many practical applications we would like to test if there is a significant
association between two features. This helps to understand the cause-and-
effect relationships. For example, is it true that smoking causes lung cancer?
Do the data confirm that drinking and driving increases the chance of a
traffic accident? We can use a chi-square test to verify
H0 : Factors A and B are independent
vs
HA : Factors A and B are dependent.
In general, we assume that there are two sets of categories A1 , . . . , Ak and
B1 , . . . , Bm where Ai ∩ Aj = ∅, Bi ∩ Bj = ∅ for any i 6= j and independence
is understood as for random variables, i.e., factors A and B are independent
if any randomly selected x belongs to categories Ai and Bj independently
of each other. Therefore, we are testing
H0 : P (x ∈ Ai ∩ Bj ) = P (x ∈ Ai )P (x ∈ Bj ) for all i, j
vs
HA : P (x ∈ Ai ∩ Bj ) 6= P (x ∈ Ai )P (x ∈ Bj ) for some i, j.

Example 74: Smokers and lung cancer


We have k = m = 2 such as A1 = set of smokers, A2 = set of non-
smokers, B1 = set of persons with lung cancer, B2 = set of persons
without lung cancer.

We next assume the following observed counts:


B1 B2 . . . Bm row total
A1 n11 n12 . . . n1m n1·
A2 n21 n22 . . . n2m n2·
... ... ... ... ... ...
Ak nk1 nk2 . . . nkm nk·
column total n·1 n·2 . . . n·m n·· = n
P P
where ni· = i nij and n·j = j nij are the row and column totals. This
table defines the observed counts Obs(i, j) = nij . We then estimate the
probabilities
P (x ∈ Ai ) and P (x ∈ Bj ) ∀i, j
as
ni· n·j
Pb(x ∈ Ai ) = and Pb(x ∈ Bj ) = , ∀i, j.
n n
If H0 is true then P (x ∈ Ai ∩ Bj ) = P (x ∈ Ai )P (x ∈ Bj ) and thus
Exp(i, j) = n · P (x ∈ Ai ∩ Bj ).

101
8.5. CHI-SQUARE TESTS

Based on the information that we have, we can only estimate


d j) = n ni· n·j = ni· n·j .
Exp(i,
n n n
Next, we compute the test statistic
 Inference II
Statistical 2 313
X k X m d j)
Obs(i, j) − Exp(i,
2
You can always checkχthat
obs =all the row totals, the d
column totals, and the whole table total
are the same for the observedi=1 the expected Exp(i,
andj=1 j) if there is a mistake, catch it
counts (so,
here).
and consider its corresponding chi square distribution. It can be shown that
(160 − 120)2 (240 − 280)2Inference
Statistical (140 −II180)2 (460 − 420)2 313
χ2obs =of freedom +
the degrees are (k − 1)(m+ − 1) due to + some linearly dependent
= 31.75.
120 280 180 420
constraints
From Table A6in with
the (2
sum of −
− 1)(2 differences.
1) = 1 d.f., we find that the P-value P < 0.001. We have a
You can always check that all the row totals, the column totals, and the whole table total
significant evidence
are the same that
for the an email
observed andhaving an attachment
the expected counts is somehow
(so, if there related to being
is a mistake, spam.
catch it
Example
Therefore,
here). 75:piece
this Internet Shopping
of information can be used in anti-spam filters. ♦
2 2 2 2
(160 − 120) (240 − 280) (140 − 180) (460 − 420)
A web χdesigner
2
obs =
suspects+ that the chance
+ for an +internet shopper to make
= 31.75.
120 280 180 420
aExample
purchase through
10.5 (Internether web site
shopping varies
on depending
different days on
of the
the day of
week). the week.
A web
From Table A6 with (2 − 1)(2 − 1) = 1 d.f., we find that the P-value P < 0.001. We have a
de-
signer
To suspects
test this that the chance
claim, she for an internet
collects data shopper to
during make
one a purchase
week, whenthrough
the her web
web site
significant evidence that an email having an attachment is somehow related to being spam.
site varies depending on the day of the week. To test this claim, she collects data during
recorded
Therefore,3758 hits.
this piece of information can be used in anti-spam filters. ♦
one week, when the web site recorded 3758 hits.
Observed Mon Tue Wed Thu Fri Sat Sun Total
ExampleNo 10.5 (Internet shopping
purchase 399 261 on different
284 263 days 393 of 531the502
week).
2633A web de-
signer suspects
Single that the chance for
purchase 119an internet
72 97shopper
51 to make
143 a145 purchase
150 through
777 her web
site variesMultiple
depending on the day39of the50week.20To test
purchases 15 this41claim,
97she 86
collects348
data during
one week,Total
when the web site recorded
557 383 3758 hits.
401 329 577 773 738 3758

Testing independence Mon Tue


Observed (i.e., probability Wed aThu
of making purchaseFri or Satmultiple
Sun Total
purchases is the
Testing
same on anyNo purchase
independence 399
(i.e.,
day of the week), 261
probability
we compute 284 263
of making
the estimated 393 531 502 2633
a purchase
expected counts, or multiple
Single purchase 119 72 97 51 143 145 150 777
purchasesMultiple
is the same
!purchases
on any
(ni·39 day
)(n·j ) 50 of the week), we compute
863. 348 estimated
the
Exp(i, j) = for20 i =151, . . .41
, 7, j97
= 1, 2,
expected counts,
Total n
557 383 401 329 577 773 738 3758
Expected
Testing independence
d j) = (i.e., i· nTue
Monnprobability
·j nofWedn·j Thu
i· making a purchase Frior multiple
Sat purchases
Sun Total
is the
No Exp(i,
purchase
same on any day of the week),
n
390.26 268.34
we
=
compute 280.96 for i
230.51 = 1, . .
404.27 .
the estimated expected counts,
, 7, j =
541.59 1, 2, 3.
517.07 2633
Single purchase
n n n
115.16 79.19 82.91 68.02 119.30 159.82 152.59 777
! 51.58 (ni· )(n·j )
We get purchases
Multiple Exp(i, j) = 35.47 37.13
for 30.47
i = 1, . . .53.43
, 7, j = 1, 71.58
2, 3. 68.34 348
Total 557 n
383 401 329 577 773 738 3758
Expected Mon Tue Wed Thu Fri Sat Sun Total
Then, the test statistic is
No purchase 390.26 268.34 280.96 230.51 404.27 541.59 517.07 2633
Single purchase 2 115.16 390.26)282.91 68.02
(399 −79.19 (86 − 68.34)
119.30 2
159.82 152.59 777
Multiple purchases χobs =51.58 35.47 37.13 + . . . +30.47 53.43 = 71.58 60.79, 68.34 348
390.26 68.34
Total 557 383 401 329 577
and it has (7 − 1)(3 − 1) = 12 degrees of freedom. From Table A6, we find that the P-value 773 738 3758
isThen,
P < 0.001, so statistic
the test indeed, there
is is significant evidence that the probability of making a single
purchase
and the ortest
multiple purchases
statistic is varies during the week. ♦
(399 − 390.26)2 (86 − 68.34)2
χ2obs = + . . . + = 60.79,
(399 390.26 2
− 390.26) − 68.34)2
(8668.34
2
Matlab χ
Demo.
and it has (7obs =
The χ 2
test for independence + .
− 1)(3 − 1) = 12 degrees of freedom. From Table. . +
takes only a few we =
A6,lines of 60.79.
findcode
thatinthe
MATLAB.
P-value
Foris Pexample,
< 0.001,here
390.26to Example
is the solution
so indeed, there is significant 10.4. that68.34
evidence the probability of making a single
purchase or multiple
− 1)(3 purchases varies during of thefreedom
week. ♦
WeX = have
[160 (7 240; 140 − 1) = 12
460]; degrees %and
Matrix get ofa very low p-value
observed counts
<Row 0.001 which means
= sum(X’)’; Col = sum(X);that weTot should
= sum(Row);reject H%0 ,Row [Link]
column is totals
significant
kMatlab
= length(Col);
evidence that the
Demo. χ=2 test
Themprobability
length(Row); of makingtakes
for independence only %a purchase
a single Dimensions
few lines of or ofmultiple
code the
in table
[Link]-
eFor= example,
zeros(size(X));
here is the solution to Example 10.4. % Expected counts
chases varies during the week.
for i=1:k; for j=1:m; e(i,j) = Row(i)*Col(j)/Tot; end; end;
X = [160
chisq 240; ^
= (X-e). 2140
./e;460]; %% Matrix
Chi-square of observed
terms counts
Row = sum(X’)’;
chistat = sum(sum(chisq)); Col = sum(X); Tot = sum(Row); % Row and column
% Chi-square statistic totals
k = length(Col);
Pvalue m = length(Row); 102
= 1-chi2cdf(chistat,(k-1)*(m-1)) %% Dimensions
P-value of the table
e = zeros(size(X)); % Expected counts
for i=1:k; for j=1:m; e(i,j) = Row(i)*Col(j)/Tot; end; end;
chisq = (X-e).^ 2./e; % Chi-square terms
chistat = sum(sum(chisq)); % Chi-square statistic
Pvalue = 1-chi2cdf(chistat,(k-1)*(m-1)) % P-value
8.5. CHI-SQUARE TESTS

In the case of small sample sizes and very unequally distributed data (among
the cells), Fisher’s exact test can be used instead to test for independence.

103
8.5. CHI-SQUARE TESTS

104
Bibliography

[1] Michael Baron. Probability and statistics for computer scientists. CRC
Press, 2013.

[2] Cambridge Mathematics. Skewed usage of skewed distribu-


tion. [Link]
distribution/, Last accessed on 2024-04-29.

[3] Radford M. Neal. Sta 247 - week 2 lecture summary. http://


[Link]/~radford/sta247.F11/[Link], Last ac-
cessed on 2020-02-11.

105

You might also like