8 Estimators

In probability theory we aim to deduce the likelihood of different outcomes based on known probability distributions. Statistics is the inverse problem. We want to infer unknown probability distributions from outcomes we observe. In this chapter we begin to explore these ideas.

8.1The Estimation Problem¶

Let’s start off by clarifying the nature of the problem we want to study.

8.1.1Definitions¶

From a pedagogical point of view, it’s best to start our study of statistics in a setting where successive observations are both independent and identically distributed. Later we’ll extend to dependent data. In an IID setting, the fundamental problem of econometrics and statistics is this:

The set $\pP$ is the universe of distributions we are willing to consider. It can be anything, including the set of all distributions on the outcome space. One of the main tasks of economic theory is to restrict $\pP$ and thereby narrow down the set of distributions we have to search over.

Let’s set out some notation:

• $Z$ is the outcome space where each observation $\boldz_n$ takes values.

• $\zdata$ denotes the sample or data set $(\boldz_1, \ldots, \boldz_N)$.

• $\Zdata := \times_{n=1}^N Z$ is the sample space in which $\zdata$ takes values.

• $\Pdata := \lL(\zdata)$ is the joint distribution of the sample.

Although calling $\Zdata$ the sample space follows standard nomenclature, we are overloading the terminology first used. (Unfortunately, statisticians and probabilists don’t always coordinate on terminology.)

As above, $P = \lL(\boldz_n)$ is the common distribution of the observations. For as long as we continue to focus on IID data, the joint distribution $\Pdata$ will be the $N$th product of $P$, as defined.

• the mean and higher moments of $P$,

• measures of dispersion, such as the variance, or properties of tails,

• the median and other quantiles, and

• $P$ itself, or the density of $P$ if it exists.

In example 8.1.2, the individual observations are scalar. In other scenarios, observations are vector-valued:

Because the marginal distribution $P$ of the observations is now multivariate, new features of the distribution come in to play, such as

• correlations across coordinates of $P$,

• the variance–covariance matrix associated with $P$, and

• parameters controlling dependence when, say, $P$ is modeled via a copula over certain marginals.

8.1.1.1Features¶

We referred above to "features" of $P$ that we might be interested in estimating. Let’s define a feature of $P$ to be an object of the form

The set $S$ is left arbitrary to accommodate all possible features. When $P$ is understood, we’ll write $\gamma(P)$ as $\gamma$. Here are some examples for univariate $P$ routinely estimated in econometric studies:

• $\gamma(P) = \int s^k P(\diff s)$, the $k$th moment of $P$.

• $\gamma(P) = \inf\setntn{s \in \RR}{P(- \infty, s] \geq 1/2}$, the median of $P$.

• $\gamma(P) = P$, when we want to estimate $P$ itself.

• $\gamma(P) =$ the density of $P$ when $P$ is absolutely continuous.

If $P$ is multivariate over $\boldz = (\boldx, y)$, then one feature of interest is the regression function $f^*(\boldx) := \EE[ y \given \boldx]$. This function is uniquely determined by $P$ (see Section 5.2.5).

8.1.1.2Parametric versus Nonparametric Classes¶

In the statistical problems listed above, it is assumed that the unknown distribution belongs to some class $\pP$. We call $\pP$ a parametric class if it can be expressed as

In other words, a class of distributions is parametric if it can be indexed by finitely many parameters. A class of distributions is called nonparametric if it is not parametric.

Many traditional methods of inference are parametric in nature: The data are assumed to be generated by an unknown element $P_{\boldtheta}$ of a parametric class $\pP$. The aim is to estimate $\boldtheta$ using the data. Once we have an estimate $\hboldtheta$ of $\boldtheta$, we can plug this estimate back into the parametric class to obtain an estimate $P_{\hboldtheta}$ of $P_{\boldtheta}$.

The concept of a feature defined in (8.1) is a generalization of the notion of a parameter or vector of parameters. In parametric settings, the feature $\gamma$ that we are most interested in estimating is the parameter vector $\boldtheta$. If we have a good estimate $\hboldtheta$ of $\boldtheta$, we can estimate any feature $\gamma = \gamma(P_{\boldtheta})$ via $\hat \gamma = \gamma(P_{\hboldtheta})$.

Following common usage, we will sometimes refer to the $\boldtheta$ associated with the $P_\boldtheta$ that generates the data as the true value of the parameter vector. However, it’s important to remember that this is an assumption, not a "truth." We assume that the data are generated by some member of a parametric class. This assumption can be completely false.

In our initial formulation of the estimation problem we considered estimation of features rather than just parameters because it is suboptimal to always restrict ourselves to parametric assumptions on $\pP$. We’ll discuss the relative merits of parametric and nonparametric estimation more in chapter 14.

8.1.2Statistics and Estimators¶

A statistic is any observable function of the data. More precisely, a statistic is any $\bB$-measurable function

that can be evaluated once the data $\zdata$ are observed. As in (8.1), the set $S$ is left arbitrary to accommodate all the possible features that we might wish to estimate. Sometimes we write $T$ as $T_N$ to emphasize dependence on the sample size. The assumption that a statistic is $\bB$-measurable is just a basic regularity condition.

An estimator is a statistic used to infer some feature $\gamma(P)$ of an unknown distribution $P$. Thus a statistic becomes an estimator when paired with and compared to a feature of the distribution. Note that there is nothing in the definition of an estimator that implies it will be a sensible estimator of the target feature, let alone a good one.

julia> using Statistics; median([1, 3, 5])
3.0
julia> median([2, 4, 6, 8])
5.0

In some texts the sample variance is defined as $\frac{1}{N-1} \sum_{n=1}^N (x_n - \bar x_N)^2$ rather than (8.2) in order to produce an unbiased estimator. However, (8.2) fits better with the theory in this text and is perhaps the more common definition.

8.1.3Empirical Distributions¶

The empirical distribution of a given $Z$-valued sample $\zdata = (\boldz_1, \ldots, \boldz_N)$ is the discrete distribution on $Z$ that puts equal probability $1/N$ on each sample point $\boldz_n$. Another way to state this is that $\hat P_N$ assigns to each Borel set $B \subset Z$ the number

This is just the fraction of the sample that falls in $B$. (Refer back to figure 4.2 if you wish to see a visual representation.) The expectation of a function $h$ with respect to $\hat P_N$ is

as follows from (5.14).

The empirical distribution is a statistic, mapping observations $\boldz_1, \ldots, \boldz_N$ into $\hat P_N = \frac{1}{N} \sum_{n=1}^N \1\{\boldz_n \in \cdot \}$, a random element of the set of all distributions on $Z$. If we think of $\boldz_1, \ldots, \boldz_N$ as independent draws from common but unknown distribution $P$, then $\hat P_N$ becomes an estimator of $P$. In particular, if the feature of $P$ we want to infer is $P$ itself, the simplest natural estimator is the empirical distribution.^[1]

With scalar data $x_1, \ldots, x_N$ we can visualize the empirical distribution by plotting its CDF. The CDF of $\hat P_N$ will be denoted in what follows by $\hat F_N$. Specializing $B$ in (8.7) to $(-\infty, s]$, we get

The CDF $\hat F_N$ is called the empirical cumulative distribution function, or ECDF, corresponding to the sample. Graphically, $\hat F_N$ is a step function, with an upward jump of $1/N$ at each data point. Figure 8.1 shows an example where draws are from the CDF $F$ of the Beta$(5,5)$ distribution. Each panel shows an observation of the ECDF, which is in turn constructed from $N = 15$ observations from $F$.

8.1.3.1Convergence¶

The empirical distribution is asymptotically an excellent estimator of $P$. If $\boldz_1, \ldots, \boldz_N$ are IID with common distribution $P$ and $\hat P_N$ is the empirical distribution, then, by the law of large numbers,

for any Borel set $B$. See, in particular, (6.8).

Specializing to the scalar case with $B := (-\infty, s]$, we have $\hat F_N(s) \toprob F(s)$ for any $s \in \RR$, where $F$ is the CDF of $P$. In fact a much stronger statement is also true. It is sometimes called the fundamental theorem of statistics, or the Glivenko–Cantelli theorem:

The norm here is the supremum norm, meaning that

Roughly speaking, this is the maximal deviation over the domain. (For a definition of the supremum, see Section 15.4.1.) In fact the convergence occurs “almost surely,” which is a stronger notion than in probability. See theorem 11.4.2 of Dudley (2002).

Figure 8.2 illustrates the convergence in (8.11). Each ECDF $\hat F_N$ is based on a single sample $x_1, \ldots, x_N$ of independent draws from $F$, where $F$ is Beta$(5,5)$. Each subplot shows 8 realizations of the ECDF at a different value of $N$. Across the subplots, the sample size is stepped up from 10 to 500.

Theorem 8.1.1 tells us that, at least in the IID setting, if we have an infinite amount of data, then in the limit we can learn the underlying distribution without having to impose any assumptions at all. Moreover, once we know the distribution $P$, we know any feature $\gamma = \gamma(P)$.

Does this offer a solution to the basic problem of estimation stated? The answer is no. In practice, we only ever have a finite amount of data. Inference is about generalization, as discussed in the introduction. There’s no need to generalize if we know the full population.

Moreover, when we still have only finitely many observations, the empirical distribution puts all of its mass on these points. In other words, it treats the sample like the unknown distribution, without acknowledging the fact that we have limited information. This can be viewed as an extreme form of overfitting. We’ll discuss related ideas in a number of contexts throughout the book.

8.1.3.2Identification¶

A class of distributions $\pP = \{P_{\boldtheta}\}$ indexed by $\boldtheta \in \Theta$ is called identifiable if the map $\boldtheta \mapsto P_{\boldtheta}$ is one-to-one (see Section 15.2) on $\Theta$.

Identifiability means that the parameter vector associated with the unknown distribution can eventually be distinguished from the data. To see this, suppose that $\{P_{\boldtheta}\}$ is identified on $\Theta$ and nature generates an infinite sequence of observations $\{\boldz_n\}$ from $P = P_{\boldtheta}$. Let $\boldtheta'$ be any other vector in $\Theta$ and let $P'= P_{\boldtheta'}$. By identifiability, there exists at least one Borel set $B$ with $P(B) \neq P'(B)$. Since the empirical distribution $\hat P_N(B)$ converges to $P(B)$, we can in the limit conclude that the data are not generated by $P'$.

8.2Estimation Principles¶

We started off with the problem of estimating a feature $\gamma = \gamma(P)$ of the unknown distribution $P$ from a sample. We’ve introduced a number of different statistics aimed at estimating various features. Where do these estimators come from and what, if anything, do they have in common?

8.2.1The Sample Analogue Principle¶

Most of the estimators defined above can be derived from a simple principle. In the modern statistical literature this is often called the plug-in method. In econometrics it goes by a variety of names, including "analogue estimation" (see, in particular, Manski (1986)). We’ll call it the sample analogue principle. The principle is:

Here $\hat P_N$ is the empirical distribution constructed from the sample. We replace the unknown distribution $P$ with the observable distribution $\hat P_N$ and then evaluate $\gamma(\hat P_N)$.

The $k$th sample moment applies the same principle to estimate the $k$th moment.

The estimators in these two examples are relatively straightforward, but the set of all applications of the sample analogue principle is very broad. As we’ll see below, estimation methods that can be obtained as special cases include least squares regression, maximum likelihood, the method of moments, and the generalized method of moments.

8.2.1.1Best Linear Prediction¶

Here’s an example that illustrates the generality of the sample analogue principle. Recall the best linear prediction problem from Section 4.1.5.1, where $\alpha$ and $\beta$ are chosen to minimize $\EE[ (y - \alpha - \beta x)^2 ]$. Letting $P$ be the distribution of $(x, y)$, we can write the problem as

There is a corresponding statistical problem when the underlying distribution $P$ is unknown. The problem is to produce the best linear predictor based only on a sample $\zdata = ((x_1, y_1), \ldots, (x_N, y_N))$ of observations from $P$. We proceed as follows: Given $\zdata$, we form the empirical distribution

Applying the sample analogue principle, we plug $\hat P_N$ into (8.12), yielding the estimator

(The term $\frac{1}{N}$ has been dropped from the right-hand side of (8.13) because it doesn’t affect the minimizers.) This is the simple linear bivariate least squares problem. The minimizers are

(See ex. 8.5.2.) Comparing with (4.23), we see that these estimators are themselves the sample analogues of the underlying features $\alpha^*$ and $\beta^*$.

8.2.1.2Limitations¶

While the sample analogue principle produces sensible estimators in many situations, there are also instances where it fails completely. Let’s look at one example.

Let $\pP$ be the set of absolutely continuous distributions on $\RR$, so that the distribution $P$ generating the data is assumed to have a density. The density of $P$ can be regarded as the feature

Here $DP :=$ the derivative of the CDF $F$ of $P$.

Let $\hat P_N$ be the empirical distribution from any given sample. What happens when we plug this empirical distribution into the right-hand side of (8.15)? In other words, what information do we get when we differentiate the function $\hat F_N$? Since $\hat F_N$ is a step function, the derivative is zero everywhere, except at a finite number of jump points where the derivative is undefined. Thus the density estimate produced by the sample analogue principle is useless. This statement remains true even when $N$ is enormous.

One way to understand this outcome is that when we wish to estimate complex objects with finite samples, we need to combine the data with some kind of regularization. Loosely speaking, this means that we penalize complexity in our search for a solution, or impose some kind of "smoothing" a priori. This idea has a long history in the field of numerical methods, where certain inverse problems are unstable, or "ill-posed," and their solution requires regularization.

Regularization is an essential element of statistics too. It says that we should not grant the empirical distribution equal status with the unknown true distribution. That is, we should regard the empirical distribution only as partial information, and seek to combine it with some form of prior information or external theory. Regularization is one application of prior information: It applies the idea that probability mass most likely falls in places other than just the sample points observed so far. We pick up the topic of regularization again in chapter 14.

8.2.2Empirical Risk Minimization¶

Let’s now look at another inductive principle that combines the sample analogue principle with additional structure. The principle is called empirical risk minimization, or ERM. The terminology and main concepts come from the machine learning literature. Nonetheless, many standard estimators in econometrics are special cases of ERM, including maximum likelihood and least squares. The ideas discussed here help to frame these tools as part of a broader set of methods.

Consider a setting where we observe an input $\boldx \in \RR^K$ to a system, followed by a scalar output $y$. Both are random variables and the joint distribution of $\boldz := (\boldx, y)$ is $P$. Our aim is to predict new output values from observed input values. We’ll do this by choosing a function $f$ such that $f(\boldx)$ is our prediction of $y$ once $\boldx$ is observed. In the machine learning literature, $f$ is called a prediction rule. In economics, $f$ is called a strategy or policy function.

Incorrect prediction incurs a loss. The size of this loss is written as $L(y, f(\boldx))$. The function $L$ is called the loss function. Common choices for the loss function include

• the quadratic loss function $L(y, f(\boldx)) = (y - f(\boldx))^2$,

• the absolute deviation loss function $L(y, f(\boldx)) = |y - f(\boldx)|$, and

• the discrete loss function $L(y, f(\boldx)) = \1\{y \neq f(\boldx)\}$.

Given loss function $L$, we consider choosing $f$ so as to minimize the prediction risk, which is defined as the expected loss

Sometimes $R(f)$ is also called the expected prediction error.^[2] The expectation in (8.16) is computed using the joint distribution $P$ of $(\boldx, y)$. Sometimes we’ll write $\EEP L(y, f(\boldx))$ to emphasize this.

In a statistical setting we are constrained in our ability to minimize risk by the fact that we cannot evaluate $\EEP$, as this requires knowledge of the joint distribution $P$. Suppose, however, that we have access to data $\boldz_1, \ldots, \boldz_N$, where each pair $\boldz_n = (\boldx_n, y_n)$ is an independent draw from $P$. To make use of the sample, we will apply the sample analogue principle, replacing $P$ in (8.16) with the empirical distribution $\hat P_N$. This produces the new objective function

This function is called the empirical risk. The empirical risk is the prediction risk evaluated under the empirical distribution of the sample rather than the true distribution $P$. Minimizing the empirical risk is called empirical risk minimization (ERM).

There is an additional step. When we choose a decision rule the problem that we actually solve is

We have restricted the domain to a class of functions $\hH$. This set of functions is called the hypothesis space. Choice of $\hH$ is structure we impose on the estimation problem.

It might seem at first pass that we should set $\hH$ to be the set of all functions $f \colon \RR^K \to \RR$. After all, if the risk-minimizing function $f^* := \argmin_f R(f)$ is not in $\hH$, as visualized in figure 8.3, then the solution to (8.19) is not equal to $f^*$, and we are making a suboptimal choice.

Actually this reasoning is false. In fact we usually want to be quite restrictive in our choice of $\hH$. This is because minimizing the empirical risk is not the same as minimizing the prediction risk. Occam’s razor comes in to play here: we are solving a complex problem on the basis of limited information. It would be a serious mistake to act as if we had unlimited information. These ideas are explored further in Section 8.2.3.

While we ended up with the same minimization problem in two different ways, the reasoning here is more natural. When we first came to this problem in (8.13) it was by applying the sample analogue principle to the problem of finding best linear predictors. Here we are simply trying to find the best predictor. At the same time we have restricted ourselves to linear approximations in recognition of the fact that we are basing our estimate on limited information. Linearity implements Occam’s razor by bounding the complexity of the learning algorithm.

8.2.2.1Quantile Regression¶

Quantile regression is used when we want to estimate a quantile of a given distribution based on a set of predictive variables. For example, quantile regression has been used to estimate quantiles of CEO compensation as a function of the market value of their firms (Koenker & Hallock (2001)). In this section we show how quantile regression can be regarded as a special case of empirical risk minimization.

To describe the idea, let $F$ be a strictly increasing CDF on $\RR$ and let $\tau \in (0, 1)$ be given. Recall from Section 4.2.6 that the $\tau$th quantile of $F$ is the $\xi$ that solves $F(\xi) = \tau$. Although it might not be immediately obvious, the $\tau$th quantile can also be defined as the solution to the optimization problem

where $y$ is a random variable with distribution $F$ and

Exercise 8.5.6 asks you to confirm that the solution to (8.23) is precisely the $\xi$ that solves $F(\xi) = \tau$.

Suppose now that we want to estimate the $\tau$th quantile of $F$ using some input variable $x$ (e.g., a CEO compensation quantile as a function of firm size). Motivated by (8.23), we can frame the search for a suitable function as a problem of minimizing the prediction risk

where $L_\tau$ is as defined above. If we employ the principle of empirical risk minimization with data set $(x_1, y_1), \ldots, (x_N, y_N)$, we get $\min_{f \in \hH} \; \sum_{n=1}^N L_\tau (y_n , f(x_n))$, where $\hH$ is the hypothesis space. When $\hH = \llL$ as defined in (8.21), the ERM problem is

This is the standard expression for the quantile regression problem.

If $\tau = 0.5$, then the objective function is proportional to $\sum_{n=1}^N | y_n - \alpha - \beta x_n |$. This case is called median regression or least absolute deviation regression.

8.2.3The Choice of Hypothesis Space¶

Let’s return to the issue of choosing the hypothesis space. We begin with a straightforward result:

In other words, we can always decrease empirical risk by increasing the hypothesis space. The reason that fact 8.2.1 holds is that we are minimizing over a larger set.

At the same time, our actual objective is to minimize the prediction risk. The prediction risk of a given predictor $f$ measures the out-of-sample fit of $f$, which is the expected performance of $f$ when confronted with new data. A statistical procedure that takes a given data set and produces a predictor $f$ with low prediction risk (relative to the benchmark $f^*$—see (8.17)) means that we have succeeded in meeting the goal we set at the start of the book: we have generalized from data.

Of course, prediction risk is unobservable. Sometimes it is tempting to use empirical risk $\Remp(f)$ as an estimator of $R(f)$. However, empirical risk, which measures in-sample fit of $f$, is a downward biased estimator of risk. In particular, as we vary the hypothesis space, prediction risk can rise even when empirical risk is falling.

Let’s illustrate this by way of an example, where empirical risk is minimized over progressively larger hypothesis spaces. The model we will consider is one that generates input-output pairs via

Here $U[-1,1]$ is the uniform distribution on the interval $[-1,1]$. Our hypothesis spaces for predicting $y$ from $x$ will be sets of polynomial functions. To fix notation, let $\pP_d$ be the set of all polynomials of degree $d$. That is,

This sequence of hypothesis spaces is increasing, in the sense that

Indeed, if $f$ is a polynomial of degree $d$, then $f$ can be represented as a polynomial of degree $d+1$ just by setting the last coefficient $c_{d+1}$ to zero. The set of linear functions $\llL$ defined in (8.21) is equal to $\pP_1$.

If we seek to predict $y$ from $x$ using quadratic loss and the set $\pP_d$ as our candidate functions, the risk minimization problem is

while the empirical risk minimization problem is

To illustrate the difference between risk and empirical risk, we first generate $N=25$ data points from the model (8.25). Taking this as our data set, we then solve (8.27) repeatedly, once for each $d$ in $1,2,\ldots,20$. The solution to the $d$th minimization problem is denoted $\hat f_d$, and is, by construction, a polynomial of degree $d$. Finally, we compare the risk $R(\hat f_d)$ and empirical risk $\Remp(\hat f_d)$.^[3] The results are in figure 8.4.

As expected, empirical risk falls monotonically with $d$. But the risk decreases slightly and then increases rapidly. For large $d$, the minimizer $\hat f_d$ of the empirical risk is associated with very high risk in the sense of large expected loss.

We can get a sense for what is happening by plotting the data and the functions. In figures 8.5–8.8, the $N$ data points are plotted alongside the function $y = \cos(\pi x)$ from the true model (8.25) in black, and fitted polynomial $\hat f_d$ in red. The function $y = \cos(\pi x)$ is the risk minimizer, and represents the ideal prediction function. In figure 8.5 we have $d=1$, and the fitted polynomial $\hat f_1$ is the linear regression line. In figures 8.6, 8.7 and 8.8 we have $d = 3$, $d=11$ and $d=14$ respectively. The fitted polynomials are $\hat f_3$, $\hat f_{11}$ and $\hat f_{14}$.

On the one hand, when $d=1$, the hypothesis space $\pP_d = \pP_1$ is small and no function in this class can effectively fit the underlying model. This is called underfitting, and is reflected in the poor fit of the red line to the black line in figure 8.5.

When $d=3$, the class of functions $\pP_d = \pP_3$ is considerably larger. Given that the data are relatively noisy, and that we only have 25 observations, the fit of the function is good (figure 8.6). If we look at the risk for $d=3$ on the right-hand side of figure 8.4, we see that it is lower than for $d=1$.

On the other hand, if we take $d$ larger, the fit to the underlying model becomes poor and the risk is high. Examining figures 8.7 and 8.8, which correspond to $d=11$ and $d=14$, we see that the fitted polynomial has been able to fit the observed data closely, passing near many of the data points. Too much emphasis has been given to this particular realization of the data. When a new input $x$ is drawn, the prediction $\hat f_d(x)$ is likely to be a poor predictor of $y$. This situation is called overfitting.

To summarize, the choice of $\hH$ is central to our ability to generalize from the data. When $\hH$ is too small, no function in $\hH$ provides a good fit to the regression function. If this is so, the prediction risk cannot be made small regardless of the data. This can be seen directly from (8.17). Any estimate $\hat f$ in $\hH$ must satisfy

Conversely, if $\hH$ is too large, then we can attain low empirical risk but the prediction risk itself is large.

Of course, in statistical applications we do not have the luxury of knowing the true data-generating process when we choose $\hH$. The best scenario is that we have firm theory that guides us to a suitable hypothesis space. The worst scenario is that we have no idea and choose blindly. Once again, the message is that statistical learning equals prior knowledge plus data.

8.3Some Parametric Methods¶

In this section we review some standard parametric estimation methods, including maximum likelihood, Bayesian estimation, and the generalized method of moments.

8.3.1Maximum Likelihood¶

One standard approach to deriving estimators in a parametric setting is the principle of maximum likelihood. Let’s begin with a simple case and then progress to more complex applications.

Suppose that we observe a single draw $x$ from distribution $\nN(\mu,1)$ where $\mu$ is unknown. The variance is known and set to $1$ for simplicity. Our task is to guess the value of $\mu$ given the observation $x$. A guess $\hat \mu$ of $\mu$ also pins down the distribution $\nN(\hat \mu, 1)$, so we could equivalently state that our job is to guess the distribution that generated $x$. In guessing this distribution, if we centered it around some number $\hat \mu$ much larger than $x$, say, then our observed data point $x$ would be an “unlikely” outcome for this distribution. See figure 8.9. The same logic would apply if we centered the density at a point much smaller than $x$.^[4]

In fact, in the absence of any additional information, the most obvious guess would be that the normal density is centered on $x$. To center the density on $x$, we just set $\hat \mu = x$.

Maximum likelihood formalizes these steps: The density of $x$ is the density $p(s ; \mu)$ of the distribution $\nN(\mu, 1)$. Plugging the observed value $x$ into this density gives

Think of $p(x ; \mu)$ as representing the probability of realizing our sample point $x$. The principle of maximum likelihood suggests that we take as our guess $\hat \mu$ of $\mu$ the value that maximizes this probability. It’s not hard to show that $\hat \mu = x$ is the maximizer. This coincides with our intuition from figure 8.9.

The same principle applies when the data are independent draws $x_1, \ldots, x_N$ from $\nN(\mu, 1)$, where $\mu$ is unknown. The joint density of the sample is the product of the marginals. Plugging the sample values into the joint density, one then maximizes the joint density with respect to $\mu$:

The maximizer $\hat \mu$ is precisely the sample mean of $x_1, \ldots, x_N$ (see ex. 8.5.4).

We can generalize these ideas in several ways. Let’s suppose now that the data $x_1,\ldots,x_N$ has joint density $p$ in the sense of (5.8). We will assume that $p = p(\cdot \,; \boldtheta)$ is a member of a parametric class $\pP$ indexed by parameter vector $\boldtheta \in \Theta$. Each choice of $\boldtheta$ pins down a particular density $p = p(\cdot \,; \boldtheta)$, but the value of $\boldtheta$ that generated the data is unknown.

In this setting, the likelihood function is $p$ evaluated at the sample $x_1,\ldots,x_N$, and regarded as a function of $\boldtheta$:

The principle of maximum likelihood tells us to estimate $\boldtheta$ by maximizing $L(\boldtheta)$ over $\boldtheta \in \Theta$. A statistic $\hat \boldtheta$ is called a maximum likelihood estimate (MLE) of $\boldtheta$ if

This is equivalent to maximizing the log likelihood function

The set of MLEs can in general be a singleton, contain multiple elements or be empty.

In the preceding discussion, $p$ was a density function, but it can be a probability mass function as well. Examples are given below.

To implement maximum likelihood estimation, we need to compute the joint distribution of the data. As we saw in (8.28), when the data points are independent this is easy because the joint density $p$ is the product of the marginals. More generally, if each $x_n$ is drawn independently from fixed arbitrary (marginal) density $p_n(\cdot \,; \boldtheta)$ on $\RR$, then

If each data point is multivariate, just replace $x_n$ with $\boldx_n$.

Maximum likelihood estimation is a much celebrated theory of estimation. MLE estimators typically have excellent asymptotic properties. But good finite sample properties are not guaranteed (more about this soon) and the attractive asymptotic theory is dependent on correct specification of the underlying parametric class. Indeed we need to bring a lot of knowledge to the table just to form the estimator. (To specify the likelihood function, we must specify the entire joint distribution of the sample.)

8.3.1.1Conditional Maximum Likelihood¶

As in Section 8.2.2, suppose that we observe inputs $\boldx_1, \ldots, \boldx_N$ to some system and corresponding outputs $y_1, \ldots, y_N$. The pairs $(\boldx_n, y_n)$ are assumed to be IID. Our aim is to estimate $\boldtheta$ in $p(y \given \boldx ; \boldtheta)$ in order to pin down the conditional density of $y$ given $\boldx$.

The principle of maximum likelihood tells us to maximize

Letting $\pi$ be the marginal density of $\boldx$, we can use the decomposition (5.27) to write

The density $\pi$ is unknown but we have not parameterized it because we aren’t trying to estimate it. We can now rewrite the log likelihood as

The second term on the right-hand side is independent of $\boldtheta$ and as such it does not affect the maximizer. Hence the MLE is

The objective function here is called the conditional log likelihood. The preceding argument tells us that, when we want to estimate parameters in the conditional density of $y$ given $\boldx$, we can maximize the conditional log likelihood directly, which is simpler and more direct than maximizing the full log likelihood.

To find the MLE in this example, we can differentiate $\ell$ to obtain first order conditions, but this does not in general lead to an analytical solution. Instead, numerical optimization of the likelihood is required. Discussion of numerical optimization is given in Section 13.2.2.

8.3.2Maximum Likelihood via ERM¶

Maximum likelihood is a special case of empirical risk minimization. To see this, suppose that we wish to learn unknown density $q$ on the basis of observations from $q$. We take our loss function to be

The meaning is this: if our guess of $q$ is $p$ and the value $x$ is realized, then our loss is $- \ln p(x)$. Thus we suffer large loss when $p$ puts small probabilities near where $x$ is realized. The corresponding risk function is

Now suppose that we observe IID draws $x_1, \ldots, x_N$ from $q$. To estimate $q$, the ERM principle indicates we should solve for

where $\pP$ is a hypothesis space of densities. To clarify the connection with maximum likelihood, take $\pP$ to be a parametric class $\{p(\cdot \,; \boldtheta)\}_{\boldtheta \in \Theta}$. Choosing our estimate $\hat p$ of $q$ now reduces to choosing an estimate $\hat \boldtheta$ of $\boldtheta$. Rewriting our optimization problem for this case, we obtain

Here $\ell$ is the log-likelihood. It follows from this expression that the ERM estimator is precisely the maximum likelihood estimator.

Incidentally, the minimizer of the true risk $R(p)$ is the unknown density $q$. To see this, let us first transform our expression for the risk function to

The term on the far right is called the entropy of $q$, and does not involve $p$. Hence minimization of the risk comes down to minimization of

This quantity is called the Kullback–Leibler (KL) deviation between $q$ and $p$. The KL deviation is possibly infinite, always nonnegative, and zero if and only if $p = q$.^[5] It follows that the unique minimizer of the risk is the true density $q$.

8.3.3The Method of Moments and GMM¶

Suppose that we wish to estimate a vector $\boldtheta$ that solves an equation of the form

In this expression, $g$ and $h$ are observable vector-valued functions (taking values in $\RR^K$, say). We cannot solve this expression because the distribution $P$ of $\boldx$ is unknown, and hence the expectation cannot be evaluated. If, however, we have observations $\boldx_1, \ldots, \boldx_N$ from $P$ then we can apply the sample analogue principle and replace the population expectation with the expectation under the empirical distribution. This yields the method of moments estimator, which is the solution $\hat \boldtheta$, if it exists, to the equation

Generalized method of moments (GMM) is a small step from method of moments. If we express (8.35) as $\EE[g(\boldtheta) - h(\boldx)] = \boldzero$, then it becomes natural to consider the more general expression

This expression is called the orthogonality condition. The generalized method of moments estimator of $\boldtheta$ is the solution $\hat \boldtheta$ to the empirical counterpart, which is

Of course, there is no guarantee that a solution will exist here, partly because the function $G$ can be nonlinear and partly because the number of equations can be greater than the number of unknowns. If the number of equations is greater, then the estimation problem is said to be overidentified.

Our study of overdetermined systems of equations in Section 3.3.2 suggests a logical way to handle the overidentified case: Minimize the norm of the left-hand side of (8.38). This leads to the expression

In practice, we usually make two adjustments to this expression. The first is trivial: we minimize the squared norm, instead of the norm, leaving the minimizer unchanged. The second is to replace the Euclidean norm $\| \cdot \|$ with a weighted norm $\| \cdot \|_W$ defined by $\| \boldx \|_W^2 = \boldx^\T \boldW \boldx$, where $\boldW$ is a positive definite weighting matrix. Positive definiteness is desired here because it shares with the Euclidean norm the property that led us to (8.39). That is, $\| \boldx \|_W = 0$ if and only if $\boldx = \boldzero$. The estimation problem is then

The weighting matrix has been written as $\hat \boldW$ because it is allowed to depend on the sample. Evidently the choice of $\hat \boldW$ affects the minimizer, and the objective in choosing this matrix is to produce an estimator that has small variance asymptotically.

8.3.4Bayesian Estimation¶

Bayesian inference adopts a rather different strategy from the methods discussed so far. The main idea is to treat parameters as unknown quantities for which we hold subjective beliefs regarding their values. These subjective beliefs are called priors. The Bayesian approach to estimation suggests that we take both data and prior knowledge into account when forming an estimate or prediction.

For the purposes of estimation, a prior can be thought of as a distribution over $\pP$, the set of distributions in play. The standard Bayesian approach is parametric, so we can specialize this further to a density over parameter space. Thus the primitives in our analysis are:

• $\boldtheta$, the parameter vector, which takes values in $\Theta \subset \RR^J$,

• $\pi$, the prior distribution, a density over $\Theta$,

• $\boldx$, the data, and

• $p(\cdot \given \boldtheta)$, the joint density of the data given $\boldtheta$.

Note that $L(\boldtheta) := p(\boldx \given \boldtheta)$ is the likelihood function.

Priors are reassessed based on evidence in the data. This process leads to an updated density over parameter space called the posterior distribution, which we represent by $\pi(\boldtheta \given \boldx)$. The posterior is obtained via an application of Bayes’ law, which leads us to

Here $p(\boldx)$ represents the unconditional density of $\boldx$ evaluated at the outcome. The term on the far right shows why this unconditional density is not listed as a primitive: we can recover it from the other primitives using the law of total probability.

The same method can be applied when the densities are replaced with probability mass functions and the integral is replaced with a sum. Here’s a standard example that mixes a density over priors with a binomial PMF for the likelihood:

Figure 8.10 shows evolution of the posterior density (8.44) in a simulation. The prior is set to Beta$(3, 5)$. The true payout probability is $\theta_0 = 0.7$. Despite the poor prior, the data shift probability mass towards $\theta_0$.

Point estimates are extracted from the posterior distribution based on some measure of central tendency, such as the mean, the median, or the mode of the posterior (i.e., the maximizer in the unimodal case).

The posterior we came up with in (8.44) belongs to the same parametric class as the prior. Priors where the parametric class is preserved under Bayesian updating for a specific likelihood function are called conjugate. Conjugacy makes application of Bayes’ law particularly simple.

In applications, the existence of a closed-form solution for the posterior is rare, and the integration over the parameter space in (8.41) has to be carried out numerically. The standard technique is Markov chain Monte Carlo. An exposition of the Metropolis–Hastings algorithm was given in Section 7.4.2. If you return to (7.38), you will notice how implementation of this algorithm only requires us to evaluate ratios of the posterior, which means that the integral term on the right-hand side of (8.41) cancels out. Because the integral can be very high dimensional, this often eliminates a huge amount of complexity.

Bayesian inference has experienced a surge of popularity in recent years. One reason is that in high dimensions, exploring the posterior via MCMC has proved to be more successful in practice than the numerical optimization required to obtain maximum likelihood estimates. Another is that Bayesian estimation provides a form of regularization that stabilizes and typically improves estimation of complex models. See Section 14.2.3. A third is that Bayesian estimation comes with an elegant, unified decision-theoretic approach to inference. Some discussion is given in Section 9.2.5.

8.4Further Reading¶

The ERM principle is a very general principle for solving statistical problems and producing estimators. For such a general method it is difficult to give a set of strong results showing that ERM produces good estimators. Indeed, there will be instances when ERM produces poor estimators, as discussed in Section 8.2.3. Having said that, some rather general consistency results have been obtained. Additional discussion can be found in Vapnik (2000).

There are many good treatments of Bayesian estimation available in the literature, including Geweke (2005), Geweke et al. (2011), and Kroese & Chan (2014).

Chernozhukov & Fernández-Val (2011) discuss quantile regression of the $\tau$th quantile in cases where $\tau$ is close to zero or one.

8.5Exercises¶

8.5.1Solutions to Selected Exercises¶