0% found this document useful (0 votes)
14 views70 pages

Harvard Statistics 211 Course Notes

Statistics 211 is a graduate-level course at Harvard focused on frequentist methods for statistical inference, aimed at first-year PhD students. The course covers topics such as point estimation, hypothesis testing, and Bayesian inference, with an emphasis on good statistical thinking. Taught by Lucas Janson, the course includes active learning components and utilizes provided notes along with supplemental textbooks.
Copyright
© All Rights Reserved
We take content rights seriously. If you suspect this is your content, claim it here.
Available Formats
Download as PDF, TXT or read online on Scribd
0% found this document useful (0 votes)
14 views70 pages

Harvard Statistics 211 Course Notes

Statistics 211 is a graduate-level course at Harvard focused on frequentist methods for statistical inference, aimed at first-year PhD students. The course covers topics such as point estimation, hypothesis testing, and Bayesian inference, with an emphasis on good statistical thinking. Taught by Lucas Janson, the course includes active learning components and utilizes provided notes along with supplemental textbooks.
Copyright
© All Rights Reserved
We take content rights seriously. If you suspect this is your content, claim it here.
Available Formats
Download as PDF, TXT or read online on Scribd

Statistics 211: Statistical Inference I

Eric K. Zhang
ekzhang@[Link]

Fall 2021

Abstract
These are notes for Harvard’s Statistics 211, a graduate-level class taught by Lucas Janson
in Fall 2021, targeted at first-year PhD students. The main focus of this class is on frequentist
methods for statistical inference, i.e., how to draw mathematical conclusions from sample data
based on likelihoods from a parametric model.
Course description: Foundations of frequentist and Bayesian inference, and decision the-
ory. Likelihood, sufficiency, and ancillarity. Point estimation, unbiasedness, maximum like-
lihood, method of moments, minimum-variance. Parametric and non-parametric hypothesis
testing, confidence intervals. Selective inference: multiple testing, familywise error rate, false
discovery rate. Bayesian inference, conjugate priors, credible intervals. Admissibility, Stein’s
phenomenon, empirical Bayes. Time permitting: post-selection inference and the bootstrap.

Contents
1 September 1st, 2021 4
1.1 Class Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.2 Statistical Philosophy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4

2 September 8th, 2021 6


2.1 Likelihood and Notation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
2.2 Sufficiency . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6

3 September 13th, 2021 9


3.1 Exponential Families and Sufficiency . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
3.2 Unbiased Estimation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
3.3 Minimal Sufficiency . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11

4 September 15th, 2021 12


4.1 Minimal Sufficiency (cont.) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
4.2 Complete Sufficiency . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
4.3 Optimal Unbiased Estimation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13

5 September 20th, 2021 14


5.1 More on Lehmann-Scheffé . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
5.2 Ancillary Statistics and Basu’s Theorem . . . . . . . . . . . . . . . . . . . . . . . . . 14

1
6 September 22nd, 2021 16
6.1 The Score Function . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
6.2 Fisher Information . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

7 September 27th, 2021 19


7.1 More on Cramér-Rao . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
7.2 Method of Moments Estimation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
7.3 Maximum Likelihood Estimation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20

8 September 29th, 2021 22


8.1 Consistency of MLE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
8.2 Asymptotic Normality of MLE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22

9 October 4th, 2021 24


9.1 The Delta Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
9.2 Hypothesis Testing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25

10 October 6th, 2021 26


10.1 The Neyman-Pearson Lemma . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
10.2 Testing Composite Hypotheses . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27

11 October 13th, 2021 28


11.1 Karlin-Rubin Test . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
11.2 Likelihood-Ratio Test . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28

12 October 18th, 2021 31


12.1 Asymptotic Hypothesis Tests . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
12.2 Nonparametric Hypothesis Tests . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32

13 October 20th, 2021 34


13.1 Sign Test and Signed-Rank Test . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
13.2 Two-Sample Tests . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35

14 October 25th, 2021 37


14.1 Tests of Independence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
14.2 Selective Inference . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38

15 October 27th, 2021 40


15.1 FWER-Controlling Procedures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
15.2 FDR-Controlling Procedures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41

16 November 1st, 2021 42


16.1 More on FDR Control . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
16.2 Confidence Intervals . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43

17 November 3rd, 2021 45


17.1 CIs by Inverting Hypothesis Tests . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
17.2 CIs by Using Pivotal Quantities . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
17.3 Criteria for Selecting CIs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46

2
18 November 8th, 2021 48
18.1 Asymptotic Confidence Intervals . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
18.2 Bayesian Inference . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
18.3 Conjugate Priors of NEFs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50

19 November 10th, 2021 52


19.1 More on Conjugacy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
19.2 The Jeffreys Prior . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
19.3 Bayesian Point and Interval Estimation . . . . . . . . . . . . . . . . . . . . . . . . . 53

20 November 15th, 2021 55


20.1 Bayesian Predictive Inference . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
20.2 Frequentist Decision Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
20.3 Bayesian Decision Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56

21 November 17th, 2021 58


21.1 Bayes Rules are Admissible . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
21.2 Admissible Rules are Bayes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59

22 November 22nd, 2021 61


22.1 Proof of the Complete Class Theorem . . . . . . . . . . . . . . . . . . . . . . . . . . 61
22.2 Admissibility of the Sample Mean . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
22.3 Least Favorable Priors . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63

23 November 29th, 2021 64


23.1 More on Minimax Decision Rules . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
23.2 Stein’s Paradox . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65

24 December 1st, 2021 67


24.1 Proof of Stein’s Theorem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
24.2 Properties of the James Stein Estimator . . . . . . . . . . . . . . . . . . . . . . . . . 67
24.3 Modern Statistical Research . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68

3
1 September 1st, 2021
This is the first lecture of the course. We will discuss logistics, an overview of the class, and a bit
of statistical philosophy.

1.1 Class Overview


Our professor is Lucas Janson, who is a professor in statistics at Harvard focusing on high-
dimensional inference and statistical machine learning problems. Our teaching fellows are Yufan
Li, Biyonka Liang, and Yash Nair. The prerequisites for this class are an undergraduate-level infer-
ence class, such as Stat 111, and some exposure to graduate-level probability (e.g., CLT, Jensen’s
inequality, Slutsky’s theorem, continuous mapping theorem, convergence).
Most of the material in this class will come from the provided notes, but we also have three
textbooks [CB21, LR06, LC06] for supplemental use. Lucas will also post modern literature that
is related to each lecture topic at the end of class. Each lecture will include an “active learning”
component where students think individually about a problem.
Now we discuss course topics. In most sections of this class, there will be some data y and
parameter θ, and our goal is to determine some estimate of θ that is as close to the true value as
possible, along with a measure of uncertainty about our estimate. The main topics are:

1. Point estimation: How can we find an estimate for θ?

2. Confidence intervals: How confident is our estimate, under Bayesian or frequentist terms?

3. Hypothesis testing: Is there significant evidence of a hypothesis θ = θ0 being false?1

4. Selective inference: How can we test multiple hypotheses simultaneously?

5. Decision making: If we assign costs to each type of error, what estimator minimizes cost?

6. Prediction: If I collect new data according to y = fθ (x), what would it look like?

Lucas emphasizes that this is a course in good statistical thinking, and although not all enrollees
are the target audience of first-year PhD students, he hopes that the class content is broadly useful
in many contexts.

1.2 Statistical Philosophy


We first discuss the difference between Bayesian and frequentist interpretations of inference. In
Bayesian inference, our model is that θ is random and y is fixed. However, in frequentist thinking,
we interpret our parameters as y being random and θ being fixed.
There are different techniques and notations for Bayesian and frequentist methods, like confi-
dence intervals versus credible intervals. However, in practice, both approaches will achieve similar
results, so we should not draw a hard distinction.2 Let’s pivot to one of Lucas’s core beliefs.

Proposition 1.1 (First Law of Statistics). The more you know about your data’s distribution, the
more you can infer about that distribution from the data.
1
This includes non-parametric hypothesis testing, such as asking if a distribution is symmetric.
2
Lucas tells us to take Stat 213 for a theoretical justification of this.

4
In other words, when you know some facts about the distribution, such as the fact that it is
Gaussian with variance 1, you can produce better inference methods for quantities like the mean
µ. If you knew nothing about the distribution, it could be the case that the mean does not exist
1
at all, as is the case with the Cauchy distribution p(x) = π(1+x 2) .

Of course, this statement is a tautology, since a more knowledgeable individual could just
pretend to know less about the data. However, it has a couple important conceptual consequences
to the way we approach statistics problems:

1. When we analyze data, first ask what we know about the distribution, then ask how to use
that knowledge to learn as much about the data as possible. This means that good statistics
practice is tied to knowledge of its domain of application.

2. Domain knowledge allows us to obtain assumptions about the distribution of data being
modeled, which are crucial to determining the type of inference method to use.

3. Bayesian and frequentist inference are different ways of encoding domain knowledge. Bayesian
inference is better at encoding assumptions about where θ will be in the parameter space,
while frequentist inference is better at problems where we do not know the prior.

In scientific literature, there are many statistical methods that are applied over and over again.
We will not learn names of domain methods in this class, but they will typically be special cases
of inference techniques we learn in this class. More importantly, we will understand where these
methods come from and how to extend them.
In statistics, methodological innovation comes from doing better with the same set of assump-
tions, or by leveraging more assumptions within a method. Neither Bayesian nor frequentist ap-
proaches are perfect, since approximations will never exactly match the truth. Therefore, we can
summarize the relevance to this class as follows:

• In Bayesian inference, there is a rigid way of specifying domain knowledge. Once we specify
the domain knowledge, we’re basically done, as we just examine the posterior. Our main
challenge is computational rather than methodological, so we will not focus on Bayesian
inference in this class. It’s too “elegant and simple” to talk about.

• In frequentist methods, this is not true, and different methods can be better at certain parts
of the parameter space. Therefore, it is often the case that complex problems do not have
optimal frequentist inference methods.

That concludes the first lecture of the course. Next week, we will continue discussing frequentist
and Bayesian philosophies and sufficiency.

5
2 September 8th, 2021
Today we discuss likelihood and relevant notation, sufficiency, and unbiased estimation (the last
topic, only if time permits).

2.1 Likelihood and Notation


We start by reiterating the difference between probability and inference problems. In probability
tasks, we are given a fixed distribution and are asked the probabilities of certain events in that dis-
tribution. However, in inference, the task is reversed: given certain observations from an unknown
distribution, find the likelihood of the distribution having certain parameters.

Definition 2.1 (Parametric model). Suppose that we have a model fθ (y) representing the prob-
ability mass of some distribution at y, which is defined on some base measure.3 We call fθ (y)
parametric if θ is finite-dimensional and nonparametric otherwise.

Now we will set up the likelihood function, which is key in statistical inference.

Definition 2.2 (Likelihood function). Given a model fθ (y) and observed data y, the likelihood
function on parameter space is defined by

L(θ) = fθ (y).

We notate the log-likelihood by `(θ) = log L(θ).

Oftentimes, our observations Y will be multi-dimensional, and we write them as a vector Y =


(Y1 , . . . , Yn ). Also, in addition to the model probability density fθ (y), we notate the cumulative
distribution function by a capital letter Fθ (y) = Prθ (Y ≤ y).

Note. We often parameterize probability and expectation by subscripting the relevant operator
with the parameter θ. For example, we could write Prθ (Y ∈ A) or Eθ [Y ].

After getting some basic notation out of the way, we can now start the first “unit” of the course,
which is about sufficiency.

2.2 Sufficiency
Sufficiency is a concept related to how useful a set of observations is for predicting the parameters
of a distribution. Essentially, a sufficient statistic is “good enough” for getting all the statistical
information about parameters from the underlying set of variables.

Definition 2.3 (Sufficiency). Given random variables Y = (Y1 , . . . , Yn ) ∼ fθ (y), we say that a
statistic T (Y ) computed from Y is sufficient for θ if for all θ1 , θ2 ∈ Θ and A ⊂ Y, we have
 
Pr Y ∈ A | T = T (Y ) = Pr Y ∈ A | T = T (Y ) .
θ1 θ2

There are a few other equivalent definitions based on notions from other fields:

• (Information Theory). The chain θ → T → Y is Markovian, i.e., Y ⊥⊥ θ | T .

• (Bayesian). The conditional distribution of θ | T is the same as θ | Y , for any prior on θ.


3
Note that this definition supports both continuous and discrete random variables.

6
• (Measure Theory). The statistic T (Y ) and variable Y generate the same σ-algebra.

Example 2.4. Suppose that Y1 , Y2 are i.i.d. ∼ Pois(λ), and T (Y ) = Y1 + Y2 . Then, the conditional
distribution is (Y1 , Y2 ) | T (Y ) ∼ Mult(T, ( 12 , 21 )). We claim that T (Y ) is sufficient for the Poisson
rate parameter λ. This is because, following the definition,

Prλ (Y1 + Y2 = t | Y1 = k) Prλ (Y1 = k)


Pr(Y1 = k, Y2 = t − k | Y1 + Y2 = t) =
λ Prλ (Y1 + Y2 = t)
λt−k e−λ λk e−λ
(t−k)! · k!
= (2λ)t e−2λ
t!
 .
t
= 2t .
k

This last expression does not depend on λ, so we have established sufficiency.

Example 2.5. If Y1 , . . . , Yn are i.i.d. ∼ Bern(p), then T (Y ) = Y1 + · · · + Yn is a sufficient statistic


for p. We can verify this property in a similar manner to the above example.

Oftentimes, we may have an intuitive suspicion that a statistic is sufficient, but this requires
algebraic verification to ensure that the property actually holds.

Exercise 2.1 (Pencil problem). You have two minutes to think about these questions.

1. If T is sufficient, and g is injective, is g(T ) sufficient?

2. If Yi ∼ Bern(pi ) for i = 1, 2, is Y1 + 2Y2 sufficient for (p1 , p2 )?

Now, here is an intuitive but nontrivial statement about probability densities of random vari-
ables with a sufficient statistic.

Theorem 2.6 (Factorization theorem). A statistic T (Y ) is sufficient for θ if and only if the joint
density fθ (y) can be written as
fθ (y) = gθ (T (y)) · h(y).
Here, the density gθ of the sufficient statistic is allowed to depend on θ, but the density h(y) is not
allowed to vary. In terms of log-likelihood, this is

`(θ) = log(gθ (T )) + log(h(y)).

Proof. This theorem holds in generality, but for simplicity of proof, here assume that Y is a discrete
random variable. Let’s first handle the easier direction, which is the “if” statement. Following the
definition of sufficiency, note that

Prθ (Y = y)
Pr(Y = y | T (Y ) = T (y)) =
θ Prθ (T (Y ) = T (y))
gθ (T (y))h(y)
=P 0 0
y 0 |T (y 0 )=T (y) g(T (y ))h(y )
h(y)
=P 0
.
y 0 h(y )

7
This does not depend on the parameter θ, so we are done. On the other hand, in the inverse
direction, we can use conditional probability and the definition of sufficiency to get

fθ (y) = Pr(Y = y)
θ
= Pr(Y = y, T (Y ) = T (y))
θ
= Pr(Y = y | T (Y ) = T (y)) · Pr(T (Y ) = T (y)).
θ θ

The former is a function free of θ, while the latter is a function of T (y), and hence we have completed
the factorization.

Note. In the above factorization, neither gθ nor h are unique, nor are they required to be actual
probability measures with physical meaning whose masses sum to 1. It helps me to intuitively
think of gθ as a measure over T and h as a probability within each quotient equivalence class of
the form {y ∈ Y | T (y) = t}. In the statement of the theorem, h is more of a conditional density
Y | T shoved into the shape of a function on all of Y .

Here’s another example of sufficiency, with applications to many fundamental distributions.

Example 2.7 (Exponential family). Let Y1 , . . . , Yn be i.i.d. ∼ fθ (y), belonging to an exponential


family
fθ (y) = exp{η(θ)T (y) − ψ(η(θ))}h(y).
Here, η is called the natural parameter of the family. This has joint density
n n
( ) n
Y X Y
fθ (yi ) = exp η(θ) T (yi ) − nψ(η(θ)) h(yi ).
i=1 i=1 i=1

We will continue this example in the next lecture.

8
3 September 13th, 2021
Today we will first start by redefining exponential families in a way that is hopefully more intuitive.
We will then talk about unbiased estimation and minimal sufficiency.

3.1 Exponential Families and Sufficiency


Exponential families are an important topic that is related to both sufficiency and this class in
general. It will show up multiple times, so we will spend some time clearly motivating and defining
it from a clean slate.

Definition 3.1 (Exponential family). let Y1 , . . . , Yn be i.i.d. ∼ fθ (y), where fθ (y) ∝ eθT (y) h(y) for
some functions T and h. Then, we define the normalizing factor ψ(θ) to be
Z ∞
eψ(θ) = eθT (y) h(y) dy.
−∞

Under this definition, we can write the density as fθ (y) = exp{θT (y) − ψ(θ)}h(y). We call fθ an
exponential family with respect to some parameter η such that θ = θ(η).

In an exponential family, the joint probability density of Y1 , . . . , Yn is


( n ) n
X Y
fθ (y) = exp θ T (yi ) − nψ(θ) h(yi ).
i=1 i=1
Pn
By Theorem 2.6, since θ only appears in the above expression once, we know that i=1 T (Yi ) is a
sufficient statistic for θ.

Note. This sufficiency is a powerful result, since exponential families are very general and include
many of the standard continuous distributions in statistics. It means that we can find a sufficient
statistic any exponential family by just adding together the values of T (Yi ).

To generalize further, when θ is a k-dimensional parameter vector, the probability density for
an exponential family is given by
 
Xk 
fθ (y) = exp θj Tj (y) − ψ(θ) h(y).
 
j=1

In the k-dimensional scenario, the analogous k-vector is a sufficient statistic:


n n n
!
X X X
T1 (Yi ), T2 (Yi ), . . . , Tk (Yi ) .
i=1 i=1 i=1

The Normal, Binomial, Gamma, and Poisson distributions are all exponential families. Note that
our analysis above has elided mention of the transformation η(θ), which turns the parameter θ that
we care about into the natural parameter η.

Example 3.2 (Normal distribution is an EF). Recall that the normal distribution N (µ, σ 2 ) has
probability density
1 (y−µ)2
fµ,σ2 (y) = √ e− 2σ2 .
2πσ 2

9
This is an exponential family with two parameters. The natural parameters are
 
2 µ 1
η(µ, σ ) = ,− .
σ 2 2σ 2

The exponential maps are T (y) = (y, y 2 ), and then the density can be written as
 
X 2 
fµ,σ2 (y) = exp ηj (µ, σ 2 )Tj (y) − ψ(η) .
 
j=1

This is called canonical form. It shows that the sum of T (y), i.e., the first and second empirical
moments of the data, are a sufficient statistic for N (µ, σ 2 ).

Generally, if we see a new distribution and are asked to find a sufficient statistic, we should try
to write that distribution as an exponential family. This will provide a clear and methodical way
to find a sufficient statistic. Furthermore, we will see later on that some exponential families have
a unique sufficient statistic governed by the natural law.
Now, let’s go over an example of sufficiency outside of the exponential family framework.

Example 3.3 (Uniform distribution sufficient statistic). Suppose that Y1 , . . . , Yn are i.i.d. ∼
Unif[0, θ]. Let Y(1) , . . . , Y(n) be the order statistics. Then, the joint uniform density is

1
fθ (y) = 1 1 .
θn {y(n) ≤θ} {y(1) ≥0}
Then, by Theorem 2.6, we immediately conclude that Y(n) is sufficient.

It actually turns out that order statistics are a more general technique.

Example 3.4 (Order statistics are sufficient). Let fθ be any density with respect to the Lebesgue
measure, parameterized by some scalar θ. Then, if Y1 , . . . , Yn are i.i.d. ∼ fθ , then the order statistics
Y(1) , . . . , Y(n) are a nontrivial sufficient statistic for θ. This is because the joint density is
n
Y n
Y
f (yi ) = f (y(i) ).
i=1 i=1

3.2 Unbiased Estimation


Estimation is a common problem in statistics, and we will place it in our inference framework.

Definition 3.5 (Unbiased estimator). For an estimated quantity g(θ), we say that an estimator
T (Y ) is unbiased if Eθ [T (Y )] = g(θ) for all values of θ.

Unbiased estimators tend to be “good” in terms of mean-squared error, since we can decompose
the error into a bias and variance term like

Eθ (T (Y ) − g(θ))2 = Eθ (T (Y ) − Eθ [T (y)])2 + (Eθ [T (Y )] − g(θ))2


   

= Varθ [T (Y )] + Bias2 .

Now, we introduce a famous theorem that connects sufficient statistics to unbiased estimation.

10
Theorem 3.6 (Rao-Blackwell). Let W (Y ) be an unbiased estimator of g(θ) and T be a sufficient
statistic for θ. Consider the estimator ψ(T ) = Eθ [W (Y ) | T ].4 Then,

1. Eθ [ψ(T )] = g(θ).

2. ψ is “better” than W , meaning that Varθ [ψ(T )] ≤ Varθ [W (Y )], with the equality case being
when ψ(T ) = W for all θ.

Proof. The first part follows directly from the law of iterated expectation, since

g(θ) = Eθ [W (Y )] = Eθ [Eθ [W (Y ) | T ]] = Eθ [ψ(T )] .

The second part follows from the law of total variance, since

Varθ [W (Y )] = Eθ [Varθ [W (Y ) | T ]] + Varθ [Eθ [W (Y ) | T ]] ≥ Varθ [ψ(T )] .

Intuitively, this theorem means that if we average the value of our unbiased estimator across
fibers of a sufficient statistic, we can essentially uniformly reduce the variance of that estimator.

3.3 Minimal Sufficiency


Recall that if a statistic T is sufficient for θ, then any injective function f (T ) is also sufficient for θ.
Some sufficient statistics convey more information than others, but we often want less information
to be encoded in our statistic, rather than more.

Definition 3.7 (Minimal sufficient statistic). A sufficient statistic T is called minimal if T is a


function of any other sufficient statistic.

We present some examples below, although we do not have time to justify them:

• If Y1 , . . . , Yn ∼ N (µ, σ 2 ), then
n n
!
X X
Yi , Yi2
i=1 i=1

is a minimal sufficient statistic. Equivalently, this means that the empirical mean and stan-
dard deviations (without the n − 1 correction) are minimal sufficient.

• If Y1 , . . . , Yn ∼ Unif[θ − 12 , θ + 12 ], then (Y(1) , Y(n) ) is minimal sufficient.

• If Y1 , . . . , Yn ∼ f , where f is any density, then (Y(1) , . . . , Y(n) ) is minimal sufficient.

In general, minimal sufficient statistics are unique if they exist, up to bijective maps. Furthermore,
they almost always exist, except for a few pathological cases.

Theorem 3.8. Let Y be an i.i.d. vector of data with Yi ∼ fθ , and let T (Y ) be a sufficient statistic
for θ. If for all pairs x, y such that ffθθ (x)
(y) is free of θ, T (x) = T (y), then T is minimal sufficient.

We’ll talk more about this theorem in the next class.

4
Notice that this is only a reasonable definition when T is a sufficient statistic, as otherwise it would depend on θ.

11
4 September 15th, 2021
Today we will finish talking about minimal sufficiency, introduce completeness, and discuss optimal
unbiased estimation.

4.1 Minimal Sufficiency (cont.)


First, we prove the theorem from the end of last lecture. It should hopefully be fairly intuitive
based on the definition of minimal sufficiency, but the formal details require some thinking.

Proof of Theorem 3.8. Let T 0 be another sufficient statistic for θ. Then, for any x and y such that
T 0 (x) = T 0 (y), we have by the factorization theorem that
fθ (x) gθ (T 0 (x))h(x) h(x)
= 0
= .
fθ (y) gθ (T (y))h(y) h(y)
This is free of θ, so by the minimal sufficiency assumption, T (x) = T (y). Therefore, we conclude
that T (x) is a function of T 0 (x) because its fibers are at least as coarse.

Intuitively, another way of thinking about minimal sufficiency is that the fibers of T exactly
encode all dependencies of the relative distribution mass on θ. Let’s provide an example of using
this theorem to prove minimal sufficiency.
Example 4.1. Let Y1 , . . . , Yn ∼ N (µ, σ 2 ). Then, using the normal distribution density,
" ! !#
fθ (x) 1 X 2 X 2 µ X X
= exp − 2 xi − yi + 2 xi − yi .
fθ (y) 2σ σ
i i i i

This quotient of probability densities is independent of our parameters (µ,P


σ) if and
P only  if the first
and second empirical moments for x and y are equal. Therefore, T (y) = y 2, y is minimal.
i i i i

4.2 Complete Sufficiency


There is a different but related special case of sufficiency that we describe now.
Definition 4.2 (Complete sufficient statistic). A sufficient statistic T (θ) for θ is called complete if
the only unbiased estimator of zero that is a function of T is the zero function. More precisely, if
for all θ,
Eθ [h(T (Y ))] = 0,
then h(T ) = 0 almost surely.
This means that any function whose expectation is independent of θ must essentially destroy
all information about T .
Example 4.3. Consider i.i.d. Bernoulli random variables. Consider T (y) = ni=1 yi , the sum of
P
these random variables, which is a sufficient statistic for the Bernoulli parameter p. It is also
complete sufficient, since if there is a function h(T ) such that Ep [h(T )] = 0 for all p ∈ (0, 1), then
n  
X n k
h(k) p (1 − p)n−k = 0.
k
k=0
p
This is an n-th degree polynomial in 1−p ∈ (0, ∞). For this to vanish at all values of p, each
coefficient must be zero, so h(k) = 0 for all k.

12
We will show very soon that complete sufficiency implies minimal sufficiency, except in patho-
logical cases. The converse does not hold.
Example 4.4. Consider Y1 , . . . , Yn to be i.i.d. ∼ Unif[θ − 12 , θ + 12 ], for some θ ∈ R. Then, as shown
earlier, T (Y ) = (Y(1) , Y(n) ) is minimal sufficient. However, T is not complete sufficient because if
h(T ) = Y(n) − Y(1) , then h(T ) | θ ∼ Beta(n − 1, 2). This means that
 
n−1
Eθ h(T ) − = 0,
n+1
so T is not complete sufficient.
Proposition 4.5. If a minimal sufficient statistic for some parameter θ exists, then any complete
sufficient statistic for θ is also minimal.
Proof. Let T and M be sufficient statistics such that T is complete and M is minimal. We would
like to prove that T is also minimal. The key observation is to consider the function

h(T ) = Eθ [T | M ] − T.

This is well-defined, since T is a function of M because M is a minimal sufficient statistic. Note


that Eθ [h(T )] = 0 by the law of iterated expectation. Since T is complete, we know that h(T ) = 0
almost surely, and the equality case of Theorem 3.6 implies that T is a function of M .

There is one more big result, which we will not prove.


Proposition 4.6. Any minimal sufficient statistic of an exponential family is complete, except in
pathological cases.
Proof. Consult the course notes for a citation about this result.

4.3 Optimal Unbiased Estimation


An important concept in statistical inference is the uniformly minimum variance unbiased estimator
(UMVUE), which is the “best” unbiased estimator for a statistic in the variance sense. Here is the
big theorem that ties everything together and motivates why we care about complete sufficiency.
Theorem 4.7 (Lehmann-Scheffé). An unbiased estimator of g(θ) that is a function of a complete
sufficient statistic is the unique UMVUE.
Proof. Let W and W̃ be two unbiased estimators of g(θ), and let T be a complete sufficient statistic
for θ. Then, φ(T ) = Eθ [W | T ] and φ̃(T ) = Eθ [W | T ]. By Rao-Blackwell (Theorem 3.6), φ and φ̃
are unbiased estimators with variance no greater than the variances of W and W̃ , respectively.
Let h(T ) = φ(T )− φ̃(T ), which means that Eθ [h(T )] = 0. By complete sufficiency, φ(T ) = φ̃(T )
almost surely. Now, assume for the sake of contradiction that W and W̃ are both UMVUEs, and
that W is a function of T , but W̃ is not a function of T . Then, W = φ, but
h i h i
Varθ [W ] = Varθ [φ] = Varθ φ̃ ≤ Varθ W̃ < Varθ [W ] .

This is a contradiction, so any unbiased estimator of g(θ) is a UMVUE. Furthermore, this is unique,
since we can just apply the Rao-Blackwell sledgehammer to arrive at a contradiction if multiple
UMVUEs exist.

13
5 September 20th, 2021
Today, we finish our discussion of the Lehmann-Scheffé theorem, then we introduce the concepts
of ancillary statistics and Basu’s theorem.

5.1 More on Lehmann-Scheffé


Recall that Proposition 4.5 implies that if the minimal and complete sufficient statistics both exist,
then they are equivalent and unique. It is very rare that the MSS does not exist and somewhat
rare that the CSS does not exist. Furthermore, no sufficient statistic in an exponential family is
complete except curved exponential families.
Example 5.1. Let Y1 , . . . , Yn ∼ Pois(λ). Then, we know that T = ni=1 Yi is a complete sufficient
P
statistic for λ. Note that Y1 is an unbiased estimator for λ, and Y1 | T ∼ Bin(T, n1 ), so
T
Eλ [Y1 | T ] =
n
is an unbiased estimator of λ that is a function of the CSS T , so it is the unique UMVUE for λ.
Note how powerful Lehmann-Scheffé was in the above example, as it made it super easy to
prove that the average of the samples is the UMVUE for the rate parameter λ, a nontrivial fact!
Example 5.2. There are some model classes for which no unbiased estimators exist, and there are
also some classes for which the UMVUE is ridiculous. Let Y ∼ Pois(λ), but let g(λ) = e−2λ . Then,
T (y) = (−1)y is a UMVUE for g(λ), even though it is clearly an absurd estimator.
Now, here’s a quick pencil exercise to test our understanding of unbiased estimators.
Exercise 5.1 (Pencil problem). Find a nontrivial model {fθ : θ ∈ Θ} such that Θ contains an
open subset of R, and there exists a biased estimator B(Y ) of the mean of Eθ [Y ] such that B(Y )
always has smaller mean-squared error than the sample mean, or UMVUE?
Proof. Yes, consider the model Y ∼ N (θ, 1), where θ ∈ Θ = [0, ∞). Then, the sample mean y has
strictly worse mean-squared error than the thresholded quantity max(0, y).

5.2 Ancillary Statistics and Basu’s Theorem


So far, we have been dealing with problems in this class where the model is already known. However,
sometimes we don’t know the model (e.g., in goodness-of-fit tests). The following definition will be
useful in cases where we want to do hypothesis testing on whether a model is satisfied.
Definition 5.3 (Ancillary statistic). A statistic A(Y ) is called ancillary for parameter θ if its
distribution does not depend on θ.
Here are some examples of ancillary statistics. Any constant is ancillary, but that is trivial.
However, the following two examples are classic and widely used in hypothesis testing.
Example 5.4 (χ2 test). If Y1 , . . . , Yn ∼ N (µ, 1), then A(y) = ni=1 (yi − y)2 ∼ χ2n−1 is ancillary.
P

Example 5.5. If Y1 , . . . , Yn ∼ Expo(λ), meaning that fλ (y) = λe−λy for y ≥ 0, then


yn
A(y) =
y1 + · · · + yn
is ancillary because we can let Ui = λYi ∼ Expo(1).

14
In general, this pattern is called a scale family, where we have a CDF of the form Pr Y ≤ y =
for σ > 0. Then, any statistic that depends on Y only through yy21 , . . . , yyn1 is ancillary.
F ( σy )

Example 5.6. Consider Yi ∼ Unif[θ − 12 , θ + 12 ], and recall that (Y(1) , Y(n) ) is a minimal sufficient
statistic. Note that the difference Y(n) − Y(1) is ancillary.

Exercise 5.2 (Pencil problem). Find an example of a model and two statistics where A1 and A2
are both ancillary, but (A1 , A2 ) is not ancillary.

Proof. Let (Xi , Yi ) ∼ N ([ 00 ], [ 1r 1r ]). Then, A1 = X ∼ N (0, 1) and A2 = Y ∼ N (0, 1) are both
ancillary, but (A1 , A2 ) is not, since its distribution depends on the covariance r.

Now, we’re ready to introduce a very elegant and useful result.

Theorem 5.7 (Basu’s theorem). If T is a complete sufficient statistic for θ and A is ancillary for
θ, then A ⊥⊥ T .

Proof. For any measurable set B, consider hB (T ) = Prθ (A ∈ B | T ) − Prθ (A ∈ B). This is a
function of T that does not depend on θ, since T is a sufficient statistic and A is ancillary. By the
law of iterated expectation,
 
Eθ [hB (T )] = ET Pr(A ∈ B | T ) − Pr(A ∈ B) = 0.
θ θ

Therefore, since T is a CSS, we conclude that hB (T ) = 0 almost surely, so T ⊥⊥ A.

What makes Basu’s theorem so useful in practice is that it allows us to find independent statistics
that are not obvious at first glance. Sometimes we might even add a variable to our model in
order to apply Basu’s theorem.

Example 5.8. If Yi ∼ N (0, σ 2 ), and Y = n1 ni=1 Yi , while S = ni=1 (Yi − Y )2 , then Y ⊥⊥ S.


P P

Proof. Consider the family of distributions N (µ, σ 2 ) parameterized by µ, with σ 2 known. Then, Y
is a complete sufficient statistic for µ, while S is ancillary. By Theorem 5.7, Y ⊥⊥ S for any µ.

Here’s the last example that we’ll discuss today.

Example 5.9. Suppose that Y1 , . . . , Yn ∼ N (0, σ 2 ). Let M be the median of Y1 , . . . , Yn . Then,

Cov(Y , M ) = Cov(Y , M − Y + Y )
= Cov(Y , M − Y ) + Cov(Y , Y ).

However, note that M − Y is ancillary for µ, so by Theorem 5.7, Y ⊥⊥ M − Y . Independence


implies zero correlation, so the covariance is zero, and we conclude that

σ2
Cov(Y , M ) = .
n
σ2
Exercise 5.3 (Pencil problem). Based on the example above, show that Cov(Y(1) , Y ) = n .

15
6 September 22nd, 2021
Today we discuss feedback, likelihood derivatives (the score function), and the Cramér-Rao lower
bound for variance of unbiased estimators.

6.1 The Score Function


So far we’ve talked about the likelihood function L(θ) = fθ (y) for given data samples y, as well as
the log-likelihood function `(θ). What happens when we take the gradient of the log-likelihood?

Definition 6.1 (Score function). Given data samples Y = (Y1 , . . . , Yn ) drawn from fθ (y), the score
function S(y, θ) is given by
∂ log fθ (y) ∂ `(θ)
S(y, θ) = = .
∂θ ∂θ
Note that the score function for multiple data points equals the sum of the score functions for each
individual data point.

Among other uses, the score function is useful for finding maximum likelihood estimators. In
particular, if the MLE for y occurs at θ̂ in the interior of parameter space Θ, and the log-likelihood
function ` is differentiable at θ̂, then it must satisfy S(y, θ̂) = 0.

Definition 6.2 (Differentiating under the integral). Given a family of densities parameterized by
θ for a random variable Y with respect to measure µ, we say that a function gθ (y) satisfies the
m-th order EDI condition if
∂m dm
Z Z
gθ (y) dy = g (y) dy.
m θ
∂θm µ µ dθ

This is a useful regularity condition, and there are various real analysis results that can be used
to prove this, such as the dominated convergence theorem. However, the proof of this result lies in
measure theory that is out of scope for this class, so we will take the liberty of assuming necessary
EDI conditions without justification.5
We will define a couple core regularity conditions on a model {fθ (y) : θ ∈ Θ}, just for conve-
nience to prove facts about these models.

• (A.1). The parameter space Θ is an open set in Rn .

• (A.2). The support of Y , i.e., {y : fθ (y) > 0}, does not depend on θ.

Proposition 6.3. For a model fθ (y) satisfying (A.1) and (A.2), if fθ is differentiable on its support
and the first-order EDI holds, then for all θ ∈ Θ,

Eθ [S(Y, θ)] = 0.

Proof. Pretty trivial argument, we just apply Definition 6.2 and notice that the total integral of
fθ (y) over all y is 1, so it does not change based on the value of θ.

Note. There are some interesting structural parallels between this fact, Eθ [S(Y, θ)] = 0, and the
maximum likelihood estimator, which satisfies S(y, θ̂) = 0.
5
Lucas says that in his research, which is statistical theory rather than probability theory, he almost always
assumes the ∞-order EDI without justification. It is very rare for this condition to break.

16
Corollary 6.3.1. If θ is 1-dimensional and the conditions in the previous proposition hold, then
2

Varθ [S(Y, θ)] = Eθ S (Y, θ) . (Analogous results hold for the covariance matrix of the score
function when θ is higher-dimensional.)

This corollary motivates the next topic of our lecture, which is a measure of how much infor-
mation a random variable Y carries about its parameter θ.

6.2 Fisher Information


Intuitively, if a random variable carries a lot of information about its parameter, then its log-
likelihood will change greatly between different values of θ. This means the score function will have
high variance, which leads to the following.

Definition 6.4 (Fisher information). The Fisher information for a model Y ∼ fθ (y) is a function
of the parameter θ given by
I(θ) = Eθ S 2 (Y, θ) .
 

Proposition 6.5. For a model fθ (y) satisfying (A.1) and (A.2), such that fθ is twice-differentiable
with respect to θ on its support and satisfies the second-order EDI,
 2 

I(θ) = − Eθ log fθ (Y ) .
∂θ2

Proof. Observe that by the chain rule and product rule,


"  #
fθ00 (Y ) fθ (Y ) 2
 0
I(θ) = Eθ − .
fθ (Y ) fθ (Y )

By moving the derivative out of the integral sign, the first term goes to zero. Meanwhile, the second
term equals the Fisher information, so we conclude.

Corollary 6.5.1. Consider any θ ∈ Θ, with two models Y1 ∼ fθ (y) and Y2 ∼ gθ (y), and Y1 ⊥ ⊥ Y2 .
If I1 is the Fisher information of Y1 , and I2 is the Fisher information of Y2 , then (Y1 , Y2 ) has
Fisher information I1 + I2 .

Now, let’s go through some examples to get a feeling for the Fisher information metric.

Example 6.6. If Y1 , . . . , Yn ∼ Pois(λ), then


 2  Y1 −λ 
∂ λ e 1
I1 (λ) = − Eλ 2
log = .
∂λ Y1 ! λ

Therefore, the Fisher information is In (λ) = nλ .

Example 6.7. For a location family where fθ (y) = f (y − θ), the Fisher information metric is a
constant function. For a scale family where fθ (y) = θ−1 f ( yθ ), a similar calculation tells us that I(θ)
is proportional to θ−2 .

Finally, we conclude with the big result about Fisher information, which is a lower bound on
how good an unbiased estimator for a parameter can be.

17
Theorem 6.8 (Cramér-Rao bound). Given a model {fθ (y) : θ ∈ Θ} such that (A.1) and (A.2)
hold, let g(θ) be a differentiable parametric function, and let T (Y ) be an unbiased estimator for
g(θ). Then, if fθ is differentiable on its support and satisfies the first-order EDI, and I(θ) > 0,

(g 0 (θ))2
Varθ [T ] ≥ .
I(θ)

A similar result holds for the case when θ is a vector, but it is an inequality between matrices.

Proof. First, note that

Covθ (S(Y, θ), T (Y )) = Eθ [S(Y, θ)T (Y )] − Eθ [S(Y, θ)] Eθ [T (Y )]


Z
∂ log fθ (y)
= T (y)fθ (y) dy
∂θ
Z
∂ fθ (y)
= T (y) dy
∂θ
Z

= T (y)fθ (y) dy
∂θ

= Eθ [T (Y )]
∂θ
= g 0 (θ).

Furthermore, as we showed previously, Varθ [S(Y, θ)] = I(θ). Therefore, by the Cauchy-Schwarz
inequality,
Covθ (S(Y, θ), T (Y ))2 ≤ Varθ [S(Y, θ)] Varθ [T (Y )] .
Substituting in our results above, we get

(g 0 (θ))2 ≤ I(θ) Varθ [T (Y )] .

The Cramér-Rao lower bound will give us a goal to aim towards in the case when we want to
find low-variance unbiased estimators. It is not always realizable, but it is still very important, and
it is most of the reason behind why we care about Fisher information.

18
7 September 27th, 2021
Today we will finish discussing the Cramér-Rao lower bound, then introduce method of moments
estimation (MOM) and maximum likelihood estimation (MLE).

7.1 More on Cramér-Rao


Recall from last lecture that the Cramér-Rao lower bound was a way to obtain a lower limit on the
variance of an unbiased estimator of a model parameter, based on the Fisher information of that
parameter. This lower bound is not always satisfied, but there are some cases where it is tight.
2
Example 7.1. If Y1 , . . . , Yn ∼ N (µ, σ 2 ), then Var Y = σn = I(µ)
1
 
.
λ 1
 
Example 7.2. If Y1 , . . . , Yn ∼ Pois(λ), then Var Y = n = I(λ) .

Therefore, the Cramér-Rao lower bound holds for the mean of Gaussian and Poisson random
variables. However, it does not hold for estimating a nonlinear function g(µ) of the mean. It turns
out that maximum likelihood estimation asymptotically achieves the Cramér-Rao lower bound
when the sample size is large, but not for any finite sample size.
The above examples were of unbiased estimators, but there is a variant of the Cramér-Rao lower
bound that also works for biased estiamtors, as long as the bias converges quickly to zero as the
sample size n increases.

Proposition 7.3 (Asymptotic Cramér-Rao). Given an estimator Tn (Y) of g(θ), if


√ d
n(Tn − g(θ)) −
→ N (0, v(θ)),
(g 0 (θ))2
then the asymptotic variance v(θ) ≥ I1 (θ) .

Proof. We omit the proof, but see §6.2 of [LC06] for more details.

Note that the convergence in distribution condition is much weaker than the original assumption
of being an unbiased estimator in Theorem 6.8, since convergence in distribution can tolerate edge
cases such as unlikely, far-out outliers.

7.2 Method of Moments Estimation


Theorizing about lower bounds is interesting, but how might we actually design estimators in prac-
tice? The method of moments (MoM) is a general strategy for obtaining estimators of parameters.
We set sample moments equal to the population moments and solve for θ.

Definition 7.4 (Moment). Given a random variable Y parameterized by θ, the r-th moment of
the distribution of Y is Eθ [Y r ], defined for positive integers r.

Definition 7.5 (Central moment). Given a random variable Y parameterized by θ, the r-th central
moment of the distribution of Y is Eθ [(Y − Eθ [Y ])r ], once again defined for positive integers r.

Sometimes, for heavy-tailed distributions like the Cauchy distribution, k-th moments may not
be defined. In these cases, we may take a transformation of our samples before computing moments,
which is known as the generalized method of moments.

19
Example 7.6 (MOM for N (µ, σ 2 )). The method of moments estimator for the normal distribution
N (µ, σ 2 ) is

µ̂MOM = Y ,
n
2 1X
σ̂MOM = (Yi − Y )2 .
n
i=1
1
Note that the empirical variance appears, but it lacks the more subtle n−1 correction factor.
Example 7.7 (MOM for Binom(k, p)). If Y1 , . . . , Yn ∼ Binom(k, p), then equating the sample
moments with the population moments yields

Y = kp,
n
1 X
Yi2 = E Y12 = kp(1 − p) + k 2 p2 .
 
n
i=1

If we solve these equations for the parameters, we get the method of moments estimators:
2
Y
k̂MOM = 1 Pn ,
Y − n i=1 (Yi − Y )2
p̂MOM = Y /k̂MOM .

Note that these estimators are sketchy, since the difference between the sample mean and raw
sample variance in the denominator could be negative, which gives strange results in that case.
But in general, finding the MLE for the binomial distribution is hard, so the MOM estimator
provides a good starting point.
Example 7.8 (MOM for exponential family). If fθ (y) belongs to an exponential family

fθ (y) = exp{η(θ)T (y) − ψ(η(θ))}h(y),

then differentiating under the integral sign (Definition 6.2) with respect to η yields
Z Z
T (y) exp{ηT (y) − ψ(y)}h(y) dy = ψ 0 (η) exp{ηT (y) − ψ(y)}h(y) dy.

This means that Eθ [T (Y )] = ψ 0 (η), so if we transform the samples of an exponential family by T ,


we get the method of moments estimator by solving:
n
1X
T (Yi ) = ψ 0 (η(θ)).
n
i=1

7.3 Maximum Likelihood Estimation


The maximum likelihood estimator (MLE) of θ is the parameter value that maximizes the likelihood
function L(θ) for the observed data. Equivalently, this also maximizes the log-likelihood `(θ).
Example 7.9 (MLE for exponential family). Once again, suppose that we have an exponential
family with parameter θ, and observe that the log-likelihood function is
n n
!
X X
`Y (θ) = T (Yi ) η(θ) − nψ(η(θ)) + log h(Yi ).
i=1 i=1

20
If we take the derivative with respect to η and set it to zero (to find local extrema), we get the
equation for η given by
n
∂` 1X
= 0 =⇒ T (Yi ) = ψ 0 (η(θ)).
∂η n
i=1

Note that this is the exact same as the equation for the MOM estimator!

Here are some basic properties of the MLE.

Proposition 7.10 (Equivariance of the MLE). If θ̂ is an MLE for θ and τ = g(θ), then τ̂ = g(θ̂)
is an MLE for τ .

Proposition 7.11 (MLE and sufficiency). If there exists a unique MLE for θ, then it is a function
of every sufficient statistic T of θ.

This implies that if the MLE is a sufficient statistic, then it is minimal sufficient. Also, if there
exists a complete sufficient statistic and the MLE is unbiased, then the MLE must be the unique
UMVUE by Theorem 4.7.

Note. The MLE is not necessarily unbiased. For example, the conjugate prior of Y1 , . . . , Yn ∼
α−1
Bern(p) is the Beta distribution Beta(α, β), which has mode α+β−2 corresponding to the MLE, but
α
it has mean α+β , which is not the same.

Finally, we state a useful theorem about maximum likelihood estimation in the context of
consistency. Roughly speaking, this means that under some mild conditions, the MLE θ̂ for a
sample of size n is guaranteed to converge to the true value of θ as n increases to infinity.

Definition 7.12 (Identifiability). A model {fθ : θ ∈ Θ} is identifiable if, for any two θ1 , θ2 ∈ Θ,
we have fθ1 (y) = fθ2 (y) almost surely, then θ1 = θ2 .

Essentially, this means that we can’t have two parameters that produce the same distribution,
since that would make it impossible to distinguish which parameter was true.

Theorem 7.13 (Consistency of MLE). Let Y1 , . . . , Yn ∼ fθ0 , and let θ̂n be the MLE with respect
to the model {fθ : θ ∈ Θ} with θ0 ∈ Θ. Then, assuming the following conditions:

(i) fθ is identifiable,

(ii) The support of fθ does not depend on θ,

(iii) Eθ0 [|log fθ (Y1 )|] < ∞ for all θ ∈ Θ, and

(iv) |Θ| < ∞,

then θ̂n exists and is unique with probability tending to 1 as n → ∞, and it is strongly consistent,
meaning that θ̂n → θ0 almost surely.

21
8 September 29th, 2021
Today we prove the consistency of MLE, discuss the asymptotic normality of MLE, and introduce
the delta method.

8.1 Consistency of MLE


First, we prove Theorem 7.13, which was our main theorem at the end of yesterday’s lecture. Recall
that this means if a statistical model is identifiable and satisfies some specific regularity conditions
(in particular, we need |Θ| < ∞), then the MLE is guaranteed to be strongly consistent as n → ∞.

Proof of Theorem 7.13. First, let `(θ) = log fθ (y) = ni=1 log fθ (yi ). Then,
P

n
1X
`n (θ) = log fθ (yi ).
n
i=1

Also, let `(θ, θ0 ) = Eθ0 [log fθ (Yi )]. By the weak law of large numbers, we know that `n (θ) → `(θ, θ0 )
in probability for all parameters θ ∈ Θ.
Second, note that `(θ, θ0 ) is uniquely maximized by θ, since KL divergence is nonnegative. This
is a result known as Gibb’s inequality, which can be verified by checking that
   
fθ (yi ) fθ (yi )
Eθ0 [log fθ (yi )] − Eθ0 [log fθ0 (yi )] = Eθ0 log ≤ Eθ0 = 0.
fθ0 (yi ) fθ0 (yi )

Combining these two facts, we conclude that Prθ0 (`n (θ) ≥ `n (θ0 )) → 0 for all θ ∈ Θ \ {θ0 }.
Therefore, the probability that θ0 uniquely maximizes the empirical likelihood function converges
almost surely to 1 as n → ∞, as desired.

Note. Unfortunately, the reasoning in the theorem above does not apply to infinite parameter
spaces, since we can’t guarantee that the likelihood maximizer exists. We need to add a couple
more conditions when |Θ| <6 ∞:

(iv) θ0 lies in the interior of Θ as a smooth manifold.

(v) fθ (y) is differentiable with respect to θ, for all θ and almost surely in y.

These conditions imply that the likelihood has a unique maximizer when the score function is zero,
with probability p → 1 as n → ∞.

8.2 Asymptotic Normality of MLE


The second result today will be a stricter set of conditions that provide a central limit theorem-like
result for the maximum likelihood estimator as n → ∞. Unlike our previous result, which only had
four regularity conditions, this one will have seven conditions!

Theorem 8.1. Let Y1 , . . . , Yn ∼ fθ0 , and assume that {fθ : θ ∈ Θ} satisfies:

(i) θ0 lies in the interior of Θ.

(ii) fθ is identifiable.

(iii) The support of fθ does not depend on θ,

22
(iv) For all θ ∈ Θ, log fθ (y) is almost surely three times differentiable in θ.

(v) There exists a function M (y), possibly dependent on θ0 , such that for all θ in a neighborhood
of θ0 ,

∂ 3 log fθ (y)
≤ M (y), and
∂θ3
Eθ0 [M (Y )] < ∞.

(vi) The score function is zero at θ0 , i.e.,


" #
∂ log fθ (Y )
Eθ0 = 0.
∂θ θ=θ0

(vii) The Fisher information satisfies


" #
∂ 2 log fθ (Y1 )
0 < I1 (θ0 ) = − Eθ0 .
∂θ2 θ=θ0

Then, there exists a consist zero θ̂n of the score function Sn , such that Sn (θ̂n ) = 0. As n → ∞,


 
d 1
n(θ̂n − θ0 ) −
→ N 0, .
I1 (θ0 )

Proof. First, we will show the existence of the consistent root θ̂n . The key in this proof is to take
a Taylor expansion of the score function Sn . A second-order Taylor expansion with Lagrange form
of the remainder tells us that
1
Sn (θ̂n ) = Sn (θ̂n ) + Sn0 (θ̂n )(θ̂n − θ0 ) + Sn00 (θ̂n∗ )(θ̂n − θ0 )2 ,
2

where θ̂n∗ is some value between θ0 and θ̂n . Roughly speaking, we can “solve” for θ̂n − θ0 to get

√ Sn (θ0 )/ n
n(θ̂n − θ0 ) = 1 1 00 ∗
.
− n Sn0 (θ0 ) − 2n Sn (θn )(θ̂n − θ0 )

The numerator of this fraction converges in distribution to N (0, I1 (θ0 )) by the central limit theorem,
and the first term − n1 Sn0 (θ0 ) of the denominator converges in probability to I1 (θ0 ) by the weak law
of large numbers.
At this point, we would almost be done by applying Slutsky’s theorem to this quotient, but
1 00 ∗
there is one more issue: what about the last quadratic term − 2n Sn (θn )(θ̂n − θ0 )? We need to
bound this last term appropriately, and that is the crucial step in making this proof rigorous. We
omit the rest of the proof in these notes, but you can finish the argument by a combination of fairly
technical steps that concludes in an application of Slutsky’s theorem.

Lucas states that this is one of the most important results that we will improve in the class.
Although only stated above for scalar θ because the notation is more convenient, it also holds in
generality for cases where θ is multi-dimensional, by extending the conditions in a straightforward
manner with linear algebra.

23
9 October 4th, 2021
Today we discuss the Delta method in asymptotics and introduce hypothesis testing.

9.1 The Delta Method


First, a remark about Theorem 8.1. Note that the normality of the MLE implies that it is unbiased,
and furthermore, the variance is approximated by (nI1 (θ0 ))−1 as n → ∞, which matches with the
Cramér-Rao lower bound, meaning that θ̂n is asymptotically optimal.
Definition 9.1 (Asymptotic relative efficiency). If Wn and Vn be two estimators of a parametric
function g(θ) such that
√ d 2
n(Wn − g(θ)) − → N (0, σW ),
√ d
n(Vn − g(θ)) − → N (0, σV2 ),
2 /σ 2 .
then the asymptotic relative efficiency (ARE) of Vn with respect to Wn is σW V

The normality of the MLE for the parameter θ̂ is a nice result, especially since we get the
asymptotic variance from the information Fisher information metric. If we apply a function to θ̂,
the result g(θ̂) is also an MLE for g(θ), but what is its distribution? It turns out that it is also
normally distributed, due to the following fact.
√ d
Theorem 9.2 (Delta method). If n(Tn − θ) − → N (0, v(θ)), then if g is continuously differentiable
at θ and g 0 (θ) 6= 0,
√ d
→ N 0, v(θ)(g 0 (θ))2 .

n(g(Tn ) − g(θ)) −
Proof. We offer a proof sketch. Take the first-order Taylor expansion of g(Tn ) to get
g(Tn ) = g(θ) + g 0 (θ∗ )(Tn − θ),
for some mean-value term specifying θ∗ between θ and Tn . Using the continuous mapping theorem

and Slutsky’s theorem, we can show that n(g(Tn ) − g(θ)) converges to the same distribution as
√ 0
ng (θ)(Tn − θ), so we conclude.
Example 9.3. If Y1 , . . . , Yn ∼ Bin(k, p), for fixed k, then the MLE for p is p̂ = Y /k. This has
asymptotic variance p(1−p) 1
n , which is consistent with the Fisher information I1 (p) = p(1−p) . Now, if
we want to estimate the odds ratio g(p) = p/(1 − p), note that g 0 (p) = 1/(1 − p)2 , so

 
d p
n(g(p̂) − g(p)) −
→ N 0, .
(1 − p)3
There are some issues with the method presented above. One is the case when we want to esti-
mate a parameter plus error bars by one standard deviation, which is hard because the asymptotic
variance is a function of the parameter, which is unknown. In this case, we can usually just plug
in our MLE for the parameter and get a result that still converges in distribution.
Alternatively, we can apply a variance-stabilizing transformation to the parameter, which is
essentially a function h(θ) p
such that the asymptotic 0 2
0
R p variance (h (θ)) /I1 (θ) is invariant as θ changes.
This means that h (θ) ∝ I(θ), so h(θ) ∼ I(θ). For the binomial parameter, this variance-

stabilizing transformation is sin−1 p.
The second issue is that the rate of convergence to normal also depends heavily on p. For
example, in a Binomial distribution, the sample mean is approximately Poisson for small values of
the parameter p, which has a heavy skew between the lengths of the two tails.

24
9.2 Hypothesis Testing
In hypothesis testing, we typically test a null hypothesis H0 against an alternative hypothesis H1 .
There are several kinds of hypotheses for a scalar parameter θ; here are some examples:

• H0 : θ = θ0 versus H1 : θ = θ1 . (Point null versus point alternative.)

• H0 : θ = θ0 versus H1 : θ > θ1 . (Point null versus one-sided alternative.)

• H0 : θ = θ0 versus H1 : θ 6= θ1 . (Point null versus two-sided alternative.)

• H0 : θ ≤ θ0 versus H1 : θ > θ1 . (Composite null versus composite alternative.)

• H0 : θ = θ0 versus H1 : θ ∈ R. (Point null within a composite alternative.)

These types of tests can all be useful in different scientific scenarios. Generally, a test partitions the
sample space into an acceptance region A, where H0 is accepted, and a rejection region R, where
H0 is rejected. Also, we will usually partition the space of a sufficient statistic T (Y), rather than
dealing with the details of Y directly.

Example 9.4 (Z-test). If Y1 , . . . , Yn ∼ N (µ, 1), then one test for H0 : µ = µ0 versus H1 : µ 6= µ0

would be to use the statistic T (Y) = Y and decision rule to reject H0 if n(Y − µ0 ) > 2.

A test is judged by two kinds of error. Type-I error (false positive) is when H0 is rejected when
it is true, while Type-II error (false negative) is when H0 is accepted despite being false. The power
of a test is the probability of rejecting H0 when it is false, and it is denoted β(θ) = Prθ (T (Y) ∈ R).

25
10 October 6th, 2021
Today, we will continue discussing hypothesis testing and introduce the notion of significance levels,
the goal of maximizing power, and offer theoretical analysis of the most powerful tests.

10.1 The Neyman-Pearson Lemma


As before, consider the same problem of hypothesis testing, with two hypotheses H0 and H1 . Let
Θ0 be the set of parameters in the null hypothesis H0 , and let Θ1 be the parameter space of
the alternative hypothesis H1 . A common goal among hypothesis tests is to guarantee that the
probability of Type-I error is bounded by a certain significance level α, meaning that
sup Pr(T (Y) ∈ R) ≤ α.
θ∈Θ0 θ

(For example, a common value in some scientific fields is α = 0.05.) The value on the left-hand
side is called the size of the test, and we generally try to achieve equality between the size and
significance level to maximize power.
Definition 10.1 (UMP). A test is called universally most powerful (UMP) for a significance level
α if it has the maximum power for all θ ∈ Θ1 , among all hypothesis tests of level α.
The following famous lemma is used to construct UMP tests.
Theorem 10.2 (Neyman-Pearson lemma). Suppose that we have i.i.d. samples Y1 , . . . , Yn ∼ fθ (y),
and suppose that we are testing the null hypothesis H0 : θ = θ0 against an alternative hypothesis
H1 : θ = θ1 , where θ0 6= θ1 . The rejection region
 
fθ1 (Y)
R= Y: ≥c
fθ0 (Y)
is the most powerful level α test if c satisfies the size-α condition
Pr(Y ∈ R) = α.
θ0

Proof. Intuitively, the way to think about this statement is that the maximum-power estimator is
simply the one that rejects when the likelihood ratio between H0 and H1 exceeds some constant
ratio. This makes sense because by adding the maximum-likelihood ratio points to the rejection
region, you get the most marginal power for a given test size. We leave the details of the proof as
an exercise; it simply involves moving around some integral inequalities in the proper way.
The NP lemma generalizes to randomized hypothesis tests, when there is no threshold c that
exactly has size α. Instead, we can interpolate between the two nearest values to α. Assume that
f (Y)
there are c1 , c2 such that fθθ1 (Y) lies in (c1 , c2 ) with probability zero assuming H0 , and
2
   
fθ1 (Y) fθ1 (Y)
Pr ≥ c1 = α1 > α > α2 = Pr ≥ c2 .
θ0 fθ2 (Y) θ0 fθ2 (Y)

Then, the maximum power level α test rejects whenever the likelihood ratio is ≥ c2 , accepts when
it is < c1 , and rejects with probability αα−α 1
1 −α2
when it is equal to c1 .
Example 10.3 (Z-tests are UMP). Consider a model of Y1 , . . . , Yn ∼ N (µ, 1). If our hypotheses
are H0 : µ = µ0 , and H1 : µ = µ1 for some µ1 > µ1 , then the likelihood ratio is monotone increasing
in Y , so the NP lemma tells us that the MP test is a simple threshold Y ≥ c for some c. Similarly,
this inequality would be reversed if µ1 < µ0 , and it can be easily generalized to the case when the
alternative hypothesis is a one-sided composite H1 : µ ≥ µ1 .

26
10.2 Testing Composite Hypotheses
We’ve seen how to construct MP tests for general point hypotheses using Theorem 10.2, but what
about UMP tests for composite alternative hypotheses? In the one-sided alternative case of N (µ, 1),
a Z-test is also UMP, as we discussed above. However, for a two-sided hypothesis H1 : µ 6= µ0 , the
UMP does not exist, since the MP test differs based on µ ∈ Θ1 = R \ {µ0 }.

Definition 10.4 (Monotone likelihood ratio). The family of distributions {fθ (y) : θ ∈ Θ} has
monotone likelihood ratio (MLR) in a statistic T (Y) if the ratio fθ2 (Y)/fθ1 (Y) can be expressed
as a function of θ1 , θ2 , T (Y), and for each θ1 < θ2 , the ratio is non-decreasing in T (Y) when at
least one of the numerator and denominator is positive.

Note. If a model has MLR in T (Y), then T is a sufficient statistic, since for any fixed θ0 ∈ Θ,

fθ (Y)
fθ (Y) = · fθ0 (Y).
fθ0 (Y)

Then, we conclude that T is sufficient by Theorem 2.6.

Example 10.5. An exponential family fθ (y) = exp{T (y)η(θ) − ψ(η(θ))}h(Y ) has MLR in the
sufficient statistic T (Y) if the natural parameter η(θ) is a non-decreasing function of θ.

Here’s the punchline. If we want to test a composite hypothesis with composite null, of the
form H0 : θ ≤ θ0 versus H1 : θ > θ0 , then the following theorem gives us a way to find a UMP.

Theorem 10.6 (Karlin-Rubin). Consider testing H0 : θ ≤ θ0 against H1 : θ > θ0 . If the family


{fθ (y) : θ ∈ Θ} has MLR in some statistic T (Y), then for any t0 , the test with rejection region
T > t0 is UMP for level α, where α = Prθ0 (T > t0 ).

Proof. See Chapter 8 of [CB21] for a proof.

27
11 October 13th, 2021
Today we continue discussing composite hypotheses and the likelihood-ratio test.

11.1 Karlin-Rubin Test


Similar to how the Neyman-Pearson lemma generalizes to randomized tests when there is no pre-
cise size-matches-power cutoff (happening when the test statistic T is discrete), Karlin-Rubin also
generalizes to randomized tests in the same way.
Example 11.1. Let Y1 , . . . , Yn ∼ Bern(p), and suppose that we want to test H0 : p ≤ p0 against
the alternative hypothesis H1 : p > p0 , using the Karlin-Rubin lemma. First, suppose that we are
doing a simple point test of p0 versus p1 > p0 , where the Neyman-Pearson test statistic would be
 T 
1 − p0 T

p1
,
p0 1 − p1
Pn
where T is the natural sufficient statistic T (Y) = i=1 Yi . Notice that the left-hand side is
increasing with p1 for any p1 > p0 , so we conclude that T has the monotone likelihood ratio
property (Definition 10.4). Therefore, the Karlin-Rubin lemma states that the most powerful test
of its size for H0 : p ≤ p0 has rejection region T ≥ c for some threshold c. This test is most powerful
for the specific level
n  
X n j
α= p (1 − p0 )n−j .
j 0
j=c+1

However, given a desired level α, there may not be a value of c that satisfies this equation exactly,
since c is a discrete variable. In this case, we can do a randomized test by finding c such that
n   n  
X n j n−j
X n j
p0 (1 − p0 ) <α< p (1 − p0 )n−j .
j j 0
j=c+1 j=c

Then, our test rejects when T > c, accepts when T < c, and when T = c, it rejects with probability

α − nj=c+1 nj pj0 (1 − p0 )n−j


P 
n c

n−c
.
c p0 (1 − p0 )

11.2 Likelihood-Ratio Test


Consider the problem of testing a null hypothesis H0 : θ ∈ Θ0 against an alternative hypothesis
H1 : θ ∈ Θ1 . Assume that Θ0 ∩ Θ1 = ∅ and Θ0 ∪ Θ1 = Θ. The likelihood-ratio test compares these
statistical models based on the ratio of their likelihoods, i.e., with statistic
supθ∈Θ L(θ)
Λ= .
supθ∈Θ0 L(θ)

Here, we reject the null hypothesis when Λ > c for some threshold c. In other words, if θ̂0 is the
likelihood maximizer for L(θ) over Θ0 , and θ̂ is the global maximizer over Θ, then the likelihood
ratio is simply
f (Y)
Λ = θ̂ .
fθ̂0 (Y)
Note that Λ ≥ 1, since θ̂ is the MLE for θ.

28
Example 11.2. Consider two parameters θ, µ > 0 and corresponding samples of i.i.d. random
variables, X1 , . . . , Xn ∼ Expo(θ) and Y1 , . . . , Ym ∼ Expo(µ). Our null hypothesis is H0 : µ = θ,
while our alternative hypothesis is H1 : µ 6= θ. This is a complex hypothesis test, as we have
multiple parameters and a difficult parameter region Θ0 , which makes it difficult to apply results
like Karlin-Rubin.
Instead, we will try to directly use a likelihood-ratio test, comparing the log-likelihoods. Observe
that the log-likelihood can be written as
Pn Pm
j=1 Xj j=1 Yj
`(θ, µ) = −n log θ − − m log µ − .
θ µ
We can immediately see that the maximum likelihood estimator over the entire parameter space is
θ̂ = X, and µ̂ = Y . Meanwhile, it can be shown that the maximum likelihood estimator in the null
hypothesis parameter space Θ0 is
Pn Pm
j=1 Xj + j=1 Yj
θ̂0 = .
n+m
Therefore, our likelihood-ratio statistic is
L(θ̂, µ̂) θ̂0n+m nn mm
Λ= = = T −n (1 − T )−m ,
L(θ̂0 , θ̂0 ) θ̂n µ̂m (n + m)n+m
Pn
Xj
j=1P
where T = Pn
Xj + m
Yj
. The log-likelihood ratio is now
j=1 j=1

log Λ = const − n log T − m log(1 − T ).


If we graph this function on T ∈ (0, 1), notice that it blows up to ∞ at either side of the interval,
n
and it achieves its minimum when T = n+m . Therefore, instead of using the likelihood-ratio test
directly, we might prefer a test of the form
 
n
R= T− >α .
n+m
To determine α, we could just simulate T , whose distribution is θ-free under the null hypothesis,
by ancillarity of the quotient of random variables in a scale family.
Note. In general, we often use the above formulation of the likelihood-ratio test when Θ is a
higher-dimensional space, and the alternative hypothesis Θ1 is a lower-dimensional subset of that
space. For instance, in the example above, µ = θ forms a one-dimensional subset of R2 .
Naturally, the next question is about the asymptotic distribution of the likelihood-ratio statistic,
since this informs our choice of the threshold α. Just like we can prove asymptotic normality of
the MLE as n → ∞, the likelihood-ratio statistic also has an asymptotic distribution.
Theorem 11.3 (Asymptotic chi-squared distribution of likelihood ratio). Consider testing H0 :
θ ∈ Θ0 versus an alternative hypothesis H1 : θ ∈ Θ \ Θ0 . Let
supθ∈Θ Ln (θ)
Λn =
supθ0 ∈Θ0 Ln (θ0 )
be the likelihood-ratio statistic for i.i.d. sampled data of size n from {fθ (y) : θ ∈ Θ}. Then, under
the same conditions as Theorem 8.1 and that the global maximum of Ln (θ) is one of the consistent
roots of the score equation, we have under H0 that
d
→ χ21 ,
2 log Λn −

29
Proof. First, observe that we can write

2 log Λn = −2(`n (θ0 ) − `n (θ̂n )).

Here, we use `n for the log-likelihood and θ̂n for the maximum likelihood estimator that is a
consistent root of the score equation, which is similar to the notation in Theorem 8.1. A second-
order Taylor series expansion of `n (θ0 ) around θ̂n yields
1
`n (θ0 ) = `n (θ̂n ) + (θ0 − θ̂n )2 `00n (θn∗ ),
2

for some θn∗ between θ0 and θ̂n . Plugging this in, we get
 
1 00 ∗
2 log Λn = −(θ0 − θ̂n )2 `00n (θn∗ ) 2
= n(θ̂n − θ0 ) − `n (θn ) .
n

Notice that the difference between the maximum likelihood estimator θ̂n and θ0 is asymptotically
√ d 1
normal by Theorem 8.1, i.e., n(θ̂n − θ0 ) −
→ N (0, I1 (θ 0)
), so by the continuous mapping theorem,

d 1
n(θ̂n − θ0 )2 −
→ χ2 .
I1 (θ0 ) 1

Finally, we can take care of the last term by using a combination of Taylor expansion and the weak
p
law of large numbers to get that − n1 `00n (θn∗ ) →
− I1 (θ0 ) under H0 , so we arrive thereafter at the desired
result by Slutsky’s theorem.

30
12 October 18th, 2021
Today, we discuss asymptotic and non-parametric hypothesis tests.

12.1 Asymptotic Hypothesis Tests


Recall that we showed the likelihood-ratio test statistic converges to an asymptotic χ2 -distribution.
There is another way of thinking about hypothesis tests in terms of how quickly certain test statistics
asymptotically converge to their limit distributions.
Definition 12.1 (Score test). Suppose that we are testing a null hypothesis H0 : θ = θ0 against
an alternative hypothesis H1 : θ 6= θ0 . The likelihood-ratio test statistic in this case would be
Λn = L(θ̂n )/L(θ0 ). As an alternative, the score statistic has asymptotic distribution
1 d
√ Sn (Y, θ0 ) −
→ N (0, I1 (θ0 )),
n
by the definition of Fisher information as the variance of the score function. Therefore, in the score
test, we reject at level α when
S(Y, θ0 )
p > z1− α2 .
nI1 (θ0 )
Definition 12.2 (Wald test). In the Wald test, observe that

 
d 1
n(θ̂n − θ0 ) −
→ N 0, .
I1 (θ0 )
Therefore, in this hypothesis test, we reject at level α when
p
nI1 (θ0 )(θ̂n − θ0 ) > z1− α2 .

Figure 1: Schematic visualization of asymptotic hypothesis tests.


Visually, we can think of these tests as in terms of a schematic like Fig. 1. Here, the Wald test
statistic is represented by w, the likelihood-ratio test statistic is visualized by log h as a quotient of
likelihoods, and the score test statistic is the slope of the tangent labeled s. Each of these is nearly
asymptotically equivalent as n → ∞ but they have different power.

31
Note. We can generalize each of these asymptotic hypothesis tests to cases when the data and
parameters are multi-dimensional. In these cases, the tests involve the Fisher information matrix,
which is required to be positive definite (in order to be invertible). The limit distributions in these
cases depend on dimensionality of the parameter sets Θ and Θ0 :

• For the LR test, we have


d
→ χ2dim(Θ)−dim(Θ0 ) .
2 log Λn −

• For the Wald test, we have


d
(θ̂ − θ0 )> I(θ0 )(θ̂ − θ0 ) −
→ χ2dim(Θ)−dim(Θ0 ) .

• For the score test, we have


d
[∇`(Y, θ0 )]> I−1 (θ0 )[∇`(Y, θ0 )] −
→ χ2dim(Θ)−dim(Θ0 ) .

12.2 Nonparametric Hypothesis Tests


Before talking about non-parametric hypothesis tests, we need to introduce a rigorous notion of
p-value, for the purposes of this course.

Definition 12.3 (p-value). For a family of hypothesis tests of size α that defines a negation region
Rα for any α ∈ [0, 1], such that Rα1 ⊆ Rα2 for any α1 ≤ α2 , the p-value for a data set X is the
smallest significance level α such that Y ∈ Rα .

Most hypothesis tests are of this form, and they reject based on a test statistic threshold
T (Y) ≥ C(α), for some increasing non-linear function C. Then, the p-value of a test statistic T
can simply be computed by taking

p = min{α : T (Y) ≥ C(α)} = C −1 (T (Y)).

In general, we know that supθ∈Θ0 Prθ (p ≤ α) = supθ∈Θ0 Prθ (Y ∈ Rα ) ≤ α, for all α, since the size
of a test is less than or equal to its significance level. This is the common way that hypothesis
tests are formulated in applied statistics. Furthermore, if the size equals the level for all α, which
is fairly common, we can also say more precisely that
d
p−
→ Unif[0, 1].

Unlike our previous hypothesis tests in this course, note that by writing p-values, we effectively
eliminate the parameter from our mathematical expressions. Here are some considerations to keep
in mind when thinking about non-parametric tests:

• Weak assumptions: Non-parametric tests can easily support very general null hypotheses.
For example, consider data distribution Y ∼ N (θ, 1). While a standard parametric test could
compare H0 : θ = 0 against H1 : θ 6= 0, non-parametric tests could also test hypotheses like
H0 : Y ∼ N (0, 1), or H0 : E [Y ] = 0.

• Power: Under parametric assumptions, we lose a bit of power in non-parametric tests com-
pared to their counterparts. However, non-parametric tests can gain relative power when the
parametric model fails.

32
• Discrete statistics: Non-parametric tests tend to have discrete test statistics.

Let’s give some concrete examples of these hypothesis tests now. Most of these examples center
around some “maximally ancillary” statistic for the null hypothesis.

Example 12.4 (Sign test). Suppose that we have random variables Y1 , . . . , Yn ∈ R sampled
i.i.d. from some data distribution. Then, the sign test evaluates the likelihood of the null hy-
pothesis H0 : Pr(Y > θ0 ) = Pr(Y < θ0 ). The test stiatistic in this case is
n  
X 1
1{Yi >θ0 } ∼ Bin n, under H0 .
2
i=1

The test statistic concentrates around the mean, so we reject when the statistic is too far away
from n2 , and our size is    
1 n
Pr Bin n, − ≥ c ≤ α.
2 2

33
13 October 20th, 2021
Today we will finish discussing the sign test, then introduce a bunch of other famous hypothesis
tests: the Wilcoxon signed-rank test, the Kolmogorv-Smirnov test, the Mann-Whitney U test, and
the permutation test.

13.1 Sign Test and Signed-Rank Test


Recall that the sign test asks whether the median of the sampled data Y is equal to θ0 . This test
has statistic distributed according to Bin(n, 12 ) under the null hypothesis, which is given by
n
X
T (Y) = 1{Yi >θ0 } .
i=1
What if Y is a discrete variable, so there may be positive probability mass on the median point
Pr(Y = θ0 ) > 0, and we would like to write a sign test for whether the median equals θ0 ? In this
case, instead of our single additive sign statistic, we could imagine describing a trinomial statistic
sgn(Yi − θ0 ) ∈ {−1, 0, 1} and test for
Pr(sgn(Yi − θ0 ) = −1) ≤ 0.5,
Pr(sgn(Yi − θ0 ) = 1) ≤ 0.5.
Alternatively, if we wanted to instead test if the two sides of θ0 were exactly balanced in mass,
there are a couple ways that we could implement such a test in practice with just one statistic:
• Ignore all samples that are equal to θ0 , testing the fact that under H0 , since
Pr(Y > θ0 | Y 6= θ0 ) = Pr(Y < θ0 | Y 6= θ0 ).

• For each sample Yi = θ0 , flip a fair coin to determine its sign, in order to break ties.
Another simple variant of the sign test can be derived by applying it to paired samples.
Definition 13.1 (Paired sample sign test). Given i.i.d. samples (X1 , Y1 ), . . . , (Xn , Yn ) ∼ D, the
paired sample sign test has null hypothesis
H0 : Pr(Yi − Xi > θ0 ) = Pr(Yi − Xi < θ0 ).
The test statistic in this case would be T (X, Y) = ni=1 1{Yi −Xi >θ0 } .
P

Next, suppose we wanted to test a null hypothesis H0 stating that the data is sampled from a
symmetrical distribution around θ0 . In other words, the variable |Yi − θ0 | is conditionally indepen-
dent from the event Yi − θ0 > 0. This can be handled with the following test.
Definition 13.2 (Wilcoxon signed-rank test). Given data Y1 , . . . , Yn ∼ D, the Wilcoxon signed-
rank test for continuous random variables has test statistic
Xn
W = sgn(Yi − θ0 )Ri ,
i=1

where Ri = |Yi −θ0 | for each i. Under the null hypothesis that D is a symmetric distribution around
θ0 , each sign is a Rademacher random variable that is independent from Ri ∈ Unif{1, . . . , n} which
are integers drawn without replacement, so
n(n + 1)(2n + 1)
E [W ] = 0, Var [W ] = 12 + · · · + n2 = .
6
Furthermore, by the central limit theorem, W is asymptotically normal under H0 .

34
Note. In the discrete case of Wilcoxon’s test, there may be cases when the ranks Ri are ambiguous
due to ties. To fix this, we can either use a randomized tie-breaking strategy (which preserves
the distribution of W but is non-deterministic), or correct for ties by setting each rank Ri to the
average of tied ranks (which is deterministic but changes the distribution of W ).

13.2 Two-Sample Tests


Now, let’s discuss two-sample hypothesis tests. Here, we have some data X1 , . . . , Xn ∼ FX and
Y1 , . . . , Ym ∼ FY , both samples i.i.d. and satisfying X ⊥⊥ Y. Our null hypothesis will be that the
two distrbutions are in fact the same, i.e., FX = FY .

Definition 13.3 (Kolmogorov-Smirnov test). The Kolmogorv-Smirnov (K-S) test tests distribu-
tional null hypotheses by comparing the cumulative distribution functions. There are two major
variants of K-S, for the one-sample and two-sample cases.

• One-sample variant: Given a distribution F0 and samples Y1 , . . . , Yn ∼ F , the one-sample


K-S test has null hypothesis H0 : F = F0 . The test statistic relies on the cumulative distri-
bution functions,
sup F0 (y) − F̂Y (y) ,
y∈R
Pn
where F̂Y (y) = n1 i=1 1{Yi ≤y} is the empirical cumulative distribution function. Notice that
this test statistic does not depend on F0 when Y has a continuous distribution, since by a
change of variables,
n
1X
sup F0 (y) − F̂Y (y) = sup ŷ − 1{F0 (Yi )≤ỹ}
y∈R ŷ∈(0,1) n
i=1
n
1X
= sup ŷ − 1{Ui ≤ỹ} ,
ŷ∈(0,1) n
i=1

where Ui = F0 (Yi ) ∼ Unif[0, 1]. Visually, we can imagine taking an empirical estimate of the
CDF for Y as a staircase function, then finding the maximum vertical deviation between this
estimate and F0 .

• Two-sample variant: In the two-sample variant, we simply compare the empirical cumu-
lative distribution functions of both variables, to get test statistic

sup F̂X (y) − F̂Y (y) .


y∈R

Here, we are comparing two staircase functions against each other. Once again, for continuous
random variables, this statistic doesn’t depend on FX and FY under the null hypothesis H0 .

Note that both of these variants have a test statistic that is asymptotically distributed according
to the maximum absolute deviation of a Brownian bridge process on [0, 1], which makes it possible

to analytically compute cutoffs for a desired confidence level α, around Θ(− log α/ n).

Although the Kolmogorv-Smirnov test is extremely general and applies in many cases, it is not
as powerful as a test that is targeted towards specific alternative distributions, if more information is
known about the data. Next, we’ll introduce a more specific test that is similar to the paired-sample
Wilcoxon signed-rank test, but it assumes independence of X and Y .

35
Definition 13.4 (Mann-Whitney U test). Assume that we have independent data X1 , . . . , Xn ∼
DX and Y1 , . . . , Ym ∼ DY . The null hypothesis H0 states that DX = DY , and the Mann-Whitney
U test statistic is
n
X n(n + 1)
U= R(Xi ) − .
2
i=1
nm(n+m+1)
Under the null hypothesis, E [U ] = nm
2 and Var [U ] = 12 . This is asymptotically normal as
n, m → ∞. Here, the alternative hypothesis for this test would be that Y stochastically dominates
X, or vice versa, in that one of their CDFs is strictly greater than the other.

We have a midterm exam next week, and review materials will be posted on the course website.

36
14 October 25th, 2021
Today we discuss the permutation test, a general family of statistical tests for independence of
paired draws from a bivariate distribution, as well as selective inference.

14.1 Tests of Independence


Suppose that we have i.i.d. random samples (X1 , Y1 ), . . . , (Xn , Yn ) drawn from an unknown bivariate
distribution. We want to test the null hypothesis H0 , that Xi ⊥⊥ Yi are independent, against the
alternative hypothesis that that they are not independent.
In the permutation test, the essential observation is that when X ⊥⊥ Y , we can permute the
samples according to some arbitrary permutation σ, taking

(X1 , . . . , Xn , Y1 , . . . , Yn ) 7−→ (Xσ(1) , . . . , Xσ(n) , Y1 , . . . , Yn ).

Notice that in the null hypothesis, the likelihood of any sample is identical to that of any of its n!
permuted variants. Therefore, for any statistic T : R2n → R, the distribution of T (X, Y) is the
same as the distribution of its permuted variants.
Definition 14.1 (Permutation test). Given random samples (X1 , Y1 ), . . . , (Xn , Yn ) from an un-
known bivariate distribution, and a test statistic T : R2n → R, we define the following for any
permutation π of {1, . . . , n}:

T π = T (Xπ(1) , . . . , Xπ(n) , Y1 , . . . , Yn ).

Then, if Πn is the set of all permutations of length n, and id ∈ Πn is the identity element, the
permutation test p-value is computed as
1 X
1{T π ≥T id } .
n!
π∈Πn

Notice that the p-value lies in the range [1/n!, 1], and it is conservative. In other words, under H0 ,
the rejection region Rα where the p ≤ α has likelihood PrH0 (Rα ) ≤ α.
Note that the above test requires computing n! statistics T π for each π ∈ Πn , which is intractable
for most values of n. Instead, we can generate estimate p-values given N ≤ n! permutations, denoted
π1 , . . . , πN , sampled uniformly from Πn . Here, the p-value is
 
N
1  Y
1+ 1{T πj ≥T id } .
N +1
j=1

In the above expression, the extra addition of 1 replaces the always-true 1{T id ≥T id } term of the full
permutation test, which makes this a conservative p-value.
Notice that the permutation test is very flexible, since T can be any function without restric-
tions. If there is some expected relationship (e.g., quadratic relationship) between X and Y in the
alternative hypothesis, T can be set to the regression coefficient of a quadratic fit, which would be
much higher in the original data than the shuffled data.
Note. Another way of viewing the permutation test is in terms of sufficient statistics. The minimal
sufficient statistic for general univariate probability distributions is simply the order statistics of
the sampled data. Therefore, we can condition on the order statistics to get an unbiased estimate
of the joint distribution under the null hypothesis.

37
14.2 Selective Inference
Now we introduce selective inference, which is the problem of soundly testing many hypotheses for
statistical significance on data. It is increasingly common in scientific fields to want to test many
hypotheses at once, but such a procedure can often lead to fallacies if care is not taken. Consider
the following scenarios:
• For a single hypothesis, we could reasonably reject a null hypothesis at significance level 5%.

• However, if we test 20 hypotheses all at this same 5% significance level, then we should expect
to get a false positive at this significance level, even if the null hypothesis H0 is true.

• Even worse, given 20000 human genes, if we did a test at the 5% level, this would generate
around 1000 false positives, which could completely drown out the signal for the few significant
genes that we actually care about.
Therefore, in selective inference, we choose to primarily study not a single hypothesis, but rather a
family of hypotheses. Given m null hypotheses H0,1 , . . . , H0,m , let H0 be the subset of those that
are true null hypotheses, and let m0 = |H0 |. (We sometimes also slightly abuse notation and let
H0 refer to the set of indices of true null hypotheses.)

Figure 2: Notation for selective inference.

Under some hypothesis test, we select some of the null hypotheses for rejection. Here, we let R
be the number of selected hypotheses, and we define quantities V , S, U , and T according to Fig. 2.
Definition 14.2 (Familywise error rate (FWER)). The familywise error rate of a selection rule is
the probability that any of the hypotheses in H0 is selected, or FWER = Pr(V > 0).
There are two kinds of control for the familywise error rate: strong control, where the error is
tracked for any arrangement of true and false null hypotheses, and weak control, where we only
bound the FWER when m = m0 . For the rest of this section, we only discuss strong control.
Definition 14.3 (False discovery rate (FDR)). We define the false discovery proportion of a selec-
tion rule to be
V V
FDP = = ,
V +S R
where 0/0 is taken to be 0. The false discovery rate is the expectation
 
V
FDR = E [FDP] = E ,
R
where 0/0 is once again taken to be 0.
Note. Under the global null hypothesis, the FDR is equal to the FWER, as both are either 0 or 1.
In general, we can observe that FWER ≥ FDR, so the FDR may be easier to control as a weaker
error rate that admits procedures with higher power.

38
Most work in selective inference focuses on taking in p-values p1 , . . . , pm and returning a selected
set of indices that controls a rate like the FWER or FDR, while maximizing power. Here’s a first
example of such a procedure.

Definition 14.4 (Bonferroni procedure). Given m p-values p1 , . . . , pm , the Bonferroni procedure


for controlling the FWER at level α rejects all hypotheses with p-values lower than α/m.

This succinctly captures the core idea that with more hypotheses, we usually need to be stricter
about p-value thresholds being low enough to reject. However, this procedure is relatively primitive,
as it doesn’t make any dependence or independence assumptions among the p-values.

Proposition 14.5. The Bonferroni procedure controls the FWER at level α.

Proof. Simply use a union bound:


X  α  m0
FWER = Pr(V > 0) ≤ Pr pk ≤ = α.
m m
k∈H0

In fact, this means that the procedure is conservative by a factor of m0 /m < 1.

39
15 October 27th, 2021
Today, we continue discussing selective inference procedures that control the FWER and FDR.

15.1 FWER-Controlling Procedures


Recall that the Bonferroni procedure rejects all p-values below α/m, where α is the FWER level
and m is the number of p-values. This does not require independent samples.

Example 15.1. Suppose that we draw m p-values independently from the uniform distribution
1 1
Unif[ m , 1]. Then, we pick one of the p-values uniformly among the m and redraw it from Unif[0, m ].
These values are now marginally Unif[0, 1], and the FWER of the Bonferroni procedure applied to
these p-values, assuming the global null hypothesis, is exactly α.

What about the case when the hypotheses are independent? We can derive a slightly stronger
procedure under this assumption, which has a larger p-value rejection threshold.

Definition 15.2 (Šidák procedure). We reject all hypotheses below αm = 1 − (1 − α)1/m .

Proposition 15.3. The Šidák procedure controls the FWER at level α when the null hypothesis
p-values are independent, and it is tight when m0 = m.

Proof. This is a simple calculation using joint probability of independent events,

FWER = Pr(V > 0) = 1 − Pr(V = 0) = 1 − (1 − α)m0 ≤ 1 − (1 − α)m .

Note that by taking a Taylor series expansion,

α α2
 
α
αm ≈ 1+ + + ··· ,
m 2 3

when m is large and α is small. This means that the p-value in the Šidák procedure is about α/2
better than the Bonferroni procedure in relative terms. For example, when α = 0.05, the Šidák
procedure has a selection threshold that is 1.025× the Bonferroni threshold.
We might be a little bit disappointed, since this procedure requires a strong independence
assumption, yet it doesn’t improve the threshold much. Indeed, the Bonferroni procedure captures
most of the selection power with a simple rule. However, the following selection rule will show that
Bonferroni can be uniformly improved without making any additional assumptions.

Definition 15.4 (Holm procedure). Let p(1) , . . . , p(m) be p-values in increasing order, and let
H0,(1) , . . . , H0,(m) be their corresponding null hypotheses. Then, the Holm procedure rejects nothing
if p(1) > α/m, and otherwise, it rejects H0,(1) , . . . , H0,(k̂) where
 
α
k̂ = min k : p(k+1) > .
m−k

This is called a step-down procedure,6 and it can be thought of as an iterative procedure that
chooses to reject each null hypothesis in turn by comparing its p-value against the current threshold,
assuming all previous null hypotheses are false.
6
The name “step-down” is confusing, since it means the opposite: going in order of increasing p-value.

40
Note that the Holm procedure is strictly less conservative than the Bonferroni procedure, since
it rejects at least as many null hypotheses.

Proposition 15.5. The Holm procedure controls the FWER at level α.

Proof. Let k0 be the sorted index of the smallest null p-value. In other words,

k0 = min{k : H0,(k) ∈ H0 }.

Then, the family-wise error rate can be bounded by


   
α α
FWER = Pr(V > 0) ≤ Pr p(k0 ) ≤ ≤ Pr min pk ≤ .
m − k0 − 1 k∈H0 m0

Finally, by a union bound, this expression is at most


   
α X α
Pr min pk ≤ ≤ Pr pk ≤ = α.
k∈H0 m0 m0
k∈H0

An equivalent way to think about the Holm method is as a non-stepwise procedure, which
applies threshold α/(m − k̂ + 1) to all p-values. Of course, you still need to use the same stepwise
method to actually compute k̂, but it’s an interesting conceptual way of thinking about it.

15.2 FDR-Controlling Procedures


The following procedure for controlling the FDR was introduced in 1995, and it has since become
one of the most cited papers in the statistical literature, ubiquitous throughout science research for
selective inference. It is also interesting from a statistical point of view.

Definition 15.6 (Benjamini-Hochberg procedure). Let q be a level for the FDR. The Benjamini-
Hochberg (BH) procedure rejects nothing if mink P(k) /k > q/m. Otherwise, it rejects the hypotheses
H0,(1) , . . . , H0,(k̂) , where
 
kq
k̂ = max k : p(k) ≤ .
m
This is called a step-up procedure, and it can be viewed as iterating in order of descending p-value
until finding one that is lower than the significance threshold.

Note. Even with the same rejection levels, the behavior of a step-up and step-down procedure
differs based on where the p-value curve intersects the threshold curve.

41
16 November 1st, 2021
Today, we go over the proofs of correctness for selective inference procedures that control the FDR,
and we discuss confidence intervals.

16.1 More on FDR Control


Recall that the Benjamini-Hochberg procedure from last lecture is a step-up procedure, which at
level q, rejects the k̂ smallest p-values, where k̂ is the largest index k such that p(k) ≤ kq/m.
Theorem 16.1. The Benjamini-Hochberg procedure controls the FDR at level q, assuming that the
p-values are independent.
Proof. Let p1 , . . . , pm be p-values, and without loss of generality, let p1 be derived from the null
hypothesis. Let R be the number of p-values rejected by the BH procedure. Observe that p1 is
rejected and R = r if and only if p1 is rejected and R̃ = r, where

R̃ = R({0, p2 , . . . , pm })

is the number of rejected p-values if p1 were replaced with 0. This means that
m
X hv i
FDR = E 1{R=r}
r
r=1
 
m
X 1 X
= E 1{pk ≤ qr } 1{R=r} 
r m
r=1 k∈H0
m
X m0 h i
= E 1{p1 ≤ qr } 1{R=r}
r m
r=1
m
X m0 h i
= E 1{p1 ≤ qr } 1{R̃=r}
r m
r=1
m
X m0  qr   
= Pr p1 ≤ Pr R̃ = r
r m
r=1
m
m0 X
≤q Pr(R̃ = r)
m
r=1
| {z }
=1
m0
=q .
m
Observe that the BH procedure is conservative by a factor of m0 /m, similar to the Bonferroni
procedure. We can also compute the variance of the FDP (not just its mean), which is
     
qm0 1 1 1
Var [FDP] = E − ≤ qE .
m 1 + R̃ m max{R, 1}
p
For example, if q = 0.1 and R ≈ 50, then Var [FDP] ≈ 0.05.

So far, we’ve seen the proof assuming independence of p-values, which is a nice simplifying
assumption. However, it turns out that the BH procedure also controls the FDR (up to a log-
factor) when the p-values are allowed to be arbitrarily dependent, which we will show now.

42
Theorem 16.2. The Benjamini-Hochberg procedure controls the FDR at level
m
!
m0 X 1 m0
q ≈q (log(m) + 0.577),
m i m
i=1

under arbitrary p-value dependence.7

Proof. Once again, we can compute that


m
X hv i
FDR = E 1{R=r}
r
r=1
 
m
X 1 X
= E 1{pk ≤ qr } 1{R=r}  .
r m
r=1 k∈H0

Now, we rewrite this inner sum as


m X
r
" #
X X 1{R=r}
FDR = E 1{ q(i−1) <p qi
r m k≤ m }
k∈H0 r=1 i=1
"m m #
X XX 1{R=r}
= E 1{ q(i−1) <p ≤ qi }
r m k m
k∈H0 i=1 r=i
m
" #
X X 1{R≥i}
= E 1{ q(i−1) <p ≤ qi }
R m k m
k∈H0 i=1
m  
X X 1 q(i − 1) qi
≤ Pr < pk ≤
i m m
k∈H0 i=1
m
!
m0 X 1
=q .
m i
i=1

In the last step, we assumed that the individual p-values were tight up to the size of the tests, but
the same bound would apply to conservative p-values by the rearrangement inequality.

16.2 Confidence Intervals


We now introduce confidence intervals, which are commonly abbreviated CIs.

Definition 16.3 (Interval estimator). Let L(Y) and U (Y) be a pair of functions of the data
Y ∼ fθ such that L(Y) ≤ U (Y) almost surely. Then, [L(Y), U (Y)] is an interval estimator of θ if
we infer that L(Y) ≤ θ ≤ U (Y).

Definition 16.4 (Coverage probability). The coverage of an interval estimator is the probability

Pr(L(Y) ≤ θ ≤ U (Y)).
θ

Our general goal is to design confidence intervals such that the coverage is large.
7
The approximate numerical value γ ≈ 0.577 is the Euler-Mascheroni constant.

43
Definition 16.5 (Confidence coefficient). The confidence coefficient of an interval estimator is
 
inf Pr(L(Y) ≤ θ ≤ U (Y)) .
θ∈Θ θ

Given these definitions, we are now ready to introduce confidence intervals.

Definition 16.6 (Confidence interval). A confidence interval of level α is an interval estimator


with confidence coefficient at least 1 − α.

Example 16.7. If Y1 , . . . , Yn ∼ Unif[0, θ] for θ > 0, and Y( n) = max1≤i≤n Yi , consider the following
two interval estimators:

• [aY(n) , bY(n) ], where 1 ≤ a < b.

• [c + Y(n) , d + Y(n) ], where 0 ≤ c < d.

In the first case, the coverage probability does not depend on θ, as we could just scale the model
appropriately. However, in the latter case, the coverage probability depends strongly on θ.

Example 16.8. If Y1 , . . . , Yn ∼ N (µ, σ 2 ), and we are interested in a 95%-confidence interval, where


α = 0.05, for µ, there are two common methods we can use:

• (z-test). If σ 2 is known, return Y ± z1−α/2 √σn .


Pn
• (t-test). If σ 2 is unknown, return Y ± tn−1,1−α/2 √sn , where s2 = i=1 (Yi − Y )/(n − 1) is the
unbiased sample variance.

The first method generally yields better results, since you have more information.

After setting the level of a confidence interval, the next question to ask is about how much
power we get. We’ll see next time that we should generally prefer confidence intervals that have
shorter length, when both have equal coverage.

44
17 November 3rd, 2021
Today, we will discuss confidence intervals in more detail, including the two most common means
of constructing confidence intervals.

17.1 CIs by Inverting Hypothesis Tests


One primary method of constructing confidence intervals is by using a hypothesis test. Suppose
that for every θ0 in the parameter space, we have a level-α test for the null hypothesis H0 : θ = θ0
versus the alternative H1 : θ 6= θ0 . Furthermore, let A(θ0 ) be the acceptance region for this level-α
test, meaning that
Pr(Y ∈ A(θ0 )) ≥ 1 − α.
θ0

We can then define the set C(Y) = {θ : Y ∈ A(θ)}. Notice that

Pr(θ ∈ C(Y)) = Pr(Y ∈ A(θ)) ≥ 1 − α.


θ θ

Therefore, if the set C(Y) turns out to be an interval (which will often happen in practice), then
we have shown that C(Y) is a (1 − α)-confidence interval for θ.
We have therefore shown that given a general hypothesis test for a point null hypothesis against
the natural alternative, we can naturally arrive at an equivalent confidence interval for that test.
Hypothesis tests are equally as powerful as confidence intervals in these cases; neither is more
informative than the other. However, there are a few limitations to confidence intervals:
• Existence of a confidence interval requires stronger assumptions than a hypothesis test. For
example, calculating a confidence interval implicitly assumes that the model generating the
data is correct, whereas you could reject the null hypothesis in a hypothesis test even if the
model was not correct.

• To generate a confidence interval using this method, you need hypothesis tests to exist for
every value of the parameter θ, which isn’t always easy.

• Hypothesis tests can answer much more general questions, including non-parametric ones.
Confidence intervals only work when we’re analyzing a parameter in a contiguous range.
Sometimes, the confidence sets C(Y) may not be contiguous or useful at all.
A classic example of getting a CI from a hypothesis test is the t-interval from Example 16.8, which
is commonly seen in high school statistics classes. You can derive this interval by applying the
method above with the t-test statistic.
Example 17.1. Let Y1 , . . . , Yn be i.i.d. ∼ θ · Expo, where θ > 0, and we are interested in obtaining
some (1 − α)-confidence interval for θ. First, we need to find a hypothesis test for H0 : θ = θ0
versus H1 : θ 6= θ0 . Since the MLE of θ is just Y , the likelihood ratio test statistic is
−n −n
Y e−n

Y −Y /θ0
Λ = −n ∝ e .
θ0 e−nY /θ0 θ0

Note that we have nY /θ0 ∼ Gamma(n), so the distribution of Λ is free of the parameter θ, and we
can find some cutoff cα such that the acceptance region of the size-α LR test is
 
Y −Y /θ0
A(θ0 ) = Y : e > cα .
θ0

45
Next, observe that the function f (x) = xe−x is unimodal with a maximum at x = 1, so the
(1) (2) (1)
acceptance region is really just an interval, and if cα ≤ cα are the two points where f (cα ) =
(2)
f (cα ) = cα , then a confidence interval for θ would be
" #
Y Y
C(Y) = (2) , (1) .
cα cα

This was a more complicated example of a confidence interval derived from a hypothesis test,
and hopefully it illustrates how the technique is applied in general.

17.2 CIs by Using Pivotal Quantities


The second method of constructing confidence intervals is using pivots. The definition of a pivot is
similar to that of an ancillary statistic, but it is also allowed to use the parameter θ.

Definition 17.2 (Pivotal quantity). A pivotal quantity, abbreviated pivot, is a function Q(Y, θ)
whose distribution does not depend on θ.

If we can find a closed set A such that Prθ (Q(Y, θ) ∈ A) ≥ 1 − α, then it follows that the set
of θ where Q(Y, θ) ∈ A is a (1 − α)-confidence set, and if Q is monotone and A is an interval, then
this is also a confidence interval.

Example 17.3. If we have samples Y1 , . . . , Yn ∼ N (µ, σ 2 ), and σ 2 is known, then the z-score of Y
is a pivot: √
(Y − µ) n
Q1 (Y, µ) = ∼ N (0, 1).
σ
Otherwise, if σ 2 is unknown, then the following is also a pivot:

(Y − µ) n
Q2 (Y, µ) = ∼ tn−1 ,
s
1 Pn
where s2 = n−1 2
j=1 (Yj − Y ) is the sample variance. Both Q1 and Q2 recover the z-interval and
t-interval from Example 16.8, respectively.

Note that we could also have used Q1 as a pivot in the case when σ 2 is unknown, but this would
produce a 2-dimensional confidence set for (µ, σ 2 ). This is not useful if we just want an interval for
µ. The advantage of the t-statistic is that it does not depend on σ 2 .

Example 17.4. In Example 17.1, we inverted the acceptance region of the likelihood ratio test to
obtain a confidence interval for the parameter θ. However, we observed that the likelihood ratio Λ
was actually invariant in distribution based on the parameter θ, so it is a pivotal quantity.
We also could have chosen to invert the pivot Y /θ0 directly, which would yield a slightly different
confidence interval for θ, based on quantiles of the gamma distribution.

17.3 Criteria for Selecting CIs


The power of the confidence interval we get from a pivot largely depends on the quality of the
pivot. Even after picking a pivot, we still need to choose some contiguous region A ⊂ R where the
quantity Q(Y, θ) contains at least 1 − α of its mass. There are a few desiderata:

• Shortest-width: Minimize the value of U (Y) − L(Y).

46
• Equal-tails: Satisfy the requirement that Prθ (θ > U (Y)) = Prθ (θ < L(Y)).
L(Y)+U (Y)
• Centered on an estimator: For some estimator θ̂, choose so that θ̂ = 2 .

It’s fairly straightforward to satisfy the equal-tails and centered on an estimator requirements by
simply looking at appropriate regions of mass 1 − α in the probability distribution of Q(Y, θ). For
the shortest-width requirement, this comes down to optimization, minimizing b − a subject to
Z b
f (x) dx ≥ 1 − α,
a

where f (x) is the probability density function of the pivot. In the case when the distribution is
unimodal, the solution to this optimization problem is just a level set A = {x : f (x) ≥ cα } of the
probability density, with cα selected so that the mass in A is at least 1 − α.

47
18 November 8th, 2021
Today we continue discussing confidence intervals and introduce Bayesian inference.

18.1 Asymptotic Confidence Intervals


Similar to how we define asymptotic level-α hypothesis tests, we can also derive methods for
constructing asymptotic (1 − α)-confidence intervals. Formally, an asymptotic CI is an interval
estimator [Ln (Y), Un (Y)], where n denotes the sample size, such that
lim inf Prθ (Ln (Y) ≤ θ ≤ Un (Y)) ≥ 1 − α.
n→∞

Analogous to last lecture, we can construct asymptotic confidence intervals by inverting asymptotic
hypothesis tests (see Section 12.1), or by finding an asymptotic pivot, either by using variance-
stabilizing transformations or other techniques.
Example 18.1. Suppose that Y1 , . . . , Yn ∼ Pois(λ), and we would like to find a confidence interval
for λ. The mean and variance are both λ, so the asymptotic distribution of the MLE Y is
√ d
n(Y − λ) −
→ N (0, λ).
√ √ d
This means that n(Y − λ)/ λ − → N (0, 1), which is an asymptotic pivot. Then, an asymptotic
(1 − α)-confidence interval is
 √ 
n(Y − λ)
C(Y) = λ : √ ≤ z1−α/2
λ
n o
2 2 2
= λ : nλ − (2nY + z1−α/2 )λ + nY ≤ 0 .

By the quadratic formula, the endpoints of this confidence interval are at


q
2
2nY + z1−α/2 2
± 4nY z1−α/2 4
+ z1−α/2
.
2n
Example 18.2. Another confidence interval of the same level can be obtained with Slutsky’s
d √ √ d
theorem and the consistency of the MLE, Y − → λ, which implies that n(Y − λ)/ Y −
→ N (0, 1).
This yields the slightly different confidence interval,
 q q 
C(Y) = λ : Y − Y /n z1−α/2 ≤ λ ≤ Y + Y /n z1−α/2 .

Example 18.3. What if we wanted to instead use √ a variance-stabilizing transformation? For the
Poisson distribution, this transformation is λ 7→ λ, from which we can use the Delta method to
get  √ !2 
√ p √  d
 
∂ λ 1
n Y − λ − → N 0, λ = N 0, .
∂λ 4

This asymptotic confidence interval has endpoints at


z1−α/2 2
p 
Y ± √ .
2 n
Each of these three examples yields a slightly different confidence interval in the finite-sample
case, but when taking the limit to asymptotic behavior as n → ∞, their corresponding confidence
regions are equivalent.

48
18.2 Bayesian Inference
In Bayesian inference, rather than carefully considering individual parameter configurations in the
set Θ, we instead imagine a distribution of parameters over the global set. In this setting, our usual
frequentist model fθ (y) is changed to a conditional distribution f (y | θ). The prior probability of
i.i.d. samples Y1 , . . . , Yn ∼ f (y | θ) is equal to the likelihood function

L(θ) = f (Y | θ).

Then, using Bayes’ rule, the posterior probability distribution is given by


L(θ)π(θ)
π(θ | Y) = R ∝ L(θ)π(θ).
Θ L(θ̃)π(θ̃) dθ̃

Here, the integral in the denominator is called a normalizing constant, and it is usually intractable
to compute this exactly due to the very difficult integral over Θ.
Example 18.4. How might we compute or estimate the value of E [θ | Y] in the Bayesian inference
setting? One way to do this is is to discretize the parameter space Θ and generate a grid of values,
then estimate the normalizing constant as a Riemann sum over this grid. This suffers from the
curse of dimensionality.
Another, often more efficient method, is to use the Markov chain Monte Carlo (MCMC) al-
gorithm or other sampling techniques to generate approximately i.i.d. samples from π(θ | Y).
Although these techniques are not exact, they usually obtain good results after a fully-polynomial
number of mixing steps. Then, an estimate for E [θ | Y] is obtained by taking the sample mean.
Although the general Bayesian inference problem is intractable (application of Bayes’ rule),
there are some distributions for which we can analytically compute the posterior distribution.
Definition 18.5 (Conjugate prior). Suppose we have a family of distributions, G = {gτ (θ) : τ ∈ T}.
Then, G is conjugate to f (y | θ) if π(θ) = gτ (θ) for some τ ∈ T means that for every y of the data
drawn from f (• | θ),
π(θ | y) = gτ 0 (θ)
for some τ 0 ∈ T.
In other words, a conjugate prior is a family of distributions for θ, such that updating your
prior on sampled data still results in a distribution of that family.
Example 18.6. If Y ∼ Bin(n, θ), where n is known, then
 
n y
f (y | θ) = θ (1 − θ)n−y .
y
Then, a conjugate prior for this distribution is π(θ) = Beta(α, β) for α, β > 0, with PDF

θα−1 (1 − θ)β−1
,
B(α, β)
where B(α, β) is the Beta function. Under this assumption on the prior, the posterior distribution
for θ has a particularly simple form:

π(θ | y) ∝ π(θ)f (y | θ) ∝ θy+α−1 (1 − θ)(n−y)+β−1 .

Therefore, the Bayesian update is simply (α, β) 7→ (y + α, n − y + β).

49
Conjugate priors are useful because they are usually already known for common distributions,
and they make it computationally efficient to apply inference techniques. However, one common
criticism of conjugate priors is that your actual prior knowledge of the parameters of a distribution
may not match any of the existing distributions gτ (θ). This limitation is somewhat reduced in
severity by taking mixtures of conjugate priors.

Proposition 18.7. If we have k prior distributions π1 (θ), . . . , πk (θ), such that each corresponding
posterior distribution
L(θ)πj (θ)
πj (θ | y) =
Cj
is computable, then if π(θ) = kj=1 wj πj (θ) for weights w1 , . . . , wk summing to 1, we have
P

Pk k
L(θ) j=1 wj πj (θ)
X wj C j
π(θ | y) = = πj (θ | y).
C C
j=1

Therefore, π(θ | y) is a mixture of πj (θ | y) for each j, with coefficients

wj C j
wj∗ = Pk .
`=1 w ` C `

Example 18.8. Given a continuous probability density function f for a scalar random variable,
supported on [−c, c], how might you approximate f as a mixture of Gaussians? One simple method
would be to divide up the interval into N bins, then add a Gaussian centered at each bin with
variance c/N and weight equal to the amount of probability mass in that bin.

18.3 Conjugate Priors of NEFs


One interesting property of natural exponential families is that they have a simple conjugate prior
in the Bayesian inference setting.

Proposition 18.9. Suppose that we have an NEF with distribution

f (y | η) = exp{ηy − ψ(η)}h(y),

where we draw n i.i.d. samples Y1 , . . . , Yn ∼ f (y, η). Our natural sufficient statistic is Y = T . The
mean and variance of the distribution are µ = E [Y | η] = ψ 0 (η) and Var [Y | η] = ψ 00 (η) = V (µ),
respectively, and our likelihood function is
( n )
X
L(η) ∝ exp n Yi − nψ(η) = exp{n(ηT − ψ(η))}.
i=1

Then, the following family of distributions, parameterized by (µ0 , r), is a conjugate prior:

π(η) ∝ exp{r(ηµ0 − ψ(η))}.

Alternatively, we can reparameterize this in terms of the mean µ = ψ 0 (η), which yields
 Z 
1 µ − µ0
π(µ) ∝ exp −r dµ .
V (µ) V (µ)

50
Proof. Applying Bayes rule by multiplying the likelihood and prior distribution density function,
we have that the posterior distribution is proportional to

π(η | Y) ∝ exp{(nT + rµ0 )η − (n + r)ψ(η)}.

From this, we can see that the prior π(η) and posterior π(η | Y) are in the same algebraic form.
The update can be effectively summarized by:

r 7−→ n + r,
rµ0 + nT
µ0 7−→ .
n+r
Furthermore, the mean of this distribution (or any distribution in the conjugate prior) is actually
E [µ] = µ0 . It is not obvious to show this, but see [DY79] for the proof. Therefore, the parameter
can be interpreted as the prior expectation for the mean.
Furthermore, we can interpret the posterior update rule as a weighted average of the prior mean
µ0 and the observed mean T = Y . The weights are determined by the relative Fisher information
of two estimates µ0 and T , treated as if they used r and n independent samples, respectively. Then,
we have E [µ | Y] = Bµ0 + (1 − B)T , where
r
V (µ)
B= r n .
V (µ) + V (µ)

Intuitively, B is the shrinkage factor, which determines how much the prior mean µ0 contributes
to the posterior distribution, compared to the new evidence Y.

51
19 November 10th, 2021
Today, we finish talking about conjugacy, then introduce the Jeffreys prior, as well as Bayesian
point and interval estimation.

19.1 More on Conjugacy


Let’s give a couple examples of the conjugate priors corresponding to natural exponential families,
as determined by application of Proposition 18.9.
Example 19.1. If Y ∼ N (µ, c), where the variance function is V (µ) = c for known c, then
Z ∞
µ − µ0 1 µ0
dµ = µ2 − µ.
−∞ V (µ) 2c c
Since this is the form of the negative log-likelihood function, up to a constant, we conclude that
the posterior distribution of µ is normal.
Example 19.2. If Y ∼ Pois(µ), then V (µ) = µ, and the corresponding integral is
Z ∞
µ − µ0
dµ = µ − µ0 log µ.
0 V (µ)
Therefore, the posterior distribution π(µ) is a gamma distribution, since the likelihood is propor-
tional to µrµ0 e−rµ ∼ 1r · Gamma(rµ0 + 1).

19.2 The Jeffreys Prior


The Jeffreys prior is a kind of “standard” prior, which does not depend on a specific parameteri-
zation. In other words, the prior is equivariant on bijective transformations of the parameter θ. It
turns out that this is enough to uniquely specify a prior.
Definition 19.3 (Jeffreys prior). Given a model parameterized by θ, with Fisher information I(θ),
the Jeffreys prior for θ is given implicitly by
p
π(θ) ∝ I(θ).

It turns out that the Jeffreys prior is equivariant, which we formalize and prove below.
Proposition 19.4. Given a model with parameter θ, the Jeffreys prior for θ is the same as the
prior that comes from reparameterizing the Jeffreys prior for β(θ) and changing variables to θ, for
any smooth bijective function β.
Proof. We can simply compute via the chain rule that
s v " v "
 2 u  2 #  2 u  2 #
dβ dβ u d` dβ u d` dβ p
π(θ) = π(β) ∝ I(β) = tEβ = tEθ · = I(θ).
dθ dθ dβ dθ dβ dθ

Example 19.5 (Jeffreys prior for a variance-stabilizing transformation). If we apply a variance-


stabilizing parameter h(θ) to our parameter such that h0 (θ) ∝ I(θ), then the Jeffreys prior on
p

h(θ) has constant likelihood at every point; it is a uniform distribution.

52
Note that in the above example, if the parameter domain is infinite, then this prior is improper
or uninformative, since it cannot be normalized to form a true distribution with nonzero mass.
However, we can still use improper priors in calculations, as long as we are careful to normalize the
posterior distribution to be proper.

Example 19.6 (Jeffreys prior for an NEF). If Y1 , . . . , Yn ∼ NEF[µ, V (µ)], then

Eµ [Yi ] = µ and Varµ [Yi ] = V (µ).

Then, observe that Y ∼ NEF[µ, V n(µ) ], and the Fisher information of Y with respect to the natural
parameter η is
 2 
d
I(η) = − Eη log f (y | η)
dη 2
= nψ 00 (η)
= nV (µ).
p
Therefore, we conclude that the Jeffreys prior is π(η) ∝ V (µ). If we wanted the prior in terms of
the mean µ instead of the natural parameter η, then we could transform it by a change of variables,
which yields the expression
−1
dη p dψ 0 (η) 1
π(µ) ∝ π(η) ∝ V (µ) =p .
dµ dη V (µ)

19.3 Bayesian Point and Interval Estimation


Let’s talk a bit more about how we might do estimation in Bayesian inference. Recall that there
are two main types of inference, which are point and interval estimation. The former, Bayesian
point estimation, is fairly simple.8

Example 19.7 (Minimum-variance Bayesian estimation). If our goal is to find an estimator θ̂ that
minimizes the mean-squared error
h i
MSE = E (θ̂ − θ)2 | Y
= E (θ − E [θ | Y ])2 + (θ̂ − E [θ | Y ])2 ,
 

then the optimal solution, is simply θ̂ = E [θ | Y ].

There are a couple other common Bayesian point estimation techniques. For example, if we
wanted to optimize for the minimum absolute error, then we would take the median of the posterior
distribution θ | Y , rather than its expectation. Likewise, if we wanted to find the highest posterior
likelihood, then we would use the mode of the posterior distribution, which is the maximum a
posteriori (MAP) estimator.
This is essentially all we will discuss about Bayesian point estimation. Hopefully you’ve seen
MAP estimation in previous statistics classes. Now, we’ll move on to interval estimation. The
Bayesian version of a confidence interval is called a credible interval.
8
Notice how Bayesian methods are much simpler than point estimation in the frequentist setting, which is why
we spend more time in this class on frequentist inference.

53
Definition 19.8 (Credible interval). An estimator C(Y ) is called a (1 − α)-credible interval if

Pr(θ ∈ C(Y )) ≥ 1 − α.

This is just a probability interval of π(θ | Y ) with mass at least 1 − α.


Similar to confidence intervals created by pivoting, the smallest credible interval for a unimodal
posterior distribution is just the level set. Also, we may prefer an equal-tailed C(Y ), or for C(Y )
to be centered on some estimator θ̂, which is analogous to the analysis in Section 17.3.
Example 19.9 (Jeffreys prior for a normal distribution). Suppose that Y | µ ∼ N (µ, σ 2 ), where
σ 2 = 100 and we observe some value of Y . Then, Y is the maximum likelihood estimator for µ, as
well as the mean, median, and mode of the posterior distribution for µ under Jeffreys prior.
If we instead used the conjugate prior N (100, 152 ) in the above example, then r = σ 2 /τ 2 =
100/225 = 0.44. The posterior distribution would be N (113.8, 69.2). Notice how the posterior
variance is equal to half the harmonic mean of σ 2 = 100 and τ 2 = 225, and the posterior mean can
be computed by taking a weighted average of r·100+1·120
r+1 .
Example 19.10 (Frequentist coverage of credible intervals). Let’s consider the frequentist coverage
of the above intervals under various priors. Essentially, the coverage probability is a function of
the parameter µ, and it specifies how likely it is that µ lies in the Bayesian credible interval.
When taking the Jeffreys prior, the Bayesian 95% credible interval is simply Y ±1.96σ, which co-
incides with the frequentist 95% confidence interval. Therefore, the coverage probability is precisely
0.95, which makes sense given that the Jeffreys prior for µ in N (µ, σ 2 ) is essentially uninformative.
However, with a conjugate prior of π ∼ N (µ0 , τ 2 ), where r = σ 2 /τ 2 , we can calculate a general
formula for the frequentist coverage of Bayesian credible intervals in our model, which is
   
r n σ
Prµ µ ∈ µ0 + Y ± z1−α/2 √
n+r n+r n+r
    
r n σ σ
= Prµ µ ∈ µ0 + µ+ √ Z ± z1−α/2 √
n+r n+r n n+r
 √  
σ n r σ
= Prµ − Z∈ (µ0 − µ) ± z1−α/2 √
n+r n+r n+r
" r #!
r n+r
= Prµ Z ∈ √ (µ − µ0 ) ± z1−α/2
σ n n
r ! r !
σ σ2 σ σ2
= Φ 2 √ (µ − µ0 ) + z1−α/2 1 + 2 − Φ 2 √ (µ − µ0 ) − z1−α/2 1 + 2 .
τ n τ n τ n τ n

If we graph this expression, it roughly follows one’s intuition for how the coverage probability should
behave based on the conjugate prior. When n is small relative to r, if our prior is roughly correct
(µ ≈ µ0 ), the frequentist coverage is higher than 1 − α. However, if |µ − µ0 | is greater than about
τ , the coverage rapidly decreases because the inaccurate prior leads our interval astray.
As n increases relative to r, the frequentist coverage of the Bayesian credible interval at any µ
converges to 1 − α, since the prior distribution gets less important relative to the evidence.
In general, as n gets large, the frequentist coverage of Bayesian credible intervals gets close to
the true coverage. This is a consequence of the Bernstein-von Mises theorem, which states that
the Bayesian credible interval for a fixed prior is asymptotically normal and converges to the Wald
confidence interval (Definition 12.2), as n → ∞.

54
20 November 15th, 2021
Today, we discuss Bayesian predictive inference and decision analysis.

20.1 Bayesian Predictive Inference


Suppose that we have observed data Y = (Y1 , . . . , Yn ), as well as some new data point Ynew . The
problem of predictive inference asks us to find the posterior predictive distribution,
Z Z
f (Ynew | Y) = f (Ynew , θ | Y) dθ = f (Ynew | θ)π(θ | Y) dθ.
Θ Θ

Assume that Ynew | θ is known, since this model is parametric. Then, the predictive inference
problem is essentially equivalent to the Bayesian update framework that we have already discussed
in class, using the data Y to find the posterior distribution π(θ) using Bayes’ rule. Predictive
intervals can also be computed using the same method.

20.2 Frequentist Decision Analysis


Decision analysis is concerned with the following choice scenario. The statistician chooses a function
δ from the sample space to the decision space D, and nature chooses a parameter θ from a data-
generating process (DGP). Then, depending on what the statistician cares about, our objective is
to minimize a loss function
L(θ, δ(Y )) ∈ R≥0 .
(This problem can be extended to multi-dimensional data, but we will consider the one-dimensional
case Y ∈ R in this section because it is simpler.) Some examples of loss functions that are commonly
used in statistics include:

• (Mean squared error). L(θ, δ(Y )) = (θ − δ(Y ))2 .

• (Mean absolute error). L(θ, δ(Y )) = |θ − δ(Y )|.

• (Asymmetric linear error). For some p, q > 0,


(
p(θ − δ(Y )) if θ ≥ δ(Y ),
L(θ, δ(Y )) =
q(δ(Y ) − θ) if θ < δ(Y ).

• (Linear exponential loss). For some small value of c,


1
L(θ, δ(Y )) = ec(θ−δ(Y )) − c(θ − δ(Y )) − 1 ≈ c2 (θ − δ(Y ))2 .
2

Now, we discuss how to analyze and solve decision analysis problems in a frequentist context. Core
to this idea is the following notion.

Definition 20.1 (Risk). The risk function for a decision analysis problem is

R(θ, δ) = Eθ [L(θ, δ(Y ))] .

55
Definition 20.2 (Domination). We say that a decision rule δ1 dominates another rule δ2 if

R(θ, δ1 ) ≤ R(θ, δ2 ) ∀θ ∈ Θ,

with strict inequality holding for at least one θ ∈ Θ. On the other hand, we say that δ1 is as good
as δ2 if neither rule dominates the other.

Definition 20.3 (Admissibility). A decision rule δ is called admissible if there is no other rule
δ 0 ∈ D that dominates it.

Roughly speaking, we can think of admissible rules as being better than every other decision
rule for at least one value of the parameter θ. These definitions have a very frequentist flavor, as we
assume no prior information about the distribution of the parameters. However, while inadmissible
rules are usually poor, not all admissible rules are necessarily good either.

Example 20.4. Consider Y1 , . . . , Yn ∼ N (θ, 1), and let our loss function be the mean squared
error. A simple and reasonable decision rule in this case might be δ(Y) = Y , which happens to be
the UMVUE. However, setting δ(Y) = 5 to be a constant is technically a valid decision rule, since
for θ = 5, this is the unique rule with the global minimum risk of R(θ, δ) = 0.

Admissibility is a simple notion that simply checks, intuitively, that δ is good somewhere. Next,
we will introduce the notion of minimaxity, which intuitively checks that δ is terrible nowhere.

Definition 20.5 (Minimax decision rule). We say that δm is minimax with respect to D if

δm ∈ argmin max R(θ, δ).


δ∈D θ∈Θ

In other words, δm minimizes the global maximum risk for an unknown value of θ.

Notice that minimax decision rules are simple to find. It simply requires computing, for each
δ ∈ D, the maximum of R(θ, δ) for all θ ∈ Θ. Then, you take δ that results in the smallest of these
maximum-risk values, and that is the minimax rule.

Proposition 20.6 (Pencil problem). If a decision rule δ ∗ with constant R(θ, δ ∗ ) = c(δ ∗ ) is admis-
sible in D, then it is minimax in D.

Proof. Assume for the sake of contradiction that there exists δ 0 ∈ D such that δ 0 has smaller
maximum risk than δ ∗ . Then,

R(θ, δ 0 ) ≤ max R(θ, δ 0 ) < max R(θ, δ ∗ ) = c(δ ∗ ) = R(θ, δ ∗ ),


θ∈Θ θ∈Θ

where the last step follows because the risk function is a constant for θ ∈ Θ. Therefore, δ 0 dominates
δ ∗ , which is a contradiction, so we conclude.

20.3 Bayesian Decision Analysis


Next, we shift gears to discuss the Bayesian formulation of decision analysis, which is fairly different
because we have a prior on the parameter distribution θ ∼ π.

Definition 20.7 (Bayes risk). The Bayes risk function is

B(π, δ) = Eπ [R(θ, δ)] = Eπ [Eθ [L(θ, δ(Y ))]] .

56
Unlike frequentist inference, where we introduced two different criteria for desirable decision
rules, the Bayesian formulation has an unambiguous criterion for a “Bayes rule” — not to be
confused with Bayes’ rule or Bayes’ theorem.

Definition 20.8 (Bayes decision rule). A decision rule δ ∈ D is called a Bayes rule with respect
to π and D if
δ π ∈ argmin B(π, δ).
δ∈D

Proposition 20.9. Suppose that Y ∼ fθ , for θ ∈ Θ, and we have a prior π(θ), decision space D,
and loss function L. If D contains a δ ∗ minimizing

Eπ [L(θ, δ(Y )) | Y ]

almost surely with respect to the marginal distribution of Y ,


Z
m(y) = Eπ [fθ (y)] = fθ (y)π(θ) dθ,
Θ

then δ ∗ is is a Bayes rule. In other words, Bayes rules are precisely the decision rules that minimize
the expected loss with respect to the posterior distribution of Y .

Proof. If all δ ∈ D have infinite Bayes risk, then all of them are Bayes rules. Otherwise, suppose
that there exists δ ∈ D with B(π, δ) < ∞. Then,

B(π, δ) = Eπ [Eδ [L(θ, δ(Y )]]


= Em [Eπ [L(θ, δ(Y ) | Y ]]
≥ Em [Eπ [L(θ, δ ∗ (Y ) | Y ]]
= B(π, δ ∗ ).

Therefore, δ ∗ has Bayes risk that is less than or equal to the Bayes risk of δ.

57
21 November 17th, 2021
Today, we illustrate a link between frequentist and Bayesian decision analysis by showing that under
some reasonable regularity conditions on the model, any Bayes rule is admissible, and conversely,
every admissible rule is Bayes.

21.1 Bayes Rules are Admissible


First, we show the basic result, that if a rule is unique and Bayes, then it is also admissible, meaning
that it is not dominated by any other decision rule.
Theorem 21.1. Any unique Bayes rule (up to equality, almost surely in m) is admissible.
Proof. Suppose for the sake of argument that there is some decision rule δ 0 such that

R(θ, δ 0 ) ≤ R(θ, δ π ),

for all θ ∈ Θ. Then, integrating this over θ ∼ π implies that

B(π, δ 0 ) = Eπ R(θ, δ 0 )
 

≤ Eπ [R(θ, δ π )]
= B(π, δ π ).

However, δ π is defined as the unique global minimizer of the Bayes risk under parameter distribution
π, so this implies that B(π, δ 0 ) = B(π, δ π ), and δ 0 = δ π almost surely.

The above argument was fairly simple, but we can also prove variants of this fact. The following
variant works for finite parameter spaces, even if for Bayes rules that are not necessarily unique.
It uses the fact that a dominating rule must be different in at least one risk value.
Theorem 21.2. Let Θ = {θ1 , . . . , θk } be a finite set of parameters, and let δ π be Bayes for π with
finite Bayes risk, where π(θj ) > 0 for all j. Then, δ π is admissible.
Proof. Suppose for the sake of argument that there exists a rule δ 0 such that

R(θj , δ 0 ) ≤ R(θj , δ π ),

for all j ∈ {1, . . . , k}, with strict equality for some j0 . Then,
k
X
0
B(π, δ ) = π(θj )R(θ, δ 0 )
j=1
X
≤ π(θj0 )R(θj0 , δ 0 ) + π(θj )R(θ, δ π )
j6=j0
X
π
< π(θj0 )R(θj0 , δ ) + π(θj )R(θ, δ π )
j6=j0
π
= B(π, δ ).

This is a contradiction, so we conclude that it is impossible.

Can we generalize this same result to continuous parameter spaces? Initially, this might seem
difficult, since the strict inequality is harder to use in an infinite parameter space. However, with
some analysis, it turns out that the answer is yes: we can, assuming extra regularity conditions.

58
Theorem 21.3. Let Θ ⊆ R be an open set, such that π(θ) has support Θ, and R(θ, δ) is continuous
in θ, for all δ ∈ D. Then, if δ π is a Bayes rule for π, it is admissible.
Proof. Suppose for the sake of argument that δ 0 is a rule such that

R(θ, δ 0 ) ≤ R(θ, δ π ), ∀θ ∈ Θ,
0 π
R(θ0 , δ ) < R(θ0 , δ ), for some θ0 ∈ Θ.

Then, let η = R(θ0 , δ π ) − R(θ0 , δ 0 ) > 0. By continuity, there exists some  > 0 such that for all θ
in an open ball A of radius  around θ0 ,
η
R(θ, δ π ) − R(θ, δ 0 ) > .
2
This region A has positive probability mass because the support of π is over all of Θ, so
η
B(π, δ π ) − B(π, δ 0 ) > π(A) > 0.
2
Thus, we have a contradiction, as we assumed that δ π is a Bayes rule.

21.2 Admissible Rules are Bayes


Next, we discuss the converse result, which is that admissible rules are Bayes. Before we can prove
anything, we’ll look a little bit at the high-dimensional geometry of admissible decision rules.

Figure 3: Geometric drawing of admissible decision rules.

As an initial example, assume that Θ = {θ1 , θ2 }, and suppose that we map each decision rule δ to
a geometric point (R(θ1 , δ), R(θ2 , δ)) ∈ R2 . Then, each point is located in the upper-right quadrant

59
of the Cartesian plane, and the convex hull of the points represents every risk configuration that can
arise as the randomized mixture of one or more decision rules. See Fig. 3 for a diagram depicting
this geometry, where each point is a rule δ ∈ D, and the lower-left boundary of the hull represents
precisely the admissible rules that are not dominated by any other rule.

Proposition 21.4. Assume that |Θ| = k, and let S ⊂ Rk be the convex hull of the points sδ =
(R(θ1 , δ), . . . , R(θk , δ)) for all δ ∈ D. Let the lower orthant of a point in Rk be defined as

Q(s) = {(x1 , xk ) ∈ Rk | x1 ≤ sq , . . . , xk ≤ sk }.

Then, δ is admissible if and only if Q(sδ ) ∩ S = {sδ }.

This proposition is fairly straightforward, and it makes sense, given that a decision rule is only
valid if it occupies an optimal point on the frontier of possibilities. Now we might ask what the
Bayes rules in the picture are. It turns out that they are precisely the same; the Bayes rules are
solutions to linear programming problems optimizing the value of some linear form (represented by
the dotted lines) on the convex set S, so they can only occur on a lower-left boundary point of the
hull. We formalize this intuition with the following theorem.

Theorem 21.5 (Complete class theorem for finite Θ). Let |Θ| = k be a set of parameters, and let
D be a set of decision rules that is closed under randomized combinations. Furthermore, assume
that R(θ, δ) is a nonnegative risk function taking finite values for all δ ∈ D and θ ∈ Θ. Then, if
δ ∗ ∈ D is admissible, it is a Bayes rule with respect to some proper prior.

We will prove this theorem in the next lecture. Note that the closure under randomized combi-
nations is an important condition, as otherwise, the set S would not be guaranteed to be convex, so
not all admissible rules on the would be accessible as the minimum of a linear form corresponding
to the supporting hyperplane.

60
22 November 22nd, 2021
Today, we first prove the complete class theorem for finite Θ as stated in last lecture, then we
discuss the admissibility of the sample mean in univariate normals.

22.1 Proof of the Complete Class Theorem


First, we pick up from where we left off last time and prove the complete class theorem.

Proof of Theorem 21.5. Once again, we use the notation where S ⊂ Rk is a convex set containing
the risk values sδ of hypotheses δ ∈ D, where

sδ = (R(θ1 , δ), . . . , R(θk , δ)).

Then, applying Proposition 21.4 to the admissible hypothesis δ ∗ , we get that Q(s∗ ) ∩ S = {s∗ },
where we abbreviate sδ∗ = s∗ for clarity. Then, if we let Q̃(s∗ ) = Q(s∗ ) \ {s∗ }, we have

Q̃(s∗ ) ∩ S = ∅.

Notice that s∗ is an extreme point of Q(s∗ ), so removing it maintains that the resultant set is still
convex. By the hyperplane separation theorem,9 there exists some nonzero normal vector w ∈ Rk
representing a linear form, such that

sup w> x ≤ inf w> s.


x∈Q(s∗ ) s∈S

By inspecting the structure of Q, clearly all coordinates of w must be nonnegative. Finally, we can
turn w into a proper prior on Θ by taking
wj
∀j : π(θj ) = .
w1 + · · · + wk

Observe that B(π, δ) ∝ w> sδ for all δ, so δ ∗ must be a Bayes rule for the prior π, as desired.

22.2 Admissibility of the Sample Mean


Before showing that the sample mean is admissible for any loss function, we will first do an technical
thought experiment to get some intuition about its properties.

Example 22.1. Suppose that Y1 , . . . , Yn ∼ N (µ, σ 2 ), and σ 2 is known. Also, suppose that our
prior is µ ∼ N (µ0 , τ 2 ). Then, our shrinkage factor is

σ 2 /n
B= .
σ 2 /n + τ2

As n → ∞, the shrinkage factor approaches B = 0, so the prior makes a smaller and smaller effect
on the final value of the posterior mean. Therefore, no matter what prior we choose, the limit of
the Bayes risk-minimizing decision rule is just the sample mean.

Motivated by this example, here is the big result that we’re excited about.
9
This is a key property of convex sets in Euclidean space and will be our sledgehammer in this proof.

61
Theorem 22.2 (Blyth’s method). Given i.i.d. Y1 , . . . , Yn ∼ N (µ, σ 2 ), with variance σ 2 known, the
sample mean Y is an admissible decision rule for µ with respect to the squared error.
Proof. Suppose for the sake of argument that Y is inadmissible, so by definition, there exists some
decision rule t(Y) that dominates Y . Therefore, for any µ,

R(µ, t) = Eµ (t(Y) − µ)2 ≤ Eµ (Y − µ)2 = R(µ, Y ),


   

with strict inequality holding for some µ0 . Without loss of generality, assume that µ0 = 0. Then,
since our risk function is continuous in µ, we can argue that there exists some  > 0 such that for
any |µ| < , we have
R(µ, Y ) − R(µ, t) ≥ .
Now, assume for the sake of simplicity that σ 2 /n = 1, as the particular value will not matter. Then,
R(µ, Y ) = 1, so we can rewrite this as

R(µ, t) = Eµ (t(Y) − µ)2 ≤ 1 − 1{µ∈(−,)} .


 

Furthermore, given a prior πτ 2 ∼ N (0, τ 2 ), we see that

τ2
Y
1 + τ2
is a Bayes rule with respect to πτ 2 for all τ 2 > 0. Hence, using the fact that Bayes rules minimize
the expectation of the posterior risk given their prior, we have

Eπτ 2 1 − 1{µ∈(−,)} ≥ Eπτ 2 Eµ (t(Y) − µ)2


    
" " 2 ##
τ2
≥ Eπτ 2 Eµ Y −µ
1 + τ2
= Eπτ 2 B 2 µ2 + (1 − B)2
 

= B 2 τ 2 + (1 − B)2
1
=1− .
1 + τ2
However, the original risk on the left-hand side of this inequality can also be written as
 
Eπτ 2 1 − 1{µ∈(−,)} = 1 − Prπτ 2 (− < µ < )
    
=1− Φ −Φ −
τ τ
!
2 Φ τ − Φ − τ


2
=1− .
τ 2/τ

Combining this with the previous inequality yields


!
22 Φ τ − Φ − τ

1
1− ≥1− .
τ 2/τ 1 + τ2

Finally, we can rearrange to get the inequality



 
!
2 Φ τ − Φ −τ τ
2 ≤ .
2/τ 1 + τ2

62
Now, we examine the limit behavior of this inequality, as our prior tends towards the Jeffreys prior
when τ 2 → ∞. Using the difference quotient on the left-hand side, we see that it converges to

22
22 Φ0 (0) = √ .

On the other hand, the right-hand side converges to zero, so we have a contradiction.

Here is a nice bonus property of the sample mean, which we get for free.

Corollary 22.2.1. If Y1 , . . . , Yn ∼ N (µ, σ 2 ), with σ 2 known, then Y is minimax for µ with respect
to the squared error.

Proof. Recall that an admissible decision rule with constant risk must also be minimax. The sample
mean Y is an admissible decision rule with constant risk σ 2 /n, invariant of µ.

22.3 Least Favorable Priors


So far, we have been looking at the relationship between Bayes rules and admissible rules. In this
section, we turn our attention to minimaxity, which is connected to least favorable priors.

Definition 22.3 (Least favorable prior). If rπ = B(π, δ π ), then a prior π ∗ is called least favorable
if rπ∗ ≥ rπ for all proper priors π. In other words,

π ∗ ∈ argmax min B(π, δ).


π δ

Theorem 22.4. If δ π is Bayes with respect to a prior π such that rπ = supθ R(θ, δ π ), then:

(i) δ π is a minimax decision rule,

(ii) δ π is the unique minimax decision rule if δ π is the unique Bayes rule for π, and

(iii) π is a least favorable prior.

Proof. To prove (i), note that any other decision rule δ has a maximum risk value over all parameters
θ that cannot be smaller than that of δ π , since

sup R(θ, δ) ≥ Eπ [R(θ, δ)] ≥ Eπ [R(θ, δ π )] = rπ = sup R(θ, δ π ).


θ θ

Furthermore, if δ π is the unique Bayes rule for π, then the second step above becomes a strict
inequality Eπ [R(θ, δ)] > Eπ [R(θ, δ π )], so (ii) follows. Finally, to prove (iii), note that for any other
prior π̃ on the parameters, we have

rπ̃ = B(π̃, δ π̃ ) ≤ B(π̃, δ π ) ≤ sup R(θ, δ π ) = rπ .


θ

What are the scenarios where we can use the above theorem? Well, one such scenario when rπ
achieves this supremum is when our risk R(θ, δ π ) is the same for all values of θ.

Corollary 22.4.1. If a Bayes rule δ π has constant risk, then it is minimax.

63
23 November 29th, 2021
Today, we discuss the admissibility and minimaxity of the sample mean (also the maximum likeli-
hood estimator) in any dimension and show the counterintuitive result that admissibility does not
generalize to higher dimensions (Stein’s paradox).

23.1 More on Minimax Decision Rules


Recall that a minimax decision rule δ ∗ minimizes the max risk over all parameters: supθ R(θ, δ ∗ ).
Furthermore, a Bayes rule for a least favorable prior is also minimax, as we showed at the end of
the last lecture.
Unfortunately, if we want to apply this theorem to prove that some decision rule δ is minimax,
we actually need to show that our estimator is a Bayes rule with respect to some proper prior, so
this doesn’t work for things like the sample mean, which is Bayes with respect to the Jeffreys prior.
To get around this, we will extend our definition to sequences of priors.

Definition 23.1 (Least favorable sequence). A sequence of priors πk is called least favorable if

lim rπk ≥ rπ .
k→∞

for all proper priors π.

Theorem 23.2 (Minimax duality). Let πk be a sequence of priors and δ a decision rule. If

sup R(θ, δ) = lim rπk ,


θ k→∞

then δ is minimax, and {πk } is a least favorable sequence of priors.

Proof. Note that the inequality version of the stated condition always holds, since if δ 0 is any
decision rule, then for any k,
sup R(θ, δ 0 ) ≥ B(πk , δ 0 ) ≥ rπk .
θ

Therefore, we have shown that δ is minimax, since it has the lowest possible value of supθ R(θ, δ).
For the other part, note that for any prior π,

rπ = Eπ [R(θ, δ π )] ≤ Eπ [R(θ, δ)] ≤ sup R(θ, δ) = lim rπk .


θ k→∞

Corollary 23.2.1. If Y ∼ N (µ, σ 2 Ik ) with σ 2 known, then µ̂ = Y is a minimax estimator for the
mean with respect to the sum of squared errors.

Proof. Consider a sequence of priors πτ 2 = N (0, τ 2 Ik ). As τ 2 → ∞, this converges to the Jeffreys


prior on the entire space, and the Bayes risk approaches kσ 2 , which is the same as the constant risk
of the UMVUE decision rule µ̂ = Y. Therefore, the sequence πτ 2 as τ 2 → ∞ is a least favorable
sequence of priors, and µ̂ = Y is a minimax estimator.

64
23.2 Stein’s Paradox
Although we just showed that the sample mean is minimax with respect to squared error, we did not
prove admissibility, which is a different criterion. Note that the sample mean is trivially admissible
in one dimension. We might also expect for it to be admissible in k dimensions, but surprisingly,
this is actually not the case for k ≥ 3!
Example 23.3 (Stein’s paradox). Consider a vector of independent random variables Y1 , . . . , Yk ∼
N (µi , 1), where our parameter vector is µ = (µ1 , . . . , µk ). The maximum likelihood estimator for
µ is the sample mean
µ̂MLE = Y.
Suppose that our loss function is the mean-squared error, so our risk is
" k #
X
R(µ, µ̂) = Eµ (µi − µ̂i )2 .
i=1

For the sample mean, we have R(µ, µ̂MLE )= k. The paradox is that the following estimator, known
as Stein’s estimator, dominates the sample mean for k ≥ 3, with risk R(µ, µ̂JS ) ≤ k:
 
JS k−2
µ̂ = 1 − Y.
kY k2
Let’s see why this paradox occurs. First, why does Theorem 22.2 (Blyth’s method) fail for
k ≥ 2, when it worked to prove admissibility of the sample mean when k = 1? The issue lies in the
last step, when we argued that πτ 2 (A) → 0 at a rate of 1/τ , which is slower than the rate-1/τ 2
convergence of the limits of Bayes risk. (Note that although Blyth’s method fails when k ≥ 2, there
are other methods that can prove admissibility for k = 2 specifically.)
To prove Stein’s theorem, which is that Stein’s estimator has risk ≤ k, we will prove two useful
intermediate lemmas, known simply as Stein’s identity and Stein’s lemma. The former is a formula
for the mean-squared error, and the latter is a step in deriving this formula.
Lemma 23.4 (Stein’s identity). Given Y ∼ N (µ, σ 2 Ik ), let
µ̂(Y) = Y + g(Y),
for any function g : Rk → Rk that is differentiable and satisfies the condition
" k #
X
Eµ |∇i gi (Y)| < ∞.
i=1
In other words, for any µ, the expectation of each diagonal entry of the Jacobian is finite. Then,
the expectation of the mean-squared error of µ̂ is
k
" #
X
Eµ kµ − µ̂k2 = kσ 2 + Eµ kg(Y)k2 + 2σ 2
 
∇i gi (Y) .
i=1
Proof. We can expand the left-hand side using linearity of expectation to get
Eµ kµ̂ − µk2 = Eµ kY + g(Y) − µk2
   
h i
= Eµ kY − µk2 + 2 Eµ (Y − µ)> g(Y) + Eµ kg(Y)k2
   
" #
Y−µ >

2 2
 
= kσ + Eµ kg(Y)k + 2σ Eµ g(Y) .
σ
The result follows immediately from application of the following lemma.

65
The following unbiased estimator for the risk is a direct corollary of Stein’s identity.

Corollary 23.4.1 (Stein’s unbiased risk estimate (SURE)). An unbiased estimate of the squared-
error risk of µ̂(Y) = Y + g(Y) is
k
X
SURE(µ̂) = kσ 2 + kg(Y)k2 + 2σ 2 ∇i gi (Y).
i=1

Let’s now prove Stein’s lemma, which will finish the argument for both results.

Lemma 23.5 (Stein’s lemma). Using the notation of the previous lemma, for any i = 1, . . . , k,
  
Y i − µi
Eµ gi (Y) = σ Eµ [∇i gi (Y)] .
σ

Proof. This is a result of integration by parts. Let φ be the density function of the multivariate
standard normal distribution N (0, Ik ), so we have
   Z    
Yi − µ i yi − µi 1 y−µ
Eµ gi (Y) = gi (y) k φ dy
σ Rk σ σ σ
 
y−u
Z
1
= −σ gi (y) k ∇i φ dy
Rk σ σ
 
y−u
Z
1
=σ ∇i gi (y) k φ dy
Rk σ σ
= σ Eµ [∇i gi (Y)] .

Finally, we use the former results to prove Stein’s theorem.

Theorem 23.6 (Stein’s theorem). Stein’s estimator has risk R(µ, µ̂JS ) ≤ k.

Proof. This follows from applying Stein’s identity to the function


k−2
g(y) = − y.
kyk2

After plugging this into Stein’s identity and doing some algebra (taking σ 2 = 1), we get

(k − 2)2
 
 JS 2

Eµ kµ̂ − µk = k − Eµ < k.
kYk2

66
24 December 1st, 2021
Today is the last lecture. We go over some details in proof of Stein’s theorem that we glossed over
last time, discuss the properties of Stein’s estimator, and finally conclude the course.

24.1 Proof of Stein’s Theorem


An astute reader might have noticed that in the previous proof of Stein’s theorem, we seem to
be missing something, since it appears as if µ̂JS even dominates Y for k = 1. The catch here
is that when we applied integration by parts in Stein’s lemma, it requires the function g to be
differentiable, so it only holds for k = 2.
Luckily, we can slightly relax the conditions to only require that for all j, g(y) is differentiable
in its j-th argument and almost surely in the remaining arguments. This criterion also holds for
k ≥ 3, but not when k = 1, which explains the difference in behavior for various dimensions.
Exercise 24.1 (Pencil problem). The risk bound in Theorem 23.6 is not in closed form, as it
involves an expectation, which make it difficult to use for analysis. Can we find a simpler upper
bound for this risk function, still lower than k, but which depends on the value of µ?
Proof. The trick is to use Jensen’s inequality to bound the value of the reciprocal of the squared
magnitude, which is
 
1 1 1 1
Eµ ≥ =  2 = .
kYk2 Eµ [kYk2 ] Eµ χk + kµk 2 k + kµk2
Therefore, the risk can be upper-bounded by
(k − 2)2 (k − 2)2
 
JS
R(µ, µ̂ ) ≤ k − Eµ ≤ k − .
kYk2 k + kµk2
This bound is best when µ = 0 and worse when kµk is large. Intuitively, the reason is that the
shrinkage method always reduces the variance of the estimator by bringing points closer together,
but it works best near the origin, where all values are being moved toward the true mean.

24.2 Properties of the James Stein Estimator


Note that Stein’s estimator µ̂JS turns out to be inadmissible, since we can adjust it slightly to
produce a dominating estimator
 
JS,+ k−2
µ̂ = 1− Y,
kYk2 +
where we use the notation (x)+ = max(x, 0). Now, we show several examples that illustrate that
Stein’s phenomenon occurs in many statistical inference situations in practice, not just the basic
case of minimizing mean-squared error on an isotropic standard normal.
Example 24.1. We can generalize Stein’s estimator for other multivariate normal distributions
Y ∼ N (µ, Σ), by adjusting the `2 norm to the quadratic form
!
k̃ − 2
µ̂JS,Σi = 1 − > Y,
Y ΣY
tr(Σ)
where we define k̃ = λmax (Σ) . This estimator also dominates the simple maximum likelihood
estimator Y with respect to mean-squared error, for all k̃ ≥ 3.

67
Example 24.2. Even when the variance σ 2 of the standard normal distribution N (0, σ 2 Ik ) is
unknown, if we can approximate it by some unbiased sample variance s2 ∼ σ 2 χ2ν /ν, then the
estimator !
ν 2
2 (k − 2) ν+2 s
µ̂JS,s = 1 − Y.
kYk2
dominates the maximum likelihood estimator Y for k ≥ 3.
Furthermore, Stein’s phenomenon holds in generality even for data distributions that are not
multivariate normal, as well as loss functions that are not as heavy-tailed as the mean-squared
error. For example, it holds for Y in dimension k ≥ 3 with respect to
L(µ, µ̂) = log 1 + kµ̂ − µk2 .


Therefore, as a somewhat sobering point, James Stein showed that multi-dimensional decision
analysis requires nontrivial adjustments to minimize risk.
Note. The philosophical message of Stein’s phenomenon is that shrinkage is sometimes desirable,
even from a purely frequentist standpoint. If you care about the loss of your estimator, you should
generally be biased towards shrinking your estimator towards zero. An application of this is in
machine learning, where ridge regression is a form of shrinkage that generally improves error.

24.3 Modern Statistical Research


That concludes the programmed material for this course. Lucas Janson concludes the class by
summarizing the topics we covered so far and providing an overview of modern statistical research.
The goal of learning the material in a class like this one is to understand the underlying moti-
vation behind standard methods in different scientific disciplines. For example, Lucas runs a free
statistical consulting service at Harvard, where most clients come with scientific data that slightly
violates some assumptions of standard methods in their discipline. By understanding the math
behind course topics like sufficient statistics and frequentist inference, we can adjust methods to fit
the specific needs of research.
On the research side, Lucas works in high-dimensional statistics, where current problems are
typically focused on larger domains to handle constantly-growing datasets. Researchers think about
many problems, including issues related to:
• More data: With bigger n, your methods become harder to execute in practice. For example,
computing the covariance matrix is an operation that scales linearly in the number of data
points and quadratically in the number of features.
• More dimensions: With bigger p, you end up with an exponential amount of parameters
in most standard statistical models, such as polynomial regression. Even computing nearest
neighbors in p dimensions is very expensive. This is known as the curse of dimensionality.
• More powerful computers: Given the increasing amount of highly available, large-scale
compute resources, experiments have been getting much larger, which necessitates the devel-
opment of statistical approaches that handle these datasets.
• Machine learning: Neural networks are not well-understood in the current statistical learn-
ing literature, but a key empirical property is that they perform shrinkage to obtain models
of high-dimensional data distributions. Otherwise, their effectiveness would not make sense,
given their lack of adversarial robustness. It’s not understood how this shrinkage happens,
but it occurs implicitly in the training methods (SGD) and architecture (Dropout).

68
• Black boxes: Many statistical learning methods of today are black-box algorithmic ap-
proaches, which perform well empirically but are not well-understood. The effectiveness of
models like CNNs therefore illustrates something about the distributions that they are trained
on, like image classification. In some sense, these are implicit assumptions, but not knowing
the distribution formally makes it hard to apply many of the methods in this class.

Finally, we will briefly introduce an approach that handles the analysis needed for these kinds of
black-box models.

Example 24.3. In conformal inference, we have an arbitrary function f (Xi ; Dy ), and we develop
a variant of the permutation test given this estimator. Using ideas directly from non-parametric
inference and hypothesis testing, we can construct a prediction interval for Ynew | Xnew , while
leveraging the benefits of a black-box machine learning algorithm.

That concludes our statistical inference course for the semester! For the undergraduates in this
class, Lucas reminds us to consider the Concurrent Masters program and suggests that we think
about writing an honors thesis in Statistics.

69
References
[CB21] George Casella and Roger L Berger. Statistical inference. Cengage Learning, 2021.

[DY79] Persi Diaconis and Donald Ylvisaker. Conjugate priors for exponential families. The Annals
of statistics, pages 269–281, 1979.

[LC06] Erich L Lehmann and George Casella. Theory of point estimation. Springer Science &
Business Media, 2006.

[LR06] Erich L Lehmann and Joseph P Romano. Testing statistical hypotheses. Springer Science
& Business Media, 2006.

70

Common questions

Powered by AI

A minimax decision rule is distinctive because it is formulated to minimize the worst-case risk over all possible parameter values. It selects a decision rule that yields the lowest possible maximum risk, unlike other rules that might optimize average performance or cater to specific parameter values, ensuring that no potential parameter configuration results in excessive risk .

Fisher information plays a critical role in determining the asymptotic distribution of test statistics under hypothesis testing. It appears as the variance in the score test, defining its asymptotic normal distribution, and influences the scaling of test statistics in the Wald and other tests, as it quantifies the amount of information that an observable variable carries about an unknown parameter .

Bayes risk in Bayesian decision analysis measures the expected loss of a decision rule averaged over the prior distribution of the parameter. It provides a criterion for evaluating decision rules in a Bayesian framework, allowing for minimization of expected loss considering both the uncertainty in the parameter estimate and the decision consequences, highlighting its practicality in actionable decision-making .

The Lehmann-Scheffé theorem establishes that an unbiased estimator of a function g(θ) that is derived from a complete sufficient statistic is the unique uniformly minimum variance unbiased estimator (UMVUE). This significance lies in identifying the most precise unbiased estimator for a parameter, ensuring minimal variance among all unbiased estimators .

To handle ties in Wilcoxon signed-rank tests with discrete data, one can either randomize tie-breaking by assigning random ranks within tied groups to maintain the distribution's fairness, or by adjusting the rank assignments to the average ranks of the tied positions, which alters the distribution but is deterministic .

The theorem uses the factorization theorem to show that for any pair of data x and y where a sufficient statistic T' satisfies T'(x) = T'(y), the ratio of their probability densities remains free of the parameter, indicating T(x) = T(y). If another statistic T satisfies this same condition, T is shown to be a function of T' due to its coarser fibering, establishing minimal sufficiency .

A Bayes decision rule is constructed by selecting a rule that minimizes the Bayes risk, the expected loss averaged over the parameter's prior distribution. This contrasts with frequentist approaches that ignore prior distributions and instead evaluate rules based on performance over repeated sampling, without incorporating prior beliefs or information in the decision-making process .

Non-parametric tests differ from parametric ones as they rely on fewer assumptions; specifically, they do not assume a specific form for the underlying distribution. This grants them flexibility in assessing broad null hypotheses and retaining power when parametric models fail. However, under correct parametric model specification, non-parametric tests typically have less power due to their general nature .

A minimal sufficient statistic is a statistic that captures the essential information needed for parameter estimation with the least data redundancy. It is unique up to bijective maps, meaning that if two minimal sufficient statistics exist, they are related by a bijective mapping, capturing the same information with respect to the parameter of interest, except in pathological cases .

Complete sufficiency is important because it ensures that no unbiased estimator of zero exists other than the zero function, thus conserving all parameter-related information. It implies minimal sufficiency since a complete statistic cannot destroy parameter information, except in pathological cases. Therefore, a complete sufficient statistic is, by necessity, also minimal .

You might also like