2 Frequentist Parameter Estimation

2.1 Overview

Parameter estimation can be used to infer unknowns about the real world (e.g, the frequency of a given disease among individuals with a certain genetic mutation) and to estimate the distribution of the data in machine learning problems.

There are two main frameworks for parameter estimation:

Frequentist methods: In the frequentists’ perspective, the true parameter value \(\theta^*\) is unknown and fixed. The estimate \(\hat\theta\) is a function of the data, which provides a single “best” estimate of \(\theta^*\). Frequentists have different methods for estimation including maximum likelihood, which we will discuss in detail, and the moment method, which finds the parameters by solving equations obtained by equating empirical moments and theoretical moments.
Bayesian methods: Parameters are considered to be random and are treated as such. The Bayesian method provides a unified approach consisting of the following steps:
1. Start with the prior distribution for the parameter
2. Collect data
3. Obtain posterior distribution by updating the prior distribution using data and Bayes’ theorem

Let us consider a problem about inferring unknown values and making decisions and use probability to solve it, using both frequentist and Bayesian views. Suppose that the probability that someone with a given allele of a gene will develop a certain disease is \(\theta\). We are interested in determining \(\theta\). Different interpretations lead to different approaches to problems. But to determine \(\theta\), both frequentists and Bayesians need data.

Data (\({\mathcal{D}}\)). Among a sample of 100 people with this allele, 2 had the disease.

A Frequentist thinks of \(\theta\) as unknown non-random parameter. She starts by asking “What is the probability of the observation as a function of \(\theta\)?” We can view each of the 100 people chosen to be an independent Bernoulli trial with probability \(\theta\). So the distribution is Binomial and the probability of the observation as a function of \(\theta\) is \[L(\theta) = \binom{100}{2}\theta^2 (1-\theta)^{98}.\] Probability of the observation as a function of the parameter is called the likelihood function. So what value for \(\theta\) makes the most sense? Since the observation has actually happened, we would expect it to have a high probability so we find \(\theta\) that maximizes the likelihood. This method is called maximum likelihood estimation, and we’ll discuss it in much more detail later. In this case, we estimate \(\theta\) to be \[\hat\theta = \arg\max_\theta L(\theta) = \frac{2}{100},\] which is a reasonable estimate. But how close is the estimate to the true value? For frequentists, this is a tricky question to answer probabilistically since the true value and the estimate are both deterministic at this point. With some clever reasoning (some would say mental gymnastics), frequentists come up with confidence intervals and confidence levels to quantify the accuracy of estimators.
A Bayesian thinks of \(\theta\) as random and assigns to it a distribution, called the prior, before seeing the data. Thinking of \(\theta\) as random is imaginative (some would say questionable) since there is no repeatable experiment and there is a single value that is true. One way to justify randomness of \(\theta\) is to think of our universe being drawn from a set of possible universes. Regardless, the Bayesian view is used widely in practice.

Our Bayesian statistician then looks at the data and updates her distribution for \(\theta\), thus obtaining the posterior distribution. Assume that before seeing the data, we believe that the distribution for \(\theta\) is uniform, i.e., \(p(\theta)\sim \mathop{\mathrm{\mathrm{Uni}}}[0,1]=\mathop{\mathrm{\mathrm{Beta}}}(1,1)\). This means that while we do not know what \(\theta\) is, we believe it is equally likely to be any value between 0 and 1. When we see the data, we can update this belief, \[p(\theta|{\mathcal{D}}) = \frac{p({\mathcal{D}}|\theta)p(\theta)}{p({\mathcal{D}})} \qquad\text{(Bayes' rule)}\] It turns out \(p(\theta|{\mathcal{D}})\sim \mathop{\mathrm{\mathrm{Beta}}}(3,99)\).

In contrast to the frequentist view, the Bayesian view is consistent and flexible: any question about \(\theta\) is answered by a computation with the posterior.

What is the source of uncertainty in this problem? It is the finite sample size. If we know the status of a very large number of people with the allele, we would know the distribution/ the value of \(\theta\).

2.2 Maximum likelihood estimation

Suppose data \(x\) is collected. We model this data as a realization of a random variable \(X\) with distribution \(p_{X}\), which has an unknown parameter \(\theta^*\). The probability of observing \(x\), assuming \(\theta\), is \(p_X(x;\theta)\). To estimate \(\theta^*\), Maximum likelihood estimation (MLE) chooses the parameter that assigns the highest probability to the data: \[\begin{equation*} {\hat\theta_{{\rm mle}} = \arg\max_\theta p_X(x;\theta)}. \end{equation*}\] The expression \(p(x;\theta)\), viewed as a function of \(\theta\), is called the likelihood; hence the name maximum likelihood estimation. As shorthand, we use \(L(\theta) = p_X(x;\theta)\) and \(\ell(\theta) = \ln L(\theta)\), where \(\ell(\theta)\) is the log-likelihood. Clearly, the value of \(\theta\) that maximizes \(L(\theta)\) is the same as the one that maximizes \(\ell(\theta)\): \[\begin{equation*} {\hat\theta_{{\rm mle}} = \arg\max_\theta \ell(\theta) = \arg\max_\theta \ln p_X(x;\theta)} \end{equation*}\]

Example 2.1 Heavy or light traffic

In this example, we attempt to show the intuition behind maximum likelihood. Suppose that a given road has heavy traffic or light traffic. We denote the probability of light traffic by \(\theta^*\). To estimate data, we count the number of times \(X\) that the road has light traffic in a period of 100 days. After collecting this data, we observe that \(X=65\). We have \[\begin{align*} p_X(65; \theta) &= \binom{100}{65}\theta^{65} (1-\theta)^{35} \end{align*}\] Let’s try a few different choices for \(\theta\), e.g., \(\theta\in\{0.2,0.4,0.6,0.8\}\), and see which one makes more sense: \[\begin{align*} p(65;\theta=0.2) &= 1.6\times 10^{-22}, \\p(65;\theta=0.4) &= 0.00000026, \\p(65;\theta=0.6) &= 0.0491, \\p(65;\theta=0.8) &= 0.00019, \end{align*}\]

In the figure below, \(p({{\boldsymbol{x}}};\theta)\) is plotted for \(\theta\in\{0.2,0.4,0.6,0.8\}\). The vertical line indicates the observation, i.e., 65 days with light traffic. Which is a more appropriate value for \(\theta\)?

If \(\theta=0.2\), the probability of 65 days with light traffic is extremely small. So observing \(x=65\) would be very unlikely, which in turn would make \(\theta=0.2\) an unreasonable guess. Among the presented choices, \(\theta=0.6\) appears the most reasonable. This reasoning suggests the following: The value of the parameter that assigns a higher probability to the observation is a better choice.

Since we are not limited to a specific set of choices, we can find the parameter that maximizes the probability of the observation. In the figure below, \(L(\theta) = p(x;\theta)\) is plotted as a function of \(\theta\). This is the likelihood function.

We can see that \(\theta=0.65\) maximizes the likelihood and hence is the maximum-likelihood estimate. We can also show this analytically. First, the likelihood is given as \[L(\theta) = p(x; \theta) = \binom{100}{65}\theta^{65} (1-\theta)^{35}.\] We usually use the log-likelihood as the function to optimize: \[\begin{equation} \ell(\theta) = \log{L(\theta)}=\log\left(\binom{100}{65}\theta^{65} (1-\theta)^{35}\right)\doteq 65\log\theta+35\log (1-\theta), \end{equation}\] where \(\doteq\) denotes equality but with ignoring additive terms that are constant in \(\theta\) (and thus do not alter the value of \(\theta\) that maximize the log-likelihood). We differentiate \(\ell(\theta)\) to find the value of \(\theta\) that maximizes \(l(\theta).\) \[\begin{equation} \frac{d\ell(\theta)}{d\theta} = \frac{65}{\theta} - \frac{35}{1-\theta} = 0 \implies 65-65\theta=35\theta\implies \hat\theta_{{\rm mle}} = \frac{65}{100}. \end{equation}\tag{2.2}\]

Note that this result is intuitive as it agrees with our observation that 65% of the days had light traffic.

A note on notation. We indicate the parameter usually as an extra argument of the probability function, e.g., \(p(x;\theta)\), separated by a semicolon to indicate that it is a parameter not a realization of a random variable. Furthermore, in general, our data is a vector, which we denote by bold symbols such as \({{\boldsymbol{x}}}\). The corresponding random variable is \({{\boldsymbol{X}}}\).

2.2.1 MLE with a finite parameter space

In many problems the parameter can only take values in a finite set \(\Theta=\{\theta_1,\dots,\theta_m\}\). In this setting, maximum likelihood is conceptually simple: evaluate the likelihood for each candidate and pick the one with the largest value \[\hat\theta_{\rm mle} \,=\,\arg\max_{\theta\in\Theta} L(\theta) \,=\,\arg\max_{\theta\in\Theta} \ell(\theta).\]

Example 2.2 N-sided die with finite candidates

An \(N\)-sided die is rolled once and we observe \(X=x\). The parameter \(N\in\Theta=\{4,6,8,12,20\}\). For each possible value of \(x\), let us find the maximum-likelihood estimate of \(N\).

Solution

Under the model of a fair die, \[p(x;N)=\begin{cases} \tfrac{1}{N}, & 1\le x\le N,\\ 0, & \text{otherwise.} \end{cases}\] Hence the likelihood is zero for any \(N<x\), and equals \(1/N\) for \(N\ge x\), so MLE selects the smallest candidate \(N\) that is still feasible: \[\hat N_{\rm mle}=\min\{N\in\Theta: N\ge x\}.\] For example, if \(x=7\) then \(\hat N_{\rm mle}=8\) and if \(x=1\) then \(\hat N_{\rm mle}=4\).

Example 2.3 Fair vs. bent coin

The coin is either fair or bent with \(p(H)=0.6\), so \(\Theta=\{0.5,0.6\}\). Consider two experiments.

Solution

One flip (\(n=1\)): If \(X=H\), compare \(0.6\) vs. \(0.5\) and choose \(\hat\theta_{\rm mle}=0.6\); if \(X=T\), compare \(0.4\) vs. \(0.5\) and choose \(\hat\theta_{\rm mle}=0.5\).
Ten flips (\(n=10\)): Let \(k\) be the number of heads. Compare the log-likelihoods using an annotated comparison: \[\ell(\theta\!=\!0.6) \underset{\theta=0.5}{\overset{\theta=0.6}{\lessgtr}} \ell(\theta\!=\!0.5) \iff k\log\tfrac{0.6}{0.5}+(10-k)\log\tfrac{0.4}{0.5} \underset{\theta=0.5}{\overset{\theta=0.6}{\lessgtr}} 0\] Equivalently, \[k \underset{\theta=0.5}{\overset{\theta=0.6}{\gtrless}} \frac{10\log(5/4)}{\log(3/2)}\approx 5.5,\] so we choose the bent coin if \(k\ge 6\); otherwise the fair coin.

We note that MLE for finite parameter spaces is closely related to hypothesis testing, which we will not study in depth. Additional concepts in hypothesis testing include likelihood ratio, type I and type II errors, trade-off between these errors, and the Neyman-Pearson lemma.

2.2.2 MLE with a continuous parameter space

Example 2.4 Parameters of the normal distribution

A device for measuring an unknown quantity \(\mu^*\) (e.g., the mass of an electron) is used \(n\) times producing values \({{\boldsymbol{Y}}}=\left(Y_1,\dotsc,Y_n\right)\). Each measurement is independent and for each \(i\) we have \(Y_i=\mu^*+Z_i\), where \(Z_i\) is the measurement noise satisfying \(Z_i\sim\mathop{\mathrm{{\mathcal{N}}}}(0,{\left( \sigma^* \right)}^2)\). Note that this implies \(Y_i\sim\mathop{\mathrm{{\mathcal{N}}}}(\mu^*,{\left( \sigma^* \right)}^2)\).

Suppose we have collected data \({{\boldsymbol{y}}}=(y_1,\dotsc,y_n)\). We consider the problem in two cases: \(\mu^*\) is unknown but \({\sigma^*}\) is known; and both \(\mu^*\) and \({\sigma^*}\) are unknown.

Solution

Known \({\sigma^*}\), unknown \(\mu^*\): We have \[\begin{align*} p_{Y_i}(y_i;\mu) &= \frac{1}{{{\sigma^*}}\sqrt{2\pi}}\exp\left(-\frac12\left(\frac{y_i-\mu}{{{\sigma^*}}}\right)^2\right)\\ L({\mu})&= p_{{{\boldsymbol{Y}}}}({{{\boldsymbol{y}}}};\mu) = \prod_{i=1}^n p_{Y_i}(y_i;\mu)\\ \ell(\mu)&= \sum_{i=1}^n \ln p_{Y_i}(y_i;\mu) = \sum_{i=1}^n \left(-\ln({{{\sigma^*}}\sqrt{2\pi}})-\frac12\left(\frac{y_i-\mu}{{{\sigma^*}}}\right)^2\right)\doteq - \frac12\sum_{i=1}^n \left(\frac{y_i-\mu}{{{\sigma^*}}}\right)^2\\ \end{align*}\] and so \[\begin{align*} \frac{d\ell}{d\mu} = \sum_{i=1}^n \frac{y_i-\mu}{{{\sigma^*}}}=0\implies \hat\mu_{{\rm mle}} = \frac{1}{n}\sum_{i=1}^n y_i={\textcolor{blue}{\bar y}}. \end{align*}\]
Unknown \({{\sigma^*}},\mu^*\): We have \[\begin{align*} \ell(\mu,{\sigma})&= \sum_{i=1}^n \left(-\ln({{\sigma}\sqrt{2\pi}})-\frac12\left(\frac{y_i-\mu}{{\sigma}}\right)^2\right)\doteq -n\ln {\sigma} - \frac12\sum_{i=1}^n \left(\frac{y_i-\mu}{{\sigma}}\right)^2\\ \end{align*}\] and so \[\begin{align*} \frac{\partial \ell}{\partial \mu} &= \sum_{i=1}^n \frac{y_i-\mu}{{\sigma}}=0,\\ \frac{\partial \ell}{\partial \sigma} &= -\frac{n}{{\sigma}}+\sum_{i=1}^n \frac{(y_i-\mu)^2}{{\sigma}^3}=0. \end{align*}\] Solving this system of equations for \(\mu\) and \({\sigma}\) yields \[\begin{align*} \hat\mu_{{\rm mle}} &= \frac{1}{n}\sum_{i=1}^n y_i={\textcolor{blue}{\bar y}},\\ \hat\sigma_{{\rm mle}}^2 &= \frac{1}{n}\sum_{i=1}^n (y_i-\bar y)^2. \end{align*}\]

2.2.3 Maximum likelihood and the closest distribution

We have described maximum likelihood as aiming to find a distribution that gives a high probability to the observed data. An alternative view relates it to the empirical distribution of the data, denoted \(p_{{\boldsymbol{x}}}\). Given \({{\boldsymbol{x}}}=\{x_1,\dotsc,x_n\}\), let \(\#_x\) denote the number of times \(x\) appears in \({{\boldsymbol{x}}}\). The empirical distribution is given as \[p_{{\boldsymbol{x}}}(x) = \frac{\#_x}{n} = \frac{1}{n}\sum_{i = 1}^n 𝟙{\left( {{\boldsymbol{x}}}= {{\boldsymbol{x}}}_i \right)}\] where \(𝟙(\cdot)\) equals 1 if the enclosed condition is true and 0 otherwise.

Now consider a parameterized family of distributions \(p_\theta\). It makes sense to choose \(\theta\) such that \(p_\theta\) is close to \(p_{{\boldsymbol{x}}}\). In other words, we would like \(\theta\) to be chosen such that \(p_\theta\) describes the observed data well. A standard way of measuring the “closeness” of \(p_\theta\) to the empirical distribution \(p_{{{\boldsymbol{x}}}}\) is relative entropy (reviewed in Chapter 1), \(D{\left( p_{{{\boldsymbol{x}}}}||p_\theta \right)}\).

It turns out the closest distribution is in fact given by maximum likelihood, i.e., \[\begin{equation} \hat\theta_{{\rm mle}} = \arg\min_{\theta} D{\left( p_{{{\boldsymbol{x}}}}||p_\theta \right)}. \end{equation}\tag{2.3}\] This fact provide further evidence for the soundness of MLE strategy. Note in particular that if there exists \(\theta\) such that \(p_\theta = p_{{\boldsymbol{x}}}\), it will be chosen by MLE. This is because relative entropy is always non-negative and equals to 0 if and only if the two distributions are the same. So choosing \(p_\theta = p_{{\boldsymbol{x}}}\), if possible, provides the smallest value for the relative entropy, i.e., 0.

Exercise 2.5 Likelihood and relative entropy

Prove (2.3).

Exercise 2.6 The other direction

Note that relative entropy is not symmetric. Instead of \(D{\left( p_{{{\boldsymbol{x}}}}||p_\theta \right)}\), we could minimize \(D{\left( p_\theta||p_{{{\boldsymbol{x}}}} \right)}\). What are the differences between the two formulations and which one is more suitable for parameter estimation?

Plug-in prediction. The estimated distribution is directly useful for prediction. In a supervised problem with features \(X\) and target \(Y\), model the joint distribution as \(p_\theta(x,y)\), estimate \(\hat\theta\) by maximum likelihood, and then predict by treating \(p_{\hat\theta}\) as the truth: for squared loss, \(\hat f(x)=\mathop{\mathrm{{\mathbb{E}}}}_{\hat\theta}[Y\mid X=x]\), and for 0-1 loss, \(\hat f(x)=\arg\max_{y} p_{\hat\theta}(y\mid x)\), the optimal rules of Chapter 1, computed under the estimated distribution. This plug-in approach is one of the two main routes from data to predictions, the other being empirical risk minimization (see the linear regression chapter).

2.3 Properties of estimators

Maximum likelihood is just one way of estimating parameters. We can choose any function of the data as the estimate. For instance, in Example 2.4, we could choose the middle (median) value among \(y_1,\dotsc,y_n\) as the estimate for \(\mu^*\). Given the fact that there are many estimators, how do we evaluate them and select one?

Clearly, we would like the estimate to be close to the true value. But stating this condition in a rigorous probabilistic way is a bit challenging in the frequentist framework. We are specifically interested in the error: \[\hat\theta({{\boldsymbol{x}}}) - \theta^*,\] where \(\hat\theta({{\boldsymbol{x}}})\) is the estimate based on data \({{\boldsymbol{x}}}\) and \(\theta^*\) is the true value¹. Evaluating \(\hat\theta({{\boldsymbol{x}}})\) is difficult because, obviously, the true value is unknown.

So instead of finding the specific error, we may try to find the probability that the true value \(\theta^*\) is within say 10% of the estimate \(\hat \theta\). But after the estimate is produced based on a given data set, the estimate is a deterministic value. For instance, in Example 2.1, the MLE is given as \(\hat\theta_{{\rm mle}}=0.65\). So questions such as “What is the probability that the difference between \(\theta^*\) and \(\hat \theta({{\boldsymbol{x}}})\) is larger than \(0.05\)?” are not meaningful because, while \(\theta^*\) is unknown, both \(\theta^*\) and \(\hat\theta({{\boldsymbol{x}}})\) are deterministic after data is collected and the estimation task is performed.

The solution to these difficulties is to study the properties of the estimator not based on a specific realization \({{\boldsymbol{x}}}\) of the data but in general, over all possible data sets that could be produced and all the resulting estimated values. We can think of the thought experiment in which many, many, data sets are collected and the estimation task is performed based on each. The estimate itself is a random variable because each time we perform the estimation task, new data samples are obtained and these are random, following a certain distribution. In other words, instead of considering a single estimate \(\hat\theta({{\boldsymbol{x}}})\) for a specific realization \({{\boldsymbol{x}}}\), we study the estimator \(\hat\theta({{\boldsymbol{X}}})\), i.e., a random variable. Then it makes sense to ask “What is the probability that the difference between \(\theta^*\) and \(\hat \theta({{\boldsymbol{X}}})\) is larger than \(0.05\)?” since \(\hat\theta({{\boldsymbol{X}}})\) is a random variable with some distribution. It may be difficult to find the distribution of \(\hat\theta({{\boldsymbol{x}}})\) and it may depend on the unknown parameter \(\theta^*\) but at least the question is meaningful. In this section, we will see some of the evaluation criteria based on this view.

A note on notation. Typically, we use \(\theta\) as the generic parameter, with \(\theta^*\) denoting its true value, according to which \({{\boldsymbol{X}}}\) is distributed. For a given data \({{\boldsymbol{x}}}\), the estimate is shown by \(\hat\theta({{\boldsymbol{x}}})\) or \(\hat\theta\). So, \(\hat\theta\) denotes both the estimator, i.e., a function that produces the estimate given the data, and the estimate; the intent should be clear from the context. Finally, we may use \(\hat\Theta=\hat\theta({{\boldsymbol{X}}})\) to denote the estimate as a random variable.

2.3.1 Bias

Bias is the expected estimation error, \[\begin{equation} \mathop{\mathrm{Bias}}(\hat\theta) = \mathop{\mathrm{{\mathbb{E}}}}[\Hat{\theta}({{\boldsymbol{X}}}) - \theta^*]=\mathop{\mathrm{{\mathbb{E}}}}[\Hat{\theta}({{\boldsymbol{X}}})]-\theta^* \end{equation}\] As discussed, the expected value is taken over the randomness in \({{\boldsymbol{X}}}\). Bias of the estimator tells us whether in general the estimator over- or under-estimates the true value. If bias is equal to 0, then the estimator is called unbiased.

Example 2.7 Example 2.1 continued

Previously, we obtained the maximum likelihood estimate for the probability \(\theta\) of having light traffic. Let us find its bias. Again we collect data over 100 days and let \(X\) denote the number of days when there is light traffic. We know that \(\hat \theta_{{\rm mle}}\), as a function of data, is given by \[\begin{equation*} \hat \theta_{{\rm mle}}(X) = \frac{X}{100}, \end{equation*}\] Note that instead of using a specific value for the number of days with light traffic, such as 65, we use a random variable \(X\) representing this quantity. Dropping the dependence on \(X\) for simplicity, the expected value of \(\hat \theta_{{\rm mle}}\) is given by \[\begin{equation*} \mathop{\mathrm{{\mathbb{E}}}}{\left[ \hat \theta_{{\rm mle}}({{\boldsymbol{X}}}) \right]} = \frac{\mathop{\mathrm{{\mathbb{E}}}}{\left[ X \right]}}{100}. \end{equation*}\] Assuming \(\theta^*\) to be the true value, the number \(X\) of days when there is light traffic follows \(\mathop{\mathrm{\mathrm{Bin}}}(100,\theta^*)\), and so \(\mathop{\mathrm{{\mathbb{E}}}}[X]=100\theta^*\). It follows that \[\begin{equation*} \mathop{\mathrm{{\mathbb{E}}}}{\left[ \hat \theta_{{\rm mle}}({{\boldsymbol{X}}}) \right]} = \frac{100\theta^*}{100} = \theta^*. \end{equation*}\] Hence, the maximum likelihood estimate is an unbiased estimator.

Example 2.8 Three estimators of a mean

Given iid data \({{\boldsymbol{y}}}=(y_1,\dotsc,y_n), n\ge 3\), with mean \(\theta^*\), let us find the bias of each of the following estimators, \[\begin{align*} \hat\theta_1({{\boldsymbol{y}}}) &= \bar y = \frac{1}{n}\sum_{i=1}^n y_i,\\ \hat\theta_2({{\boldsymbol{y}}}) &= y_1,\\ \hat\theta_3({{\boldsymbol{y}}}) &= \frac{2y_2+y_3}{3}.\\ \end{align*}\] Let \(Y_i\) be the random variable corresponding to observation \(y_i\) and \(\bar Y = \sum_{i=1}^n Y\). We have \[\begin{align*} \mathop{\mathrm{{\mathbb{E}}}}\hat\theta_1({{\boldsymbol{Y}}}) &= \mathop{\mathrm{{\mathbb{E}}}}\bar Y = \frac{1}{n}\sum_{i=1}^n \mathop{\mathrm{{\mathbb{E}}}}Y_i= \frac{1}{n}\sum_{i=1}^n \theta^* = \theta^*,\\ \mathop{\mathrm{{\mathbb{E}}}}\hat\theta_2({{\boldsymbol{Y}}}) &= \mathop{\mathrm{{\mathbb{E}}}}Y_1 = \theta^*,\\ \mathop{\mathrm{{\mathbb{E}}}}\hat\theta_3({{\boldsymbol{Y}}}) &= \mathop{\mathrm{{\mathbb{E}}}}\left[\frac{2Y_2+Y_3}{3}\right]= \frac{2\mathop{\mathrm{{\mathbb{E}}}}Y_2+\mathop{\mathrm{{\mathbb{E}}}}Y_3}{3}=\theta^*. \end{align*}\] So all of these estimators are unbiased.

Example 2.9 Sample mean and variance

Given \(n\) samples \({{\boldsymbol{y}}}=(y_1,\dotsc,y_n)\) from a distribution with mean \(\mu^*\) and variance \({\left( \sigma^* \right)}^2\), are the estimators \[\hat \mu= \bar y = \frac{1}{n}\sum_{i=1}^ny_i, \quad \hat{\sigma}^2 = \frac{1}{n}\sum_{i=1}^n(y_i-\bar y)^2\] for the mean and variance, respectively, unbiased?

For \(\hat\mu\), we have \[\mathop{\mathrm{{\mathbb{E}}}}[\hat\mu({{\boldsymbol{Y}}})] = \mathop{\mathrm{{\mathbb{E}}}}[\bar Y] = \mathop{\mathrm{{\mathbb{E}}}}\left[\frac1n\sum_{i=1}^n Y_i\right] = \frac1n\sum_{i=1}^n\mathop{\mathrm{{\mathbb{E}}}}[Y_i] = \frac1n n \mathop{\mathrm{{\mathbb{E}}}}[Y_1] = \mu^*\] and so the estimator for the mean is unbiased. We can show (how?) that \[\mathop{\mathrm{{\mathbb{E}}}}\left[\hat\sigma^2({{\boldsymbol{Y}}})\right] =\frac{n-1}{n}{\left( \sigma^* \right)}^2\] and the bias of estimating \({\left( \sigma^* \right)}^2\) is \[\begin{equation*} \mathop{\mathrm{{\mathbb{E}}}}\left[\hat\sigma^2({{\boldsymbol{Y}}})\right]-{\left( \sigma^* \right)}^2 = -\frac{1}{n}{\left( \sigma^* \right)}^2. \end{equation*}\] Based on this, we can create an unbiased estimator for the variance as \[\hat{\sigma}_{u}^2({{\boldsymbol{y}}}) = \frac{1}{n-1}\sum_{i=1}^n(y_i-\bar y)^2.\]

Example 2.10 How many balls in the urn?

[1, Example 2.8.2] An urn has \(m^*\) balls, numbered \(1, 2, ..., m^*\). Suppose however that \(m^*\) is unknown to us. We pick one random ball from the urn and the number on the ball is \(y\). We estimate \(m^*\) using maximum likelihood. First, let \(Y\) be the random variable corresponding to observation \(y\), with distribution \(p_Y(y;m^*)\). We have \[\begin{equation*} p_Y(y;m) = \begin{cases} \frac{1}{m} & y \leq m,\\ 0& y > m. \end{cases} \end{equation*}\] and thus \[\begin{equation*} L(m) = \begin{cases} \frac{1}{m} & m\ge y,\\ 0& m<y. \end{cases} \end{equation*}\] Hence, \(L(m)\) is maximized by choosing \(m(y) =y\) and so \(\hat m_{{\rm mle}}=y\). To find the bias of \(\hat m_{{\rm mle}}\), \[\begin{align*} \mathop{\mathrm{{\mathbb{E}}}}[\hat m_{{\rm mle}}(Y)] &= \mathop{\mathrm{{\mathbb{E}}}}[Y] = \sum_{i=1}^{m^*} i\cdot \frac{1}{m^*}= \frac{m^*+1}{2},\\ \mathop{\mathrm{Bias}}(\hat m_{{\rm mle}}) &= \frac{m^*+1}{2} - m^* = -\frac{m^*-1}{2}, \end{align*}\] which means that the ML estimator tends to underestimates \(m^*\) by almost a factor of 2.

Example 2.11 Linear unbiased estimator

Can we design an unbiased estimator for Example 2.10? There are many options, but for simplicity we may choose an estimator that is linear in the data, in particular, one of the form \[\begin{equation*} \hat m_L(y) = a y + b. \end{equation*}\] We find \(a\) and \(b\) such that \(\hat m_L\) is unbiased. We have \[\mathop{\mathrm{{\mathbb{E}}}}[\hat m_L(Y)] = a\mathop{\mathrm{{\mathbb{E}}}}Y+b = a\frac{m^*+1}{2}+b.\] Setting this equal to \(m^*\) (equality should hold for any \(m^*\)) yields \(a=2\) and \(b=-1\), i.e., \[\hat m_L(y) = 2y-1.\]

Example 2.12 Survival of humanity (!)

The human species will eventually die out. We use two methods to estimate the total number of humans \(m\) who will ever live. Let humans be enumerated by birth order as \(h_1, h_2, ..., h_y, ..., h_m\), where \(h_1\) represents Adam, \(h_2\) represents Eve, \(h_y\) represents you, and \(h_n\) represents the last human to live. Assuming that your birth order \(y\) is random, the problem is similar to estimating the number of balls in an urn in Example 2.10.

Assuming that 100 billion humans have been born so far, we have \(\hat m_{{\rm mle}}=100\) billion and \(\hat m_{L}=200\) billion. The ML estimate predicts that the end is here. Further, assuming that there will be 140 million births each year, the unbiased estimator predicts the end of humanity to occur in around 700 years.

2.3.2 Mean squared error and variance

Example 2.13 Unbiased but useless

Consider an unbiased estimator \(\hat\theta\) and define \(\hat\theta'=\hat\theta + W\), where \(W\) is a zero-mean random variable with a large variance. Now, \(\hat\theta'\) is unbiased, similar to \(\hat\theta\), but it is not a good estimator (regardless of how good \(\hat\theta\) is). So clearly, being unbiased alone is not sufficient to ensure that an estimator is “good.”

For an estimator \(\hat\theta\), where the random variable describing data is denoted by \({{\boldsymbol{X}}}\), the mean squared error (MSE) is defined as \[\begin{equation*} \mathop{\mathrm{MSE}}(\hat \theta) = \mathop{\mathrm{{\mathbb{E}}}}\left[\left(\hat\theta({{\boldsymbol{X}}}) - \theta^*\right)^2\right]. \end{equation*}\] The smaller the MSE, the more accurate the estimator.

Let \(\hat\Theta = \hat \theta({{\boldsymbol{X}}})\). Note that \[\begin{align*} \mathop{\mathrm{MSE}}(\hat \theta) &=\mathop{\mathrm{{\mathbb{E}}}}\left[\left(\hat{\Theta}-\theta^*\right)^{2}\right] \\ & =\mathop{\mathrm{{\mathbb{E}}}}\left[\left(\left(\hat{\Theta}-\mathop{\mathrm{{\mathbb{E}}}}\hat{\Theta}\right)+\left(\mathop{\mathrm{{\mathbb{E}}}}\hat\Theta-\theta^*\right)\right)^{2}\right]\\ & =\mathop{\mathrm{{\mathbb{E}}}}\left[\left(\hat{\Theta}-\mathop{\mathrm{{\mathbb{E}}}}\hat{\Theta}\right)^2\right] +\left(\mathop{\mathrm{{\mathbb{E}}}}\hat\Theta-\theta^*\right)^{2}+2\mathop{\mathrm{{\mathbb{E}}}}\left[\left(\hat{\Theta}-\mathop{\mathrm{{\mathbb{E}}}}\hat{\Theta}\right)\right]\left(\mathop{\mathrm{{\mathbb{E}}}}\hat\Theta-\theta^*\right)\\ & =\mathop{\mathrm{{\mathbb{E}}}}\left[\left(\hat{\Theta}-\mathop{\mathrm{{\mathbb{E}}}}\hat{\Theta}\right)^2\right] +\left(\mathop{\mathrm{{\mathbb{E}}}}\hat\Theta-\theta^*\right)^{2}, \end{align*}\] where, the third equality uses the fact that \(\mathop{\mathrm{{\mathbb{E}}}}\hat\Theta-\theta^*\) is a deterministic constant and the fourth equality the fact that \(\mathop{\mathrm{{\mathbb{E}}}}\left[\left(\hat{\Theta}-\mathop{\mathrm{{\mathbb{E}}}}\hat{\Theta}\right)\right]=0\). Hence, \[\mathop{\mathrm{MSE}}(\hat\theta)=\mathop{\mathrm{Var}}(\hat\theta)+(\mathop{\mathrm{Bias}}(\hat\theta))^2.\]

For unbiased estimators, the variance is an important quantity since it is equal to the MSE.

Example 2.14 Example 2.1 re-revisit

We saw in Example 2.7 that the maximum likelihood estimate for the probability of traffic \(\theta^*\) is unbiased. Now, let us find its variance. Again, we write \(\hat \theta_{{\rm mle}}({{\boldsymbol{X}}}) = \frac{X}{100}\) and \[\begin{equation*} \mathop{\mathrm{Var}}(\hat\theta_{{\rm mle}}) = \frac{\mathop{\mathrm{Var}}(X)}{100^2} = \frac{\theta^*(1-\theta^*)}{100}, \end{equation*}\] where \(X\) is the number of days without traffic, which follows \(\mathop{\mathrm{\mathrm{Bin}}}(100,\theta^*)\) with variance \(100\theta^*(1-\theta^*)\). As we can see, the variance (hence, MSE) increases as the true value of \(\theta^*\) approaches \(1/2\), i.e., every data point contains more uncertainty. Furthermore, we can extend this result to the more general case where we collect data for \(n\) days. By the same argument, we get \[\begin{equation*} \mathop{\mathrm{MSE}}(\hat\theta_{{\rm mle}}) = \mathop{\mathrm{Var}}(\hat\theta_{{\rm mle}}) = \frac{\theta^*(1-\theta^*)}{n}. \end{equation*}\]

Example 2.15 The sample mean is unbiased

Consider data \({{\boldsymbol{y}}}={{\left( y_1,...,y_n \right)}}\), where the corresponding random variables \(Y_i\) are iid with distribution \(\mathop{\mathrm{{\mathcal{N}}}}{{\left( \mu,\sigma^2 \right)}}\). The ML estimator for the mean \(\mu\) is \(\hat\theta_{{\rm mle}}({{\boldsymbol{y}}})=\bar y=\frac{1}{n}\sum_{i=1}^n y_i\) is unbiased. We have \[\begin{equation*} \mathop{\mathrm{MSE}}(\hat\theta_{{\rm mle}})=\mathop{\mathrm{Var}}(\bar Y) = \frac{\sigma^2}{n}. \end{equation*}\]

Note that as \(n\) increases, the MSE decreases and the estimate becomes more accurate, as would be expected. This property is studied next.

Exercise 2.16 MSE of the mean estimators

For the estimators in Example 2.8, find the MSE, assuming the variance is \({\left( \sigma^* \right)}^2\).

Exercise 2.17 Bias-variance trade-off

Given iid data \({{\boldsymbol{y}}}=(y_1,\dotsc,y_n), n\ge 3\), with mean \(\theta^*\) and variance \(\sigma^2\), show that the MSE of \[\begin{align*} \hat\theta_1 &= ay_1,\\ \hat\theta_n &= a\bar y = \frac{a}{n}\sum_{i=1}^n y_i, \end{align*}\] for some constant \(a\in{\mathbb{R}}\) is given as \[\begin{align*} \mathop{\mathrm{MSE}}(\hat\theta_1) &= (a-1)^2{{\left( \theta^* \right)}}^2+a^2\sigma^2,\\ \mathop{\mathrm{MSE}}(\hat\theta_n) &= (a-1)^2{{\left( \theta^* \right)}}^2+a^2\sigma^2/n. \end{align*}\] What is a good value for \(a\)? Does anything other than \(a=1\) make sense? The components of the MSE are given in the plots below for \(\hat\theta_1\) and \(\hat\theta_n\) with \(n=10\), for \(\theta^*=0.5,\sigma^2 = 0.1\). A trade-off between the bias and variance is evident. Why is it not feasible to design an estimator by optimizing for \(a\)? What is the difference between estimation based on little data (\(\hat\theta_1\)) and a lot of data (\(\hat\theta_n,n=10\))?

2.3.3 Consistency

Consider an estimator \(\hat\theta_n({{\boldsymbol{x}}})\) based on \(n\) samples \({{\boldsymbol{x}}}=(x_1,\dotsc,x_n)\). Let \({{\boldsymbol{X}}}=(X_1,\dotsc,X_n)\) be the random variables that describe the \(n\) data samples and let \(\hat\Theta_n=\hat\theta_n({{\boldsymbol{X}}})\) be the random variable that corresponds to the estimate. The estimator \(\hat\theta_n\) is said to be consistent if \(\hat \Theta_n \to \theta^*\) as \(n\to \infty\). More precisely, for all \(\epsilon>0\), we need \[\lim_{n\to\infty}{\Pr}(|\hat\Theta_n-\theta^*|\ge\epsilon)=0.\] That is, \(\hat\Theta_n\) must converge in probability to \(\theta^*\) (the mode of convergence named in the probability review chapter): the estimator is accurate if the size of the data is large.

Example 2.18 Both estimators with more data

The ML and linear estimators described in Examples 2.10 and 2.11 are very different for a single data point. But how do they behave if we have a lot of data. First we need to define these for \(n\) data samples. Suppose that we take \(n\) samples from the urn with replacement, resulting in \({{\boldsymbol{y}}}=(y_1,y_2, \dotsc, y_n)\). Define \[\bar y=\frac{1}{n}\sum_{i=1}^{n} y_i.\] To extend the linear estimator to \(n\) data points, we can choose \[\hat m_{L,n} = 2\bar y -1.\] For the ML estimator, we have (why?) \[\begin{align*} \hat m_{{\rm mle},n}&=\max_i y_i. \end{align*}\] Both of these, although they look very different, are consistent and converge to \(m^*\) as \(n\to\infty\).

As \(n\to \infty\), by LLN, \(\bar Y\) converges to the mean of the distribution, i.e., \(\mathop{\mathrm{{\mathbb{E}}}}[ Y_1]=\frac{m^*+1}{2}\). Hence, \(\hat m_{L,n}\to 2\cdot\frac{m^*+1}{2}-1=m^*\).
For the ML estimator, as \(n\to\infty\), at some point, we will pick the ball numbered \(m^*\) and so we will eventually have \(\hat m_{{\rm mle}}=m^*\).

Given the two estimators, the bad news is that the estimators disagree significantly for small data. However, as the size of the sample data increases, the two estimators agree.

2.4 The Cramer-Rao lower bound (graduate only)

For an unbiased estimator, the MSE is equal to the variance, and thus the variance represents the accuracy of the estimator. This leads to the following question: For a given distribution of data, what is the smallest possible variance of an unbiased estimator?

The accuracy of estimating a parameter \(\theta\) depends on how strongly the distribution of the data \({{\boldsymbol{X}}}\) depends on \({\theta}\). If the dependence is strong, i.e., for values of \({\theta}\) other than the true value \(\theta^*\), the probability of the observed data falls sharply, then we may expect to find \({\theta^*}\) with accuracy. On the other hand, if the dependence is week, then it will be difficult to find \({\theta^*}\) with precision. These two cases are shown in Fig. 2.1.

**Figure 2.1:** The likelihood function on the left demonstrates strong dependence on \({\theta}\) compared to the one on the right.

Let the data be encoded as a vector \({{\boldsymbol{X}}}\), whose distribution is given by \(p\) with parameter \(\theta^*\). Assuming \({{\boldsymbol{X}}}={{\boldsymbol{x}}}\), the log-likelihood is \(p({{\boldsymbol{x}}};\theta)\). The sharpness of the log-likelihood \(\ell({\theta})\) at the true value \(\theta^*\) can be quantified as \[\begin{equation} \left.-\frac{\partial^2 \ell({\theta})}{\partial \theta^2}\right|_{\theta=\theta^*}=\left.-\frac{\partial^2 \ln{p({{\boldsymbol{x}}};\theta)} }{\partial \theta^2}\right|_{\theta=\theta^*}. \end{equation}\tag{2.5}\] Given the randomness of the data \({{\boldsymbol{X}}}\), the above quantity is random, \[\left.-\frac{\partial^2 \ln{p({{\boldsymbol{X}}};\theta)} }{\partial \theta^2}\right|_{\theta=\theta^*}\] So to average over the data, we define \[\begin{equation*} {\mathcal{I}}(\theta^*)=-\mathop{\mathrm{{\mathbb{E}}}}\left[\left.\frac{\partial^2 \ln{p({{\boldsymbol{X}}};\theta)} }{\partial \theta^2}\right|_{\theta=\theta^*} \right]=-\int \left.\frac{\partial^2 \ln{p({{\boldsymbol{x}}};\theta)} }{\partial \theta^2}\right|_{\theta=\theta^*}p({{\boldsymbol{x}}}; \theta^*)d{{\boldsymbol{x}}}, \end{equation*}\] which is called the Fisher information.

The following theorem provides a lower bound on the variance, which is referred to as the Cramer-Rao lower bound (CRLB).

Theorem 2.19 Cramer-Rao lower bound

Given that the log-likelihood \(\ell(\theta)\) satisfies certain regularity conditions , the variance of any unbiased estimator \(\hat \theta\) of \(\theta^*\) satisfies \[\begin{equation*} \mathop{\mathrm{Var}}(\hat\theta)\geq \frac{1}{{\mathcal{I}}(\theta^*)}. \end{equation*}\]

If an estimator achieves the CRLB, i.e., \(\mathop{\mathrm{Var}}(\hat\theta)=1/{\mathcal{I}}(\theta^*)\), then it is called efficient.

As a special case, consider when we have \(n\) iid data points, and denote the estimator based on this data as \(\hat\theta_n\). Denote the Fisher information based on \(n\) data points as \(I_n({\theta^*})\) and based on one data point as \(I_1({\theta^*})=I({\theta^*})\). Since the Fisher information is additive (Why? Hint: definition), we have \(I_n({\theta^*})= n I({\theta^*})\). Thus, the variance of an unbiased estimator \(\hat \theta_n\) based on \(n\) independent observations satisfies \[\begin{equation} \mathop{\mathrm{Var}}(\hat \theta_n)\geq \frac{1}{n{\mathcal{I}}(\theta^*)}. \end{equation}\]

Example 2.20 Fisher information for a Gaussian

In Example 2.4, where we estimated the mean \(\mu^*\) of a Gaussian distribution with known \({\sigma}^2\) based on \(n\) iid samples \(y_1,\dotsc,y_n\), the log-likelihood, ignoring constant terms, was given as \[\ell(\mu) \doteq -\sum_{i=1}^{n}\frac{(y_i-\mu)^2}{2\sigma^2}.\] And, \[\begin{equation} \frac{\partial \ell(\mu)}{\partial \mu}=\frac{1}{\sigma^2}\sum_{i=1}^{n}(y_i-\mu). \end{equation}\] Observe that \[\frac{\partial^2 \ell(\mu)}{ \partial \mu^2}=-\frac{n}{\sigma^2}\implies I(\mu^*) = -\mathop{\mathrm{{\mathbb{E}}}}\left[\frac{\partial^2 \ell(\mu^*)}{ \partial \mu^2}\right] =\frac{n}{\sigma^2}.\] Based on the CRLB, the variance of the estimator satisfies \[\mathop{\mathrm{Var}}(\hat \mu)\geq \frac{\sigma^2}{n}.\] The variance of the estimator is \(\mathop{\mathrm{Var}}(\hat \mu)=\frac{\sigma^2}{n}\). Hence, the ML estimator is efficient in this case.

2.5 Asymptotic normality of the MLE

As shown before, the maximum-likelihood estimator is not necessarily unbiased. However, if we have a large amount of data, under some regularity conditions, the ML estimator \(\hat{\Theta}_n\) based on \(n\) iid data points satisfies \[\sqrt{n}(\hat \Theta_n-\theta^*)\to {\mathcal{N}}(0, I^{-1}(\theta^*)).\] So for large data, \(\hat{\Theta}_{n}\) is nearly normally distributed with mean \(\theta^*\) (hence unbiased) and variance \(I^{-1}({\theta^*})/n\) (efficient).

While we stated the CRLB and the asymptotic normality of the MLE for scalar parameters, almost identical results also hold for a vector of parameters.

References

[1] B. Hajek, Random Processes for Engineers, 2014. http://hajek.ece.illinois.edu/Papers/randomprocJuly14.pdf

Note the slight abuse of notation: sometimes \(\theta\) is used as the generic parameter, e.g., as the argument of the likelihood function, and sometimes as the true value of the parameter. The distinction should be clear from the context↩︎