Lecture 4 Likelihood and the Method of Moments

In the previous two lectures we have introduced the problem of estimation, but all the estimators we have considered so far were just guesses. In this lecture we introduce the likelihood function, the most important object in this module: the concepts of sufficiency (in lecture 8), Fisher information (in lecture 11), maximum likelihood (in lecture 5) and the Bayesian approach (in lecture 26) are all based on this function. We also discuss the method of moments, a simple method for constructing an estimator from a model. We will see the weaknesses of this approach which will lead to the introduction of the maximum likelihood estimator in the next lecture.

4.1 The Likelihood Function

From lecture 1 we know that a parametric model specifies a joint density for the data, for each $\theta\in\Theta$ . For an i.i.d. sample this is the product of the marginal densities. The idea of this section is to consider the converse problem, i.e. to consider, once the data $x_{1},\dots,x_{n}$ are fixed, how well a given value of $\theta$ explains the data.

Definition 4.1.

Let $x_{1},\dots,x_{n}$ be data from a parametric model $\bigl{(}f(x;\theta)\!\mathrel{\big{|}}\!\theta\in\Theta\bigr{)}$ . Then the likelihood function is the function $L\colon\Theta\to[0,\infty)$ , given by

L(\theta)=\prod_{i=1}^{n}f(x_{i};\theta),

and the log-likelihood is $\ell(\theta)=\log L(\theta)$ , given by

\ell(\theta)=\sum_{i=1}^{n}\log f(x_{i};\theta).

While the formula for $L$ is the same as the joint density of the sample, the interpretation of this function has changed: while the joint density is a function of the data with the parameter fixed, the likelihood is a function of the parameter with the data fixed. To emphasise that the data are used to compute the likelihood, we write $L(\theta;x_{1},\dots,x_{n})$ , and if we substitute the random variables $X_{i}$ for the numbers $x_{i}$ , the value of $L(\theta)$ is a random quantity. Both the numerical values and the random quantities can be used freely. For discrete data, $L(\theta)$ equals $\mathbb{P}_{\theta}(X_{1}=x_{1},\dots,X_{n}=x_{n})$ , i.e. the probability that the model with parameter $\theta$ generates the observed data.

Example 4.2.

Assume we have observed $7$ heads in $n=10$ tosses of a coin with unknown heads probability $\theta$ . Using the Bernoulli weights from table 1.1, we find the likelihood to be

L(\theta)=\prod_{i=1}^{10}\theta^{x_{i}}(1-\theta)^{1-x_{i}}=\theta^{7}(1-% \theta)^{3},\qquad\theta\in(0,1).

For the fair coin case we have $L(0.5)=0.5^{10}\approx 0.00098$ and $L(0.7)=0.7^{7}\cdot 0.3^{3}\approx 0.00222$ . Both values are very small, since any sequence of ten tosses is unlikely, but the ratio of these values is more informative: $\theta=0.7$ is $2.3$ times better at explaining the data than the fair coin model is, and this is the kind of information we are interested in.

Warning.

The likelihood is not a probability distribution over $\theta$ : the parameter is an unknown constant, statements like “the probability that $\theta=0.7$ is small” make no sense, and $L$ does not integrate to one over $\Theta$ . We will learn more about this in lecture 26, where we will see how the Bayesian approach can be used to turn the parameter into a random variable.

In practice we almost always use the log-likelihood: the likelihood is a product of $n$ factors, typically each less than one, so for realistic sample sizes it is astronomically small. In contrast, $\ell$ is a sum of moderate size and it is much easier to take derivatives of a sum than of a product. Since the logarithm is strictly increasing, the functions $L$ and $\ell$ order the parameter values in the same way, so no information is lost when we use the logarithm.

Example 4.3.

For an i.i.d. sample $x_{1},\dots,x_{n}$ from the exponential model with density $f(x;\theta)=\theta e^{-\theta x}$ for $x\geq 0$ , the likelihood is

L(\theta)=\prod_{i=1}^{n}\theta e^{-\theta x_{i}}=\theta^{n}\exp\Bigl{(}-% \theta\sum_{i=1}^{n}x_{i}\Bigr{)}.

Thus, the log-likelihood is

\ell(\theta)=n\log\theta-\theta\sum_{i=1}^{n}x_{i},\qquad\theta>0.

The product in this expression has been replaced by a smooth function of $\theta$ , which uses only two summaries of the data, $n$ and $\sum_{i}x_{i}$ . We will learn more about this approach to summarising data in lecture 8.

4.2 The Score Function

If the log-likelihood function is differentiable, then the derivative is sometimes called the score function. In this section we only use the score function to find the maximum of the likelihood. (We will learn the deeper meaning of the score function, as the carrier of information in the sample, in lecture 11.)

Definition 4.4.

Let $\ell(\theta)$ be the log-likelihood function and assume that it is differentiable at $\theta$ . Then the score function is given by

\ell^{\prime}(\theta)=\frac{\mathrm{d}}{\mathrm{d}\theta}\,\ell(\theta)=\sum_{% i=1}^{n}\frac{\partial}{\partial\theta}\log f(x_{i};\theta).

The score is the sum of $n$ terms, one for each observation, and if we assume that the data are random, these terms are i.i.d. This additive structure allows us to apply the law of large numbers and the central limit theorem to the likelihood function later in the module. Since the score is the derivative of the log-likelihood, a zero of the score is a good starting point to find the parameter value which gives the best explanation for the data.

Example 4.5.

Continuing from example 4.3, the score for the exponential sample is given by

\ell^{\prime}(\theta)=\frac{n}{\theta}-\sum_{i=1}^{n}x_{i},\qquad\theta>0.

The score is positive for $\theta<1/\bar{x}$ and negative for $\theta>1/\bar{x}$ , where $\bar{x}=(1/n)\sum_{i}x_{i}$ . Thus the log-likelihood increases up to $1/\bar{x}$ and then decreases, so $\theta=1/\bar{x}$ is the best explanation for the data. The idea of using the maximiser of the likelihood to estimate parameters is the basis of maximum likelihood estimation, as discussed in lecture 5.

4.3 The Method of Moments

Before we consider how to maximise the likelihood, we will discuss a much simpler method, the oldest general purpose method for constructing estimators. This method is based on the fact that the model often gives the moments of the observations as functions of the model parameters $\theta$ , and that the data can be used to compute the empirical moments. The idea of the method is then to find an estimator such that the two sets of moments coincide.

Let

\mu_{j}(\theta)=\mathbb{E}_{\theta}\bigl{(}X_{1}^{j}\bigr{)},\qquad m_{j}=% \frac{1}{n}\sum_{i=1}^{n}X_{i}^{j}

for $j=1,2,\dots$ . Then $\mu_{j}(\theta)$ is the $j$ -th population moment, a deterministic function of $\theta$ , and $m_{j}$ is the corresponding sample moment, a statistic. By the law of large numbers (theorem A.14 in appendix A), we have $m_{j}\xrightarrow{\ \mathrm{p}\ }\mu_{j}(\theta_{0})$ : for large sample size the sample moments are close to the population moments corresponding to the true parameter value.

Definition 4.6.

Let $\theta=(\theta_{1},\dots,\theta_{k})$ be a parameter with $k$ components. Then the method of moments estimator $\hat{\theta}$ is found by solving the system of $k$ equations

\mu_{j}(\hat{\theta})=m_{j},\qquad j=1,\dots,k,

for $\hat{\theta}\in\Theta$ , if this system of equations has a unique solution.

The recipe uses as many moment equations as there are unknown parameter components, and the required moments can usually be found in tables like tables A.1 and A.2. We illustrate the method for three of our running examples.

Example 4.7.

For the exponential model with rate $\theta$ , we have $\mu_{1}(\theta)=1/\theta$ and the only moment equation is $1/\hat{\theta}=\bar{X}$ . Thus, the method of moments estimator is given by

\hat{\theta}=\frac{1}{\bar{X}}.

Coincidentally, this is the same as the zero of the score from example 4.5.

Example 4.8.

For the normal model with $\theta=(\mu,\sigma^{2})$ the parameter has two components. Thus we need to match the first two moments, $\mu_{1}(\theta)=\mu$ and $\mu_{2}(\theta)=\sigma^{2}+\mu^{2}$ . We can solve these two moment equations (see exercise 4.3) to find

\hat{\mu}=\bar{X}\qquad\text{and}\qquad\hat{\sigma}^{2}=\frac{1}{n}\sum_{i=1}^% {n}(X_{i}-\bar{X})^{2}.

Note the divisor $n$ : since we are matching moments without considering the bias, the method gives the biased plug-in variance estimator from lecture 2, instead of the unbiased sample variance $S^{2}$ . In this example, from the mean squared error comparison in lecture 2, this is not a problem.

Example 4.9.

Let $X_{1},\dots,X_{n}\sim\text{Gamma}(\alpha,\beta)$ be i.i.d. where the shape $\alpha$ and the rate $\beta$ are unknown. Then $\theta=(\alpha,\beta)$ . From table A.2 we know $\mathbb{E}_{\theta}(X_{1})=\alpha/\beta$ and $\mathop{\mathrm{Var}}\nolimits_{\theta}(X_{1})=\alpha/\beta^{2}$ . Instead of matching $\mu_{1}$ and $\mu_{2}$ we can equivalently match the mean and the variance: if we write $V=m_{2}-m_{1}^{2}=(1/n)\sum_{i}(X_{i}-\bar{X})^{2}$ for the plug-in variance, we get the equations

\frac{\hat{\alpha}}{\hat{\beta}}=\bar{X}\qquad\text{and}\qquad\frac{\hat{% \alpha}}{\hat{\beta}^{2}}=V.

Dividing the first equation by the second we find $\hat{\beta}=\bar{X}/V$ and thus

\hat{\alpha}=\frac{\bar{X}^{2}}{V}\qquad\text{and}\qquad\hat{\beta}=\frac{\bar% {X}}{V}.

This example shows the method in a good light: for the gamma distribution the maximum likelihood equations from lecture 5 have no closed-form solution, whereas the method of moments finds explicit formulas in two lines.

We now weigh the merits of the method. On the credit side, it is quick, often gives explicit formulas where the likelihood does not, and is consistent under weak conditions: the sample moments converge by the law of large numbers, and if the solution of the moment equations depends continuously on the moments, the continuous mapping theorem (theorem A.17) carries the convergence over to the estimator.

The method has three main weaknesses. First, the estimates can be contradicted by the data: for the uniform model $X_{1},\dots,X_{n}\sim\text{Uniform}(0,\theta)$ we have $\mu_{1}(\theta)=\theta/2$ and thus $\hat{\theta}=2\bar{X}$ , but for the dataset $x_{1}=0.1$ , $x_{2}=0.2$ , $x_{3}=0.9$ this gives $\hat{\theta}=2\cdot 0.4=0.8$ , an estimate where the observation $x_{3}=0.9$ could never have occurred. This problem is immediately detected by the likelihood, since $L(0.8)=0$ , but the moment equations do not detect this problem. Secondly, the choice of moments is arbitrary: if we had used the second moment instead of the first, we would in general have obtained a different estimator. Finally, by taking only a few moments, much of the information in the data is lost. A method which makes use of all the information in the data should be more successful. This is the aim of the next few lectures.

Remark.

In some cases, both approaches coincide: for the exponential model the moment estimator $1/\bar{X}$ is also the maximiser of the likelihood (from example 4.5), the same is true for the Poisson distribution (exercise 4.1) and for the normal distribution the two methods coincide (lecture 5). In contrast, for the uniform model considered above, the maximum likelihood estimator is $\max_{i}X_{i}$ instead of $2\bar{X}$ and has a very different mean squared error (as shown in exercises 2.2 and 4.4).

4.4 Identifiability

Both recipes of this lecture, and the whole estimation programme, assume that the parameter can be recovered from the distribution of the data: if two parameter values lead to the same distribution, no data can distinguish between them and no good estimator for $\theta$ can exist. This property of a model deserves a name.

Definition 4.10.

A parametric model $\bigl{(}\mathbb{P}_{\theta}\!\mathrel{\big{|}}\!\theta\in\Theta\bigr{)}$ is identifiable, if different parameter values lead to different distributions, i.e. if we have $\mathbb{P}_{\theta_{1}}\neq\mathbb{P}_{\theta_{2}}$ for all $\theta_{1},\theta_{2}\in\Theta$ with $\theta_{1}\neq\theta_{2}$ .

All the families of models listed in tables 1.1 and 1.2 are identifiable. Non-identifiability usually occurs in models which are over-parametrised, as the following example shows.

Example 4.11.

Assume that we want to model a measurement which is subject to two additive effects. We can assume that $X_{1},\dots,X_{n}\sim N(\alpha+\beta,1)$ are i.i.d. with parameter $\theta=(\alpha,\beta)\in\mathbb{R}^{2}$ . Then the distribution of the data depends on $\theta$ only via the sum $\alpha+\beta$ . Thus, the parameters $(1,2)$ and $(0,3)$ give rise to the same distribution and the model is not identifiable. Both recipes from this lecture fail in this case: the likelihood is constant on the lines $\alpha+\beta=\text{const}$ and the moment equations only contain the unknowns in the combination $\alpha+\beta$ . To solve the problem we can either reparametrise the model using $\mu=\alpha+\beta$ or we can restrict the parameter space, for example by setting $\beta=0$ .

Identifiability is a question about the parametrisation, not about the data or the sample size. A more subtle version of the problem is the case of the mixture models from lecture 19, where the models are only identifiable up to relabelling the mixture components. For the rest of the module, we will assume that the models are identifiable, unless we explicitly state otherwise.

Summary.

•

The likelihood $L(\theta)=\prod_{i}f(x_{i};\theta)$ is the joint density of the observed data, as a function of $\theta$ with the data fixed, and is not a probability distribution over $\theta$ .
•

The score $\ell^{\prime}(\theta)$ , the derivative of the log-likelihood function $\ell(\theta)=\sum_{i}\log f(x_{i};\theta)$ , finds the parameter value which best explains the data. The score returns in lecture 11.
•

The method of moments solves the equations $\mu_{j}(\hat{\theta})=m_{j}$ for $j=1,\dots,k$ . The method is fast, explicit and consistent, but it can return impossible estimates and does not make good use of the data. The maximum likelihood method from lecture 5 is better.
•

A model is identifiable, if different parameter values lead to different distributions. If the model is not identifiable, no estimator can work and we need to change the parametrisation.

Exercise 4.1.

Let $X_{1},\dots,X_{n}$ be an i.i.d. sample from a Poisson distribution with parameter $\theta>0$ and let $x_{1},\dots,x_{n}$ be the observed data.

1.

Write down the likelihood $L(\theta)$ and show that the log-likelihood is $\ell(\theta)=\log\theta\sum_{i}x_{i}-n\theta-\sum_{i}\log(x_{i}!)$ .
2.

Compute the score $\ell^{\prime}(\theta)$ and find the value of $\theta$ where the score equals zero.
3.

Compute the method of moments estimator for $\theta$ and compare your result to your answer in the previous part.

Exercise 4.2.

Let $X_{1},\dots,X_{n}$ be an i.i.d. sample from the geometric distribution with parameter $p\in(0,1)$ , i.e. from the distribution with $f(x;p)=(1-p)^{x-1}p$ for $x\in\{1,2,\dots\}$ and $\mathbb{E}_{p}(X_{1})=1/p$ . (There is some disagreement among authors whether the geometric distribution counts the number of trials until the first success, inclusive of this trial, starting at $1$ as in the present exercise, or the number of failures before the first success, starting at $0$ . In this module we always count trials.)

1.

Using the method of moments, determine an estimator $\hat{p}$ for $p$ .
2.

Show that the value of $\hat{p}$ always lies in the interval $(0,1]$ , i.e. that the estimate can never be contradicted by the data (this is in contrast to the example of the uniform distribution in the text).

Exercise 4.3.

Let $X_{1},\dots,X_{n}\sim N(\mu,\sigma^{2})$ i.i.d. with parameter $\theta=(\mu,\sigma^{2})$ , as in example 4.8.

1.

Using the relation $\mathop{\mathrm{Var}}\nolimits(X_{1})=\mathbb{E}(X_{1}^{2})-\mathbb{E}(X_{1})^% {2}$ , show that the first two population moments are $\mu_{1}(\theta)=\mu$ and $\mu_{2}(\theta)=\sigma^{2}+\mu^{2}$ .
2.

Write down the two moment equations and solve them to verify the estimators $\hat{\mu}=\bar{X}$ and $\hat{\sigma}^{2}=(1/n)\sum_{i}(X_{i}-\bar{X})^{2}$ from example 4.8.

Exercise 4.4.

Let $X_{1},\dots,X_{n}\sim\text{Uniform}(0,\theta)$ i.i.d. and let $\hat{\theta}=2\bar{X}$ be the method of moments estimator as given in the text.

1.

Show that $\hat{\theta}$ is unbiased and that $\mathop{\mathrm{MSE}}\nolimits(\hat{\theta})=\theta^{2}/(3n)$ .
2.

In exercise 2.2 we have seen that $\mathop{\mathrm{MSE}}\nolimits(\max_{i}X_{i})=2\theta^{2}/\bigl{(}(n+1)(n+2)% \bigr{)}$ for the biased estimator $\max_{i}X_{i}$ . Show that $\max_{i}X_{i}$ has strictly smaller mean squared error for all $n\geq 3$ and comment on your result.

Exercise 4.5.

Consider the model $X_{1},\dots,X_{n}\sim N(0,\theta^{2})$ i.i.d. for the parameter $\theta\in\Theta=\mathbb{R}\setminus\{0\}$ .

1.

Show that the model is not identifiable.
2.

Find a smaller parameter space $\Theta^{\prime}\subseteq\Theta$ , where the model is identifiable.
3.

On your space $\Theta^{\prime}$ , determine the method of moments estimator for $\theta$ , using the second moment.