Course Notes Book
Course Notes Book
Sofya Titarenko
2
Contents
1 Probability Distributions 9
1.1 Random variables and distributions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
1.2 Relationship between distributions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
1.3 Joint, Marginal, and Conditional Distributions. Convolution. . . . . . . . . . . . . . . . 12
2 Moments 19
2.1 Raw Moments and Central Moments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
2.2 Sample Moments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
2.3 Covariance and Correlation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
4 Random Vectors 37
4.1 Mean and covariance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
4.2 Bivariate normal distribution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
6 Limit Theorems 55
6.1 Markov’s Inequality . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
6.2 Chebyshev’s Inequality . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
6.3 Sequence of Random Variables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
6.4 Weak Law of Large Numbers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
6.5 Central Limit Theorem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
6.5.1 The Theorem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
6.5.2 Interpretation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
6.5.3 Examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68
7 Estimation Methods 75
7.1 Basic Concepts . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75
7.2 Methods of Moments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
7.2.1 Examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
3
4 CONTENTS
12 Appendix 203
12.1 Vectors . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 203
12.2 Matrices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 203
12.3 Eigenvalues . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 208
12.4 Matrix calculus . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 208
12.5 Real Analysis Reminders . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 209
12.6 Probability Reminders . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 210
12.7 Statistical Methods Reminders . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 221
Introduction
Welcome to MATH2701!
This module introduces key concepts in probability and statistical inference. Please refer to the table
below for the approximate layout of lectures, workshops, and lecture notes, and see the Assessment
Information section for detailed information on grading and deadlines.
Before each workshop, you are expected to attempt the problems from the Exercise sheets;
before the workshop, we will hold a vote to decide which problems to explain-the most popular will be
discussed.
All solutions will be unlocked at the end of the week in which the workshop occurs. Note
that some workshop problems require R, while others are pen-and-paper exercises.
The core literature is listed below, and additional references will be provided at the end of each
chapter to support your learning.
Core Literature
6 CONTENTS
1. Rice, J. A. (2007). Mathematical statistics and data analysis (3rd international ed.). Brooks/Cole.
3. Verzani, J. (2014). Using R for introductory statistics (2nd ed.). CRC Press.
CONTENTS 7
Probability Distributions
Overview
Gamma function
The gamma function sometimes called the factorial function. The gamma function was first investigated
in 1729 by Leonhard Euler (1707-1783),
√ aPSwiss mathematician who also introduced the notations e for
the exponential function, i for −1 and for summation. Euler gave the definition
Z ∞
Γ(z) = xz−1 e−x dx.
0
The notation Γ(z) was introduced in 1814 by Adrien-Marie Legendre (1752-1833), a French mathe-
matician who was one of the people responsible for the development of the method of least squares.
Legendre referred to the above integral definition of Γ(z) as the Eulerian integral of the second kind.
√
It can be shown that Γ( 12 ) = π.
Γ(z) is a generalises factorial function that satisfies
Gamma distribution
If X ∼ Gamma(α, λ), then the expected value, variance and the PDF are
α α
E[X] = , Var[X] = 2 ,
λ λ
2 2 α(1 + α)
E[X ] = Var[X] + E[X] = ,
λ2
λα α−1 −λx
fX (x) = x e , x ≥ 0.
Γ(α)
9
10 CHAPTER 1. PROBABILITY DISTRIBUTIONS
F-distribution
If X ∼ F (d1 , d2 ) ((d1 , d2 ) > 0 are degrees of freedom), then the expected value, variance, and the PDF
are:
d2
E[X] = , for d2 > 2,
d2 − 2
2d22 (d1 + d2 − 2)
Var[X] = , for d2 > 4,
d1 (d2 − 2)2 (d2 − 4)
s
(d1 x)d1 dd21
(d1 x + d2 )d1 +d2
fX (x) = , x ≥ 0 (x > 0 if d1 = 1)) .
xB d21 , d22
It can be derived from the standard normal by summing the squares of independent standard normal
variables:
k
X
X= Zi2 ∼ χ2k , where Zi ∼ N (0, 1)
i=1
(X1 /d1 )
F = ∼ F (d1 , d2 ), where X1 ∼ χ2d1 , X2 ∼ χ2d2
(X2 /d2 )
It is commonly used in ANOVA and regression model comparison to test whether group variances
differ significantly.
The Normal distribution can also be used to approximate discrete distributions, which is incredibly
useful for simplifying calculations for large samples.
12 CHAPTER 1. PROBABILITY DISTRIBUTIONS
Optional
Binomial Approximation by Normal
The Binomial distribution X ∼ Binomial(n, p), can be approximated by the Normal distribution
as
approx.
N µ = np, σ 2 = np(1 − p)
X ∼ Binomial(n, p) ∼
The approximation is generally considered good when both np and n(1 − p) are sufficiently large. A
common rule of thumb is np > 5 and n(1 − p) > 5.
For greater accuracy, especially when working with discrete values, the continuity correction is
often applied.
The Poisson distribution X ∼ Poisson(λ), which models the number of events in a fixed interval,
can be approximated by the Normal distribution when the rate parameter λ is large.
approx.
N µ = λ, σ 2 = λ
X ∼ Poisson(λ) ∼
The approximation is generally good for λ > 20, and becomes excellent for larger values of λ.
In many experiments, we measure two or more random variables simultaneously. For example, we
might record a person’s height (X) and weight (Y), or the daily temperature (X) and energy
consumption (Y). To understand the relationship between these variables, we must study their joint
behavior.
The complete joint behavior of two random variables, X and Y , is captured by their joint cumulative
distribution function.
Definition 1.3.1 (Joint Cumulative Distribution Function). The joint cumulative distribu-
tion function of random variables X and Y is defined as:
FX,Y (x, y) = P (X ≤ x, Y ≤ y)
This definition applies universally, regardless of whether X and Y are discrete, continuous, or a mixture
of both.
Discrete Random Variables When X and Y are discrete random variables, their joint distribution
is described by a joint probability mass function.
1.3. JOINT, MARGINAL, AND CONDITIONAL DISTRIBUTIONS. CONVOLUTION. 13
Definition 1.3.2 (Joint Probability Mass Function). For discrete random variables X and Y ,
their joint probability mass function is:
pX,Y (x, y) = P (X = x, Y = y)
Properties:
The individual distribution of X, ignoring Y (or Y , ignoring X), is called the marginal PMF.
Definition 1.3.3 (Marginal PMF). The marginal probability mass function of X is:
X
pX (x) = P (X = x) = pX,Y (x, y)
all y
P
Similarly, pY (y) = all x pX,Y (x, y).
Theorem 1.3.4 (Independence of Discrete Random Variables). Two discrete random variables
X and Y are independent if and only if their joint PMF factors into the product of their marginals
for all x and y:
pX,Y (x, y) = pX (x)pY (y) ∀x, y
Definition 1.3.5 (Conditional PMF). The conditional probability mass function of Y , given
that X = x (provided pX (x) > 0), is:
pX,Y (x, y)
pY |X (y|x) = P (Y = y|X = x) =
pX (x)
Continuous Random Variables For two continuous random variables X and Y , their joint distri-
bution is defined by a joint probability density function.
Definition 1.3.6 (Joint Probability Density Function). The joint probability density func-
tion of continuous random variables X and Y is a function fX,Y (x, y) with the properties:
Definition 1.3.7 (Joint Cumulative Density Function). The joint cumulative distribution
function for continuous random variables is:
Z x Z y
FX,Y (x, y) = P (X ≤ x, Y ≤ y) = fX,Y (u, v)dvdu
−∞ −∞
The individual or marginal PDF is found by ”integrating out” the other variable.
Definition 1.3.8 (Marginal PDF). The marginal probability density function of X is:
Z ∞
fX (x) = fX,Y (x, y)dy
−∞
R∞
Similarly, fY (y) = −∞ fX,Y (x, y)dx.
fX,Y (x, y)
fY |X (y|x) =
fX (x)
Summary Table 1.1 provides a concise comparison of key concepts for discrete and continuous joint
distributions and the notation used in this module.
On April 15, 1912, the ocean liner Titanic collided with an iceberg and sank. Out of the 2201
passengers on board, 1490 died.
1.3. JOINT, MARGINAL, AND CONDITIONAL DISTRIBUTIONS. CONVOLUTION. 15
Table 1.2: Frequency distribution of Titanic passengers by class and survival status.
The question as to whether passenger class was related to survival has been discussed extensively.
If a passenger is selected at random, what is the chance that he/she is a non-crew member who
survived?
There are two random variables in play here:
(
1 First class
0 if passenger survived 2 Second class
X= , Y =
1 if passenger died
3 Third class
4 Crew
We define the joint PMF of X and Y as:
Table 1.3: Joint probability distribution of survival status and passenger class.
To find the probability that a randomly selected passenger is a non-crew member who
survived, we sum the probabilities for first, second, and third class survivors:
Finally, one maybe interested to calculate the conditional probabilities. Using the definition
of the conditional probability, we obtain
So we see that probability of survival across passenger class is decreasing, although we should be
careful while making remarks like this and consider other factors present.
A certain process for producing an industrial chemical yields a product that contains two main
types of impurities. Let X denote the proportion of impurities of Type I and Y denote the
proportion of impurities of Type II. Suppose that the joint density of X and Y can be modeled
as,
2(1 − x) 0 ≤ x ≤ 1, 0 ≤ y ≤ 1
f (x, y) =
0 otherwise
Compute P (0 ≤ X ≤ 0.5, 0.4 ≤ Y ≤ 0.7), fX (x) and fY (y).
Using the integration
R 0.5 R 0.7
P (0 ≤ X ≤ 0.5, 0.4 ≤ Y ≤ 0.7) = 0 0.4 2(1 − x)dydx
R 0.5 R 0.5
= 0 (0.7 − 0.4) × 2(1 − x)dx = 0.3 × 0 2(1 − x)dx
0.5
= 0.3 × −(1 − x)2 0 = 0.3 × −(1 − 0.5)2 + (1 − 0)2
Finding fX (x)
R∞ R∞
fX (x) = −∞ f (x, y)dy = −∞ 2(1 − x)1{0≤x≤1,0≤y≤1} dy
R1
= 0 2(1 − x)1{0≤x≤1} dy
= 2(1 − x)1{0≤x≤1}
hence:
2(1 − x) if 0 ≤ x ≤ 1
fX (x) =
0 otherwise
Finding fY (y)
R∞ R∞
fY (y) = −∞ f (x, y)dx = −∞ 2(1 − x)1{0≤x≤1,0≤y≤1} dx
R1 1
−(1 − x)2 0 = 1{0≤y≤1} × 1
= 0 2(1 − x)1{0≤y≤1} dx = 1{0≤y≤1}
1.3. JOINT, MARGINAL, AND CONDITIONAL DISTRIBUTIONS. CONVOLUTION. 17
hence:
1 if 0 ≤ y ≤ 1
fY (y) =
0 otherwise.
Convolution
Definitions Assume X and Y are discrete and independent r.v. with the PMFs pX (x) and
pY (y). Let Z = X + Y .
X
pZ (z) = pX (x)pY (z − x),
x
X
pZ (z) = pY (y)pX (z − y).
y
Examples
Example (Poisson Distribution)
λx e−λ µy e−µ
pX (x) = , x = 0, 1, 2, . . . , and pY (y) = , y = 0, 1, 2, . . . .
x! y!
As X and Y are discrete independent random variables, then the convolution formula becomes
u
X X λu−v e−λ µv e−µ
P(U = u) = pX (u − v)pY (v) = ·
v
(u − v)! v!
v=0
u
−(λ+µ)
X λu−v µv
=e .
(u − v)!v!
v=0
n
X n!
The binomial formula (a + b)n = an−r br shows that
r!(n − r)!
r=0
(λ + µ)u e−(λ+µ)
P(U = u) = , u = 0, 1, 2, . . . .
u!
Hence X + Y ∼ Poisson(λ + µ).
18 CHAPTER 1. PROBABILITY DISTRIBUTIONS
If X and Y are independent and identically distributed geometric(θ) random variables, what is
the distribution of U = X + Y ?
Here X ∼ geometric(θ) and Y ∼ geometric(θ) with probability functions
and
pY (y) = θ(1 − θ)y , y = 0, 1, 2, . . . .
As X and Y are discrete independent random variables, then the convolution formula becomes
X u
X
P (U = u) = pX (u − v)pY (v) = θ(1 − θ)u−v θ(1 − θ)v
v v=0
u
X
= θ2 (1 − θ)u = (u + 1)θ2 (1 − θ)u .
v=0
Moments
For the case of a discrete random variable with probability distribution p(x)
X
µ0r = E[X r ] = xr p(x).
x
We call µ0r = E(X r ) the rth moment of X about the origin or raw.
Definition 2.1.2. The rth moment about the mean (or rth central moment) of a random variable
X, denoted by µr , is the expected value of (X − µX )r :
µr = E[(X − µX )r ].
Note that:
µ1 = E[(X − µX )] = 0,
µ2 = E[(X − µX )2 ] = Var[X].
A continuous random variable X has probability density function fX (x) = cx(1 − x) for 0 < x < 1.
Obtain:
R 1 c, E(X) and Var(X).
Since 0 fX (x) dx = 1, it follows that
Z 1 1
1 2 1 3 1
cx(1 − x) dx = cx − cx = c,
0 2 3 0 6
19
20 CHAPTER 2. MOMENTS
+∞ 1
3 4 1 1
Z Z
2 3
E[X] = xfX (x) dx = 6x (1 − x) dx = 2x − x = ,
−∞ 0 2 0 2
Z +∞ Z 1 1
3 4 6 5 3
E[X 2 ] = x2 fX (x) dx = 6x3 (1 − x) dx = x − x = .
−∞ 0 2 5 0 10
∞
λα α−1 −λx
Z
fX (x) = x e , x > 0 =⇒ E[X] = xfX (x) dx,
Γ(α) 0
Z ∞
λα α−1 −λx
E[X] = x x e dx,
0 Γ(α)
Z ∞ α Z ∞
λ α−1 −λx Γ(α)
x e dx = 1 =⇒ xα−1 e−λx dx = α ,
0 Γ(α) λ
Z ∞ 0
Γ(α + 1)
xα+1−1 e−λx dx = ,
0 λα+1
λα Γ(α + 1) λα αΓ(α) α
E[X] = · α+1
= · = ,
Γ(α) λ Γ(α) λα+1 λ
Z ∞ Z ∞ α
λ
E[X 2 ] = x2 fX (x) dx = xα+1 e−λx dx,
0 0 Γ(α)
λα Γ(α + 2) λα α(α + 1)Γ(α) α(α + 1)
E[X 2 ] = · α+2
= · α+2
= ,
Γ(α) λ Γ(α) λ λ2
α(α + 1) α2 α
Var[X] = E[X 2 ] − E[X]2 = 2
− 2 = 2.
λ λ λ
Definition 2.1.3. For a random variable X with mean µ and standard deviation σ, the population
skewness is defined as:
E[(X − µ)3 ]
γ1 = .
σ3
Kurtosis measures the “tailedness” of a distribution - how heavy or light the tails are compared to a
normal distribution.
2.2. SAMPLE MOMENTS 21
(2) they have the same distribution, i.e., FX1 (x) = FX2 (x) = . . . = FXn (x), ∀x ∈ R.
Definition 2.2.2. If X1 , . . . , Xn are i.i.d. random variables, the rth sample moment (also called the
raw moment or moment about the origin) is defined as:
n
1X r
m0r = µ̂0r = Xi .
n
i=1
Definition 2.2.3 (Custom Example: Sample First and Second Moments). Given the dataset
{1, 2, 2, 2}, derive m0r for r = 1, 2.
Indeed,
n
1X r
m0r = µ̂0r = Xi .
n
i=1
Note that:
n
1X
m1 = (Xi − X) = 0.
n
i=1
is related to skewness:
m3
Skewness = .
(m2 )3/2
The fourth central moment:
n
1X
m4 = (Xi − X)4
n
i=1
is related to kurtosis:
m4
Kurtosis = .
(m2 )2
A random sample X1 , X2 , . . . , Xn has sample mean X. Prove that: m3 = m03 − 3m02 m01 + 2(m01 )3 .
Indeed, we know that:
n
1X
m3 = (Xi − X)3
n
i=1
and
n
1X r
m0r = Xi .
n
i=1
Hence:
n
1X
m3 = (Xi − X)3
n
i=1
n
1 X 3 2 3
= Xi − 3XXi2 + 3X Xi − X
n
i=1
n n 2 n n
1 X 3 3X X 2 3X X 1X 3
= Xi − Xi + Xi − X
n n n n
i=1 i=1 i=1 i=1
3 3
= m03 − 0
3Xm2 + 3X − X
= m03 − 3m02 m01 + 2(m01 )3 since X = m01 .
1. m1 - location,
2. m2 - spread,
2.3. COVARIANCE AND CORRELATION 23
3. m3 - skewness,
4. m4 - kurtosis.
Different software packages may use slightly different formulas for sample skewness and kurtosis,
especially when correcting for bias in small samples.
In R, using the e1071::skewness() and e1071::kurtosis() functions with type = 1 gives the
classical moment-based definitions, which align with the formulas above.
library ( e1071 )
skewness ( data , type = 1)
kurtosis ( data , type = 1)
Notice, that sample moments, especially of higher order (e.g., skewness and kurtosis), are highly
sensitive to outliers. This will be discussed later.
The sample moments can serve as estimators for the population moment µr and µ0r . However, some
sample moments are biased estimators. For example, to estimate the population variance, we use
the estimator
n
2 1 X
σ̂ = (Xi − X)2 .
n−1
i=1
It is an unbiased estimator, whereas the second central moment m2 (with denominator n) is biased.
This will be discussed later in the module.
Theorem 2.3.3. Suppose Xi is a sequence of random variables with joint pdf f (x)
"N # N
X X X
Var ai Xi = a2i Var[Xi ] + 2 ai aj Cov[Xi , Xj ].
i=1 i=1 1≤i<j≤N
24 CHAPTER 2. MOMENTS
Theorem 2.3.4. For the random variables (X, Y, Z) and constant a the following holds:
1. Cov[aX, Y ] = a Cov[X, Y ],
4. Cov[a, X] = 0.
Cov[X, Y ] = 0.
Proof.
We use
Var[U ] = Var[X − 2Y + Z]
= 12 Var[X] + (−2)2 Var[Y ] + 12 Var[Z]
+ 2 · 1 · (−2) Cov[X, Y ] + 2 · 1 · 1 · Cov[X, Z] + 2 · 1 · (−2) Cov[Y, Z]
= 2 + 4 · 4 + 6 − 4 · 1 − 4 · 1 = 16.
Example (Covariance)
• Obtain Cov[X; Y ]
Y
X 1 2 3
-1 0.2 0.1 0.1
0 0.1 0.1 0.1
+1 0.1 0.0 0.2
X
E[X] = x pX (x) = −0.4 + 0 + 0.3 = −0.1,
x
X
E[Y ] = y pY (y) = 0.4 + 0.4 + 1.2 = 2.0,
y
XX
E[XY ] = xy p(x, y)
x y
p(1, y) p(1, y)
P(Y = y|X = 1) = = .
pX (1) 0.3
Thus
0.1 1
P(Y = 1|X = 1) = = ,
0.3 3
0.0
P(Y = 2|X = 1) = = 0,
0.3
0.2 2
P(Y = 3|X = 1) = = .
0.3 3
d −1
fY (y) = fX (g −1 (y)) g (y) for y ∈ g(supp(fX )).
dy
−1 1
= 1 − elog(y θ)
= 1 − y − θ =⇒
1 1
fY (y) = FY0 (y) = y −( θ +1) .
θ
27
28 CHAPTER 3. FUNCTION OF A RANDOM VARIABLE
dx1 dx2
fY (y) = fX (x1 (y)) + fX (x2 (y))
dy dy
dxi 1
where = 0 for i = 1, 2.
dy g (xi (y))
√ dx √ dx
fU (u) = fX (x = − u) √
+ fX (x = + u) √
du x=− u du x=+ u
1 1 −1 1 1 1 1
= − u 2 + + u− 2 = √ , y ∈ (0, 1].
2 2 2 2 2 u
Definition 3.1.3 (Mixed Transformation: Two-to-One and One-to-One). Let Y = g(X) and
partition supp(X) into regions where g has constant monotonicity:
−1
1. A: g is strictly monotonic with inverse gA
2. B: For each y ∈ g(B), the equation g(x) = y has exactly two solutions x1 (y), x2 (y) ∈ B
Then:
2
−1 d −1 X 1
fY (y) = 1{y∈g(A)} fX (gA (y)) gA (y) + 1{y∈g(B)} fX (xi (y)) 0
dy g (xi (y))
i=1
3.1. UNIVARIATE CASE 29
The determinant
∂x ∂x
∂(x, y) ∂u ∂v
J= = ∂y ∂y
∂(u, v) ∂u ∂v
is referred to as a Jacobian (Carl Jacobi (1804–1851) was a German mathematician who worked
extensively with determinants) and the area |Rxy | can be shown to satisfy |Rxy | ≈ |J|∆u∆v.
Indeed, let point A map to (u, v) and B map to (u + ∆u, v) so A has x-coordinate x(u, v) while B
has x-coordinate x(u + ∆u, v) where x = x(u, v) and y = y(u, v) denotes the inverse transformation
from the (u, v) plane to the (x, y) plane. Since
∂y
y(u + ∆u, v) − y(u, v) ≈ ∆u.
∂u
3.2. BIVARIATE CASE 31
−−→
Hence the vector AB satisfies
−−→ ∂x ∂y
AB = ∆u i + ∆u j
∂u ∂u
−−→ −→ ∂x ∂y ∂x ∂y ∂(x, y)
|AB × AC| ≈ − ∆u∆v = ∆u∆v,
∂u ∂v ∂v ∂u ∂(u, v)
∂(x, y) v u
= = v,
∂(u, v) 0 1
while
∂(u, v) 1/y −x/y 2
= = 1/y = 1/v.
∂(x, y) 0 1
Note that SY , the support of (Y1 , Y2 ), is usually found by considering the image of SX under the
transformation Y1 , Y2 . Meaning that ∀(x1 , x2 ) ∈ SX we find (y1 , y2 ) ∈ SY :
x1 = v1 (y1 , y2 ), x2 = v2 (y1 , y2 ).
32 CHAPTER 3. FUNCTION OF A RANDOM VARIABLE
Note
This exercise sheet covers material from Sections 1, 2 and 3 of the course. Please
attempt all questions and prepare your solutions before the Week 3 workshop.
Your solutions will not be assessed. Before the workshop, you will be invited to vote for the
questions you would like to see discussed. The tutor will select the most popular questions for
discussion during the workshop. Complete solutions for all questions will be available afterwards
(by the end of the week). The solutions will be available only in HTML version of the notes.
Exercise 3.1
Let Y1 and Y2 have joint density:
Exercise 3.2
Let X ∼ Uniform(a, b). Using the definition of raw moments, compute:
Exercise 3.3
Let X ∼ Exponential(λ) with probability density function (pdf)
35
36 CHAPTER 3. FUNCTION OF A RANDOM VARIABLE
Exercise 3.4
Let random variables X and Y satisfy Cov(X, Y ) = −8. Define new variables via linear
transformations:
U = 3X − 4, V = 2Y + 7.
Compute Cov(U, V ).
Exercise 3.5
Let X ∼ N (0, 1). If Y = X 2 , find the pdf of Y .
Exercise 3.6
If X has a standard normal distribution and Y , independent of X, has probability density function
1 2
(
ye− 2 y if y > 0,
fY (y) =
0 if y ≤ 0,
What is E[U ]?
Exercise 3.7
Let X and Y be independent random variables, each uniformly distributed on the interval (0, a),
where a > 0. Define the transformed variables
X
U = X 2 + Y 2, V = .
Y
Find the marginal probability density function of the random variable V .
Hint: Consider using the inverse transformation
r
U
X = V Y, Y = 2
.
V +1
Apply the Jacobian inversion rule
∂(x, y) 1
=
∂(u, v) ∂(u, v)
∂(x, y)
Random Vectors
Properties
Theorem 4.1.1. Suppose Y is a random vector with mean µ = E[Y ] and covariance matrix Σ =
Cov[Y ]. Let C be a constant matrix and b a constant vector. Then:
Proof
Proof. Expectation.
37
38 CHAPTER 4. RANDOM VECTORS
By linearity of expectation:
E[CY + b] = E[CY ] + E[b].
Claim: E[CY ] = C E[Y ]. Write Y = (Y1 , . . . , Yn )> and C = (cij ) of size m × n. Then:
n
X
(CY )i = cij Yj .
j=1
Thus:
E[CY ] = C E[Y ] = Cµ.
E[b] = b.
Combining:
E[CY + b] = Cµ + b.
Covariance Matrix.
The covariance matrix of a random vector X is:
Let X = CY + b. Then:
Thus:
Cov[X] = E C(Y − µ)(Y − µ)> C > .
Note (Reminder)
−1 1 1 −ρ
Σ = 2 .
σ (1 − ρ2 ) −ρ 1
1 1 −ρ x1
x> Σ−1 x = x1 x2 .
σ 2 (1 − ρ2 ) −ρ 1 x2
First compute:
1 −ρ x1 x1 − ρx2
= .
−ρ 1 x2 −ρx1 + x2
Then:
x1 − ρx2
= x1 (x1 − ρx2 ) + x2 (−ρx1 + x2 ) = x21 + x22 − 2ρx1 x2 .
x1 x2
−ρx1 + x2
And,
x21 + x22 − 2ρx1 x2
x> Σ−1 x = .
σ 2 (1 − ρ2 )
Finally, we have
2
x1 + x22 − 2ρx1 x2
1
fX (x1 , x2 ) = exp − .
2σ 2 (1 − ρ2 )
p
2πσ 2 1 − ρ2
The Figure below shows the bivariate normal distribution with σ = 2, ρ = 0.5.
4.2. BIVARIATE NORMAL DISTRIBUTION 41
σ12
µ1 ρσ1 σ2
µ= , Σ= ,
µ2 ρσ1 σ2 σ22
The expression above can be derived in the similar way as we did in the previous example.
The Figure below shows the bivariate normal distribution with µ1 = µ2 = 0, σ1 = 1, σ2 = 2,
ρ = 0.5.
42 CHAPTER 4. RANDOM VECTORS
Figure 4.3: PDF for the bivariate normal distribution with µ1 = µ2 = 0, σ1 = 1, σ2 = 2, ρ = 0.5
Real-life example
The Old Faithful geyser dataset is a real-world example of bivariate normality (for an eruption type).
This famous dataset contains records of eruption durations and waiting times between eruptions for the
Old Faithful geyser in Yellowstone National Park.
While the complete dataset exhibits a clear bimodal structure due to two distinct eruption patterns,
we isolate the subset of long eruptions-defined as those lasting 3.5 minutes or longer-to reveal an excellent
approximation of bivariate normal distribution.
Below is the plot showing the contours of the bivariate normal distribution and its image in the 3D
space. The normality can be verified by Mardia test (this is not the topic for this module, I refer to the
further reading all interested Mardia, K. V., Kent, J. T., & Taylor, C. C. (2024). Multivariate Analysis
(2nd ed.). John Wiley & Sons. ISBN: 9781118738023)
4.2. BIVARIATE NORMAL DISTRIBUTION 43
Important properties
We will finish the bivariate normal distribution section by stating the following important properties.
Cor[X1 , X2 ] = 0,
Note (Important)
Cor[X1 , X2 ] = 0,
The above can be proved using moment generating function. We will look into this later, as a
part of the exercise.
44 CHAPTER 4. RANDOM VECTORS
Chapter 5
5.1 Definition
Definition 5.1.1. Given is a random variable X, the function M (t) = E[etX ] is defined moment-
generating function (MGF) of X, if the expectation is defined.
If X is a continuous r.v. with pdf f (x):
Z ∞
tX
M (t) = E[e ] = etx f (x)dx
−∞
45
46 CHAPTER 5. MOMENT GENERATION FUNCTION
MX (t) = E etX
1 2
= et + e2t .
3 3
That is well defined for ∀t ∈ R.
5.2 Properties
The next lemma helps us to find moments of a random variable.
Lemma 5.2.1. If the MGF exists in an open interval containing zero, then
∂ n M (t)
= E[X n ]
∂nt 0
Proof
Proof. Recall the Taylor Expansion of etX at 0,
t2 x2 t3 x3
etX = 1 + tx + + + ···
2! 3!
t2 t3
E[etX ] = 1 + tE[X] + E[X 2 ] + E[X 3 ] + · · ·
2! 3!
Differentiate once:
∂M (t) 2t
= 0 + E[X] + E[X 2 ] + · · ·
∂t 2!
0
M (0) = 0 + E[X] + 0 + 0 · · ·
Differentiate n times:
∂ n M (t) n × n − 1 × · · · 2 × t0 E[X n ]
= 0 + 0 + 0 + · · · +
∂nt n!
n!tE[X n+1 ] n!E[X n ] n!tE[X n+1 ]
+ + ··· = + + ···
(n + 1)! n! (n + 1)!
Suppose Z ∼ N(0, 1). Find its MGF and calculate the moments using it.
Using the definition: Z ∞
1 2
tX
E[e ] = √ etx e−x /2 dx
2π −∞
tx − 12 x2 = − 12 (x − t)2 − t2
5.2. PROPERTIES 47
Z ∞
1 t2 2 /2 t2
E[e tX
]= √ e2 e−(x−t) dx = e 2
2π −∞
Using the Lemma 5.2.1:
0 2 /2
M (0) = E[X] = et t|0 = 0
00 t2 /2
M (0) = E[X 2 ] = e (t2 + 1)|0 = 1
000 2 /2
M (0) = E[X 3 ] = et t(t2 + 3)|0 = 0
0000 2 /2
M (0) = E[X 4 ] = et (t4 + 6t2 + 3)|0 = 3
2 /2
M 5 (0) = E[X 5 ] = et t(t4 + 10t2 + 15)|0 = 0
2 /2
M 6 (0) = E[X 6 ] = et (t6 + 15t4 + 45t2 + 15)|0 = 15
Suppose X ∼ Poisson(λ), find its MGF and calculate the moments using it.
We know that if X ∼ Poisson(λ), then
λk e−λ
P(X = k) = for k = 0, 1, 2, 3, . . .
k!
and E[X] = λ, Var[X] = λ.
By definition:
∞
tX
X λk −λ
E[e ]= etk e
k!
k=0
∞
X (λet )k
= e−λ
k!
k=0
−λ λet
=e e
λ(et −1)
=e
t)
= e−λ(1−e
E[X] = λ
E[X 2 ] = λ2 + λ
Var[X] = E[X 2 ] − E[X]2 = λ
48 CHAPTER 5. MOMENT GENERATION FUNCTION
Suppose X ∼ Gamma(α, λ), find its MGF and calculate the moments using it.
By definition:
Z ∞
tX λα α−1 −λx
E[e ] = etk x e dx
0 Γ(α)
∞
λα
Z
= xα−1 ex(t−λ) dx
Γ(α) 0
It is easy to see that the previous integral converges for t < λ. Note that for a Gamma(α, λ − t)
it stands: Z ∞
(λ − t)α α−1 −(λ−t)x
x e dx = 1
0 Γ(α)
It follows that: Z ∞
Γ(α)
= xα−1 e(t−λ)x dx
(λ − t)α 0
and hence: α
λα Γ(α)
λ
M (t) = =
Γ(α) (λ − t)α λ−t
If we differentiate M (t) and evaluate the derivatives in t = 0 we get:
0 α
M (0) = E[X] =
λ
00 α(α + 1)
M (0) = E[X 2 ] =
λ2
As a consequence we find that:
α
Var[X] = E[X 2 ] − E[X]2 =
λ2
Suppose X ∼ Exponential(λ), find its MGF and calculate the moments using it.
By definition: Z ∞ Z ∞
tk −λx
tX
E[e ] = e λe dx = λe−(λ−t)x dx
0 0
the latter integral converges for t < λ.
It follows that:
λ −(λ−t)x ∞
E[etX ] = − e
λ−t 0
hence:
λ
M (t) =
λ−t
5.2. PROPERTIES 49
and
0 1
M (0) = E[X] =
λ
00 2
M (0) = E[X 2 ] =
λ2
1
Var(X) = E[X 2 ] − E[X]2 =
λ2
Note that this result could have been indirectly obtained observing that an exponential distribution
of parameter λ is a special case of a Gamma(α, λ), with the parameter α = 1.
Theorem 5.2.2. Suppose X and Y are random variables whose MGFs exist and are equal in some
open interval containing 0: MX (t) = MY (t) for all t ∈ (−δ, δ), δ > 0. Then X and Y have the same
distribution, i.e., FX (u) = FY (u) for all u ∈ R.
Then
N
Y
MY (t) = MXi (t)
i=1
50 CHAPTER 5. MOMENT GENERATION FUNCTION
Proof
Proof.
PN
MY (t) = E[etY ] = E[et i=1 Xi
]
tX1 +tX2 +···tXN
= E[e ] = E[etX1 ]E[etX2 ] · · · E[etXN ]
N
Y N
Y
tXi
= E[e ]= MXi (t),
i=1 i=1
Lemma 5.2.4. If X has MGF MX (t) and Y = a + bX, then Y has the MGF MY (t) = eat MX (bt).
Proof
Proof. By definition of MGF we have:
and so:
n n
E[Sn ] = , Var[Sn ] = .
λ λ2
The moment generating function of Xi is
λ
MXi (t) =
λ−t
so that the moment generating function of Sn = X1 + · · · + Xn is
n
λ
MSn (t) = .
λ−t
We recognize that this is the MGF of a Gamma random variable and hence Sn ∼ Gamma(n, λ).
X = X1 + X2 + · · · + Xn ,
Thus:
n
Y
MX (t) = MP Xi (t) = MXi (t) = (pet + 1 − p)n (since the Xi ’s are i.i.d.)
i=1
If X is a discrete random vector with joint probability mass function (pmf) p(x1 , x2 , . . . , xn ):
h > i X X
MX (t) = E e t X = ··· e t1 x1 +···+tn xn p(x1 , . . . , xn )
x1 xn
52 CHAPTER 5. MOMENT GENERATION FUNCTION
∂2M
X (t)
For example, E[Xi Xj ] = ∂ti ∂t j t=0
.
2. Independence: The components X1 , X2 , . . . , Xn of the random vector are independent if and
only if the joint MGF (assuming it exists in neighborhood of zero) factors into the product of the
marginal MGFs:
MX (t) = MX1 (t1 ) · MX2 (t2 ) · · · MXn (tn )
This is a very powerful result for proving independence.
3. Distribution: The MGF uniquely determines the joint distribution of the random vector (when
it exists in neighborhood of zero). If two random vectors have the same MGF, they have the same joint
distribution.
4. Transformations: The MGF is particularly useful for finding the distribution of linear transfor-
mations. If Y = AX + b for some matrix A and vector b, then the MGF of Y is:
>b
MY (t) = e t MX (A> t)
Finally, we will finish the section by demonstrating how the properties above help us to prove the
following lemma.
Lemma 5.3.1. If a n-dimensional random vector X ∼ N (µ, Σ) then for a linear transformation
Y = AX + b (where A is an m × n matrix of fixed constants and b is an m × 1 vector of constants) Y
comes from the multivariate normal distribution N (Aµ + b, AΣAT ).
Proof
Proof. We are given that X ∼ N (µ, Σ). The Moment-Generating Function (MGF) of X is
defined as:
h
>
i
> 1 >
MX (t) = E exp t X = exp t µ + t Σt
2
n
for t ∈ R . This is the defining property of the multivariate normal distribution.
Consider the linear transformation Y = AX + b, where A is an m × n matrix and b is an m × 1
vector. The MGF of Y is, by definition:
h i
MY (s) = E exp s> Y
for s ∈ Rm .
Substitute the expression for Y :
h i h i
MY (s) = E exp s> (AX + b) = E exp s> AX + s> b
Using the property of exponentials, exp(a + b) = exp(a) exp(b), and noting that s> b is a scalar
constant with respect to the expectation, we can factor it out:
h i
MY (s) = exp s> b · E exp (A> s)> X
5.3. EXTENSION FOR A RANDOM VECTOR 53
Notice that the expectation E exp (A> s)> X is precisely the MGF of X evaluated at the
vector t = A> s: h i
E exp (A> s)> X = MX (A> s)
This is the MGF of a multivariate normal distribution with mean vector µY = Aµ + b and
covariance matrix ΣY = AΣA> .
Therefore, by the uniqueness of moment-generating functions, we conclude:
Y = AX + b ∼ N (Aµ + b, AΣA> )
σ12
>
ρσ1 σ2 t1
t Σt = t1 t2 .
ρσ1 σ2 σ22 t2
Multiplying out:
t> Σt = σ12 t21 + σ22 t22 + 2ρσ1 σ2 t1 t2 .
Substitute into the general formula:
1 2 2 2 2
MX1 ,X2 (t1 , t2 ) = exp µ1 t1 + µ2 t2 + σ t + σ2 t2 + 2ρσ1 σ2 t1 t2 .
2 1 1
54 CHAPTER 5. MOMENT GENERATION FUNCTION
Chapter 6
Limit Theorems
In this chapter we will mainly deal with the behavior of the sum of random variables, as the number of
summands becomes large. Many statistical methods rely on large samples. Hence, the properties are
very important to know.
E[X]
P(X ≥ a) ≤
a
Proof
Proof. For a > 0,
Z ∞
E[X] = xf (x)dx
0
Z a Z ∞
= xf (x)dx + xf (x)dx
0 a
Because X ≥ 0,
Z ∞ Z ∞
E[X] ≥ xf (x)dx ≥ af (x)dx = aP(X ≥ a)
a a
E[X]
≥ P(X ≥ a)
a
Suppose a company models the time (in hours) a server stays online before needing a reboot as a
non-negative random variable X, with an expected uptime of:
We want to estimate the probability that the server stays online for at least 500 hours.
Applying Markov’s Inequality
55
56 CHAPTER 6. LIMIT THEOREMS
σ2
P(|X − µ| ≥ k) ≤
k2
Proof
Proof. Define the non-negative r.v. random variable Y = (X − µ)2 . Set a = k 2 and apply
Markov’s inequality:
E[Y ]
P(Y ≥ k 2 ) ≤
k2
E[(X − µ)2 ]
P((X − µ)2 ≥ k 2 ) ≤
k2
σ 2
P((X − µ)2 ≥ k 2 ) ≤ 2
k
Further we know that, if (X − µ)2 ≥ k 2 , we get |X − µ| ≥ k. Thus:
σ2
P(|X − µ| ≥ k) ≤
k2
Pafnuti Chebyshev (1821-1894) was a Russian mathematician. Many variants of the spelling of
his name exist, Pafnut-i/y (T/Ts)-cheb-i/y-sh/ch/sch-e-v/ff/w.
6.2. CHEBYSHEV’S INEQUALITY 57
Suppose the time it takes to complete a data processing task is modeled by a random variable X,
with:
• Mean µ = 10 minutes
We want to estimate the probability that the task takes more than 16 minutes (i.e., 6 minutes
above the mean).
Applying Chebyshev’s Inequality
Chebyshev’s inequality states:
σ2
P(|X − µ| ≥ k) ≤ 2
k
Let k = 6. Then:
22 4
P(|X − 10| ≥ 6) ≤ 2 = = 0.111
6 36
Interpretation
Even without knowing the exact distribution of X, Chebyshev’s inequality tells us:
“The probability that the task takes less than 4 minutes or more than 16 minutes is at
most 11.1%.”
Note
Why Chebyshev’s inequality?
It can be useful for setting performance guarantees or bounding risk when only the
mean and variance are known.
Suppose X ∼ N(µ, σ 2 ). Calculate P(|X − µ| ≥ kσ) for a number of different values of k and
compare the exact results with the bound given by Chebyshev’s inequality.
We can write that
P(|X − µ| ≥ kσ) = P(X − µ ≥ kσ) + P(X − µ ≤ −kσ)
= P(X ≥ µ + kσ) + P(X ≤ µ − kσ).
But for the normal distribution we evaluate cumulative probabilities by standardising. Thus
X −µ X −µ
P(|X − µ| ≥ kσ) = P ≥k +P ≤ −k
σ σ
= P(Z ≥ k) + P(Z ≤ −k)
58 CHAPTER 6. LIMIT THEOREMS
X −µ
where Z = ∼ N(0, 1) with Φ(z) = P(Z ≤ z). Thus
σ
P(|X − µ| ≥ kσ) = (1 − P(Z ≤ k)) + P(Z ≤ −k)
= (1 − Φ(k)) + Φ(−k)
= (1 − Φ(k)) + (1 − Φ(k)) = 2(1 − Φ(k)).
Alternatively it is 2Φ(−k).
k exact prob chebyshev bound ratio
0.25 8.026e-01 1.00000 1.246
0.50 6.171e-01 1.00000 1.621
0.75 4.533e-01 1.00000 2.206
1.00 3.173e-01 1.00000 3.151
1.50 1.336e-01 0.44444 3.326
2.00 4.550e-02 0.25000 5.494
2.50 1.242e-02 0.16000 12.883
3.00 2.700e-03 0.11111 41.155
3.50 4.653e-04 0.08163 175.457
4.00 6.334e-05 0.06250 986.700
4.50 6.795e-06 0.04938 7267.138
5.00 5.733e-07 0.04000 69771.116
x
2k k 2k
Mn = max(X1 , X2 , . . . , Xn )
Consider the scaled random variable:
Yn = n(1 − Mn )
60 CHAPTER 6. LIMIT THEOREMS
x x
P (Yn ≤ x) = P (n(1 − Mn ) ≤ x) = P 1 − Mn ≤ = P Mn ≥ 1 −
n n
x x n
= 1 − P Mn ≤ 1 − =1− 1− .
n n
Using the well-known limit:
a n
lim 1+ = ea
n→∞ n
with a = −x, we obtain:
x n
lim P (Yn ≤ x) = 1 − lim 1− = 1 − e−x
n→∞ n→∞ n
The function F (x) = 1 − e−x for x ≥ 0 is precisely the CDF of an Exponential(1) random variable.
Therefore:
d
Yn −
→ Exponential(1)
Lemma 6.4.1. Suppose X1 , X2 , . . . , Xn is a random sample from a distribution with mean µ and
Var[Xi ] = σ 2 . Then, for all ε > 0,
X1 + X2 + . . . + Xn
P − µ ≥ ε → 0 as n → ∞
n
Proof
σ2
Proof. We will prove that P X −µ ≥ε ≤ .
It is easy to show that:
nε2
Pn
E[X1 + X2 + · · · + Xn ] E[Xi ]
= i=1 =µ
n n
Further,
" P 2 #
n
i=1 Xi −µ Var[(X1 + X2 + · · · + Xn )
E =
n n2
Pn
Var[(Xi ) σ2
= i=1 2 =
n n
Apply Chebyshev’s Inequality:
σ2
X1 + X2 + . . . + Xn
P −µ ≥ε ≤ 2
n nε
The Weak Law of Large Numbers (WLLN) tells us that the sample average of independent and
identically distributed (i.i.d.) random variables converges in probability to their expected value.
Let X1 , X2 , . . . be i.i.d. random variables with Xi ∼ N(0, 10). Suppose we want to guarantee
that we have at most a 0.01 probability of the sample mean being more than 0.1 away from the
true mean µ = 0.
Chebyshev’s inequality provides a general bound that works for any distribution with finite
variance:
σ2
P Xn − µ ≥ ε ≤ 2
nε
62 CHAPTER 6. LIMIT THEOREMS
σ2
= 0.01
nε2
10
= 0.01
n(0.1)2
10
= 0.01
0.01n
1000
= 0.01
n
1000
n= = 100, 000
0.01
So, using Chebyshev’s inequality, we find that n = 100, 000 guarantees our requirement.
Since we know the Xi are normally distributed, we can use the exact sampling distribution:
10
X n ∼ N 0,
n
We want:
P |X n | ≥ 0.1 ≤ 0.01
Standardizing: ! r
0.1 n
P |X n | ≥ 0.1 = 2Φ − p = 2Φ −0.1
10/n 10
Setting this equal to 0.01:
r
n
2Φ −0.1 = 0.01
10
r
n
Φ −0.1 = 0.005
10
r
n
−0.1 = z0.005 ≈ −2.576
10
r
n
= 25.76
10
n = 10 × (25.76)2 ≈ 6, 635
So, using the exact normal distribution, we only need n ≈ 6, 635. Compare this with 100, 000
bound!
6.4. WEAK LAW OF LARGE NUMBERS 63
Recall, that the WLLN states that for independent and identically distributed random variables
X1 , X2 , . . . with finite mean µ and finite variance σ 2 , the sample mean X n converges in probability
to µ. Formally, for any ε > 0:
The plots below (upper panel) show how mean of samples taken from the Exponential(1) and
Uniform(0, 1) converges to the true means, µexp = 1 and µuni = 0.5 as the sample size increases.
The plots on the lower panel show the values of the P (|X n − µ| ≥ ε) converging to 0 as the
sample size increases.
Observe, how convergence happens quicker for the uniform distribution. Indeed, similar to the
example above, calculations give us (using conservative Chebyshev’s bound, assuming ε = 0.01
σ2
and nε 2 = 0.01)
σ2 1
Exponential: µ = 1, σ 2 = 1 =⇒ 2
= ≤ 0.01
nε n(0.01)2
=⇒ n ≥ 1, 000, 000,
1 σ 2 1/12
Uniform: µ = 0.5, σ2 = =⇒ ≤ 0.01
12 nε2 n(0.01)2
=⇒ n ≥ 83, 333.
Z x 2
Sn − nµ 1 t
lim P √ ≤ x = Φ(x) = √ exp − dt
n→∞ σ n 2π −∞ 2
Proof
Proof. Proof plan: We will follow the following plan:
1. Use the fact that convergence of MGFs implies convergence in distribution (MGF con-
vergence theorem).
2. Show that the MGF of the normalized sum converges to the MGF of the standard normal
distribution.
Note
MGF convergence theorem.
Let Fn be a sequence of cumulative distribution functions with corresponding moment
generating functions Mn (t), and let F be a CDF with MGF M (t). If lim Mn (t) = M (t)
n→∞
for all t in some neighborhood of 0, then Fn (x) F (x) for all x where F is continuous.
Note
The limit exponential theorem.
If lim an = a, then
n→∞ an n
lim 1+ = ea
n→∞ n
Note
Property of MGF.
If M (t) is the moment generating function of a random variable X, then M (0) = 1.
Let X1 , . . . , Xn be i.i.d. random variables with E[X] = 0, variance σ 2 , and MGF MX (t).
Define:
• Sn = ni=1 Xi
P
S√n
• Zn = σ n
Then: n
n t
MSn (t) = (MX (t)) , MZn (t) = MX √
σ n
Using Taylor expansion around t = 0:
0 1 00
MX (s) = MX (0) + sMX (0) + s2 MX (0) + εs
2
where εs /s2 → 0 as s → 0.
Since E[X] = 0, Var(X) = σ 2 , and using property of MGF, we have:
σ2 2
MX (s) = 1 + 0 + s + εs
2
t
Let s = √
σ n
, then as n → ∞, s → 0.
2 !n n
σ2 t2
t 1 1
MZn (t) = 1+ √ +o = 1+ +o
2 σ n n 2n n
t2
Recall the limit
n exponential theorem. Substituting an = 2 into above and using
limn→∞ 1 + ann = ea , we get
2
lim MZn (t) = et /2
n→∞
This is the MGF of the standard normal distribution N (0, 1), so by the MGF convergence
theorem, Zn N (0, 1), that is convergence in distribution. This proves the CLT.
66 CHAPTER 6. LIMIT THEOREMS
6.5.2 Interpretation
The Central Limit Theorem (CLT) deals with the long-run behaviour of the sample mean as n grows.
In general, the colloquial result is that everything becomes normal eventually.
Consider i.i.d. random variables X1 , X2 , . . . , Xn each with mean µ and variance σ 2 .
Note (Important)
It doesn’t necessarily have to be a Normal distribution! For example, if we had that each X was
distributed uniform(0, 1), then µ = 1/2 and σ 2 = 1/12.
Again, define X n as the sample mean of X. We can write this out as:
X1 + . . . + Xn
Xn =
n
As discussed above, we know that the sample mean is itself a random variable, and we know that it
approaches the true mean in the long run by the Law of Large Numbers (LLN). However, are we able
to nail down a specific distribution for this random variable as we approach the long run, not just the
value that it converges to? A good place to start is to find the mean and variance; these parameters
won’t tell us what the distribution is, but they will be useful once we determine the distribution. To
find the expectation and variance, we can just brute force our calculations. First, for the expectation,
we take the expectation of both sides:
X1 + . . . + Xn
E[X n ] = E
n
By linearity:
X1 Xn
=E + ... + E
n n
Since n is a known constant, we can factor it out of the expectation:
1 1
= E[X1 ] + . . . + E[Xn ]
n n
Now, we are left with the expectation, or mean, of each X. Do we know these values? Well, recall above
that, by the set-up of the problem, each X is a random variable with mean µ. That is, the expectation
of each X is µ. We get:
µ µ
= + ... +
n n
We have n of these terms, so they sum to:
=µ
Hence, we get that E[X n ] = µ. Think about this result: it says that the average of the sample mean
is equal to µ, where µ is the average of each of the random variables that make up the sample mean.
This is intuitive; the sample mean should have an average of µ (you could say that X n is unbiased for
µ, since it has expectation µ; this is a concept that you will explore more in detail). Let’s now turn to
the Variance:
X1 + ... + Xn
Var[X n ] = Var .
n
We know that the X terms are independent, so the variance of the sum is the sum of the variances:
X1 Xn
= Var + ... + Var .
n n
6.5. CENTRAL LIMIT THEOREM 67
σ2
=
n
Consider this result for a moment. First, consider when n = 1. In this case, the sample mean is just
2
X1 (since, by definition, we would have X11 = X1 ). The variance that we calculated above, σn , comes
out to σ 2 when n = 1, which makes sense, since this is just the variance of X1 . Next, consider what
happens to this variance as n grows; it gets smaller and smaller, since n is in the denominator. Does
this make sense? As n grows, we are essentially adding up more and more random variables in our
sample mean calculation. It makes sense, then, that the overall sample mean will have less variance;
among other things, adding up more random variables means that the effect of outlier random variables
is lessened (i.e., if we observe an extremely large value for X1 , it is mediated by the sheer number of
random variables).
2
So, we found that the sample mean has mean µ and variance σn , where µ is the mean of each
underlying random variable, σ 2 is the variance of each underlying random variable, and n is the total
number of random variables. Now that we have the parameters, we are ready for the main result of the
CLT.
Note (Important)
The CLT states that, for large n, the distribution of the sample mean approaches a Normal
distribution.
This is an extremely powerful result, because it holds no matter what the distribution of the
underlying random variables (i.e., the X’s) is.
We know that Normal random variables are governed by the mean and variance (i.e., these are the
two parameters), and we already found the mean and variance of X n , so we can say:
d σ2
X n → N (µ, )
n
d
Where → means converges in distribution; it’s implied here that this convergence takes place as n,
or the number of underlying random variables, grows.
Think about this distribution as n gets extremely large. The mean, µ, will be unaffected, but the
variance will be close to 0, so the distribution will essentially be a constant (specifically, the constant µ
with no variance). This makes sense: if we take an extremely high number of draws from a distribution,
we should get that this sample mean is at the true mean, with very little variance. It’s also the result
we saw from the LLN, which said that the sample mean approaches a constant: as n grows here, we
approach a variance of 0, which essentially means we have a constant (since constants have variance 0).
The CLT just describes the distribution on the way to the LLN convergence.
Hopefully this brings some clarity to the statement ”everything becomes Normal”: taking the sum
of i.i.d. random variables (we worked with the sample mean here, but the sample mean is just the sum
68 CHAPTER 6. LIMIT THEOREMS
divided by a constant n), regardless of the underlying distribution of the random variables, yields a
Normal distribution.
Let us highly this once more.
Note
The distribution of an average tends to be Normal, even when the distribution from which the
average is computed is decidedly non-Normal (except when the moments of the parent
distribution do not exist)!
Note
This normal distribution will have the same mean as the parent distribution, and, variance equal
to the variance of the parent divided by the sample size.
Note
The distribution of the phenomenon under study does not have to be Normal, but its average
will be.
6.5.3 Examples
µ = 1/λ
σ 2 = 1/λ2
λ
MX (t) =
λ−t
and also: n
λ
MP Xi (t) =
i λ−t
and P
Xi
− λ1 Xi √
i
P
n λ X
Zn = 1 = √i − n=a Xi + b.
√
λ n
n
i
It follows that:
√
n
λ
MZn (t) = E[etZn ] = ebt MPi Xi (at) = e− nt
√
λ − λt/ n
√ t −n
= e− nt
1− √ .
n
6.5. CENTRAL LIMIT THEOREM 69
Let’s consider the log and remember that: log(1 + x) ' x − 21 x2 + 13 x3 + ... we get
√
t
log(MZn (t)) = − nt − n log 1 − √
n
√ t2 t3
t
' − nt − n − √ − − √
n 2n 3n n
1 t3
' t2 + √
2 3 n
1 2
It is easy to see that as n → ∞, then log(MZn (t)) → 12 t2 and (MZn (t)) → e 2 t hence:
d
Zn → N(0, 1)
We use an exponential distribution with rate parameter λ = 2. In this case, the population
mean is µ = λ1 = 0.5 and the population variance is σ 2 = λ12 = 0.25.
Central Limit Theorem says that sample means X n = n1 ni=1 Xi converge to the normal
P
distribution with the theoretical mean E[X n ] = µ = 0.5 and theoretical variance Var(X n ) =
σ2 0.25
n = 100 = 0.0025 as n grows.
The simulations below compare the theoretical results with the simulations for sample size n = 100
and 2000 sample means.
Observe, how the distribution of means is approaching to the expected normal distribution shape
with the number of samples increasing.
70 CHAPTER 6. LIMIT THEOREMS
6.5. CENTRAL LIMIT THEOREM 71
72 CHAPTER 6. LIMIT THEOREMS
Exercise Sheet
Note
This exercise sheet covers material from Sections 4, 5 and 6 of the course. Please
attempt all questions and prepare your solutions before the Week 5 workshop.
Your solutions will not be assessed. Before the workshop, you will be invited to vote for the
questions you would like to see discussed. The tutor will select the most popular questions for
discussion during the workshop. Complete solutions for all questions will be available afterwards
(by the end of the week). The solutions will be available only in HTML version of the notes.
Exercise 6.1
Suppose X has a bivariate normal distribution satisfying X ∼ N (µ, Σ) where
X1 2 4 5
X= , µ= , Σ= .
X2 3 5 8
What are the values of E[X1 ], E[X2 ], Var[X1 ], Var[X2 ], Cov(X1 , X2 ), and Cor(X1 , X2 )?
Exercise 6.2
Prove Lemma 4.2.1.
Hint: use the expression for the bivariate normal distribution (see Example 4.2)).
Exercise 6.3
If X1 and X2 are independent N (µ, σ 2 ) variables, find the joint distribution of U1 and U2 where
U1 1 1 X1
= .
U2 1 −1 X2
Exercise 6.4
If X ∼ Exponential(λ), use the moment generating function (MGF) to obtain the distribu-
tion of U = aX, where a > 0 is a known constant. Using R functions plot() and rexp() plot
the empirical and theoretical MGF for U (choose a = 2 and λ = 1).
73
74 CHAPTER 6. LIMIT THEOREMS
λ
MX (t) = E[etX ] = , t < λ.
λ−t
Exercise 6.5
If X and Y are independent Exponential(λ) random variables, use the moment generating
function (MGF) to obtain the distribution of Z = X + Y . Use R to illustrate, by comparing
the distributions coming from Z and from X + Y (use λ = 2). Use R functions rgamma(), rexp()
and hist() for the illustration.
Reminder: The MGF of an exponential random variable X ∼ Exponential(λ) is
λ
MX (t) = E[etX ] = , t < λ.
λ−t
Exercise 6.6
Recall the Note 4.2: show that if the marginal distributions of X1 and X2 are normal and
Cor(X1 , X2 ) = 0, the variables are not necessarily independent.
Use the following counterexample and moment generating functions for the Normal Dis-
tribution
1 − 1 (x21 +x22 ) 4 2 2 − 1 (x21 +x22 )
fX1 ,X2 (x1 , x2 ) = e 2 1 + x1 x2 1 − x1 x2 e 2 .
2π 9
to verify that:
X1 ∼ N (0, 1), X2 ∼ N (0, 1), Cor(X1 , X2 ) = 0.
while the variables X1 and X2 are not independent.
Exercise 6.7
Write R code to illustrate the CLT for Xi ∼ Beta(0.2, 2). Use the R function rbeta() to sample
from the distribution.
Exercise 6.8
Credit to Gary Ernest Davis, ”EXPLORATORY DATA ANALYSIS with R”.
The dataset Danish Fire Insurance Claims 734 59. These data describe large fire insurance claims
in Denmark from Thursday 3rd January 1980 until Monday 31st December 1990. The data are
contained in a numeric vector.
1. Explore the dataset: Load the Danish Fire Insurance Claims data, compute key statistics
(mean, variance, skewness, kurtosis), and visualize it with a histogram and boxplot. Summa-
rize what these measures reveal about the distribution. You can use the R package e1071
to calculate skewness and kurtosis and R functions hist() and boxplot() for boxplots
and histograms.
2. Apply Chebyshev’s inequality P (|X − µ| ≥ kσ) ≤ k12 to estimate the proportion of
data within kσ of the mean. For k = 4.5 calculate the bound, find the corresponding range,
and compare it to the actual data. Discuss whether the bound is tight and what this tells
you about the dataset.
Chapter 7
Estimation Methods
Definition 7.1.1. An estimator θbn of θ is a function of the observed random variables which (we
hope) forms a useful approximation to the parameter, once the data are observed:
θbn = g(X1 , X2 , . . . , Xn ).
Note that θbn can depend only on the observed random variables, and not on any unknown parameters.
Note (Important)
The estimator is a function of random variables, so is itself a random variable, with a distribution,
mean, and variance, etc.
Suppose we have observed some particular values for the random variables, so X = (X1 =
x1 , X2 = x2 , . . . , Xn = xn ). Then θ̂(x) becomes a fixed value and is the estimate of θ for the data
x = (x1 , x2 , . . . , xn ).
Recall, that there are several methods for obtaining point estimates:
1. Method of Moments: This technique involves equating sample moments to theoretical moments
of the distribution and solving for the parameters.
2. Maximum Likelihood Estimation (MLE): This method finds the parameter values that
maximize the likelihood function, i.e., the probability of observing the given data under the
assumed model.
3. Bayesian Estimation: In this approach, parameters are treated as random variables with prior
distributions. The estimate is derived from the posterior distribution, which combines the prior
information with the observed data using Bayes’ theorem.
75
76 CHAPTER 7. ESTIMATION METHODS
Estimators can be compared using several key properties that reflect their quality and reliability.
Bias measures the difference between the expected value of the estimator and the true parameter; an
unbiased estimator has zero bias. Variance captures the spread of the estimator’s values across different
samples, with lower variance indicating more stable estimates (efficiency). The mean squared error
(MSE) combines both bias and variance, and is often used as a comprehensive criterion. Consistency
refers to an estimator’s ability to converge to the true parameter value as the sample size increases.
Asymptotic unbiasedness refers to estimators whose bias approaches zero as the sample size grows.
Lastly, invariance means that if an estimator is applied to a transformed parameter (e.g., a function of
the original), the resulting estimate should correspond to the transformation of the original estimator—
this property is especially important in maximum likelihood estimation. These criteria help in selecting
the most appropriate estimator for a given statistical task.
Some of the properties you have met before, so they are suggested as part of the revision.
7.2.1 Examples
X ∼ Gamma(α, λ)
Now, instead of the central second moment n1 (Xi − X)2 we can use raw moment: 1
Xi2 =
P P
2 n
E[X 2 ] = Var[X] + E[X]2 = λα2 + αλ = α(α+1))
λ2
.
We obtain:
α α(α + 1))
X= , m02 =
λ λ2
2
X X
=⇒ α
b= 2 , λb=
0
.
m02 − X m2 − X
Suppose X1 , . . . , Xn form a random sample from a geometric(θ) distribution. Obtain the method
of moments estimator for θ (θbM M ).
P (X = x) = θ(1 − θ)x , x = 0, 1, 2, 3, . . . =⇒
X∞ ∞
X
µ0 ≡ µ = xθ(1 − θx ) = θ x(1 − θ)x
x=0 x=0
∞
" ∞ #!
X
x−1 d X
x
= θ(1 − θ) x(1 − θ) = θ(1 − θ) − (1 − θ)
dθ
x=0 x=0
d 1 1−θ
= θ(1 − θ) − = .
dθ θ θ
∞
1 X
We used that Taylor’s expansion at 0: = xn
1−x
n=0
78 CHAPTER 7. ESTIMATION METHODS
We therefore get
m01 ≡ µ01 , X ≡ (1 − θ)/θ =⇒ θb = 1/(1 + X).
The value of θ which maximises the likelihood L(θ; x) is called the maximum likelihood estimate
of θ. In the lectures we often will use the notation θ̂M LE as opposed to θ̂M M for the estimation done by
the Method of Moments.
Let X ∼ Poisson(λ) and suppose we have observed x = 1. What is the maximum likelihood
estimate for λ?
For one observation we have
λx e−λ
L(λ; x) = ,
x!
l(λ; x) = log(λx e−λ ) − log(x!) = x log(λ) − λ + c,
∂l(λ; x) x
= − 1 = 0 =⇒ λ̂ = x = 1,
∂λ λ
∂ 2 l(λ; x) x
2
= − 2 < 0,
∂λ λ
Now, let Xi ∼ Poisson(λ) and suppose we have observed x = {1, 5, 3}. What is the maximum
likelihood estimate and the method of moments estimate for λ?
7.4. BAYESIAN ESTIMATION 79
3
Y λxi e−λ
L(λ; x) = ,
xi !
i=1
3
X X
l(λ; x) = (xi log(λ) − λ − log(xi !)) = xi log(λ) − 3λ + c,
i=1 i
P
∂l(λ, x) xi
= i − 3 = 0 =⇒ λ̂ = x̄ = 3,
∂λ λ
∂ 2 l(λ; x)
P
xi
= − i2 < 0,
∂λ2 λ
Prior
Suppose we want to learn about k unknown parameters θ = (θ1 , θ2 , . . . , θk ). If there is only one
parameter, we simplify the notation to θ.
We express our uncertainty about θ using a prior density π(θ). This prior reflects our beliefs before
observing any data. If we were to guess a value for θ, we might report the mean of the prior. Note that
the prior is subjective-your belief about θ may differ from mine.
Likelihood
Suppose we have observed n data points x = (x1 , x2 , . . . , xn ). For a single observation, we use the
notation x.
The likelihood function tells us how plausible different values of θ are, given the observed data x. It
encodes our assumptions about the data-generating process and links the data to the parameters.
80 CHAPTER 7. ESTIMATION METHODS
Importantly, the likelihood is a function of θ for fixed data x, and not a probability density function.
Common notations include L(θ; x) and p(x | θ).
Sometimes, data is denoted by y instead of x.
Predictive Distribution
The marginal distribution of the data is given by:
Z Z
π(x) = π(x, θ) dθ = π(x | θ)π(θ) dθ.
Θ Θ
Posterior
Once data x is observed, we update our beliefs using Bayes’ theorem:
Here, π(θ | x) is the posterior density, representing our updated beliefs about θ after seeing the
data.
Example (Binomial Distribution)
Consider the experiment with the coin. You have no idea if the coin is fair. Therefore, you can
choose the uniform prior distribution for the probability of observing heads. You tossed it n
times and observed y heads. Using Bayes theorem, you would like to assess if the coin is fair,
that is, to obtain the posterior distribution for the probability of observing heads θ.
Note (Assumptions)
We assume that all tosses are independent and have only two outcomes.
Note
We have used that the PDF for the Binomial distribution is
Recall, that for the Uniform(a, b), the expected value is E[θ] = 12 (a + b) and the variance is
1
Var[θ] = (b − a)2 12 . Therefore, in this example, the prior mean is 12 , and the prior variance is 12
1
.
α
Recall, that for the Beta distribution E[θ] = α+β .
y+1
Since in our model posterior follows Beta(α = y + 1, β = n − y + 1), the posterior mean is n+2 .
y
It is a compromise between the prior mean and the sample proportion, n , where clearly the
prior mean has a smaller and smaller role as the size of the data sample increases. This is a
very general feature of Bayesian inference: the posterior distribution is centered at a point that
represents a compromise between the prior information and the data, and the compromise is
controlled to a greater extent by the data as the sample size increases.
Obviously, the solution can be updated for the general case of prior θ ∼ Beta(α, β):
Below is the plot showing the prior and posterior distributions for the case of observing y = 14
successes out of n = 20 trials. I have used prior Beta(0.5, 0.5).
82 CHAPTER 7. ESTIMATION METHODS
The plot below shows how the sample size influence the posterior. Indeed, with small sample
sizes, the posterior is heavily influenced by the prior. As sample size increases, the posterior
converges toward the sample proportion.
7.4. BAYESIAN ESTIMATION 83
The property that the posterior distribution follows the same parametric form as the prior distribution
is called conjugacy. We say that the beta prior distribution is a conjugate family for the binomial
likelihood. The main reason to choose a conjugate family is mathematical convenience. Indeed this will
ensure that the posterior distribution follows a known parametric form.
Note (Important)
Beta distribution is a conjugate prior for binomial model. Indeed, we showed above that if
an observation y ∼ Binomial(θ, n), and the prior for θ is θ ∼ Beta(α, β), the posterior is
θ | y ∼ Beta(α + y, β + n − y).
In practice it could happen that the conjugate prior is not a good fit to the specific example we are
analysing, it may be necessary to use a more realistic prior distribution.
We can use the produced posterior distribution π(θ | x) to obtain a point estimate for the unknown
parameter θ. Expected value is one of the methods to obtain a point estimate.
We also call it Bayes estimate of θ and use notation θ̂Bayes .
More generally, the Bayes estimate of the unknown parameter θ can be obtained as
Z
θ̂Bayes = E[θ|x] = θπ(θ|x) dθ.
Θ
In the example above we have obtained posterior for Beta-Binomial model θ | y ∼ Beta(α +
y, β + n − y).
α
Using that E[X] for X ∼ Beta(α, β) is E[X] = α+β , we can write the Bayes estimate for θ | y:
α0 α+y
θbBayes = E[θ | y] = 0 0
= .
α +β α+β+n
Let us consider an example of Gamma-Poisson Bayesian model applied to the real life dataset
[Link] It can be shown that Gamma distribution is a conjugate for Poisson,
therefore the posterior is also Gamma.
For the two decades after 1st January 1900 there were x1 = 4 and x2 = 3 mining disasters during
1910s and 1920s respectively. The details of the disasters are summarized in the table below.
2. Prior. Prior information: number of disasters in five decades before 1900 were 4, 11, 2, 11,
11. As a prior distribution for θ we choose Gamma(α, λ) = Gamma(3.86, 0.495). Here the
exact values α and λ were estimated using the method of moments. There are other
ways to choose α and λ for the prior. One method (equivalence sample size, not covered
in this module) would suggest λ = 5 (number Rof prior observations) and α = 39 (sum of
∞
the data). A vague improper prior (integral 0 π(θ) dθ doesn’t converge) would suggest
λ = 0 and α = 1 (corresponds to the so called ”flat” prior, a horizontal line). In choosing
λ = 0.495 and α = 3.86 we are effectively treating the observed prior data as observed
values of θ (mean number of disasters per decade) and use a method of moments to fit in
the prior distribution.
Suppose n observations are available from a geometric(θ) distribution. Use a beta distribution
prior to obtain the posterior distribution of θ.
Since X ∼ geometric(θ), P(X = x) = (1 − θ)x θ, x = 0, 1, 2, . . . so that for n observations
x = (x1 , x2 , . . . , xn ) the likelihood is
θa−1 (1 − θ)b−1
π(θ) = for 0 < θ < 1.
B(a, b)
7.5. UNBIASEDNESS, MSE, CONSISTENCY AND EFFICIENCY (REVISION) 85
π(θ|x) ∝ p(x|θ)π(θ)
∝ (1 − θ)nx̄ θn · θa−1 (1 − θ)b−1
∝ θa+n−1 (1 − θ)b+nx̄−1 .
θa+n−1 (1 − θ)b+nx̄−1
π(θ|x) = , 0 < θ < 1.
B(a + n, b + nx̄)
Definition 7.5.2 (Mean Squared Error). For a given estimator θ̂ = g(X), we define its Mean
Squared Error as
MSE[θ̂] = E[(θ̂ − θ)2 ],
which is a measure of accuracy.
Definition 7.5.5 (Efficiency). θbA is said to be more efficient than θbB if (given they are unbiased)
h i h i
Var θbA < Var θbB , ∀θ ∈ Θ.
86 CHAPTER 7. ESTIMATION METHODS
Bias(θ̂n ) = E(θ̂n ) − θ → 0, as n → ∞.
n
!
1X 2 1
E (Xi − X n ) =E
n n
i=1
n
!
1X
(Xi − µ)2 − 2(Xi − µ)(X − µ) + (X − µ)2
=E
n
i=1
n n
!
1 X 2 X
=E (Xi − µ)2 − (X − µ) (Xi − µ) + (X − µ)2 ,
n n
i=1 i=1
n n
1X 1X
notice: X −µ= Xi − µ = (
n n
i=1 i=1
n
1 X 2
=⇒ E (Xi − µ)2 − (X − µ
n n
i=1
n n
! !
1X 1X
(Xi − µ)2 − (X − µ)2 (Xi − µ)2 − E (X − µ)2 ,
=E =E
n n
i=1 i=1
2
σ 1
= σ 2 − E (X − µ)2 = σ 2 − σ2 < σ2 ,
= 1−
n n
2 2
recall that E (X − µ)2 = Var X = σn . We proved that the estimator σ̂ 2 = n1 ni=1 Xi − X
P
is biased.
However,
σ2
1
Bias = 1 − σ2 − σ2 = − → 0, as n → ∞.
n n
2
We showed that the estimator σ̂ 2 = n1 ni=1 Xi − X is asymptotically unbiased.
P
Proof
Proof. If g is one-to-one, then
L(θ; x) = L g −1 (g(θ)) ; x
If g is many-to-one, then θb which maximises L(θ; x) still corresponds to g θb , so g θb still
corresponds to the maximum of L(θ; x).
Then
∞ λ
X λk X λk
P(X > λ) = exp(−λ) =1− exp(−λ) ≡ g(λ),
k! k!
k=λ+1 k=0
λ̂ n
d = g(λ̂) = 1 −
X λ̂k 1X
g(λ) exp(−λ̂) , λ̂ = xi
k! n
k=0 i=1
The log-likelihood is
!
X X
l(θ; x) = xi log θ + n− xi log(1 − θ)
i i
and P P
∂l(θ; x) i xi n − i xi
= −
∂θ θ 1−θ
so it is easily shown that the MLE of θ is θbM LE = X.
Using Lemma 7.7.1
b − θ)
gb = g θb = θ(1 b =X 1−X .
∂ 2 l(θ; x)
J(θ) = −
∂θ2
l(θ; x) ≈ l(θ;
b x) + (θ − θ) b x) + 1 (θ − θ)
b · l0 (θ; b 2 · l00 (θ;
b x)
2
Since θb maximizes the log-likelihood, the first derivative vanishes:
l0 (θ;
b x) = 0
b x) − 1 (θ − θ)
l(θ; x) ≈ l(θ; b 2 J(θ)
b
2
where J(θ)b = −l00 (θ;
b x) is the observed information.
This quadratic form resembles a parabola centered at θ.
b The curvature J(θ)
b tells us how concentrated
the likelihood is:
7.8. FISHER’S INFORMATION 89
• Large J(θ)
b → steep parabola → high certainty
• Small J(θ)
b → flat parabola → low certainty.
Definition 7.8.2 (Fisher Information). Let X be a random variable with probability density
function (PDF) fX (x; θ) (or probability mass function (PMF)). Let l(θ; X) is the likelihood, where X is
a random sample before the data has been observed. Then the Fisher information In (θ) is defined as
2
∂ l(θ; X)
In (θ) = E − .
∂θ2
In (θ) = ni(θ)
where
∂2
i(θ) = E − 2 log f (Xi ; θ) ,
∂θ
the single observation Fisher information.
It can be shown that:
" #
∂l(θ; X) 2
2
∂ l(θ; X)
In (θ) = E − =E
∂θ2 ∂θ
b ' 1
Var(θ)
In (θ)
90 CHAPTER 7. ESTIMATION METHODS
Note (Important)
Remember that the precision to which we can estimate θ is limited by the Fisher infor-
mation of the likelihood function.
7.8.2 Examples
Univariate Distribution
giving log-likelihood
X X
l(θ; X) = log L(θ; X) = Xi log θ − nθ − log Xi !
i i
Thus
∂ 2 l(θ; X)
P
i Xi
J(θ; X) = − = .
∂θ2 θ2
To find I(θ), we need E (Xi ) = θ and
1 X n
I(θ) = 2
E (Xi ) = .
θ θ
i
Notice, above we used: E(Xi ) = θ (the proof can be obtained using MGF MX (t) = exp(θ(et − 1))
using that MX0 (0) = E[X]).
Finally,
b ' θ.
Var(θ)
n
Optional
Multivariate Case
For simplicity of notation in this section we omit X in l(θ; X) and L(θ; X).
If θ is a (p × 1) vector of parameters, then I(θ) and J(θ) are (p × p) matrices.
∂ 2 l(θ) ∂ 2 l(θ)
{J(θ)}rs = − and {I(θ)}rs = E − .
∂θr ∂θs ∂θr ∂θs
so
n n 1 X
l µ, σ 2 = − log 2π − log σ 2 − 2 (Xi − µ)2 .
2 2 2σ
i
92 CHAPTER 7. ESTIMATION METHODS
and
∂l 1 X
= 2 (Xi − µ) ,
∂µ σ
i
∂l n 1 X
2
=− 2 + 4 (Xi − µ)2 ,
∂σ 2σ 2σ
i
∂2l n
2
= − 2,
∂µ σ
2
∂ l 1 X
= − (Xi − µ) ,
∂µ∂σ 2 σ4
i
∂2l n 1 X
2 = 4− 6 (Xi − µ)2 .
∂ (σ 2 ) 2σ σ
i
n 1 P
(Xi − µ)
J(µ, σ ) = 1 P σ
2 2 σ4 i
n .
1 P 2
(Xi − µ) (Xi − µ) − 4
σ4 i σ6 i 2σ
To find I(µ, σ 2 ), use
E (Xi ) = µ,
h i
Var (Xi ) = E (Xi − µ)2 = σ 2 ,
so that n
2
0
I(µ, σ 2 ) = E J(µ, σ 2 ) = σ
n .
0
2σ 4
λα α−1 −λx
fX (x) = x e ,
Γ(α)
λα α−1 −λx
L(α, λ; x) = x e ,
Γ(α)
l(α, λ; x) = − log Γ(α) + α log(λ) + (α − 1) log x − λx,
!
∂l(α,λ;x)
∂α −ψ(α) + log(λ) + log x
∂l(α,λ;x) = α ,
∂λ λ −x
∂ 2 l(α,λ;x) ∂ 2 l(α,λ;x)
!
−ψ 0 (α) λ1
∂α2 ∂α∂λ = 1
∂ 2 l(α,λ;x) ∂ 2 l(α,λ;x)
λ − λα2
∂α∂λ ∂λ2
0
ψ (α) − λ1
I(α, λ) = .
− λ1 α
λ2
Chapter 8
8.1 Robustness
In statistical analysis, the concept of robustness refers to the property of a statistical method or estimator
to be relatively unaffected by small deviations from the model’s underlying assumptions. While many
classical statistical techniques are derived under ideal conditions (e.g., data drawn from a perfectly
normal distribution), real-world data often violate these assumptions. A robust procedure is specifically
designed to be resistant to such violations, providing reliable and stable results even when the data is
“messy.”
The most common and impactful violation is the presence of outliers (data points that deviate
strongly from the rest of the sample). These values can arise from measurement error, data entry
mistakes, or simply from the natural variation of a process. Non-robust methods can be badly influenced
by a single outlier, leading to misleading conclusions. Therefore, robustness is, in essence, a measure of
an estimator’s insensitivity to outliers and other small-scale model failures.
The difference between robust and non-robust estimators can be understood by how they weight
the data. Non-robust estimators are often highly efficient under perfect assumptions but give equal or
significant weight to all data points, including extreme outliers. Robust estimators, in contrast, are
designed to automatically down-weight or ignore extreme values, focusing the analysis on the core
structure of the data.
Classical non-robust estimators are typically moment-based and are highly sensitive to outliers
because they incorporate the value of every data point directly. In contrast, many common robust
estimators are quantile-based, as quantiles determine their value from the order of the data, not the
extreme values. However, the field of robust statistics also provides powerful moment-based estimators
(e.g., trimmed means) that achieve robustness by down-weighting or ignoring extreme data points.
Let us introduce some definitions formally.
8.1.1 Definitions
Definition 8.1.1. Given a sample of size n, we can consider the ordered data (order statistics):
The pth quantile (0 ≤ p ≤ 1), denoted qp , is a value below which approximately a fraction p of the
observations fall.
So, the median q0.5 is a value such that roughly half of the observations are less than it and half are
greater. In practice, calculating a specific quantile from data requires an interpolation rule to handle
values between data points. There are several (at least 9) different conventions for this though the
93
94 CHAPTER 8. DATA AND STATISTICAL MODELLING
differences are often minor. Beware that different statistical programs use different defaults (see example
in R documentation).
The median is a robust “competitor” to the mean as a measure of the location or central tendency
of a dataset.
The lower and upper quartiles are given by q0.25 and q0.75 , respectively.
Definition 8.1.2. The interquartile range (IQR) is q0.75 − q0.25 and the semi-interquartile range
is (q0.75 − q0.25 )/2.
These measure the variability of the middle 50% of the data and are robust “competitors” to the
standard deviation as measures of the scale or spread of the data.
Definition 8.1.3. The range is the difference between the largest and smallest observations x(n) − x(1) .
Definition 8.1.4. The trimmed mean is defined as the arithmetic mean of the remaining data, after
we discard the lowest 100α% and the highest 100α% of the data
x([nα]+1) + . . . + x(n−[nα])
x̄α = .
n − 2[nα]
The code to calculating quantiles, mean, median, variance, standard deviation, skewness, kurtosis
and range in R is
# library to calculate skewness and kurtosis
library ( e1071 )
mysample <-c (1:10)
mean ( mysample )
sd ( mysample )
var ( mysample )
skewness ( mysample , type = 1)
kurtosis ( mysample , type = 1)
range ( mysample )
quantile ( mysample )
One of the most common ways (Tukey’s rule, used by default in R function boxplot()) to define
the outliers is to evaluate if the point is located outside of (q0.25 − 1.5 · IQR, q0.75 + 1.5 · IQR).
Sometimes, a wider range is used (q0.25 − 3 · IQR, (q0.75 + 3 · IQR.
8.1.2 Examples
Example (Robust estimators)
Consider the data: 18, 18, 19, 22, 22, 23, 25, 25, 25, 28, 30, 31, 33, 37, 41. Find mean, trimmed
mean, variance, standard deviation, skewness, kurtosis and quantiles. Next, suppose that the
dataset is now altered to include an additional observation with value 100. Compute the mean,
skewness, kurtosis, trimmed mean, median, variance, standard deviation and quantiles and reflect
on the result.
The mean (and first moment about zero) is
1 397
x̄ = m1 (0) = (18 + 18 + · · · + 41) = = 26.5
15 15
8.1. ROBUSTNESS 95
Many statistical methods and theories were developed to analyze normally distributed data. However,
in real life, we often observe data that don’t come from normal distributions. Sometimes, it is possible
to transform them into normal. Below, there are two examples of such data and transformations.
Figure 8.1: Comparison of original and log-transformed Ozone data with normal density curves
The Shapiro-Wilk test [5] was introduced in 1965 by Samuel Shapiro and Martin Wilk in their
landmark paper “An Analysis of Variance Test for Normality (Complete Samples)” published in
Biometrika. It was designed as a powerful statistical test to assess whether a sample comes from
a normally distributed population. Unlike earlier methods, Shapiro and Wilk’s approach uses the
correlation between the ordered sample values and their expected normal scores, making it highly
sensitive to departures from normality, especially for small sample sizes. Today, it remains one of the
most widely used normality tests in statistics and is implemented in most statistical software, including
R and Python.
The USArrests dataset provides statistics on violent crime arrests per 100,000 residents across
the 50 U.S. states in 1973, including variables for murder, assault, rape, and the percentage of
the urban population. The Assault variable exhibits right skewness, which can be mitigated by
98 CHAPTER 8. DATA AND STATISTICAL MODELLING
applying a square root transformation. This transformation helps bring the distribution closer to
normal, improving the suitability for parametric analysis.
Figure 8.2: Comparison of original and square root transformed Assault data with normal density
curves
Note
In the analysis of real-world data, the assumptions of classical parametric statistics are often
ideals rather than realities. Data is frequently contaminated with outliers, skewed, or exhibits
heavy tails, violating the assumptions of normality and homoscedasticity. This is where robust
methods and data transformations become so useful. Robust statistics, such as the median
and interquartile range, provide stable and reliable estimates of location and scale even in the
presence of outliers, making sure that a handful of unusual data points do not distort the entire
narrative of the analysis. Transformations like (e.g., logarithm or square root) can stabilize
variance, improve symmetry, and make data more appropriate for the use of both classical and
8.3. STATISTICAL MODELLING 99
robust techniques.
M = {f (x; θ) | θ ∈ Θ}
where f (x; θ) is the joint probability density (or mass) function; θ is the parameter vector. It is a
collection of the specific, fixed numbers that define a particular distribution within our model family.
Example (Women heights)
Consider modeling the height of adult women in a city. We might assume that heights are
normally distributed.
We then assume, that the model family F is Normal distribution family N (µ, σ 2 ); the pa-
rameter vector θ = (µ, σ 2 ); the parameter space Θ is the set of all possible values these
parameters can take. For this example, the mean can be any real number, and the variance must
be positive (Θ = R × (0, ∞)).
Therefore, the model M is the set of all Normal distributions we could get by plugging in different
values for (µ, σ 2 ), or mathematically
2. Efficiency (related to variance) quantifies the precision or spread of the sampling distribution.
3. Consistency describes the asymptotic behavior of the entire sampling distribution, confirming
whether uncertainty vanishes with more data.
1. Sampling Uncertainty: This is the uncertainty about θ coming from the fact that we work with
a finite sample. It is the uncertainty captured by the sampling distribution of T (X). Probabilistic
statements like confidence intervals are designed to quantify this: a 95% confidence interval is
an interval that, under repeated sampling, would contain the true θ 95% of the time. The null
hypothesis uses the sampling uncertainty to evaluate a specific claim. It asks: “Assuming the
null hypothesis is true, is the observed statistic T (X) compatible with the sampling variability we
would expect to see?” The p-value that results is the probability, under the null hypothesis, of
obtaining a result as extreme as or more extreme than the one observed. A small p-value indicates
that the observed data is unlikely under the assumed model of the null hypothesis, given the
expected sampling uncertainty.
2. Structural Uncertainty (or Model Uncertainty): This uncertainty is about the model
specification M itself. It questions whether the assumed family of distributions F is an adequate
approximation of the true DGP. This includes uncertainty about:
While sampling uncertainty is quantified by the model (e.g., with standard errors), structural
uncertainty must be assessed through model criticism and diagnostics.
ei = Xi − X̂i
If the model is well-specified, the residuals should behave in a manner consistent with the model’s
assumptions (e.g., appearing random and homoscedastic). Systematic patterns in the residuals provide
evidence of structural uncertainty.
Beyond residual analysis, a variety of formal methods are used for assessment and refinement:
1. Goodness-of-Fit Tests: A null hypothesis can be formulated that the data follows the proposed
model, H0 : Model M is adequate. A low p-value from a test statistic (e.g., a chi-square statistic)
leads to rejecting the null, forcing a search for a better model. Examples are: checking assumptions
on normality (e.g. Q-Q plots) or assumptions on equal variance.
8.3. STATISTICAL MODELLING 101
2. Comparing Nested Models: In frameworks like linear regression and ANOVA, a complex
model M1 can be compared to a simpler, nested model M0 (where M0 is a special case of M1 ).
We test the null hypothesis H0 : The additional parameters in M1 are zero. Rejecting H0 with
an F-test justifies refining the model from M0 to M1 by adding those parameters. Examples are:
checking an assumption on linearity.
These diagnostic procedures provide a principled path for model refinement, such as transforming
variables, adding interaction terms, or changing the underlying probability distribution. In the following
chapters, we will see how this general framework-data gathering, model formulation, quantification of
sampling uncertainty, and diagnostic assessment of structural uncertainty-is implemented in specific,
powerful models like linear regression and ANOVA.
102 CHAPTER 8. DATA AND STATISTICAL MODELLING
Exercise Sheet
Note
This exercise sheet covers material from Sections 7 and 8 of the course. Please
attempt all questions and prepare your solutions before the Week 7 workshop.
Your solutions will not be assessed. Before the workshop, you will be invited to vote for the
questions you would like to see discussed. The tutor will select the most popular questions for
discussion during the workshop. Complete solutions for all questions will be available afterwards
(by the end of the week). The solutions will be available only in HTML version of the notes.
Exercise 8.1
An event E has an unknown probability θ of occurring on each of a sequence of independent
trials, where θ is believed to be small. To obtain information about θ, trials are carried out until
E occurs for the first time. There are x consecutive trials without event E occurring before it
occurs for the first time. As a prior distribution for θ is used the density
Exercise 8.2
Suppose to have n observations X1 = x1 , . . . , Xn = xn of a random sample from a Poisson
distribution with mean θ (θ > 0) and probability function
θx e−θ
P (X = x; θ) = , x = 0, 1, 2, . . . .
x!
Consider using a prior density for θ which is a gamma(α, λ) distribution with α and λ known.
1. Show that the posterior density of θ is also a gamma distribution and obtain its parameters.
(a) What does the mean of the posterior distribution of θ then equal?
(b) What does the mode of the posterior distribution then equal?
(c) Comment briefly on your result.
103
104 CHAPTER 8. DATA AND STATISTICAL MODELLING
Exercise 8.3
Suppose a random variable Y | φ ∼ Exponential(φ) and φ ∼ Uniform(1, 2).
What is the density of φ | X = 0.5?
Exercise 8.4
Let y be the number of heads in n spins of a coin, whose probability of heads is θ.
1. If your prior distribution for θ is uniform on the range [0, 1], derive your prior predictive
distribution for y,
Z 1
P (y = k) = P (y = k|θ)dθ,
0
for each k = 0, 1, . . . , n.
2. Suppose you assign a Beta(α, β) prior distribution for θ, and then you observe y heads out
of n spins. It can be shown algebraically that your posterior mean of θ always lies between
your prior mean, α+βα
, and the observed relative frequency of heads, ny . This is provided as
food for thought and you do not need to show it.
3. Show that, if the prior distribution on θ is uniform, the posterior variance of θ is always
less than the prior variance.
Exercise 8.5
Let y be the number of heads in n spins of a coin, whose probability of heads is θ.
Give an example of a Beta(α, β) prior distribution and data y, n, in which the posterior variance
of θ is higher than the prior variance.
Hint: you can use R to find such an example.
Exercise 8.6
i.i.d.
Suppose Xi ∼ N(0, σ 2 ) with σ unknown. Find the σ̂M
2 2 2
LE , Var(σ̂ ) and In (σ ).
Exercise 8.7
Suppose a random sample X1 , X2 , . . . , Xn of n values is taken from a distribution with probability
density function fX (x; θ) = θ(1 − x)θ−1 with 0 < θ < 1. Obtain the maximum likelihood estimate
for θ. Find the Fisher’s Information In (θ).
8.3. STATISTICAL MODELLING 105
Exercise 8.8
Suppose a single value X is available from a geometric(θ) distribution with 0 < θ < 1. Obtain
the maximum likelihood estimate for θ and obtain its mean. Find its Fisher’s information.
106 CHAPTER 8. DATA AND STATISTICAL MODELLING
Chapter 9
Hypothesis testing
The main content of this chapter relies on Prof Simon Myers’ teaching material ([Link]
[Link]/~myers/stats_materials.html) and [3].
A question associated with the data might be formulated in terms of a hypothesis. In particular,
we have a so-called null hypothesis which refers to some basic premise which to we will adhere unless
evidence from the data causes us to abandon it.
In a clinical treatment trial data may be collected to compare two treatments (old versus new).
The null hypothesis is likely to be: no difference between treatments.
The alternative hypothesis might be:
In general we are often in a position to specify the form of the PMF, p(x; θ), say, or the PDF, f (x; θ),
but there is doubt about the value of the parameter θ.
All that is known is that θ is some element of a specified parameter space Θ. We assume that the
null hypothesis of interest specifies that θ belongs to a subset Θ0 ⊆ Θ. Thus, the null hypothesis is true
if θ ∈ Θ0 and false otherwise.
As a convention, we denote the complement of Θ0 in Θ by Θ1 = Θ \ Θ0 . We call the original
hypothesis that θ ∈ Θ0 the null hypothesis, denoted by H0 . The hypothesis that θ ∈ Θ1 is referred
to as the alternative hypothesis, denoted by H1 .
107
108 CHAPTER 9. HYPOTHESIS TESTING
Definition 9.1.1. A hypothesis test is conducted using a test statistic (X1 , X2 , . . . , Xn ) whose distri-
bution is known under the null hypothesis H0 , and is used to consider the likely truth of the null
hypothesis as opposed to a stated alternative hypothesis H1 . With hypothesis test we also associate
the critical region C1 ⊂ Rn , that is connected with rejection of the hypothesis H0 .
Definition 9.1.2. Type I error is the error of rejecting the null hypothesis when it is, in fact, true.
Definition 9.1.3. Type II error is the error of not rejecting the null hypothesis when it is, in fact,
false.
Definition 9.1.4. The p-value (or significance level or size) is the probability of the test statistic
taking a value, in the light of the alternative hypothesis, at least as extreme as its observed value. It is
calculated under the assumption that the test statistic has the distribution which it would have if the
null hypothesis were true.
If the alternative hypothesis is two-sided it will usually be the case that extreme values occur in two
disjoint regions, referring to two tails of the distribution under the null hypothesis.
9.1. BASIC DEFINITIONS 109
2. H0 : X ∼ N (0, 1) (The population follows a normal distribution with mean 0 and variance
1.)
Rules are developed in order to decide whether a so-called null hypothesis H0 or an alternative
hypothesis H1 should be considered as being true.
Below we consider to-called Neyman-Pearson lemma.
sup {L(θ; x) : θ ∈ Θ0 }
λ(x) = , x ∈ RnX .
sup {L(θ; x) : θ ∈ Θ}
In other words, the supremum is the least upper bound of the set of values {f (x) : x ∈ A}. It is the
smallest real number that is greater than or equal to every value of f (x) on A.
If the supremum is attained at some point x0 ∈ A, i.e., f (x0 ) = supx∈A f (x), then the supremum is
equal to the maximum of f on A.
Note (Intuition)
For a given observation x, we compare its maximum likelihood under the null hypothesis to its
maximum likelihood over the entire parameter space. Since the numerator is constrained to Θ0 ,
the ratio λ(x) is always less than or equal to 1. A small value of λ(x) suggests that x is unlikely
under H0 , providing evidence against the null hypothesis.
C1 = {x : λ(x) ≤ k} ,
sup P (λ(X) ≤ k | θ) = α.
θ∈Θ0
In the special case where H0 is a simple hypothesis, i.e., Θ0 = {θ 0 }, the condition simplifies to:
P (λ(X) ≤ k | θ 0 ) = α.
To determine the threshold k, we must examine the cumulative distribution function (CDF) of the
random variable λ(X), where the sample X has joint probability density function:
fX (x; θ 0 ).
In many practical situations, the null hypothesis is simple, meaning it specifies a single value of the
parameter, while the alternative hypothesis is composite, allowing a range of values.
For example:
H0 : θ = θ 0 versus H1 : θ ∈ Θ1 ,
where Θ1 is a subset of the parameter space excluding θ 0 .
In this case, the likelihood ratio simplifies to:
L(θ 0 ; x)
λ(x) = .
sup {L(θ; x) : θ ∈ Θ1 ∪ {θ 0 }}
The critical region is still of the form:
C1 = {x : λ(x) ≤ k} ,
P (λ(X) ≤ k | θ 0 ) = α.
9.1. BASIC DEFINITIONS 111
C1 = {x ∈ Rn : −2 log λ(x) ≥ c} .
112 CHAPTER 9. HYPOTHESIS TESTING
Define:
Λ(x) = −2 log λ(x) = 2 l(θ;
b x) − l(θ0 ; x) ,
where l(θ; x) is the log-likelihood function and θb is the maximum likelihood estimator (MLE) and θ = θ0
under null hypothesis. The critical region becomes:
C1 = {x ∈ Rn : Λ(x) ≥ c} ,
H0 : θ = θ0 vs H1 : θ 6= θ0 ,
gives θb = X.
Thus:
X
Λ = 2n θ0 − X + X log .
θ0
Here dim(Θ) = 1 and dim(Θ0 ) = 0, so degrees of freedom = 1. Critical value at α = 0.05:
χ21 (0.05) = 3.84.
Note (Notice)
In R the function qchisq() the default flag [Link]=TRUE. This means that you would
need either change the flag or enter qchisq(0.95,1).
Critical region:
X
C1 = x : 2n θ0 − X + X log ≥ 3.84 .
θ0
9.1. BASIC DEFINITIONS 113
Note
Recall, in the general notation:
H0 : θ = θ 0 versus H1 : θ ∈ Θ1 ,
2. L1 = sup {L(θ; x) : θ ∈ Θ} is the maximum likelihood over the full parameter space
Θ = {θ > 0}
Likelihood Function
For the exponential(θ) distribution, the likelihood function is:
n
!
X
n
L(θ; x) = θ exp −θ xi = θn exp(−nθx)
i=1
This represents how well the data fit the model when we assume θ = 1.
L1 : Maximized likelihood under H0 ∪ H1
Under the union of hypotheses, θ > 0 can be any positive value. The maximum likelihood
estimator is:
1
θ̂ =
x
Substituting into the likelihood function:
n
1 1 1
L1 = L(θ̂) = exp −n · · x = n exp(−n)
x x x
This represents the best possible fit to the data when θ is freely estimated.
9.1. BASIC DEFINITIONS 115
L0 exp(−nx)
λ(x) = = 1 = xn exp(n(1 − x))
L1 x n exp(−n)
Examples
Example (Critical region)
ZZ Z1 Z1
P((X1 , X2 ) ∈ C1 |H0 ) = 2x1 · 2x2 dx1 dx2 = 4x1 x2 dx1 dx2 =
C1 k k/x1
Z1
k2
2x1 12 − 2 dx1 = 1 − k 2 + 2k 2 log k
x1
k
We can see that the function is strictly descreasing and is continuous. Moreover,
lim = 1, lim = 0.
k→0 k→1
1 − β = P((X1 , X2 ) ∈ C1 |H1 )
Z1 Z1
= 9x21 x22 dx1 dx2 = 1 − k 3 + 3k 3 log k = 0.10
k k/x1
H0 : θ = θ0 against H1 : θ > θ0 .
1 Pn
where x = n i=1 xi .
Define:
L0 = L(θ0 ; x) = θ0n e−nθ0 x , L1 = sup{L(θ; x) : θ ∈ Θ1 ∪ {θ0 }}.
The likelihood ratio is:
L0 L(θ0 ; x)
λ(x) = = .
L1 sup{L(θ; x) : θ ≥ θ0 }
To compute L1 , note that the log-likelihood is:
Differentiating:
∂l(θ; x) n
= − nx,
∂θ θ
which is zero at:
1
θb = .
x
Behavior of L(θ; x):
2. If 1/x < θ0 , then the unconstrained maximum is outside the region. Since L(θ; x) decreases
for θ > 1/x, the supremum occurs at the smallest allowed value, θ = θ0 . Thus:
Combine cases: (
x−n e−n , if 1/x ≥ θ0 ,
sup L(θ; x) = n −nθ x
θ≥θ0 θ0 e 0 , if 1/x < θ0 .
Since λ(x) is non-decreasing in x, the critical region of the likelihood ratio test is:
n
( )
X
C1 = x : xi ≤ c or equivalently x ≤ c/n.
i=1
Proof
Proof. Let us consider the case H1 : µ 6= µ0 .
We have
Θ = (µ, σ 2 ) : µ ∈ R, σ 2 ∈ R+
Θ0 = (µ, σ 2 ) : µ = µ0 , σ 2 ∈ R+
and
1 1
f (x; µ, σ 2 ) = √ exp − 2 (x − µ)2 , x ∈ R.
2πσ 2 2σ
120 CHAPTER 9. HYPOTHESIS TESTING
Since
n
n 1 X
l(µ0 , σ 2 ; x) = − log(2πσ 2 ) − 2 (xi − µ0 )2
2 2σ
i=1
and
n
∂l n 1 X
= − + (xi − µ0 )2 ,
∂σ 2 2σ 2 2σ 4
i=1
which is zero when
n
1X
σ2 = (xi − µ0 )2
n
i=1
we conclude that
n
!−n/2
2π X
sup L(µ0 , σ 2 ; x) = (xi − µ0 )2 e−n/2 .
n
i=1
For the denominator, we already know from previous examples that the MLE of µ is x, so
n
!−n/2
2π X
sup L(µ, σ 2 ; x) = (xi − x)2 e−n/2
n
i=1
and Pn −n/2
i=1 (xi − µ0 )2
λ(x) = Pn 2
.
i=1 (xi − x)
This may be written in a more convenient form. Note that
n
X n
X
(xi − µ0 )2 = ((xi − x) + (x − µ0 ))2
i=1 i=1
n
X
= (xi − x)2 + n(x − µ0 )2
i=1
so that −n/2
n(x − µ0 )2
λ(x) = 1 + Pn 2
.
i=1 (xi − x)
The critical region is
C1 = {x : λ(x) ≤ k}
so it follows that H0 is to be rejected when the value of
|x − µ0 |
pPn
2
i=1 (xi − x)
Now
X − µ0
√ ∼ t(n − 1)
s/ n
where
n
1 X
s2 = (Xi − X)2 .
n−1
i=1
which is the standard form of the two-sided t-test for a single sample.
Lemma 9.2.3. Suppose a random sample Xi ∼ N(µ, σ 2 ) size of n. Suppose µ is unknown. We wish to
test H0 : σ = σ0 against H1 : σ 6= σ0 (or one-sided alternatives). Here the likelihood ratio is a function
of the test statistic
n
(n − 1)s2 2 1 X
, s = (Xi − X)2
σ02 n−1
i=1
Next we consider two populations, one N(µX , σX 2 )-distributed and the other N(µ , σ 2 )-distributed,
Y Y
from which respectively the following independent samples are drawn: X = (X1 , . . . , Xm ) and Y =
m
1 Pm 1 Pn 2 1 P
(Y1 , . . . , Yn ). We will use the notation: X = m i Xi , Y = n i Yi , sX = m−1 (Xi − X)2 and
i=1
n
1
s2Y X)2 .
P
= n−1 (Xi −
i=1
s2Y
,
s2X
n−1
which, under H0 is Fm−1 -distributed.
Y −X −∆
q ,
σX2 /m + σ 2 /n
Y
Lemma 9.2.6. We assume that σX and σY are equal, but unknown, denoting the unknown common
value of σX and σY by σ. We would like to test H0 : µY − µX = ∆ against H1 : µY − µX 6= ∆ (or
one-sided alternatives), where ∆ ∈ R. The test statistic is
Y −X −∆ (m − 1)s2X + (n − 1)s2Y
p , s2p = .
sp 1/m + 1/n m+n−2
H0 : θ = 0 against H1 : θ 6= 0.
Assuming the data to be normally distributed (for the sake of this example), the mean difference
is x = −4 and the estimated standard deviation is s = 10.744. Under H0 we obtain a t-statistic
of √
√ x −4 8
t= n = = −1.053.
s 10.744
The t-statistic has 7 degrees of freedom for a p-value of 0.327. The significance level is set at 0.05
and H0 is rejected if t ≤ −2.365 or t ≥ 2.365 .
9.2. HYPOTHESIS FOR NORMALLY DISTRIBUTED POPULATIONS 123
If the null hypothesis does not completely specify the distribution of the Xi ’s, it is called
composite. If the null hypothesis were refined to state that the distribution was Poisson with a
specified mean, then it would be simple. The alternative hypothesis does not completely specify
the distribution, so it is composite.
We now develop a generalized likelihood ratio test for the goodness of fit of a model for multino-
mial cell probabilities. The multinomial distribution generalizes the binomial distribution: instead
of one event with probability p, we have k categories (types of events). Each category i has probability
pi and observed count xi , with total trials n satisfying ki=1 xi = n. More details about the multinomial
P
distribution can be found in Wikipedia: Multinomial distribution.
Under the model, the vector of cell probabilities p = (p1 , . . . , pk ) is described by the null hypothesis
H0 , which specifies that:
p = p(θ), θ ∈ Θ0 ,
where θ is a parameter that may be unknown.
Suppose we fit Poisson probabilities depending on an unknown parameter to the cell counts in a
table. We want to judge the plausibility of this model relative to an alternative H1 . Under H0 ∪ H1 the
9.3. GOODNESS-OF-FIT TEST 125
cell probabilities are free except for the constraints that they are non-negative and sum to 1. If there
are k cells, then:
k
X
Θ = {(p1 , . . . , pk ) : pi ≥ 0, pi = 1}.
i=1
The likelihood under the multinomial model is:
k
Y pi (θ)xi
L(θ; x) = n! ,
xi !
i=1
where xi are the observed counts in the k cells. By definition and properties of the maximum likelihood
estimate, this likelihood is maximized when θb is the MLE of θ. The corresponding probabilities are
denoted by pi (θ),
b so:
k b xi
Y pi (θ)
sup L(θ; x) = n! .
θ∈Θ0 xi !
i=1
To find supθ∈Θ L(θ; x), we solve the optimization problem under the constraint ki=1 pi = 1 using
P
Lagrange multipliers:
∂
l − γ(g(p) − 1) = 0,
∂pi
Pk
where l = log L and g(p) = i=1 pi . This gives:
xi X X xi
−γ =0 ⇒ xi = γ pi ⇒ γ = n, and pbi = .
pi n
Thus, the likelihood ratio statistic is:
supθ∈Θ0 L(θ) h i
Λ = −2 log = 2 sup l(θ) − sup l(θ) .
supθ∈Θ L(θ) θ∈Θ θ∈Θ0
Compute:
k k k k
" # !
x xi Oi
i
X X X X
Λ=2 xi log − xi log pi (θ)
b =2 xi log =2 Oi log ,
n npi (θ)
b Ei
i=1 i=1 i=1 i=1
It can be shown that Pearson’s statistic and the likelihood ratio statistic are asymptotically equivalent
under H0 (see [4]).
126 CHAPTER 9. HYPOTHESIS TESTING
9.4 Examples
Example (Genotypes)
Under the Hardy-Weinberg equilibrium model, the genotype probabilities are given by:
p2 + 2pq + q 2 = 1,
where p2 is the dominant homozygous probability (AA), 2pq is the heterozygous probability (Aa)
and q 2 is the recessive homozygous probability (aa).
We will use the notation: θ = p and 1 − θ = q.
Suppose x = {x1 , x2 , x3 } are observed frequencies. Recall, in the multinomial model
k
X k
X
log-Likelihood: l(θ; x) = xi log pi (θ) + log n! − log xi ! =
i=1 i=1
3
X
log n! − log xi ! + 2x1 log(1 − θ) + x2 log(2θ(1 − θ)) + 2x3 log θ =
i=1
3
X
log n! − log xi ! + (2x1 + x2 ) log(1 − θ) + (2x3 + x2 ) log θ + x2 log 2,
i=1
∂l(θ, x) 2x1 + x2 2x3 + x2
Maximising: =− + = 0 =⇒
∂θ 1−θ θ
2x3 + x2 2x3 + x2
θb = = , n = x1 + x2 + x3 .
2(x1 + x2 + x3 ) 2n
In a sample from the Chinese population of Hong Kong in 1937, blood types occurred with the
following frequencies, where M and N are erythrocyte antigens. Below is the table with the
frequencies for the blood types. In this MN blood group system, both alleles are codominant -
there is no recessive or dominant allele. We will use θ frequency of allele M and (1 − θ) for allele
N.
M MN N Total
Frequency 342 500 187 1029
H0 : The data fit into the Hardy-Weinberg Law where the frequencies are described by the
probabilities (1 − θ)2 , 2θ(1 − θ), θ2 .
H1 : The data don’t fit into this model.
Recall that MLE for θ is (according to Hardy-Weinberg Law)
2x3 + x2 2x3 + x2
θb = = , n = x1 + x2 + x3 ,
2(x1 + x2 + x3 ) 2n
θb = 0.4247, n = 1029 =⇒
Expected: x1 = n(1 − θ)2 , x2 = 2nθ(1 − θ), x3 = nθ2
The calculations give us
M MN N Total
Observed 342 500 187 1029
Expected 340.6 502.8 185.6
9.4. EXAMPLES 127
We can see that the result is the same as the one obtained using the likelihood ratio test statistic.
Below are the data on the number of flying bomb hits recorded in each of 576 small areas of
1 2
4 km in the south of London during World War II.
i 0 1 2 3 4 5 ≥6
πi (θ̂) 0.3949 0.3669 0.1704 0.0528 0.0123 0.0023 0.0004
and hence
k
!
X ni
Λ=2 ni log = 1.4995.
i=1 nπi (θ)
b
Though the p-value become smaller we still can not reject the null hypothesis.
Using the Pearson’s chi-squared approach we obtain p ≈ 0.796.
(See the code in the HTML version of the notes).
Note
Notice, that in R we can also use function [Link](). The function performs Pearson’s
chi-squared test. Be careful, the function cannot automatically adjust for parameter
estimation from the same data.
Exercise Sheet
Note
This exercise sheet covers material from Section 9 of the course. Please attempt all
questions and prepare your solutions before the Week 9 workshop.
Your solutions will not be assessed. Before the workshop, you will be invited to vote for the
questions you would like to see discussed. The tutor will select the most popular questions for
discussion during the workshop. Complete solutions for all questions will be available afterwards
(by the end of the week). The solutions will be available only in HTML version of the notes.
Exercise 9.1
A sample X1 , X2 , . . . , Xn is drawn from a N(µ, σ 2 )-distribution (with known µ) is order to test
H0 : σ 2 = σ02 against H1 : σ 2 6= σ02 . Perform the test.
Exercise 9.2
Suppose that deaths in crashes involving a particular make of car have been at an average rate of
6 per week and that the company has introduced air bags. They want to use the figures over
the next year (i.e. 52 weeks) to test their effectiveness. Assume the data are from a Poisson
distribution with mean µ. The company plans to test
H0 : µ = 6 against H1 : µ < 6,
C1 = {y : y ≤ k} .
Find the critical region C1 = {y : y ≤ k} for the test at significance level α = 0.05 using the
Poisson/normal approximation. Find the power of the test.
Exercise 9.3
Suppose you throw a coin; the probability that a head will be observed is θ. You wish to test
H0 : θ = 12 against H1 : θ = 34 . TO do so, you draw a sample X1 , X2 , where Xi , i = 1, 2 is the
number of throws needed before for the first time a heard is observed. Now, Xi is geometrically
distributed with parameter θ. More specific, under H0 the Xi have a discrete probability density
129
130 CHAPTER 9. HYPOTHESIS TESTING
n 1 n−1 3
given by p0 (xi = n) = (1 − θ)n−1 θ = 21 and p1 (xi = n) = (1 − θ)n−1 θ =
4 4.
Suppose the critical region of size α is
C1 = {(n1 , n2 ) ∈ N × N : 2 ≤ n1 + n2 ≤ c}.
1. Show that C1 is, among all critical regions of size α, the one of maximal power.
4. For the critical region formed by sample (1, 1), (2, 1), (3, 1), (1, 3), find significance α and
power 1 − β.
Exercise 9.4
Starting from an exponentially distributed population with parameter θ, we wish to test H0 : θ = 1
against H1 : θ 6= 1.
Recall, that it can be shown that the λ(x) = (x̄e1−x̄ )n (see the Example for the exponential
distribution 9.1.1)
(1) (2) (1)
1. Show that the critical region C1 = C1 ∪ C1 , where C1 = {x ∈ Rn+ : x̄ ≤ c1 } and
(2)
C1 = {x ∈ Rn+ : x̄ ≥ c1 }.
1
2. Now consider sample X1 , X2 and c1 = 2 and c2 = 32 . Determine the level of significance α.
Exercise 9.5
Suppose X1 , . . . , Xn is a random sample from a Poisson(θ1 ) distribution, and Y1 , . . . , Yn is an
independent random sample from a Poisson(θ2 ) distribution. Derive an expression for the
likelihood ratio λ(x) for H0 : θ1 = θ2 , versus H1 : θ1 6= θ2 .
Exercise 9.6
The samhda (package UsingR) data set contains information about health behavior for school-age
children. For example, the variable [Link] measures how often a child smoked in the previous
month. There are seven levels: a 1 means he smoked every day and a 7 means not at all. Values
98 and 99 indicate missing data. See ?samhda for a description.
Perform the goodness of fit test to investigate whether the sample proportions are statistically
different from the probabilities:
Exercise 9.7
A die is rolled 100 times and yields these frequencies
9.4. EXAMPLES 131
side 1 2 3 4 5 6
count 13 17 9 17 18 26
Is this a fair die? Answer using a significance test with H0 : π = 1/6 for each i and H1 : π 6= 1/6
for at the least one i.
132 CHAPTER 9. HYPOTHESIS TESTING
Chapter 10
Linear Regression
Note
In this chapter, we will study the multiple linear regression. First, we will cover its mathematical
formulation, parameter estimation using the least squares method, and techniques for statistical
inference. Next, we will learn how to check model assumptions through diagnostic plots, assess
model fit, address multicollinearity, and apply methods for model selection. All concepts will be
illustrated with practical examples in R.
Y = Xβ + ε,
where the components are defined as follows
Y1 1 x1,1 x1,2 · · · x1,p
Y2 1 x2,1 x2,2 · · · x2,p
Y = . ∈ Rn , X = . .. ∈ R
n×(p+1)
,
.. .. ..
.. .. . . . .
Yn 1 xn,1 xn,2 · · · xn,p
β0
β1 ε1
ε2
β = β2 ∈ Rp+1 , ε = . ∈ Rn .
.. ..
.
εn
βp
Let us discuss the ingredients of the above representation in details.
In the above, the X is fixed and known n × (p + 1) matrix and is called the design matrix.
The first column in the design matrix is all ones, representing the intercept term. Each subsequent
column j (for j = 1, . . . , p) contains the observed values of the j-th predictor variable.
133
134 CHAPTER 10. LINEAR REGRESSION
The vector β is the parameter vector containing the regression coefficients (fixed, but unknown).
You can often find in literature, that β0 is the call an intercept term. The coefficient βj (for j = 1, . . . , p)
is the coefficient for the j-th predictor variable.
Finally, ε is error vector representing the unobserved random error terms.
We further assume that the errors are spherical. This means the errors are uncorrelated and share
a constant variance (homoscedasticity):
Cov(ε) = σ 2 In ,
ε ∼ N (0, σ 2 In ).
This implies that the errors are independent and identically distributed normal random variables.
Note
The general matrix formulation contains the simple linear regression model as a special case.
When there is only one predictor (p = 1), the design matrix reduces to
1 x1
1 x2
X = . . = (1n x),
.
. . .
1 xn
where x = (x1 , . . . , xn )> , and the parameter vector is β = (β0 , β1 )> . The model equation then
simplifies to the familiar form Yi = β0 + β1 xi + εi for i = 1, . . . , n. All subsequent general results
for estimation and inference, when applied to this simple case, will yield the well-known formulas
for the slope and intercept.
Given the assumptions above, we can also formally prove the following lemma
10.1. MULTIPLE LINEAR REGRESSION 135
Lemma 10.1.1. For the linear model 10.1.1 the following holds
Y ∼ N (Xβ, σ 2 In ).
The proof follows from the lemma 5.3.1. Recall, that the lemma 5.3.1 states that if a random vector
follows multivariate normal distribution X ∼ N (µ, Σ), then any affine transformation AX + b is also
multivariate normal
AX + b ∼ N (Aµ + b, AΣA> ).
Applying to our case (Y = In ε + Xβ ≡ Aε + b)
where the matrix H = X(X> X)−1 X> is commonly called the hat matrix (because it ”puts the hat
on Y”).
The hat matrix possess a set of nice and useful properties. I will list them here without a proof.
2. H is idempotent, i.e. H 2 = H.
The difference between the observed response vector and the fitted values gives the vector of
estimated residuals:
residual vector. It is a random vector that serves as the observable analogue of the true error.
Once we observe the data y, we calculate its realization ε̂ = y − Xβ̂, which we call the vector
of estimated residuals.
The least squares method finds the estimate β̂ that minimizes ε> ε.
Theorem 10.2.1 (Least Squares). For the linear model 10.1.1, assume that the matrix X> X ∈
R(p+1)×(p+1) is invertible. Then the function r(β) = ε> ε, representing the residual sum of squares,
takes its minimum at the vector β̂ ∈ Rp+1 given by
Proof
Proof. Using the model equation ε = Y − Xβ, we can write the residual sum of squares as
n
X
r(β) = ε2i = ε> ε
i=1
>
= (Y − Xβ) (Y − Xβ)
= Y Y − Y Xβ − (Xβ)> Y + (Xβ)> (Xβ).
> >
Using standard linear algebra rules, we find that Y> Xβ = (Xβ)> Y = β > X> Y and
(Xβ)> (Xβ) = β > X> Xβ. Thus we obtain
r(β) = Y> Y − 2β > X> Y + β > X> Xβ.
∂
To find the minimum of this function, we set all partial derivatives ∂β i
r(β) equal to 0. Examining
each term:
1. First term: Y> Y does not depend on β, so ∂β ∂
i
Y> Y = 0 for all i; 2. Second term:
∂ > > ∂ Pp > >
∂βi β X Y = ∂βi j=0 βj (X Y)j = (X Y)i 3. Third term:
p
∂ > > ∂ X
β X Xβ = βj (X> X)j,k βk
∂βi ∂βi
j,k=0
p
X
(X> X)i,k βk = 2 (X> X)β i .
=2
k=0
To confirm this is a minimum, we examine the Hessian matrix. The second partial derivatives
are:
∂
r(β) = 2(X> X)i,j
∂βi ∂βj
so the Hessian matrix is H = 2X> X. Since X has full column rank, X> X is positive definite,
confirming that β̂ is the unique minimum.
Proof
Proof. Substituting the model Y = Xβ + ε into the estimator:
Lemma 10.2.3. Under the normality assumption ε ∼ N (0, σ 2 In ), the estimator follows:
β̂ ∼ N β, σ 2 (X> X)−1 .
Proof
Proof. Since β̂ = β + (X> X)−1 X> ε is a linear transformation of the normal random vector ε,
it follows a multivariate normal distribution. We have E[β̂] = β and the covariance matrix is:
Note
We have used the following properties
1. lemma 5.3.1
>
2. (X> X)−1 = (X> X)−1 or, in words, (X> X)−1 is symmetric
138 CHAPTER 10. LINEAR REGRESSION
1 >
Lemma 10.3.1 (Cochran). 1. ε Hε ∼ χ2 (p + 1),
σ2
1 >
2. ε (In − H)ε ∼ χ2 (n − p − 1),
σ2
3. Hε and (In − H)ε are independent.
1 >
Lemma 10.3.2. The estimator σ̂ 2 = n−p−1 ε̂ ε̂ is unbiased, so:
E(σ̂ 2 ) = σ 2 .
Proof
Proof. We first note that:
(n − p − 1)σ̂ 2 = (Y − Ŷ)> (Y − Ŷ) = (Y − HY)> (Y − HY)
= Y> (In − H)> (In − H)Y = Y> (In − H)Y.
Note
We used
H> = H, H2 = H,
(In − H)> (In − H) = I2n − H> In − In H + H> H = In − H.
Note
We used
HX = X(X> X)−1 X> X = X, (In − H)X = X − HX = 0
By lemma 10.3.1:
1 1
2
(n − p − 1)σ̂ 2 = 2 ε> (In − H)ε ∼ χ2 (n − p − 1).
σ σ
Since the expectation of a χ2 (ν) distribution equals ν (see appendix 12.6), we have:
1
(n − p − 1)E(σ̂ 2 ) = n − p − 1,
σ2
and thus:
E(σ̂ 2 ) = σ 2 .
This completes the proof.
Independence of Estimators
Lemma 10.3.3 (Independence). Under the normality assumption, the estimators β̂ and σ̂ 2 are
conditionally independent given X.
This result follows from the fact that β̂ is a linear function of Y while σ̂ 2 is a quadratic function
involving the residual vector, and these turn out to be independent under normality. This independence
is essential for constructing proper test statistics for hypothesis testing. I leave this theorem without a
proof here.
10.4 Inference
10.4.1 T statistic for individual coefficients
We have seen that β̂ ∼ N (β, σ 2 C), where C = (X> X)−1 . Restricting this to a single coefficient, we find
β̂i ∼ N βi , σ 2 Cii ,
(10.2)
since the diagonal elements of the covariance matrix contain the variances of the elements of a random
vector. In practice we will not know the value of σ 2 , so we have to estimate this from data, using the
estimator
n
1 X
σ̂ 2 = (Yi − Ŷi )2
n−p−1
i=1
from equation 10.3. As a first application of Cochran’s theorem we showed in equation 10.3 that
1
(n − p − 1)σ̂ 2 ∼ χ2 (n − p − 1).
σ2
Note that in the equations above, we index the rows and columns of C using i, j ∈ {0, 1, . . . , p}, i.e.
the first row and column are using the index 0 each. This is to match the convention for the components
of β = (β0 , β1 , . . . , βp ).
140 CHAPTER 10. LINEAR REGRESSION
We now construct a quantity T which measures the distance between the estimated value β̂i and
the unknown true value βi :
β̂i − βi
T =√ 2 . (10.3)
σ̂ Cii
While there are many ways to measure this distance, the T constructed here has two main advantages:
• The value of T can be computed from the given data, without any reference to unknown quantities
(except for βi ).
• Below, we will be able to find the distribution of T . This will allow us to use T to construct
confidence intervals and statistical tests.
Lemma 10.4.1. Assume that the data follows the model Y = Xβ + ε. Then T ∼ t(n − p − 1), i.e. T
follows a t-distribution with n − p − 1 degrees of freedom (see appendix 12.6).
Proof
Proof. We have
√
(β̂i − βi )/ Cii
T = √
σ̂ 2
√
(β̂i − βi )/ σ 2 Cii
= p
σ̂ 2 /σ 2
√
(β̂i − βi )/ σ 2 Cii
=p
((n − p − 1)σ̂ 2 /σ 2 )/(n − p − 1)
Z
=: p ,
Y /(n − p − 1)
√
where Z = (β̂i − βi )/ σ 2 Cii ∼ N (0, 1) and Y = (n − p − 1)σ̂ 2 /σ 2 ∼ χ2 (n − p − 1) are independent,
by lemma 10.3.3. Thus, T ∼ t(n − p − 1) as required.
√
The quantity σ 2 Cii is sometimes called the standard error of the estimator β̂i , denoted by
SE(β̂i ).
where tn−p−1 (α/2) is the (1 − α/2)-quantile of the t(n − p − 1)-distribution. Rewriting this expression
as a condition on β̂i instead of on T gives a confidence interval for βi .
Proof
Proof. We have to show that P (βi ∈ [U, V ]) ≥ 1 − α. We have
βi ∈ [U, V ]
p
⇐⇒ β̂i − βi ≤ σ̂ 2 Cii tn−p−1 (α/2)
β̂i − βi
⇐⇒ √ 2 ≤ tn−p−1 (α/2)
σ̂ Cii
⇐⇒ T ∈ −tn−p−1 (α/2), +tn−p−1 (α/2)
For the case of simple linear regression we know that Y = β0 + β1 X + ε, or in matrix form
Y = Xβ + ε. We know that C = (X> X)−1 and the design matrix is:
1 x1
1 x2
X = . .
.. ..
1 xn
1 xn
Determinant:
X X 2 X
det(X> X) = n x2i − xi =n (xi − x)2
Inverse:
P 2 P P 2 P
Pxi − xi = P Pxi − xi
1 1
(X> X)−1 =
det(X> X) − xi n n (xi − x)2 − xi n
Note
(xi − x)2 = x2i − n(x)2 .
P P
In the above, we used
We obtain: s s
P 2
xi 1
SE(β̂0 ) = σ̂ , SE(β̂1 ) = σ̂ .
n (xi − x)2 (xi − x)2
P P
142 CHAPTER 10. LINEAR REGRESSION
where ε∗ is independent of ε.
The natural point predictor is the estimator:
Note
In R we use the function predict() to predict a value ŷ ∗ for the new set of observations x∗
3. Var(Ŷ ∗ ) = Var((x∗ )> β̂) = (x∗ )> Var(β̂)x∗ and Var(β̂) = σ 2 (X> X)−1
Therefore:
Var(e∗ ) = Var(Y ∗ ) + Var(Ŷ ∗ ) = σ 2 + σ 2 (x∗ )> (X> X)−1 x∗
h i
Var(e∗ ) = σ 2 1 + (x∗ )> (X> X)−1 x∗
Y ∗ − Ŷ ∗
T = p
σ̂ 1 + (x∗ )> (X> X)−1 x∗
This statistic follows a Student’s t-distribution with n − p − 1 degrees of freedom because:
10.4. INFERENCE 143
• Y ∗ is independent of σ̂ 2
A 100(1 − α)% prediction interval for Y ∗ is given by:
q
Ŷ ∗ ± tα/2,n−p−1 · σ̂ 1 + (x∗ )> (X> X)−1 x∗
More explicitly: q
(x∗ )> β̂ ± tα/2,n−p · σ̂ 1 + (x∗ )> (X> X)−1 x∗
The prediction interval accounts for both estimation uncertainty in β̂ and inherent variability
of the new observation Y ∗ around its mean. The term (x∗ )> (X> X)−1 x∗ measures the leverage of x∗ :
points farther from the data center have wider intervals. The additional ”1” under the square root
makes prediction intervals always wider than confidence intervals for the mean response.
It is interesting to note that as n → ∞, the interval approaches Ŷ ∗ ± zα/2 · σ, since estimation
uncertainty vanishes but observation variability remains.
Lemma 10.4.3. The test which rejects H0 if and only if |T | > tn−p−1 (α/2) has significance level α.
Proof
Proof. We have to show that the probability of type I errors (i.e. of wrongly rejecting H0 when
it is true) is less than or equal to α. Assume that H0 is true. Then we have βi = b and thus the
144 CHAPTER 10. LINEAR REGRESSION
T defined in this section coincides with the expression from equation 10.3. From lemma 10.4.1
we know that T ∼ t(n − p − 1). Thus we have
P (type I error) = P |T | > tn−p−1 (α/2)
= P T < −tn−p−1 (α/2) + P T > tn−p−1 (α/2)
= 2P T > tn−p−1 (α/2)
= 2P T > tn−p−1 (α/2)
α
=2
2
= α.
As usual with statistical tests, one needs to be extremely careful when performing several tests on
the same data. In particular, it would be unwise to test more than one component of β using this
procedure for the same data. Instead, in the next section we will consider how to perform tests for
several components of β simultaneously. This can be done using F-statistic.
Y = Xβ + ε, ε ∼ N (0, σ 2 In )
n
X
Total Sum of Squares (SST) : SStotal = (Yi − Y )2 = Y> (I − n1 J)Y
i=1
n
X
Regression Sum of Squares (SSR) : SSreg = (Ŷi − Y )2
i=1
Xn
Residual Sum of Squares (SSE) : SSres = (Yi − Ŷi )2 = Y> (I − H)Y
i=1
The ANOVA identity is not just an algebraic curiosity; it gives the foundation for constructing the
F-test. Let us prove this identity for the case of simple linear regression and understand what does it
mean.
10.4. INFERENCE 145
Proof
Proof. Simple linear regression case
We have
X X X X X X
(Ŷi − Y )ε̂i = Ŷi ε̂i − Y ε̂i = β̂0 ε̂i + β̂1 xi ε̂i − Y ε̂i .
i i i i i
P P
It can be shown that ε̂i = 0 and xi ε̂i = 0. Indeed, writing the least squares
!
X
2
min εi ,
β0 ,β1
i
!
∂ X X
(Yi − β0 − β1 xi )2 = −2 (Yi − β0 − β1 xi ) = 0
∂β0
i i
X
=⇒ ε̂i = 0,
i
!
∂ X X
(Yi − β0 − β1 xi )2 = −2 xi (Yi − β0 − β1 xi ) = 0
∂β1
i i
X
=⇒ xi ε̂i = 0.
i
P
Therefore, (Ŷi − Y )ε̂i = 0 and SStotal = SSreg + SSres .
i
The identity tells us that the total variability in the data SStotal is partitioned into a component
explained by the model SSreg and and an unexplained, residual component SSres .
Let us define the residual mean square (for simple linear regression) to be
SSres
M Sres = .
n−2
It can be shown that
X (n − 2)M Sres
E [M Sres ] = E[ ε̂2 ] = σ 2 , ∼ χ2 (n − 2).
σ2
i
This means that M Sres is an unbiased estimator of the error variance, regardless of whether H0 is
true.
SS
Let us also define the regression mean square to be M Sreg = 1reg (for the case of simple linear
regression). It can be shown that
X
E[M Sreg ] = σ 2 + β12 (xi − x)2 .
i
146 CHAPTER 10. LINEAR REGRESSION
From this we can conclude that under null hypothesis β1 = 0, the ratio M Sreg /M Sres is close to 1.
Under H1 , M Sreg estimates a larger quantity, so the ratio is greater than 1.
The F-statistic formalizes this comparison.
As SSreg and SSres are independent we construct F statistic as
SSreg /1
∼ F (1, n − 2).
SSres /(n − 2)
F ∼ F(p, n−p−1)
where F1−α (p, n − p − 1) is the (1 − α)-quantile of the F -distribution with (p, n − p − 1) degrees of
freedom.
• Normally distributed errors: εi ∼ N (0, σ 2 ) for all i. In particular, the variance of the εi does not
depend on i.
Sometimes, plotting the data can reveal when these assumptions are violated. Here we discuss
different plot types which can be used for this purpose.
To demonstrate the use of residual plots, we first produce such a plot for simulated data, where
we know that the modelling assumptions are satisfied.
n <- 50
x <- seq(0, 10, [Link] = n)
y <- 2 + 0.5 * x + rnorm(n)
m <- lm(y x)
plot(m, which=1, pch = 19, main="modelling assumptions satisfied")
We can see that the residuals are clustered around zero as expected.
The red line is a local polynomial regression fitted through the residuals. It shows the trend
in the residuals as the fitted values change. The prefect horizontal line would mean residuals
have no systematic pattern with fitted values. In practice it rarely happen though, even with
normally distributed residuals. You can observe slight wiggles due to random variation. What
we are looking for is “fairly straight” line.
Next we consider a residual plot where the modelling assumptions are known to be violated. Here
we generate data using a quadratic dependency (instead of linear) of the output on the input:
n <- 50
x <- seq(0, 10, [Link] = n)
y <- 2 + 0.5 * x^2 + rnorm(n)
m <- lm(y ~ x)
plot(m, which=1, pch = 19, main="modelling assumptions violated")
148 CHAPTER 10. LINEAR REGRESSION
It is clear that the residuals are not centred around the horizontal zero line, but form a parabola
instead. Notice, how the red line becomes a very distinct parabola, showing a clear violation of
linearity assumptions.
The structure of the x-values can affect the visual appearance of this type of plot, even if the
modelling assumptions are exactly satisfied. For example, if the x-values ”thin out” towards the edges
of the range, the ”band” of residuals may appear narrower there, just because there are fewer chances
for large residuals to occur near the edges.
The following example shows a residual plot for real data, where it is unknown whether or not the
modelling assumptions are satisfied.
The stackloss dataset is a classic dataset used in regression analysis. It contains data from an
industrial process for producing nitric acid. It includes 21 observations with four variables:
We can see that the residuals are clustered around zero as expected. The red line shows some
nonlinear behaviour, however the magnitude is not high and the pattern is not obvious to naked
eye. Therefore, we make conclusion that the plot does not indicate any problems with model fit.
A Q-Q plot of the residuals can be used to assess whether the residuals are normally distributed. These
plots plot quantiles of the data against quantiles of a normal distribution. If the residual are normally
distributed, the points in a Q-Q plot will approximately lie on a straight line. The R command to
produce Q-Q plots for the normal distribution is qqnorm().
We first show an example with simulated data, where the modelling assumptions are known to
be satisfied.
150 CHAPTER 10. LINEAR REGRESSION
n <- 50
x <- seq(0, 10, [Link] = n)
y <- 2 + 0.5 * x + rnorm(n)
m <- lm(y x)
plot(m, which=2, pch = 19, main="modelling assumptions satisfied")
Next we consider an example where the errors are exponentially instead of normally distributed.
n <- 50
x <- seq(0, 10, [Link] = n)
eps <- rexp(n) - 1 # we subtract 1, so that the errors still have mean zero
y <- 2 + x + eps
m <- lm(y x)
plot(m, which=2, pch = 19, main="modelling assumptions violated")
10.5. DIAGNOSTIC PLOTS 151
Now the points of the Q-Q plot lie on a curved line, indicating that the residuals are not normally
distributed.
Next we consider an example (stackloss dataset) where the errors are exponentially instead of
normally distributed.
m <- lm([Link] ., data = stackloss)
plot(m, which=2, pch = 19, main="stackloss data")
152 CHAPTER 10. LINEAR REGRESSION
The points in this Q-Q plot lie more or less on a straight line, so it seems a reasonable to assume
that the errors were normally distributed.
1. A scatter plot of the response y against each of the explanatory variables xj for j ∈ {1, . . . , p}.
Great care is needed when interpreting such plots: even if the full model is linear, the dependency
of the output on individual inputs may appear non-linear. This can be caused by non-linear
structures in the input data (which don’t violate the modelling assumptions).
2. A scatter plot of ε̂ against each of the explanatory variables xj for j ∈ {1, . . . , p}. Again, care is
needed when trying to interpret such plots.
3. A plot of xj against xk for j 6= k can also be useful. For example, if two or more regressors are
highly correlated, we say that multicollinearity is present. When this occurs, the least squares
estimate β̂ becomes numerically unstable. (We will learn more about this effect in Chapter 10.8.)
4. Data are often collected in time order. Even if time is not an explanatory variable, a plot of y
against time can be of interest. It can reveal serial correlation in the data. Similarly, plotting the
residuals against time.
plot(stackloss$[Link], stackloss$[Link],
xlab = "water temperature", ylab = "Stack loss")
Clearly the stack loss increases with water temperature, but the relationship seems not to be
linear. It may make sense to include the square of the water temperature as an additional input
variable in the model.
SSreg SSres
R2 = =1−
SStotal SStotal
We can easily compute the R2 value for the stackloss dataset manually:
m <- lm([Link] ~ ., data = stackloss)
[Link] <- 1 - var(resid(m)) / var(stackloss$[Link])
[Link]
[1] 0.9135869
We can also find this value near the bottom of the output of summary(m):
summary(m)
One problem with the R2 value is, that it always increases when another input variable is added to
the model. Thus, the R2 value cannot be used to compare model fit for models with different numbers
of variables. An attempt to compensate for this effect is the adjusted R2 value:
The adjusted R2 value for the stackloss dataset can be found as follows:
n <- nrow(stackloss)
p <- ncol(stackloss) - 1
[Link] <- 1 - (n - 1) / (n - p - 1) * (1 - [Link])
[Link]
[1] 0.8983257
We can also find this value near the bottom of the output of summary(m), listed as Adjusted
R-squared.
using the ANOVA identity notation for the residual sum of squares SSres .
Recall, that the adjusted R2 can be written as
2 SSres /(n − p − 1)
Radj =1− ,
SStotal /(n − 1)
10.6. MEASURES FOR MODEL FIT 155
SStotal
or, if we introduce s2Y = n−1 , as
2 σ̂ 2
Radj =1− .
s2Y
[1] 10.51556
The
√ output means the typical deviation of actual stack loss from predicted stack loss is about
10.52 ≈ 3.24 units. If we are to choose among the alternative models, we prefer the models
with the smaller value.
We can also find this value (actually, the standard deviation σ) near the bottom of the output of
summary(m), listed as Residual standard error.
n
(i) 2
X
PRESS := yi − ŷi .
i=1
(i)
In the comparison, ŷi is used instead of ŷi so that in models with too many parameters, overfitting
does not result in overoptimistic estimates. ”Good” models have small PRESS.
The following lemma helps to compute PRESS without having to fit n separate models, which would
be expensive (we take it here without a proof).
[1] 280.2634
We could use this number to compare PRESS values across different models to select the one
with best predictive performance.
Notice, PRESS is not directly available in the summary() of an lm model, but can be easily computed
as shown above.
AIC alone has no absolute meaning - it’s used for comparing models. The AIC() function in R
automatically computes this value from an lm object.
Definition 10.7.1. Two linear regression models are considered nested if one model (the reduced
or restricted model) is a special case of the other (the full model). This means that all the predictor
variables in the reduced model are also contained in the full model. The full model contains the reduced
model’s predictors plus at least one additional predictor.
The reduced model is nested within the full model because we can obtain it from the full model
by imposing the restrictions β2 = 0 and β3 = 0.
The primary method for comparing nested models is the Partial F-Test. This test assesses whether
the additional predictors in the full model provide a statistically significant improvement in explaining
the variation in the response variable.
Null Hypothesis (H0 ): The coefficients for the extra predictors in the full model are all zero. In
our example, H0 : β2 = β3 = 0.
Alternative Hypothesis (H1 ): At least one of the extra coefficients is not zero.
The test statistic for the Partial F-Test is based on the difference in the Residual Sum of Squares
(RSS) between the two models.
F statistic can be constructed in the following way
(Yi − Ŷ )2 is the residual sum of squares for the reduced model (with k + 1
P
where SSres (reduced) =
i
parameters), SSres (full) is the residual sum of squares for the full model (with p + 1 parameters, where
p > k) and n is the total number of observations.
Under the null hypothesis, this F-statistic follows an F-distribution with (p − k, n − p − 1) degrees
of freedom. A large F-value (and a corresponding small p-value) leads to the rejection of the null
hypothesis, providing evidence that the full model is superior.
For illustrating of the partial F-test we will use the MLBattend dataset.
158 CHAPTER 10. LINEAR REGRESSION
The [Major League Baseball attendance] data contains records from 1969 to 2000, including
variables such as attendance, year, runs scored, wins, and games behind. The dataset is built in
R package UsingR.
For the baseball dataset, we want to test if adding runs scored and games behind to a model
that already contains year and wins significantly improves the prediction of attendance.
Nested Models
We test hypothesis
H0 : β3 = β4 = 0, H1 : At least one of β3 , β4 6= 0.
3. AIC (Akaike Information Criterion) and BIC (Bayesian Information Criterion, not covered in
this module)
Now we describe algorithms that use these criteria to select models efficiently.
Forward selection starts with the intercept-only model y = β0 + ε, then iteratively adds variables. At
each step, for each variable not yet in the model, we compute the t−statistic for testing whether that
variable’s coefficient is zero. We add the variable (if any) with the largest |t|− value provided this is
above a chosen threshold (commonly |t| > 2), corresponding to approximate 5
Backward selection starts with the full model containing all p predictors, then iteratively removes
variables. At each step, we compute the t− statistic Tj for testing H0 : βj = 0 for each variable currently
in the model. We remove the variable with the smallest |Tj |, provided this is below a chosen threshold
(commonly |t| < 2). The process continues until all remaining variables have |t|− values above the
threshold.
Bidirectional Selection
Bidirectional (or hybrid) methods combine forward and backward selection. At each step, we consider
both adding the most significant variable not in the model and removing the least significant variable
currently in the model. Typically, different significance thresholds are used for entry and removal to
avoid cycling.
Note
Remark. While computationally efficient, stepwise procedures have important limitations: the
resulting test statistics and confidence intervals are biased (since variable selection and testing
use the same data), the methods are prone to overfitting, and small changes in the data can lead
to very different selected models. These procedures should be used primarily for exploratory
analysis, and any final model should ideally be validated on independent data.
We will illustrate the approaches on the MLBattend dataset discussed in the previous example.
library(UsingR)
data(MLBattend)
# Fit different models
model1 <- lm(attendance ~ wins, data = MLBattend)
model2 <- lm(attendance ~ wins + [Link], data = MLBattend)
model3 <- lm(attendance ~ wins + [Link] + [Link], data = MLBattend)
model4 <- lm(attendance ~ wins + [Link] + [Link] + year, data = MLBattend)
# Create a comparison table
model metrics <- [Link](
Model = c(" wins", " wins + runs", " wins + runs + behind", "Full"),
R squared = c(summary(model1)$[Link], summary(model2)$[Link],
summary(model3)$[Link], summary(model4)$[Link]),
Adj R squared = c(summary(model1)$[Link], summary(model2)$[Link],
summary(model3)$[Link], summary(model4)$[Link])
)
print(model metrics)
Interpretation: Based on the R2 and Adjusted R2 values, each additional predictor improves
the model’s ability to explain variance in MLB attendance, but with diminishing returns. The
largest improvement comes from adding [Link] to the model with wins (a 6.8% point
increase in Adjusted R2 ). Adding [Link] provides a further 3.1% point gain, while
including year offers only a minor 1.5% point improvement. However, even the full model explains
only about 32% of the variance, leaving substantial variation in attendance unexplained.
Stepwise Forward Selection
We start with an empty model and let R add variables one by one based on significance.
# Start with the model containing only the intercept
null model <- lm(attendance ~ 1, data = MLBattend)
# Perform forward selection (using AIC as the criterion)
forward model <- step(null model,
scope = list(lower = null model, upper = model4),
direction = "forward", trace = 0) # trace=0 hides detailed output
# See the final model selected
summary(forward model)
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -7.412e+06 4.962e+05 -14.940 < 2e-16 ***
[Link] 3.169e+03 4.005e+02 7.914 1.23e-14 ***
[Link] -1.165e+04 1.462e+03 -7.969 7.76e-15 ***
year 3.721e+03 2.476e+02 15.026 < 2e-16 ***
---
Signif. codes: 0 ’***’ 0.001 ’**’ 0.01 ’*’ 0.05 ’.’ 0.1 ’ ’ 1
Residual standard error: 75580 on 693 degrees of freedom
Multiple R-squared: 0.3233, Adjusted R-squared: 0.3203
F-statistic: 110.3 on 3 and 693 DF, p-value: < 2.2e-16
Interpretation: The stepwise forward selection algorithm selected a final model with
[Link], [Link], and year as predictors, notably excluding wins. All three included
variables are highly statistically significant (p < 0.001). The model explains approximately 32.3%
of the variance in MLB attendance. The positive coefficient for year indicates a general upward
10.7. MODEL SELECTION 161
trend in attendance over time, while the positive coefficient for [Link] suggests that higher-
scoring teams attract larger crowds. Interestingly, the negative coefficient for [Link]
shows that teams further from first place have lower attendance, even after accounting for their
offensive production.
Stepwise Backward Selection
We start with the full model and let R remove the least significant variables.
# Perform backward selection on the full model
backward model <- step(model4, direction = "backward", trace = 0)
# See the final model selected
summary(backward model)
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -7.412e+06 4.962e+05 -14.940 < 2e-16 ***
[Link] 3.169e+03 4.005e+02 7.914 1.23e-14 ***
[Link] -1.165e+04 1.462e+03 -7.969 7.76e-15 ***
year 3.721e+03 2.476e+02 15.026 < 2e-16 ***
---
Signif. codes: 0 ’***’ 0.001 ’**’ 0.01 ’*’ 0.05 ’.’ 0.1 ’ ’ 1
Residual standard error: 75580 on 693 degrees of freedom
Multiple R-squared: 0.3233, Adjusted R-squared: 0.3203
F-statistic: 110.3 on 3 and 693 DF, p-value: < 2.2e-16
Interpretation: The stepwise backward selection converged on the exact same model as the
forward selection: attendance [Link] + [Link] + year. This consistency is
reassuring and suggests that this is plausible model structure for this dataset. The algorithm,
starting from the full model that included wins, systematically removed it, indicating that wins
does not provide statistically significant explanatory power for attendance beyond what is already
captured by the other three variables- likely reflecting its correlation with [Link]. This
final model explains 32.3% of the variance in attendance with three highly significant predictors.
Bidirectional (Stepwise) Selection
This is often the most robust automatic method, as it checks at each step if any variable should
be added or removed.
# Perform stepwise (bidirectional) selection
stepwise model <- step(null model,
scope = list(lower = null model, upper = model4),
direction = "both", trace = 0)
# See the final model selected
summary(stepwise model)
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -7.412e+06 4.962e+05 -14.940 < 2e-16 ***
[Link] 3.169e+03 4.005e+02 7.914 1.23e-14 ***
[Link] -1.165e+04 1.462e+03 -7.969 7.76e-15 ***
year 3.721e+03 2.476e+02 15.026 < 2e-16 ***
---
Signif. codes: 0 ’***’ 0.001 ’**’ 0.01 ’*’ 0.05 ’.’ 0.1 ’ ’ 1
Residual standard error: 75580 on 693 degrees of freedom
Multiple R-squared: 0.3233, Adjusted R-squared: 0.3203
F-statistic: 110.3 on 3 and 693 DF, p-value: < 2.2e-16
Interpretation: The bidirectional selection algorithm also converged on the identical model:
162 CHAPTER 10. LINEAR REGRESSION
10.8 Multicollinearity
Definition 10.8.1. In linear regression, Multicollinearity denotes the situation when the columns
of the design matrix X are (approximately) linearly dependent.
The columns of X are linearly dependent, if and only if there is a v ∈ Rp+1 , v 6= 0 such that Xv = 0.
The columns are approximately linearly dependent, if there is a non-zero vector v such that Xv ≈ 0.
y x1 x2
2 1 1
4 2 2
6 3 3
In this case, an exact match with no residuals can be achieved as yi = xi1 + xi2 , i.e. for β1 = β2 = 1.
But this solution is not unique: we could just as well write y = 2x1 or y = 2x2 . Any choice of β1
and β2 with β1 + β2 = 2 will lead to zero residuals.
10.8. MULTICOLLINEARITY 163
The problem in the example above occurs, since the two input columns in the data are identical.
The same problem, in a less obvious way, would occur in datasets with more inputs when one column of
X can be written as a linear combination of other columns. This is a problem of the given input data,
not of the responsees or of the statistical model.
In the case where there is only approximate linear dependency between the columns of X, the
inverse (X> X)−1 exists and an estimate β̂ can be computed, but there will be huge uncertainties in
some of the estimated coefficients. We illustrate this effect using a numerical example.
Example (Multicollinearity Affects CI)
Here we simulate two datasets with 100 observations each. In both cases, the true relationship is
y = 2x1 + 3x2 + ε, where ε is small random noise.
The key difference lies in the predictors. In the multicollinear case, x2 is generated to be
almost identical to x1 , making them highly correlated. In the orthogonal case, x1 and x2 are
independent standard normal variables with zero correlation. We then fit linear regression models
to both datasets and examine their 95% confidence intervals.
# Compare multicollinear vs orthogonal predictors
[Link](20211101)
n <- 100
# Case 1: Multicollinear predictors (highly correlated)
x1 multi <- rnorm(n, 1:n, 0.1)
x2 multi <- rnorm(n, x1 multi, 0.1) # Almost identical to x1 multi
y multi <- 2*x1 multi + 3*x2 multi + rnorm(n, 0, 0.5)
# Case 2: Orthogonal predictors (uncorrelated)
x1 ortho <- rnorm(n)
x2 ortho <- rnorm(n)
y ortho <- 2*x1 ortho + 3*x2 ortho + rnorm(n, 0, 0.5)
# Fit models
m multi <- lm(y multi ~ x1 multi + x2 multi)
m ortho <- lm(y ortho ~ x1 ortho + x2 ortho)
# Compare confidence intervals
print(confint(m multi))
print(confint(m ortho))
2.5 % 97.5 %
(Intercept) -13.3354478 12.745565
x1 multi -6.1371039 7.148548
x2 multi -1.6203284 9.248167
2.5 % 97.5 %
(Intercept) -0.1037912 0.08872902
x1 ortho 1.9582176 2.03828037
x2 ortho 2.9677943 3.05317868
We can observe that for the multicollinear predictors, the intervals for both coefficients are
wide, indicating substantial uncertainty about the exact coefficient values. In contrast, for the
orthogonal predictors, the intervals are narrow and precise.
This demonstrates the effect of multicollinearity: when predictors are highly correlated, the
regression model struggles to distinguish their individual contributions.
1. In the case of exact multicollinearity, X> X is not invertible. If the columns of X are approximately
linearly dependent, then (X> X)−1 does exist, but small changes of X lead to large changes of
(X> X)−1 and numerical computation of the inverse is strongly affected by rounding errors.
2. If the columns of X are approximately linearly dependent, the computed value of the estimator
β̂ is strongly affected by small changes to the system: the noise in the model strongly affects β̂,
164 CHAPTER 10. LINEAR REGRESSION
leaving out a single observation may lead to large changes in β̂ and computation of β̂ is sensitive
to rounding errors.
While multicollinearity can make the regression coefficients ambiguous, the outputs y are not affected.
Predictions made using a model where multicollinearity is present are still reliable.
Note (Revision)
As a rule of thumb, if κ(X) < 10 there are no significant problems with multicollinearity, and if
κ(X) > 30 the regression problem suffers severe problems with multicollinearity. In the case where
the columns of X are exactly linearly dependent, we have σmin (X) = 0 and κ(X) = ∞.
In R, the condition number can be computed using the function kappa(). The function can either
be called as kappa(X, exact = TRUE), where X is the design matrix, or as kappa(m, exact = TRUE)
where m is the object returned by lm(). If the optional argument exact = TRUE is omitted, only an
approximate result is returned (using a faster algorithm).
Example (Multicollinearity)
To inspect these data, we use a pair scatter plot of the input variables, leaving out the responses
in column 1:
pairs(data[,-1])
10.8. MULTICOLLINEARITY 165
While there are visible patterns, there are no clear signs of a linear dependency between any two
columns. To see whether there are problems with multicollinearity, we fit a linear model and
determine the condition number.
m <- lm(y ~ ., data = data)
kappa(m, exact = TRUE)
[1] 148.7497
Since the condition number is larger than 30, we can say that the data shows severe multicollinear-
ity.
Note
To find out which variables are involved, we can use singular value decomposition. The
exact procedure, however, is not covered in this module.
10.8.3 Mitigations
There are various ways problems resulting from multicollinearity can be addressed:
• Sometimes columns are linearly dependent, because a redundant input is present. Removing one
of the input variables in a group of linearly dependent inputs can solve this problem.
• The inputs can be transformed by choosing new variables as functions of the original inputs. For
example, in example 10.8.1 one could try the input variables x̃1 = (x1 +x2 )/2 and x̃2 = (x1 −x2 )/2.
166 CHAPTER 10. LINEAR REGRESSION
• If n is small, the problem may be resolved by getting more data. In particular, for n < p + 1, we
always have strict multicollinearity, so in this case we definitely need more data.
• Alternative estimation approaches may be used. The basic idea is to sacrifice the requirement
that β̂ is unbiased, in exchange for a large variance reduction. Example of this approach is ridge
regression. It is not covered in this module.
10.9 Examples
10.9.1 Stackloss data
Fitting the model
Below is the code that shows how to fit in the model [Link] = [Link] + [Link] +
[Link]..
m <- lm ( stack . loss ~ . , data = stackloss )
summary ( m )
# Coefficients :
# Estimate Std . Error t value Pr ( >| t |)
# ( Intercept ) -39.9197 11.8960 -3.355 0.00375 * *
# Air . Flow 0.7156 0.1349 5.307 5.8 e -05 * * *
# Water . Temp 1.2953 0.3680 3.520 0.00263 * *
# Acid . Conc . -0.1521 0.1563 -0.973 0.34405
# ---
# Signif . codes : 0 * * * 0.001 * * 0.01 * 0.05 0.1 1
#
# Residual standard error : 3.243 on 17 degrees of freedom
# Multiple R - squared : 0.9136 , Adjusted R - squared : 0.8983
# F - statistic : 59.9 on 3 and 17 DF , p - value : 3.016 e -09
# [1] 10.51556
The square root of this number, so the estimated standard deviation of the εi is shown as Residual
standard error in the summary output above. We check that we get the same result:
10.9. EXAMPLES 167
# [1] 3.242771
This result is also listed as the Residual standard error near the bottom of the summary(m)
output, above, or by
sigma ( m )
These values are also listed in the Std. Error column of the summary(m) output, or, by entering
summary ( m ) $ coefficients [ , " Std . Error " ]
Hypothesis tests
Let us now test the hypothesis H0 : βi = 0. The test statistic for this case is the following (recall 10.4
for the case b = 0):
T <- coef ( m ) / se
T
These values are also listed in the t value column of the summary(m) output, or, by entering
summary ( m ) $ coefficients [ , " t value " ]
Before we can perform the test, we need to choose α and to find the corresponding critical value:
alpha <- 0.05
t <- qt (1 - alpha / 2 , n - p - 1)
t
# [1] 2.109816
Note
Notice: In mathematical notation, we use tn−p−1 (α/2) to represent the critical value where
the right-tail probability is α/2 (see Lemma 10.4.3). However, R’s qt() function calculates the
quantile where the cumulative probability to the left is the given value. We use qt(1 - alpha/2,
n-p-1) or qt(alpha/2, n-p-1, [Link] = FALSE)
Using the critical value t we can decided whether H0 should be accepted or rejected. For example,
looking at the intercept β0 , we find |T0 | = | − 3.3557234| > 2.109816 = tn−p−1 (α/2) and thus we can
reject the hypothesis H0 : β0 = 0. This means that the intercept is significantly different from 0.
168 CHAPTER 10. LINEAR REGRESSION
Confidence Intervals
Using the quantile t we can also get confidence intervals.√Here we only show the confidence interval for
the intercept β0 (recall Lemma 10.4.2 and that we call σ 2 Cii ≡ SE(β̂i ) the standard error of β̂i ):
c ( coef ( m ) [1] - se [1] * t , coef ( m ) [1] + se [1] * t )
# 2.5 % 97.5 %
# ( Intercept ) -65.01893 -14.82047
# Air . Flow 0.43185 0.99943
# Water . Temp 0.51908 2.07158
# Acid . Conc . -0.48205 0.17779
# Coefficients :
# Estimate Std . Error t value Pr ( >| t |)
# ( Intercept ) 3.25377 0.15014 21.672 < 2e -16 ***
# TPSA -0.00471 0.00107 -4.408 1.12 e -05 ***
# SAacc 0.01524 0.00377 4.041 5.64 e -05 ***
# H050 0.22060 0.05238 4.211 2.72 e -05 ***
# MLOGP 0.18639 0.02644 7.050 2.54 e -12 ***
# RDCHI 1.29531 0.36796 3.520 0.000449 ***
# GATS1p -0.57131 0.15388 -3.713 0.000214 ***
# nN 0.14604 0.07190 2.031 0.042467 *
# C040 -0.09101 0.05749 -1.583 0.113710
# ---
# Signif . codes : 0 * * * 0.001 * * 0.01 * 0.05 0.1 1
#
# Residual standard error : 0.833 on 986 degrees of freedom
# Multiple R - squared : 0.1457 , Adjusted R - squared : 0.1387
# F - statistic : 21.01 on 8 and 986 DF , p - value : < 2.2 e -16
10.9. EXAMPLES 169
Confidence interval
Let us find a 95% confidence interval for GATS1p.
We can obtain the result using only the quantile t and numbers shows in the summary(m) output.
We focus on the row corresponding to GATS1p in the output shown above:
GATS1p -0.571313 0.153882 -3.713 0.000227 ***
√The first number in this row is the centre of the confidence interval, the second number corresponds
to σ̂ 2 Cii :
alpha <- 0.05
n <- nrow ( stackloss )
p <- ncol ( stackloss ) - 1
t <- qt (1 - alpha / 2 , n - p - 1)
cat ( " [ " , -0.571313 - 0.153882 *t , " , " , -0.571313 + 0.153882 *t , " ]\ n " , sep = " " )
# [ -0.8960051 , -0.2466209]
# Coefficients :
# Estimate Std . Error t value Pr ( >| t |)
# ( Intercept ) 1.7970 0.3252 5.527 0.0004 * *
# x1 0.4921 0.4613 1.067 0.3144
# x2 -0.4272 0.3040 -1.405 0.1923
# x3 -0.2270 0.3077 -0.738 0.4783
# x4 -0.1608 0.3158 -0.509 0.6224
# ---
# Signif . codes : 0 * * * 0.001 * * 0.01 * 0.05 0.1 1
#
170 CHAPTER 10. LINEAR REGRESSION
From the R output it is clear that the intercept is non-zero, but none of the remaining regression
coefficients are significantly different from zero when testing at 5% level. This poses the question
whether the inputs have any effect on the output at all. To answer this question, we test the hypothesis
H0 : β1 = β2 = β3 = β4 = 0,
omitting only β0 in the hypothesis. We use the method from section 10.4.5 to perform this test.
The information required to perform this F test is contained in the last row of the summary(m)
output:
F-statistic: 4.007 on 4 and 9 DF, p-value: 0.03891
Here we can see that the test statistic is 4.007 (a rounded version of the value we found above)
and the degreen of freedom for the F -test are k = 4 and n − p − 1 = 9. The ”p-value” listed is the
smallest significance level at which the test still rejects H0 . We see that H0 is rejected at significance
level 0.03891, and thus also at α = 0.05.
Interlude: Understanding the lm()
Output
We have seen that the lm() function returns an object which contains a lot of information about the
fitted model. As you have seen before, we can inspect this information by using the summary() function.
The aim of this section is to summarise how the output of summary() relates to the mathematical
expressions we have considered in the previous sections.
We will use the stackloss dataset as an example:
m <- lm ( stack . loss ~ . , data = stackloss )
summary ( m )
Listing 10.1: Fitting a linear model to the stackloss dataset
Figure 10.9: Output of summary() for a linear model fitted to the stackloss dataset. Different parts of
the output are marked with letters A to I.
171
172 CHAPTER 10. LINEAR REGRESSION
Here I have marked different parts of the output using red shaded boxes and using the letters A to I.
We discuss each of the marked sections in turn:
• The first section, part A, contains summary statistics for the fitted residuals ε̂1 , . . . , ε̂n . The values
shown are the minimum, first quartile, median, third quartile, and maximum of the residuals.
This is the same information we can get using the command ‘summary(resid(m))‘. The mean is
omitted from A, since it always equals zero.
• Column D shows the t-test statistic for the coefficients. The values are computed using equation
10.4.
• Column E replicates the information from column D in different form, showing p-values instead
of the test statistics.
• The field F shows the estimated standard deviation σ̂. This is computed as the square root of σ̂ 2
from equation 10.3.
• Field G shows the value n − p − 1. This is the number of degrees of freedom in lemma 10.4.1.
The value is needed when performing hypothesis tests and computing confidence intervals for
individual coefficients.
• Field H shows the R2 value and adjusted R2 value. These are computed using definitions 10.6.1
and 10.6.2.
• Field I shows the F -test statistic for testing the hypothesis H0 : β1 = · · · = βp = 0 (omitting the
coefficient β0 for the intercept). This value can be used to test the hypothesis that the inputs
have no effect on the output. The degrees of freedom shown are p and n − p − 1 from lemma.
Finally, below is the table that summarises the different ways to find inference:
10.9. EXAMPLES 173
Table 10.1: Summary of R commands for extracting inference statistics from linear model objects.
174 CHAPTER 10. LINEAR REGRESSION
Chapter 11
One-way ANOVA
The overall F-test in multiple linear regression has an important special case: when all predictors
are categorical, we obtain the Analysis of Variance (ANOVA) framework. While regression
typically concerns continuous predictors, ANOVA addresses whether group means differ significantly.
Mathematically, one-way ANOVA is a special instance of the multiple linear regression model where the
design matrix X consists solely of dummy variables encoding group membership.
Example (Drug tests)
1. Group 1: Placebo
k
X
Yij = µ + τj + εij , i = 1, . . . , nj , j = 1, . . . , k, s.t.c. nj τj = 0. (11.1)
j=1
where Yij is the continuous response for observation i in group j, µ is the overall population mean,
τj is the treatment effect for group j and εij is the random error term.
The following assumptions are met
2. εij ∼ N (0, σ 2 )
175
176 CHAPTER 11. ONE-WAY ANOVA
H0 : τ1 = τ2 = · · · = τk = 0 vs H1 : τj 6= 0 for some j
This tests whether all group means are equal versus at least one group mean differs.
Example (Drug test, ANOVA context)
In the context of the example above, µ + τ1 is the mean reduction for Placebo, µ + τ2 is the mean
reduction for the low does and µ + τ3 is the mean reduction of high dose.
Under ANOVA model we test
Note
Note: ANOVA can be extended to more complex designs, such as two-way ANOVA, which
considers two factors simultaneously (e.g., drug dose and gender). In this module, we focus
exclusively on the one-way ANOVA case.
i.i.d.
Lemma 11.2.1 (Sample Variance from a Normal Distribution). If X1 , . . . , Xm ∼ N (µ, σ 2 ),
then:
1 Pm
1. X and s2 = m−1 i=1 (Xi − X)2 are independent
(m−1)s2
2. σ2
∼ χ2m−1
k
X
Uj ∼ χ2ν1 +ν2 +···+νk .
j=1
Finally, recall the ANOVA identity introduced earlier (see Lemma 10.4.4). For one-way ANOVA
problem we, therefore, can write:
SStotal = SSreg + SSres .
However, the components are calculated now using the information about the group, so:
Pn
• SStotal = i=1 (Yi − Y •• )2 is the total sum of squares
11.2. DISTRIBUTIONS OF SUMS OF SQUARES 177
Pk
• SSreg = j=1 nj (Y j• − Y •• )2 is the regression (treatment) sum of squares
Pk Pnj
• SSres = j=1 i=1 (Yij − Y j• )2 is the residual (error) sum of squares
In these lectures, we use regression notation for ANOVA, where SST otal = SSReg + SSRes . This
is equivalent to the traditional one-way ANOVA notation SST otal = SSBetween + SSW ithin , where
SSBetween measures variability between groups and SSW ithin measures variability within groups. In
our regression framework, SSReg (variation explained by the model) corresponds to SSBetween , and
SSRes (unexplained variation) corresponds to SSW ithin . We adopt this notation for the sake of being
consistent across the module.
The identity and the lemmas above will be used to prove the following statement.
Lemma 11.2.3 (Distribution and Expectation of the Error Sum of Squares). For the model
11.1, the residual sum of squares SSres follows a scaled chi-square distribution:
SSres
∼ χ2n−k ,
σ2
Pk
where n = j=1 nj . Furthermore, the Mean Square Error (MSE) M Sres is an unbiased estimator for
σ2:
SSres
E[M Sres ] = E = σ2.
n−k
Proof
Proof. For a fixed group j, the model implies that observations Y1j , . . . , Ynj j form a random
sample from N (µj , σ 2 ), where µj = µ + τj . Applying lemma 11.2.1 to this group, the scaled sum
of squared deviations from the group mean follows a chi-square distribution:
Pnj
(Yij − Y j• )2
Uj = i=1 ∼ χ2nj −1 .
σ2
The total residual sum of squares is the sum of these within-group variations across all k groups:
nj
k X
X
SSres = (Yij − Y j• )2 .
j=1 i=1
Scaling by σ 2 gives:
k
SSres X
= Uj ,
σ2
j=1
where Uj ∼ χ2nj −1 .
Since the k groups are independent, the random variables U1 , . . . , Uk are independent. Applying
lemma 11.2.2, their sum is also chi-square, with degrees of freedom equal to the sum of the
individual degrees of freedom:
k
SSres X
∼ χ2ν , ν= (nj − 1) = n − k.
σ2
j=1
SSres
Thus, σ2
∼ χ2n−k .
178 CHAPTER 11. ONE-WAY ANOVA
Rearranging:
E[SSres ] = (n − k)σ 2 .
Consequently:
SSres
E[M Sres ] = E = σ2,
n−k
which completes the proof.
Lemma 11.2.4 (Distribution and Expectation of the Regression Sum of Squares). For the
model 11.1, under the null hypothesis H0 : τ1 = τ2 = · · · = τk = 0, the scaled regression sum of squares
follows a chi-square distribution:
SSreg
∼ χ2k−1 .
σ2
Furthermore, its expectation is:
k
X
2
E[SSreg ] = (k − 1)σ + nj τj2 ,
j=1
which implies:
SSreg
1. E = k − 1 under H0
σ2
SSreg
2. E = σ 2 under H0
k−1
Proof
Proof.
Part 1: Expectation of SSreg
We use the computational formula:
k
X
SSreg = nj (Y j• − Y •• )2 .
j=1
σ2 σ2
E[Y j• ] = µ + τj , Var(Y j• ) = , E[Y •• ] = µ, Var(Y •• ) = .
nj n
11.2. DISTRIBUTIONS OF SUMS OF SQUARES 179
k
1 X
Cov(Y j• , Y •• ) = nm Cov(Y j• , Y m• ).
n
m=1
σ2
Groups are independent, so Cov(Y j• , Y m• ) = 0 for m 6= j, and equals Var(Y j• ) = nj when
m = j. Thus:
nj · σ 2 /nj σ2
Cov(Y j• , Y •• ) = = .
n n
Substituting back:
σ2 σ2
E[(Y j• − Y •• )2 ] = τj2 + − .
nj n
Therefore:
k
σ2 σ2
X
E[SSreg ] = nj τj2 + −
nj n
j=1
k k k
X X X nj σ 2
= nj τj2 + σ2 −
n
j=1 j=1 j=1
k
X
= nj τj2 + kσ 2 − σ 2
j=1
k
X
= (k − 1)σ 2 + nj τj2 .
j=1
Optional
Below is the sketch of the prove. Some statements will be discussed in the exercises.
σ2
Y j• ∼ N µ + τj , .
nj
√
2. Define Xj = nj Y j• for j = 1, . . . , k. Then the {Xj } are independent and, under H0
(τj = 0),
√
nj µ, σ 2 .
Xj ∼ N
and for ` = 2, . . . , k,
k
X
Z` = a`j Xj ,
j=1
k
X
a`j amj = δ`m , `, m = 1, . . . , k,
j=1
√ √
with a1j = nj / n (the coefficients for Z1 ).
5. Since the {Z` } are linear combinations of independent normals, they are jointly normal.
Uncorrelatedness implies independence. Hence, under H0 ,
√
Z1 ∼ N ( n µ, σ 2 ), Z` ∼ N (0, σ 2 ) for ` = 2, . . . , k.
(The term Z12 corresponds to the overall mean Y •• and is not part of the regression
sum of squares.)
7. Therefore, under H0 ,
k 2
SSreg X Z`
= ∼ χ2k−1 ,
σ2 σ
`=2
Yij = µ + τj + εij , i = 1, . . . , nj , j = 1, . . . , k,
i.i.d.
with the constraint kj=1 nj τj = 0, and where εij ∼ N (0, σ 2 ).
P
The parameters {µ, τ1 , . . . , τk } are estimated using the method of least squares. This involves
finding the values that minimize the sum of squared errors:
nj
k X nj
k X
X X
S(µ, τ1 , . . . , τk ) = ε2ij = (Yij − µ − τj )2 .
j=1 i=1 j=1 i=1
To find the minimizers, we take partial derivatives and set them to zero.
Differentiate with respect to µ:
k nj
∂S XX
= −2 (Yij − µ − τj ).
∂µ
j=1 i=1
Y•• − N µ̂ = 0 ⇒ µ̂ = Y •• . (11.2)
Substituting µ̂ = Y •• , we obtain:
τ̂j = Y j• − Y •• . (11.3)
The effect for group j is estimated by the difference between the group mean and the overall mean.
P The individual parameter estimates µ̂ and τ̂j depend on the specific constraint applied (e.g.,
nj τj = 0 or τ1 = 0). However, scientifically meaningful quantities derived from them are often
estimable functions—their estimates are unique and do not depend on the constraint.
Important examples include:
182 CHAPTER 11. ONE-WAY ANOVA
+ τj = µ̂ + τ̂j = Y •• + (Y j• − Y •• ) = Y j• .
µ\
− τm = τ̂j − τ̂m = (Y j• − Y •• ) − (Y m• − Y •• ) = Y j• − Y m• .
τj\
These estimates are intuitive and remain the same regardless of the coding scheme used for the factor
levels.
The error variance σ 2 is estimated using the Mean Square Error (MSE), which is the residual sum
of squares scaled by its degrees of freedom. As established in Lemma 11.2.3, this provides an unbiased
estimator:
SSres
σ̂ 2 = M Sres = , where E[M Sres ] = σ 2 .
n−k
2. Within each group j, Y j• and within-group variance s2j are independent (by lemma 11.2.1)
Different types of display panel for a chemical production process were tested to see how quickly
simulated emergencies could be brought under control. The time (in seconds) for three panel
types were:
Panel type A B C
21 24 40
27 21 36
24 18 35
26 19 32
X
2
yij =212 + 272 + 242 + 262 + 242 + 212 + 182 + 192
ij
X 2
y••
2
SStotal = yij −
N
ij
3232 104329
= 9269 − = 9269 −
12 12
= 9269 − 8694.08 = 574.92
k 2 2
X yj• y••
SSreg = −
nj N
j=1
982
+ 822 + 1432 3232
= −
4 12
9604 + 6724 + 20449
= − 8694.08
4
36777
= − 8694.08 = 9194.25 − 8694.08 = 500.17
4
184 CHAPTER 11. ONE-WAY ANOVA
df(SSreg ) = k − 1 = 3 − 1 = 2
df(SSres ) = n − k = 12 − 3 = 9
df(SStotal ) = n − 1 = 12 − 1 = 11
SSreg 500.17
M Sreg = = = 250.085
df(SSreg ) 2
SSres 74.75
M Sres = = = 8.306
df(SSres ) 9
Finally, F-statistic is
M Sreg 250.085
F = = = 30.11
M Sres 8.306
It is convenient present the results in ANOVA Table
Source df SS MS F p-value
Group (regression) 2 500.17 250.085 30.11 < 0.001
Residual 9 74.75 8.306
Total, SStotal 11 574.92
The calculated F-statistic of 30.11 follows an F-distribution with (2, 9) degrees of freedom under
the null hypothesis. Comparing this to critical value F0.05 (2, 9) = 4.26, we have F = 30.11 > 4.26.
We also have p < 0.001. We reject the null hypothesis that all panel types have the same mean
response time. There is strong evidence that at least one display panel type differs significantly
from the others in terms of emergency response time. The data suggest that Panel C (mean =
35.75s) is significantly slower than Panels A (24.5s) and B (20.5s) in bringing emergencies under
control.
Since F = 30.11 > 16.39, we have p < 0.001. We reject the null hypothesis that all panel types
have the same mean response time. There is strong evidence that at least one display panel type
differs significantly from the others in terms of emergency response time.
The data suggest that Panel C (mean = 35.75s) is significantly slower than Panels A (24.5s) and
B (20.5s) in bringing emergencies under control.
11.4. F-TEST AND ITS DISTRIBUTION 185
where β0 represents the mean of the reference group (group k) and βj represents the difference
between group j.
The ANOVA null hypothesis H0 : τ1 = · · · = τk = 0 becomes H0 : β1 = β2 = · · · = βk−1 = 0
Example (Drug test, MLR view)
Now in the context of linear regression, we create dummy variables: D1 = 1 if low dose and
D1 = 0 if otherwise and D2 = 1 if high does and D2 = 0 otherwise. Reference group is placebo.
The model becomes
Yi = β0 + β1 D1,i + β2 D2,i + εi ,
where β represents the mean reduction for the placebo group, β1 represents the difference (low
does-placebo) and β2 represents different (high dose-placebo).
We test
H0 : β1 = β2 = 0 vs H1 : At least one βj 6= 0.
SSR /2
We construct F statistic F = SSE /42 ∼ F (2, 42) under H0 .
SSreg df = k − 1 df = k − 1
Pk Pnj Pn
Unexplained Variation SSwithin = j=1 i=1 (Yij − Y j• )2 i=1 (Yi − Ŷi )
2
SSres df = n − k df = n − k
SS /(k−1)
Test Statistic F = SS between /(k−1)
SSwithin /(n−k)
reg
F = SSres /(n−k)
∼ F (k − 1, n − k) ∼ F (k − 1, n − k)
Null Hypothesis H0 : τ1 = τ2 = · · · = τk = 0 H0 : β1 = β2 = · · · = βk−1 = 0
(All group effects are zero) (All dummy coefficients are zero)
186 CHAPTER 11. ONE-WAY ANOVA
After fitting an ANOVA model, it’s crucial to examine the residuals to check model assumptions. The
residuals are defined as:
1. Residuals vs. Fitted Values: Checks for constant variance and linearity assumptions
While other diagnostic tools are available (Scale-Location plots for homoscedasticity, Residuals vs.
Order for time-dependent patterns, and Leverage plots with Cook’s distance for identifying influential
observations), the primary focus in this module will be on the first two plots for fundamental assumption
checking.
A plot of the residuals ε̂ij against the fitted values ŷij is a fundamental diagnostic tool. In the one-way
ANOVA model 11.1 the fitted value for an observation in group j is simply the group mean ŷij = y j• .
Therefore, all observations from the same group have the same fitted value on this plot.
If the modelling assumptions are correct, the errors εij are independent, normally distributed with
mean zero, and constant variance σ 2 . The residuals are estimates of these errors and should reflect
these properties. A key assumption to check is homoscedasticity—that the error variance is the same
across all groups.
To demonstrate a residual plot where the ANOVA assumptions are satisfied, we simulate data
for three groups with equal means and constant variance.
[Link](123) # For reproducibility
k <- 3 # Number of groups
n j <- 50 # Observations per group
n <- k * n j
# Simulate data: same mean (10), same variance (1) for all groups
group <- factor(rep(1:k, each = n j))
y <- 10 + rnorm(n)
# Perform one-way ANOVA
m <- aov(y group)
plot(m, which = 1, pch=19, main = "ANOVA Assumptions Satisfied\n(Constant Variance)")
11.5. ANOVA ASSUMPTIONS 187
We can see that the residuals are randomly scattered around zero, and the vertical spread is
approximately the same for all three groups (which correspond to the three distinct fitted values).
This pattern suggests the constant variance assumption is reasonable.
Now we consider a case where the constant variance assumption is violated. We simulate data
where the group means are the same, but the variances differ.
[Link](456)
k <- 3
n j <- 50
n <- k * n j
group <- factor(rep(1:k, each = n j))
# Simulate data with different variances per group
y <- numeric(n)
y[group == 1] <- 10 + rnorm(n j, sd = 0.2) # Low variance
y[group == 2] <- 10 + rnorm(n j, sd = 1.0) # Medium variance
y[group == 3] <- 10 + rnorm(n j, sd = 5.0) # High variance
m <- aov(y group)
plot(m, which = 1, pch=19, main = "ANOVA Assumption Violated\n(Non-Constant Variance)")
188 CHAPTER 11. ONE-WAY ANOVA
The “fanning” pattern is clear: the spread of the residuals is different for different groups. In this
case, the group with the highest fitted value also has the largest variance. This heteroscedasticity
violates a key ANOVA assumption.
The structure of the fitted values in a one-way ANOVA—clustered at the group means—simplifies
the interpretation of this plot. The primary question is whether the vertical spread of the residuals is
roughly equal across these clusters.
We can create a residual plot for the caloric intake example from Verzani
# Enter the data from Verzani’s Example 12.1
may <- c(2166, 1568, 2233, 1882, 2019)
sep <- c(2279, 2075, 2111, 2009, 1793)
dec <- c(2226, 2154, 2583, 2010, 2190)
# Create a data frame in long format
calories <- stack(list(may = may, sep = sep, dec = dec))
names(calories) <- c("values", "month")
# Perform one-way ANOVA
m cal <- aov(values month, data = calories)
# Create the base plot
plot(m cal, which = 1, pch = 19)
# Add group labels to each point
with(calories,
text(fitted(m cal), residuals(m cal),
labels = month,
pos = 4, # position: 1=below, 2=left, 3=above, 4=right
cex = 0.8, # character size
col = "blue"))
11.5. ANOVA ASSUMPTIONS 189
Based on the residual summaries, the spread of residuals for different groups looks fairly similar.
There is an indication that there may be problems with violating homogeneity of variances.
However, it appears to be not significant (checking using formal tests, such as Bartlett’s test.)
A Q-Q plot of the residuals can be used to assess whether the errors in an ANOVA model are normally
distributed. These plots compare the quantiles of the residuals against the quantiles of a normal
distribution. If the normality assumption is satisfied, the points should approximately follow a straight
line. For ANOVA models, we examine the pooled residuals across all groups.
We first show an example with simulated data where the ANOVA assumptions are known to be
satisfied.
[Link](123)
k <- 3
n j <- 50
group <- factor(rep(1:k, each = n j))
y <- 10 + rnorm(k * n j) # Normal errors
m <- aov(y group)
plot(m, which=2, pch = 19, main = "ANOVA Assumptions Satisfied\n(Normal Errors)")
190 CHAPTER 11. ONE-WAY ANOVA
As expected, the points lie approximately on a straight line, indicating that the normality
assumption is reasonable.
Next we consider an example where the errors follow an exponential distribution instead of a
normal distribution.
[Link](456)
k <- 3
n j <- 50
group <- factor(rep(1:k, each = n j))
eps <- rexp(k * n j) - 1 # Exponential errors with mean zero
y <- 10 + eps
m <- aov(y group)
plot(m, which=2, pch = 19, main ="ANOVA Assumptions Violated\n(Non-Normal Errors)")
11.5. ANOVA ASSUMPTIONS 191
The points show systematic curvature, particularly in the tails, indicating that the residuals are
not normally distributed. This violation could affect the validity of the F-test.
Now we examine the Q-Q plot for the caloric intake data from Verzani.
may <- c(2166, 1568, 2233, 1882, 2019)
sep <- c(2279, 2075, 2111, 2009, 1793)
dec <- c(2226, 2154, 2583, 2010, 2190)
calories <- stack(list(may = may, sep = sep, dec = dec))
names(calories) <- c("values", "month")
m cal <- aov(values month, data = calories)
plot(m cal, which=2, pch=19, main ="Q-Q Plot: Caloric Intake by Month")
# Add Shapiro-Wilk test for formal assessment
[Link](resid(m cal))
192 CHAPTER 11. ONE-WAY ANOVA
The Q-Q plot shows minor deviations at the tails, but these are comparable to what we would
expect from a true normal distribution with this sample size. The Shapiro-Wilk test (W = 0.977,
p = 0.945), given here without details, provides no evidence against the normality assumption.
This demonstrates an important point: with small samples (n=15 total), Q-Q plots often show
apparent ”deviations” that are actually within the range of normal sampling variability.
Let σj2 be the variance within treatment group j. Bartlett’s test evaluates:
Pk
j=1 (nj − 1)−1 − (n − k)−1
D =1+ .
3(k − 1)
Note
Bartlett’s test is sensitive to non-normality and should only be used after verifying the normality
assumption. A rule of thumb: if one Sj2 is more than 4 times another Sm 2 , this indicates potential
unequal variance.
Let us demonstrate how the test works in R for the case of Verzani Calories dataset. We can do
this since the Q-Q plot doesn’t show evidence for violation of normality.
[Link](values month, data = calories)$[Link]
Fligner-Killeen Test
This test is robust to non-normality. Let Ỹj be the median of observations in group j. The test examines
the ranks of absolute deviations dij = |Yij − Ỹj |. Under H0 , the ranks should be randomly distributed
across groups.
The Fligner-Killeen test is implemented in statistical software and is recommended when normality
is in doubt. In this module, we will use the R function [Link]().
1. Transformations: Apply transformations to the response variable (e.g., logarithm, square root,
reciprocal) to stabilize variance or improve normality.
Note
This exercise sheet covers material from Sections 10 and 11 of the course. Please
attempt all questions and prepare your solutions before the Week 11 workshop.
Your solutions will not be assessed. Before the workshop, you will be invited to vote for the
questions you would like to see discussed. The tutor will select the most popular questions for
discussion during the workshop. Complete solutions for all questions will be available afterwards
(by the end of the week). The solutions will be available only in HTML version of the notes.
Exercise 11.1
Let X be the design matrix of a model including an intercept, and let H = X(X > X)−1 X > ∈ Rn×n
be the hat matrix. Finally, let 1 = (1, 1, . . . , 1) ∈ Rn . Show that H1 = 1.
Exercise 11.2
For the stackloss dataset built into R, predict a value for [Link] when the inputs are
[Link] = 60, [Link] = 21 and [Link] = 87.
Exercise 11.3
Assume that X ∼ N (0, 1), Y ∼ χ2 (2) and Z ∼ χ2 (3) are independent. What are the distributions
of the following random variables?
Use the following rules: (1) the sum of squares of d independent, standard normally distributed
random variables follows a χ2 (d)-distribution. (2) If V1 ∼ χ2 (d1 ) and V2 ∼ χ2 (d2 ) are independent,
then (V1 /d1 )/(V2 /d2 ) ∼ F (d1 , d2 ).
• X2 + Y + Z
2X 2
• Y
Exercise 11.4
Consider the following R commands:
195
196 CHAPTER 11. ONE-WAY ANOVA
Either using this output, or using R to further inspect the built-in stackloss dataset, find a
99%-confidence interval for the parameter β[Link]. .
Exercise 11.5
Consider the following dataset. Our aim is to predict y from the variables x1 , x2 and x3 .
3. Perform a hypothesis test, at the 95%-level, for the hypothesis H0 : β2 = 0 with alternative
H1 : β2 6= 0.
Exercise 11.6
In this question we consider four different datasets, given by inputs xi and responses yi for
i ∈ {1, 2, 3, 4}.
1. Based on the following plots, discuss model fit of each model. Describe all relevant features
of the plots. Describe any problems with the model fit you find.
par(mfrow=c(2,2))
m1 <- lm(y1 ~ x1)
plot(fitted(m1), resid(m1), xlab="fitted values", ylab="residuals")
qqnorm(resid(m1), main=NULL)
m2 <- lm(y2 ~ x2)
plot(fitted(m2), resid(m2), xlab="fitted values", ylab="residuals")
qqnorm(resid(m2), main=NULL)
Exercise 11.7
The Prediction Error Sum of Squares (PRESS) is defined as
n
(i) 2
X
PRESS = yi − ŷi .
i=1
Using lemma 10.6.3, or otherwise, determine the PRESS value for the stackloss dataset built
into R (using [Link] as the response variable).
Exercise 11.8
This exercise refers to some parts of the proof of Lemma 11.2.4 (part II).
√
We construct Xj = nj Y j• for j = 1, . . . , k (see Lemma 11.2.4).
Define an orthonormal transformation from the basis {Xj } to {Z` } where:
k
1 √ √ √ X
Z1 = √ ( n1 X1 + n2 X2 + · · · + nk Xk ) , n= nj ,
n
j=1
and for ` = 2, . . . , k:
k
X
Z` = a`j Xj ,
j=1
where the coefficients a`j are chosen so that the transformation is orthonormal, i.e.,
11.5. ANOVA ASSUMPTIONS 199
k
X
a`j amj = δ`m , `, m = 1, . . . , k,
j=1
√ √
with a1j = nj / n.
Show that
(a) First sub-question:
Exercise 11.9
A study was conducted to investigate the effect of different diets on blood coagulation in animals.
Twenty-four animals were randomly assigned to one of four dietary regimes (A, B, C, D), and
their blood coagulation time (in seconds) was recorded.
Data: The observed coagulation times are as follows:
2. εij ∼ N (0, σ 2 )
1. Is there a statistically significant difference in the mean blood coagulation times among the
four diets?
2. Formally test the null hypothesis that all diet effects are zero against the alternative that
at least one is not
H0 : τA = τB = τC = τD = 0 H1 : At least one τj 6= 0.
Exercise 11.10
Using the input from the Exercise 11.9 estimate {µ, τA , τB , τC , τD } for the ANOVA model (11.1)
Exercise 11.11
For the data set Cars93 (MASS) perform a one-way analysis of variance of [Link] for each
level of DriveTrain. Does the data support the null hypothesis of equal population means?
Verify the ANOVA assumptions.
Exercise 11.12
A factory produces batches of dyestuff (chemical used in the dying of fabrics). There was concern
over the consistency of the strength of the dyestuff. To measure the strength of the dyestuff,
squares of fabric were dyed under controlled conditions. These squares were compared to a
standard by a technician, who scored them numerically (perfect match to the standard scored
100%).
This is a laborious process, so rather than sample every batch, t = 6 batches were selected at
random and studied. For each batch, a large number of squares of fabric were dyed and 6 of
them chosen at random. The resulting 36 samples were mixed to avoid bias and scored. The
percentage strengths were as follows (with group sizes and sums already calculated).
I II III IV V VI
94.5 89.0 88.5 100.0 91.5 98.5
93.0 90.0 93.5 99.0 93.0 100.0
91.0 92.5 93.5 100.0 90.0 98.0
89.0 88.5 88.0 98.0 92.5 100.0
96.5 91.5 92.5 95.0 89.0 96.5
88.0 91.5 91.5 97.5 91.0 98.0
ni 6 6 6 6 6 6
yi• 552.0 543.0 547.5 589.5 547.0 591.0
2 = 316019.
P
The other summary statistics we need are n = 36, y•• = 3370 and i,j yij
From the boxplots below, it looks like the strength of dyestuff could well be varying from batch
to batch.
11.5. ANOVA ASSUMPTIONS 201
Appendix
12.1 Vectors
We write v ∈ Rd if v = (v1 , . . . , vd ) for numbers v1 , . . . , vd ∈ R. We say that v is a d-dimensional vector,
and Rd is the d-dimensional Euclidean space. Vectors are often graphically represented as ”column
vectors”:
v1
v2
v = . .
..
vd
d
X
>
u v= ui vi . (12.1)
i=1
Note that the two vectors must have the same length for the inner product to exist.
Using this notation, the Euclidean length of a vector v can be written as
v
u d
uX √
kvk = t vi2 = v> v.
i=1
Vectors v1 , . . . , vn are said to be orthogonal, if vi> vj = 0 for all i 6= j. The vectors are said to be
orthonormal, if they are orthogonal and satisfy kvi k = 1 for all i.
12.2 Matrices
We write A ∈ Rm×n if
a1,1 ... a1,n
a2,1 ... a2,n
A= . .. ,
.. ..
. .
am,1 . . . am,n
where ai,j , sometimes also written as aij are numbers for i ∈ {1, . . . , m} and j ∈ {1, . . . , n}.
203
204 CHAPTER 12. APPENDIX
Transpose
If A ∈ Rm×n , then the transpose of A is the matrix A> ∈ Rn×m , with (A> )ij = aji for all i ∈ {1, . . . , n}
and j ∈ {1, . . . , m}. Graphically, this can be written as
a1,1 a2,1 . . . am,1
A> = ... .. .. .. ,
. . .
a1,n a2,n . . . am,n
Matrix-matrix addition
The matrices must have the same order.
C = A + B = B + A has elements cij = aij + bij ∀i, j.
Example (Matrix-matrix addition)
a11 a12 b11 b12 a11 + b11 a12 + b12
For A = and B = , A+B= .
a21 a22 b21 b22 a21 + b21 a22 + b22
Matrix-vector Product
If A ∈ Rm×n and v ∈ Rn , then Av ∈ Rm is the vector with
n
X
(Av)i = aij vj
j=1
for all j ∈ {1, . . . , n}. Going one step further, this notation also motivates the expression u> v in
equation 12.1.
Example (Row vector-matrix product)
> b11 b12
Let a = (a1 , a2 ) and B = . Then,
b21 b22
> b11 b12
a B = (a1 a2 ) = (a1 b11 + a2 b21 a1 b12 + a2 b22 ),
b21 b22
and
b11 b12 a1 b11 a1 + b12 a2
Ba = = .
b21 b22 a2 b21 a1 + b22 a2
12.2. MATRICES 205
Matrix-matrix Product
m
X
(AB)ik = aij bjk
j=1
for all i ∈ {1, . . . , `} and k ∈ {1, . . . , n}. This is called the matrix product of A and B. Note that A
and B must have compatible shapes for the product to exist.
yn yn i
Note (Properties)
• The matrix product is associative: if A, B and C are matrices with shapes such that AB
and BC exist, then we have A(BC) = (AB)C. It does not matter in which order we
perform the matrix products here.
• The matrix product is transitive: if A, B and C have the correct shapes, we have A(B+C) =
AB + AC.
• The matrix product is not commutative: if AB exists, in general A and B don’t have the
correct shapes for BA to also exist, and even if BA exists, in general we have AB 6= BA.
Identity matrix
1 0
For example I = .
0 1
Note (Properties)
Diagonal matrix
Definition 12.2.3. A matrix D ∈ Rn×n is called a diagonal matrix if all its off-diagonal entries are
zero:
dij = 0 for all i 6= j.
Graphically:
d11 0 · · · 0
0 d22 · · · 0
D= . .
.. .. ..
.. . . .
0 0 ··· dnn
d1 0
For example D = .
0 d2
Determinant of a matrix
Definition 12.2.4. Let A ∈ Rn×n be a square matrix. The determinant of A, denoted by det(A)
or |A|, is defined as:
X n
Y
det(A) = sgn(σ) ai,σ(i) ,
σ∈Sn i=1
where:
1. Sn is the set of all permutations of {1, 2, . . . , n}, 2. sgn(σ) is the sign of the permutation σ (+1
for even permutations, −1 for odd permutations), 3. ai,σ(i) are the entries of A.
Graphically, for n = 2:
a11 a12
det = a11 a22 − a12 a21 .
a21 a22
12.2. MATRICES 207
Note
The determinant is a scalar value that characterizes a square matrix. It is zero if and only if the
matrix is singular (non-invertible).
Rank
Definition 12.2.5. The rank of a matrix A ∈ Rm×n is the dimension of the image space of A:
Rank(A) = dim Av v ∈ Rn }.
The rank can also be characterised as the largest number of linearly independent columns of
A, or the largest number of linearly independent rows of A. Thus, for A ∈ Rm×n we always have
Rank(A) ≤ min(m, n). The matrix A is said to have ”full rank” if Rank(A) = min(m, n).
Trace
Definition 12.2.6. The trace of a matrix A ∈ Rn×n is given by
n
X
tr(A) = aii .
i=1
Properties:
• tr(A + B) = tr(A) + tr(B)
• tr(A> ) = tr(A)
• tr(ABC) = tr(BCA) = tr(CAB). The individual matrices A, B, C don’t need to be square for
this relation to hold, but the relation only holds for cyclic permutations as shown. In general
tr(ABC) 6= tr(ACB).
Matrix Inverse
If A is a square matrix and if there is a matrix B such that AB = I, then A is called invertible and
the matrix B is called the inverse of A, denoted by A−1 = B.
Some important properties of the inverse:
• The inverse, if it exists, is unique.
• Left-inverse and right-inverse for matrices are the same: A−1 A = I holds if and only if AA−1 = I.
• If A is symmetric and invertible, then A−1 is also symmetric. This is true because A(A−1 )> =
(A−1 A)> = I> = I and thus (A−1 )> is an inverse of A. Since the inverse is unique, (A−1 )> = A−1 .
Theorem 12.2.7. Let A ∈ Rn×n . Then the following statements are equivalent:
1. A is invertible
Orthogonal Matrices
• If U is orthogonal, the inverse and the transpose are the same: U> = U−1 .
• We have kUxk2 = x> U> Ux = x> x = kxk2 . Thus, multiplying a vector x by an orthogonal
matrix U does not change its length.
x> Ax > 0
x> Ax ≥ 0
for all x ∈ Rn .
Idempotent Matrices
12.3 Eigenvalues
Definition 12.3.1. Let A ∈ Rn×n be a square matrix and λ ∈ R. The number λ is called an
eigenvalue of A, if there exists a vector v 6= 0 such that Ax = λx. Any such vector x is called an
eigenvector of A with eigenvalue λ.
While there are very many results about eigenvectors and eigenvalues in Linear Algebra, here we
will only use a small number of these results.
We summarise what we need for this module: If A is idempotent and x is an eigenvector with
eigenvalue λ, then we have λx = Ax = A2 x = λAx = λ2 x. Thus we have λ2 = λ. This shows that the
only eigenvalues possible for idempotent matrices are 0 and 1.
Theorem 12.3.2. Let A ∈ Rn×n be symmetric. Then there is an orthogonal matrix U such that
D = UAU> is diagonal. The diagonal elements of D are the eigenvalues of A and the rows of U are
corresponding eigenvectors.
d >
Example (Differentiation dx c x)
x1 X
Let c> = (c1 , c2 ) and x = so c> x = ci xi = c1 x1 + c2 x2 . Then
x2
i
! !
d > d
d >
c x= dx1 c x = dx1 (c1 x1 + c2 x2 )
=
c1
= c.
d > d c2
dx dx2 c x dx2 (c1 x1 + c2 x2 )
d >
Example (Differentiation dx x x)
x1 X
If x = , then x> x = x2i and
x2
i
! !
d > d 2 + x22 )
d >
x x= dx1 x x = dx1 (x1 =
2x1
= 2x.
d > d 2 + x22 ) 2x2
dx dx2 x x dx2 (x1
Definitions
We start by revising the definition of the sequence of real numbers and sequence of functions.
{ai }∞
i=1 = {a1 , a2 , a3 , . . . , an , . . .}
{fi }∞
i=1 = {f1 , f2 , f3 , . . . , fn , . . .}
Now let us clarify the concept of a sequence of independent and identically distributed
random variables.
In any probability model, we have a sample space S and a probability measure P. For simplicity,
assume the sample space consists of a finite number of elements:
S = {s1 , s2 , · · · , sk }
A random variable X is a mapping that assigns a real number to each outcome si :
X(si ) = xi , for i = 1, 2, · · · , k
For a sequence of random variables X1 , X2 , . . ., each Xj is a function from S to R:
Xj (si ) = xji , for i = 1, 2, · · · , k
Thus, a sequence of random variables is a sequence of functions Xi : S → R. Given a sequence of
random variables X1 , X2 , . . . , Xn and a corresponding sequence of sampled data x1 , x2 , . . . , xn , each
observation xi is a realization of the random variable Xi .
210 CHAPTER 12. APPENDIX
Consider the experiment of tossing a fair coin. Each toss is independent, and the outcomes are
”head” or ”tail”, coded as 0 and 1, respectively.
The sample space S contains all possible sequences of these outcomes.
If the experiment continues indefinitely, we model it with a sequence of random variables
X1 , X2 , . . .. Each Xi takes the value 0 (head) or 1 (tail). All variables are independent and
identically distributed (i.i.d.), and follow a Bernoulli distribution:
i.i.d. 1
Xi ∼ Bernoulli
2
Convergence
Let us now revise the definitions of convergence for sequences of numbers and functions.
Definition 12.5.3. A sequence {ai }∞ i=1 converges to a real number A if, for every ε > 0, there exists a
positive integer nε such that for all n ≥ nε , we have:
|an − A| < ε
Definition 12.5.4. Suppose fi : X → R for all i. The sequence {fi }∞ i=1 converges pointwise to a
function f if, for all x ∈ X and ε > 0, there exists an nε such that for all n ≥ nε :
often denoted by µ.
Variance of X is
Var[X] = E[(X − µ)2 ] = E[X 2 ] − µ2 ,
where E[X 2 ] = 2 P (X
P
xx = x).
Discrete distributions
Bernoulli trials
A random variable X is said to be a Bernoulli random variable with parameter p, shown as X ∼
Bernoulli(p), if its Probability Mass Function (PMF) is given by
12.6. PROBABILITY REMINDERS 211
p for x = 1
PX (x) = 1−p for x = 0
0 otherwise,
Binomial distribution
n x
If X ∼ Bin(n, θ), then P (X = x) = θ (1 − θ)n−x for x = 0, 1, 2, . . . , n.
x
The expected value is E[X] = nθ and variance is Var[X] = nθ(1 − θ).
As an example, X might represent the number of heads in n tosses of a coin.
Poisson distribution
µx e−µ
If X ∼ Poisson(µ), then P (X = x) = for x = 0, 1, 2, 3, . . ..
x!
The expected value is E[X] = µ and variance is Var[X] = µ.
212 CHAPTER 12. APPENDIX
Geometric distribution
1−θ 1−θ
The expected value is E[X] = and variance is Var[X] = .
θ θ2
As an example, X might represent the number of tails observed before the first head occurs in a
coin tossing experiment. This is the same as the number of tosses of the coin before the first head!
Z x0
FX (x0 ) = P (X ≤ x0 ) = f (x)dx,
−∞
Z b
P (a < X ≤ b) = fX (x)dx = FX (b) − FX (a).
a
Mean of X is
Z ∞
E[X] = xfX (x)dx,
−∞
often denoted µ.
214 CHAPTER 12. APPENDIX
Z ∞
E[X 2 ] = x2 fX (x)dx.
−∞
Continuous distributions
1 1 2
φ(z) = √ e− 2 z , −∞ < z < ∞
2π
Normal distribution
1 (x−µ)2
fX (x) = √ e− 2σ 2 , −∞ < x < ∞
2πσ 2
X −µ
Z= ∼ N(0, 1) =⇒
σ
X −µ x−µ x−µ x−µ
FX (x) = P (X ≤ x) = P ≤ =P Z≤ =Φ
σ σ σ σ
12.6. PROBABILITY REMINDERS 215
Exponential distribution
If X ∼ exponential(λ), then the expected value, variance, PDF and CDF are
1 1
E[X] = , Var[X] = 2 ,
λ λ
−λx
fX (x) = λe , x > 0, λ > 0,
Z x h ix
FX (x) = λe−λx dx = −e−λx = 1 − e−λx , x > 0.
0 0
216 CHAPTER 12. APPENDIX
Uniform distribution
If X ∼ uniform(0, 1), then the expected value, variance, PDF and CDF are
1 1
E[X] = , Var[X] = ,
2 12
1 if 0 ≤ x ≤ 1,
fX (x) =
0 otherwise,
Z x 1 if x > 1,
FX (x) = fX (x)dx = x if 0 ≤ x ≤ 1,
0
0 if x < 0.
More generally a uniform distribution might be defined for −∞ < a < b < ∞
12.6. PROBABILITY REMINDERS 217
a+b (b − a)2
E[X] = , Var[X] = ,
2 1
12
fX (x) = b−a if a ≤ x ≤ b,
0 otherwise,
Z x 1 if x > b,
x−a
FX (x) = fX (x)dx = if a ≤ x ≤ b,
0 b−a
0 if x < a.
Beta function
The beta function is defined as
Γ(α)Γ(β)
B(α, β) = .
Γ(α + β)
This was first studied by Euler and Legendre and arose from work on what Legendre referred to as
the Eulerian integral of the first kind, namely
Z 1
B(α, β) = xα−1 (1 − x)β−1 dx.
0
Beta distribution
if X ∼ Beta(α, β), then the expected value, variance and the PDF are
α αβ
E[X] = , Var[X] = 2
,
α+β (α + β) (α + β + 1)
Γ(α + β) α−1
fX (x) = x (1 − x)β−1 , 0 ≤ x ≤ 1, α > 0, β > 0,
Γ(α)Γ(β)
1
fX (x) = xα−1 (1 − x)β−1 , 0 ≤ x ≤ 1, α > 0, β>0 .
B(α, β)
Chi-squared Distribution
If X ∼ χ2k , then the expected value, variance, and the PDF are:
The t-distribution
Definition 12.6.1. Let Z ∼ N (0, 1) and Y ∼ χ2 (ν) be independent. Then the distribution of
Z
T =p (12.3)
Y /ν
is called the t-distribution with ν degrees of freedom. This distribution is denoted by t(ν).
• The R command pt(x, nu) gives the value Φν (x) of the CDF of the t(ν)-distribution.
• The R command qt(alpha, nu) can be used to obtain the α-quantile of the t(ν)-distribution.
Theorem 12.6.2. For a random variables X and constants (a, b) the following holds
E[aX] = aE[X],
E[X + b] = E[X] + b.
Proof
Proof.
Z Z
E[aX] = axfX (x)dx = a xfX (x)dx = aE[X],
R R
Z
E[X + b] = (x + b)fX (x)dx
R
Z Z
= xfX (x)dx + bfX (x)dx
R R
Z Z
= xfX (x)dx + b fX (x)dx = E[X] + b.
R R
Theorem 12.6.3. For the random variables X and Y the following holds
Proof
Proof.
ZZ
E[X + Y ] = (x + y)fXY (x, y)dxdy
2
ZRZ ZZ
= xfXY (x, y)dxdy + yfXY (x, y)dxdy
R2 R2
Z Z Z Z
= xdx fXY (x, y)dy + ydy fXY (x, y)dx
R R R R
Z Z
= xfX (x)dx + yfY (y)dy = E[X] + E[Y ].
12.7. STATISTICAL METHODS REMINDERS 221
Proof
Proof.
Var[X] = E (X − E[X])2
= E X 2 − 2E[X]X + E[X]2
Theorem 12.6.5. For a random variables X and constants (a, b) the following holds
Proof
Proof.
Definition 12.7.2. We define sample covariance and correlation as following. For paired data
(xi , yi ):
n
1 X sxy
sxy = (xi − x)(yi − y), rxy = q ,
n−1 s2 s2
i=1 x y
1. sxx = s2x .
σ2 x2 σ 2 1
, where x2 = x2i .
P
2. Var(β̂) = , Var(α̂) = n
(n − 1)s2x (n − 1)s2x
4. σ̂ 2 is unbiased for σ 2 .
Chapter 13
Assessment Information
Assessment Summary
• Each set will be available for practice for two weeks only; the first set will be made available
at the beginning of Week 3.
• The problems may come from the following topics: moments, univariate transformation,
bivariate transformation, Bivariate Normal, MGF, Limits Theorems, Estimations, CLT, Fisher
Information, Null hypothesis, Multiple LR.
Note
Note: Practising these quizzes will mirror the style, and problem types in the OLTA, helping
you build fluency and confidence. You will unlimited attempts during the 2 weeks gap.
223
224 CHAPTER 13. ASSESSMENT INFORMATION
Exam Structure
The exam will be divided into two sections:
• Section Bank of questions (located in the HTML version of the notes) contains a bank of
questions (without answers). Some of these questions will appear in Section A of the
exam. Use the bank to guide your revision.
• A Section B template will be provided to to help you understand the exam structure for the
Section B.
Preparation Tips
• Start early: Use the Week 3-11 window to complete all formative sets while they are available.
• Use the question bank: Focus on patterns and concepts that recur in the bank, some items
will appear in Section A.
• Practise with the template: Once released in Week 6, use the template to practise full
solutions for Section B.
• Extensions: If you require an extension, you must apply via the University’s Mitigating Cir-
cumstances process. The maximum extension available for this assessment is 7 days (until 14
May).
• Late submissions without an extension: A 3-day grace period applies immediately after
the deadline. Late penalties are applied as follows:
After 72 hours (3 days), submissions will be recorded as not submitted unless an extension has
been granted.
• Late submissions with an approved extension: Students who are granted an extension also
receive the same 3-day late-submission grace period, with the same penalties applied.
• Feedback release: Feedback and solutions will be released after 17 May, once all possible
submission windows (including extensions and grace periods) have closed.
226 CHAPTER 13. ASSESSMENT INFORMATION
Bibliography
[1] D. De Angelis, W. R. Gilks, and N. E. Day. “Bayesian projection of the acquired immune
deficiency syndrome epidemic”. In: Journal of the Royal Statistical Society: Series C (Applied
Statistics) 47.4 (1998), pp. 449–498. doi: https : / / doi . org / 10 . 1111 / 1467 - 9876 . 00123.
eprint: [Link] 9876.00123. url:
[Link]
[2] Niels Ehlers and Mogens Norn. “History of Acta Ophthalmologica in the first 50 years”. In: Acta
Ophthalmologica 90.1 (2012), pp. 3–9. doi: [Link]
x. eprint: [Link]
url: [Link]
[3] Wiebe R. Pestman and Ivo B. Alberink. Mathematical statistics. Berlin: Walter de Gruyter, 1998.
[4] John A. Rice. Mathematical statistics and data analysis. eng. Third, international edition. Belmont,
Calif: Brooks/Cole, 2007. isbn: 9781473758063.
[5] S. S. SHAPIRO and M. B. WILK. “An analysis of variance test for normality (complete samples)†”.
In: Biometrika 52.3-4 (Dec. 1965), pp. 591–611. issn: 0006-3444. doi: 10.1093/biomet/52.3-
4.591. eprint: [Link]
[Link]. url: [Link]
227