Harvard Statistics 211 Course Notes
Harvard Statistics 211 Course Notes
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
1
6 September 22nd, 2021 16
6.1 The Score Function . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
6.2 Fisher Information . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
2
18 November 8th, 2021 48
18.1 Asymptotic Confidence Intervals . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
18.2 Bayesian Inference . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
18.3 Conjugate Priors of NEFs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
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.
2. Confidence intervals: How confident is our estimate, under Bayesian or frequentist terms?
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.
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).
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).
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:
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,
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.
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
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 .
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.
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 θ = θ(η).
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
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
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
= 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
The second part follows from the law of total variance, since
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.
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.
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.
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.
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
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 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.
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.
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.
θ θ
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.
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 µ.
Cov(Y , M ) = Cov(Y , M − Y + Y )
= Cov(Y , M − Y ) + Cov(Y , Y ).
σ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.
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.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 θ.
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
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.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.
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
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).
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.
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.
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
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.
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!
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,
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.
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 ∞:
(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 → ∞.
(ii) fθ is identifiable.
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 )] < ∞.
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.
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:
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.
(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)
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.
27
11 October 13th, 2021
Today we continue discussing composite hypotheses and the likelihood-ratio test.
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
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
29
Proof. First, observe that we can write
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.
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 :
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
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.
• 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 ).
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.
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
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.
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.)
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.
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.
39
15 October 27th, 2021
Today, we continue discussing selective inference procedures that control the FWER and FDR.
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.
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.
α α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.
Proof. Let k0 be the sorted index of the smallest null p-value. In other words,
k0 = min{k : H0,(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.
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.
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
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.
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)) .
θ∈Θ θ
Example 16.7. If Y1 , . . . , Yn ∼ Unif[0, θ] for θ > 0, and Y( n) = max1≤i≤n Yi , consider the following
two interval estimators:
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 θ.
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.
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.
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.
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.
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 .
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
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 | θ).
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:
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
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.
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:
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
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.
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θ
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.
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 (µ)
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 ,
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 − α.
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.
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.
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
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
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,
where the last step follows because the risk function is a constant for θ ∈ Θ. Therefore, δ 0 dominates
δ ∗ , which is a contradiction, so we conclude.
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 ]
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,
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.
R(θ, δ 0 ) ≤ R(θ, δ π ),
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(π, δ ).
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.
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 }.
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.
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
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
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.
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 µ,
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
τ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
= 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/τ
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) = √ .
2π
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 µ.
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,
Theorem 22.4. If δ π is Bayes with respect to a prior π such that rπ = supθ R(θ, δ π ), then:
(ii) δ π is the unique minimax decision rule if δ π is the unique Bayes rule for π, and
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
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
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 θ.
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).
Definition 23.1 (Least favorable sequence). A sequence of priors πk is called least favorable if
lim rπk ≥ rπ .
k→∞
Theorem 23.2 (Minimax duality). Let πk be a sequence of priors and δ a decision rule. If
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 π,
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.
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)] .
Theorem 23.6 (Stein’s theorem). Stein’s estimator has risk R(µ, µ̂JS ) ≤ k.
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.
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.
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
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 .