18.05 S22 Class 10: Maximum Likelihood Estimates

Maximum Likelihood Estimates Class 10, 18.05 1 Learning Goals 1. Be able to define the likelihood function for a parametric model given data. 2. Be able to compute the maximum likelihood estimate of unknown parameter(s). 2 Introduction Suppose we know we have data consisting of values 𝑥1 , … , 𝑥𝑛 drawn from an exponential distribution. The question remains: which exponential distribution?! We have casually referred to the exponential distribution or the binomial distribution or the normal distribution. In fact the exponential distribution exp(𝜆) is not a single distribution but rather a one-parameter family of distributions. Each value of 𝜆 defines a different distribution in the family, with pdf 𝑓𝜆 (𝑥) = 𝜆𝑒−𝜆𝑥 on [0, ∞). Similarly, a binomial distribution bin(𝑛, 𝑝) is determined by the two parameters 𝑛 and 𝑝, and a normal distribution 𝑁 (𝜇, 𝜎2 ) is determined by the two parameters 𝜇 and 𝜎2 (or equivalently, 𝜇 and 𝜎). Parameterized families of distributions are often called parametric distributions or parametric models. We are often faced with the situation of having random data which we know (or believe) is drawn from a parametric model, whose parameters we do not know. For example, in an election between two candidates, polling data constitutes draws from a Bernoulli(𝑝) distribution with unknown parameter 𝑝. In this case we would like to use the data to estimate the value of the parameter 𝑝, as the latter predicts the result of the election. Similarly, assuming gestational length follows a normal distribution, we would like to use the data of the gestational lengths from a random sample of pregnancies to draw inferences about the values of the parameters 𝜇 and 𝜎2 . Our focus so far has been on computing the probability of data arising from a parametric model with known parameters. Statistical inference flips this on its head: we will estimate the probability of parameters given a parametric model and observed data drawn from it. In the coming weeks we will see how parameter values are naturally viewed as hypotheses, so we are in fact estimating the probability of various hypotheses given the data. 3 Maximum Likelihood Estimates There are many methods for estimating unknown parameters from data. We will first consider the maximum likelihood estimate (MLE), which answers the question: For which parameter value does the observed data have the biggest probability? The MLE is an example of a point estimate because it gives a single value for the unknown parameter (later our estimates will involve intervals and probabilities). Two advantages of 1 18.05 Class 10, Maximum Likelihood Estimates , Spring 2022 2 the MLE are that it is often easy to compute and that it agrees with our intuition in simple examples. We will explain the MLE through a series of examples. Example 1. A coin is flipped 100 times. Given that there were 55 heads, find the maximum likelihood estimate for the probability 𝑝 of heads on a single toss. Before actually solving the problem, let’s establish some notation and terms. We can think of counting the number of heads in 100 tosses as an experiment. For a given value of 𝑝, the probability of getting 55 heads in this experiment is the binomial probability 𝑃 (55 heads) = ( 100 55 )𝑝 (1 − 𝑝)45 . 55 The probability of getting 55 heads depends on the value of 𝑝, so let’s include 𝑝 in by using the notation of conditional probability: 𝑃 (55 heads | 𝑝) = ( 100 55 )𝑝 (1 − 𝑝)45 . 55 You should read 𝑃 (55 heads | 𝑝) as: ‘the probability of 55 heads given 𝑝,’ or more precisely as ‘the probability of 55 heads given that the probability of heads on a single toss is 𝑝.’ Here are some standard terms we will use as we do statistics. • Experiment: Flip the coin 100 times and count the number of heads. • Data: The data is the result of the experiment. In this case it is ‘55 heads’. • Parameter(s) of interest: We are interested in the value of the unknown parameter 𝑝. • Likelihood, or likelihood function: this is 𝑃 (data | 𝑝). Note it is a function of both the data and the parameter 𝑝. In this case the likelihood is 𝑃 (55 heads | 𝑝) = ( 100 55 )𝑝 (1 − 𝑝)45 . 55 Notes: 1. The likelihood 𝑃 (data | 𝑝) changes as the parameter of interest 𝑝 changes. 2. Look carefully at the definition. One typical source of confusion is to mistake the likelihood 𝑃 (data | 𝑝) for 𝑃 (𝑝 | data). We know from our earlier work with Bayes’ theorem that 𝑃 (data | 𝑝) and 𝑃 (𝑝 | data) are usually very different. Definition: Given data the maximum likelihood estimate (MLE) for the parameter 𝑝 is the value of 𝑝 that maximizes the likelihood 𝑃 (data | 𝑝). That is, the MLE is the value of 𝑝 for which the data is most likely. Solution: For the problem at hand, we saw above that the likelihood 𝑃 (55 heads | 𝑝) = ( 100 55 )𝑝 (1 − 𝑝)45 . 55 18.05 Class 10, Maximum Likelihood Estimates , Spring 2022 3 We’ll use the notation 𝑝̂ for the MLE. We use calculus to find it by taking the derivative of the likelihood function and setting it to 0. 𝑑 100 𝑃 (data |𝑝) = ( )(55𝑝54 (1 − 𝑝)45 − 45𝑝55 (1 − 𝑝)44 ) = 0. 𝑑𝑝 55 Solving this for 𝑝 we get 55𝑝54 (1 − 𝑝)45 = 45𝑝55 (1 − 𝑝)44 55(1 − 𝑝) = 45𝑝 55 = 100𝑝 the MLE is 𝑝̂ = 0.55 Note: 1. The MLE for 𝑝 turned out to be exactly the fraction of heads we saw in our data. 2. The MLE is computed from the data. That is, it is a statistic. 3. Officially we need to check that this critical point is actually the maximum. We could use the second derivative test. Another way is to notice that we are interested only in 0 ≤ 𝑝 ≤ 1; that the probability is bigger than zero for 0 < 𝑝 < 1; and that the probability is equal to zero for 𝑝 = 0 and for 𝑝 = 1. From these facts it follows that the critical point must be the unique maximum. 3.1 Log likelihood If is often easier to work with the natural log of the likelihood function. For short this is simply called the log likelihood. Since ln(𝑥) is an increasing function, the maxima of the likelihood and log likelihood coincide. Example 2. Redo the previous example using log likelihood. 55 45 Solution: We had the likelihood 𝑃 (55 heads | 𝑝) = (100 55 )𝑝 (1 − 𝑝) . Therefore the log likelihood is 100 ln(𝑃 (55 heads | 𝑝) = ln (( )) + 55 ln(𝑝) + 45 ln(1 − 𝑝). 55 Maximizing likelihood is the same as maximizing log likelihood. We check that calculus gives us the same answer as before: 𝑑 𝑑 100 [ln (( )) + 55 ln(𝑝) + 45 ln(1 − 𝑝)] (log likelihood) = 𝑑𝑝 𝑑𝑝 55 55 45 = − =0 𝑝 1−𝑝 ⇒ 55(1 − 𝑝) = 45𝑝 ⇒ 𝑝̂ = 0.55 18.05 Class 10, Maximum Likelihood Estimates , Spring 2022 3.2 4 Maximum likelihood for continuous distributions For continuous distributions, we use the probability density function to define the likelihood. We show this in a few examples. In the next section we explain how this is analogous to what we did in the discrete case. Example 3. Light bulbs Suppose that the lifetime of Badger brand light bulbs is modeled by an exponential distribution with (unknown) parameter 𝜆. We test 5 bulbs and find they have lifetimes of 2, 3, 1, 3, and 4 years, respectively. What is the MLE for 𝜆? Solution: We need to be careful with our notation. With five different values it is best to use subscripts. Let 𝑋𝑖 be the lifetime of the 𝑖th bulb and let 𝑥𝑖 be the value 𝑋𝑖 takes. Then each 𝑋𝑖 has pdf 𝑓𝑋𝑖 (𝑥𝑖 ) = 𝜆e−𝜆𝑥𝑖 . We assume the lifetimes of the bulbs are independent, so the joint pdf is the product of the individual densities: 𝑓(𝑥1 , 𝑥2 , 𝑥3 , 𝑥4 , 𝑥5 | 𝜆) = (𝜆e−𝜆𝑥1 )(𝜆e−𝜆𝑥2 )(𝜆e−𝜆𝑥3 )(𝜆e−𝜆𝑥4 )(𝜆e−𝜆𝑥5 ) = 𝜆5 e−𝜆(𝑥1 +𝑥2 +𝑥3 +𝑥4 +𝑥5 ) . Note that we write this as a conditional density, since it depends on 𝜆. Viewing the data as fixed and 𝜆 as variable, this density is the likelihood function. Our data had values 𝑥1 = 2, 𝑥2 = 3, 𝑥3 = 1, 𝑥4 = 3, 𝑥5 = 4. So the likelihood and log likelihood functions with this data are 𝑓(2, 3, 1, 3, 4 | 𝜆) = 𝜆5 e−13𝜆 , ln(𝑓(2, 3, 1, 3, 4 | 𝜆) = 5 ln(𝜆) − 13𝜆 Finally we use calculus to find the MLE: 5 5 𝑑 . (log likelihood) = − 13 = 0 ⇒ 𝜆̂ = 𝑑𝜆 𝜆 13 Note: 1. In this example we used an uppercase letter for a random variable and the corresponding lowercase letter for the value it takes. This will be our usual practice. 2. The MLE for 𝜆 turned out to be the reciprocal of the sample mean 𝑥,̄ so 𝑋 ∼ exp(𝜆)̂ satisfies 𝐸[𝑋] = 𝑥.̄ The following example illustrates how we can use the method of maximum likelihood to estimate multiple parameters at once. Example 4. Normal distributions Suppose the data 𝑥1 , 𝑥2 , … , 𝑥𝑛 is drawn from a N(𝜇, 𝜎2 ) distribution, where 𝜇 and 𝜎 are unknown. Find the maximum likelihood estimate for the pair (𝜇, 𝜎2 ). Solution: Let’s be precise and phrase this in terms of random variables and densities. Let uppercase 𝑋1 , … , 𝑋𝑛 be i.i.d. N(𝜇, 𝜎2 ) random variables, and let lowercase 𝑥𝑖 be the value 𝑋𝑖 takes. The density for each 𝑋𝑖 is 𝑓𝑋𝑖 (𝑥𝑖 ) = √ (𝑥𝑖 −𝜇)2 1 e− 2𝜎2 . 2𝜋 𝜎 Since the 𝑋𝑖 are independent their joint pdf is the product of the individual pdf’s: 𝑛 𝑛 (𝑥𝑖 −𝜇)2 1 𝑓(𝑥1 , … , 𝑥𝑛 | 𝜇, 𝜎) = ( √ ) e− ∑𝑖=1 2𝜎2 . 2𝜋 𝜎 18.05 Class 10, Maximum Likelihood Estimates , Spring 2022 5 For the fixed data 𝑥1 , … , 𝑥𝑛 , the likelihood and log likelihood are 𝑛 𝑓(𝑥1 , … , 𝑥𝑛 |𝜇, 𝜎) = ( √ 𝑛 (𝑥𝑖 −𝜇)2 1 ) e− ∑𝑖=1 2𝜎2 , 2𝜋 𝜎 𝑛 √ (𝑥 − 𝜇)2 ln(𝑓(𝑥1 , … , 𝑥𝑛 |𝜇, 𝜎)) = −𝑛 ln( 2𝜋)−𝑛 ln(𝜎)−∑ 𝑖 2 . 2𝜎 𝑖=1 Since ln(𝑓(𝑥1 , … , 𝑥𝑛 |𝜇, 𝜎)) is a function of the two variables 𝜇, 𝜎 we use partial derivatives to find the MLE. The easy value to find is 𝜇:̂ 𝑛 𝑛 𝑛 ∑ 𝑥𝑖 𝜕𝑓(𝑥1 , … , 𝑥𝑛 |𝜇, 𝜎) (𝑥 − 𝜇) =∑ 𝑖 2 = 0 ⇒ ∑ 𝑥𝑖 = 𝑛𝜇 ⇒ 𝜇̂ = 𝑖=1 = 𝑥. 𝜕𝜇 𝜎 𝑛 𝑖=1 𝑖=1 To find 𝜎̂ we differentiate and solve for 𝜎: 𝑛 𝑛 ∑ (𝑥𝑖 − 𝜇)2 𝜕𝑓(𝑥1 , … , 𝑥𝑛 |𝜇, 𝜎) 𝑛 (𝑥 − 𝜇)2 =− +∑ 𝑖 3 = 0 ⇒ 𝜎̂ 2 = 𝑖=1 . 𝜕𝜎 𝜎 𝑖=1 𝜎 𝑛 We already know 𝜇̂ = 𝑥, so we use that as the value for 𝜇 in the formula for 𝜎.̂ We get the maximum likelihood estimates 𝜇̂ = the mean of the data =𝑥 𝑛 𝑛 1 1 𝜎̂ 2 = ∑ (𝑥𝑖 − 𝜇)̂ 2 = ∑ (𝑥𝑖 − 𝑥)2 = the unadjusted variance of the data. 𝑛 𝑛 𝑖=1 𝑖=1 𝑛 (Later we will learn that the sample variance is ∑𝑖=1 (𝑥𝑖 − 𝜇)̂ 2 .) 𝑛−1 Example 5. Uniform distributions Suppose our data 𝑥1 , … 𝑥𝑛 are independently drawn from a uniform distribution 𝑈 (𝑎, 𝑏). Find the MLE for 𝑎 and 𝑏. Solution: This example is different from the previous ones in that we won’t use calculus 1 to find the MLE. The density for 𝑈 (𝑎, 𝑏) is 𝑏−𝑎 on [𝑎, 𝑏]. Therefore our likelihood function is 𝑛 ( 1 ) if all 𝑥𝑖 are in the interval [𝑎, 𝑏] 𝑓(𝑥1 , … , 𝑥𝑛 | 𝑎, 𝑏) = { 𝑏−𝑎 0 otherwise. This is maximized by making 𝑏 − 𝑎 as small as possible. The only restriction is that the interval [𝑎, 𝑏] must include all the data. Thus the MLE for the pair (𝑎, 𝑏) is 𝑎̂ = min(𝑥1 , … , 𝑥𝑛 ) 𝑏̂ = max(𝑥1 , … , 𝑥𝑛 ). Example 6. Capture/recapture method The capture/recapture method is a way to estimate the size of a population in the wild. The method assumes that each animal in the population is equally likely to be captured by a trap. Suppose 10 animals are captured, tagged and released. A few months later, 20 animals are captured, examined, and released. 4 of these 20 are found to be tagged. Estimate the size of the wild population using the MLE for the probability that a wild animal is tagged. 18.05 Class 10, Maximum Likelihood Estimates , Spring 2022 6 Solution: Our unknown parameter 𝑛 is the number of animals in the wild. Our data is that 4 out of 20 recaptured animals were tagged (and that there are 10 tagged animals). The likelihood function is 𝑃 (data | 𝑛 animals) = 10 (𝑛−10 16 )( 4 ) 𝑛 (20 ) (The numerator is the number of ways to choose 16 animals from among the 𝑛−10 untagged ones times the number of was to choose 4 out of the 10 tagged animals. The denominator is the number of ways to choose 20 animals from the entire population of 𝑛.) We can use R to compute that the likelihood function is maximized when 𝑛 = 50. This should make some sense. It says our best estimate is that the fraction of all animals that are tagged is 10/50 which equals the fraction of recaptured animals which are tagged. Example 7. Hardy-Weinberg. Suppose that a particular gene occurs as one of two alleles (𝐴 and 𝑎), where allele 𝐴 has frequency 𝜃 in the population. That is, a random copy of the gene is 𝐴 with probability 𝜃 and 𝑎 with probability 1 − 𝜃. Since a diploid genotype consists of two genes, the probability of each genotype is given by: genotype probability AA 𝜃2 Aa 2𝜃(1 − 𝜃) aa (1 − 𝜃)2 Suppose we test a random sample of people and find that 𝑘1 are 𝐴𝐴, 𝑘2 are 𝐴𝑎, and 𝑘3 are 𝑎𝑎. Find the MLE of 𝜃. Solution: The likelihood function is given by 𝑃 (𝑘1 , 𝑘2 , 𝑘3 | 𝜃) = ( 𝑘 𝑘 + 𝑘3 𝑘1 + 𝑘 2 + 𝑘 3 )( 3 )𝜃2𝑘1 (2𝜃(1 − 𝜃))𝑘2 (1 − 𝜃)2𝑘3 . )( 2 𝑘3 𝑘2 𝑘1 So the log likelihood is given by constant + 2𝑘1 ln(𝜃) + 𝑘2 ln(𝜃) + 𝑘2 ln(1 − 𝜃) + 2𝑘3 ln(1 − 𝜃) We set the derivative equal to zero: 2𝑘1 + 𝑘2 𝑘2 + 2𝑘3 − =0 𝜃 1−𝜃 Solving for 𝜃, we find the MLE is 𝜃̂ = 2𝑘1 + 𝑘2 , 2𝑘1 + 2𝑘2 + 2𝑘3 which is simply the fraction of 𝐴 alleles among all the genes in the sampled population. 4 Why we use the density to find the MLE for continuous distributions The idea for the maximum likelihood estimate is to find the value of the parameter(s) for which the data has the highest probability. In this section we ’ll see that we’re doing this 18.05 Class 10, Maximum Likelihood Estimates , Spring 2022 7 is really what we are doing with the densities. We will do this by considering a smaller version of the light bulb example. Example 8. Suppose we have two light bulbs whose lifetimes follow an exponential(𝜆) distribution. Suppose also that we independently measure their lifetimes and get data 𝑥1 = 2 years and 𝑥2 = 3 years. Find the value of 𝜆 that maximizes the probability of this data. Solution: The main paradox to deal with is that for a continuous distribution the probability of a single value, say 𝑥1 = 2, is zero. We resolve this paradox by remembering that a single measurement really means a range of values, e.g. in this example we might check the light bulb once a day. So the data 𝑥1 = 2 years really means 𝑥1 is somewhere in a range of 1 day around 2 years. If the range is small we call it 𝑑𝑥1 . The probability that 𝑋1 is in the range is approximated by 𝑓𝑋1 (𝑥1 |𝜆) 𝑑𝑥1 . This is illustrated in the figure below. The data value 𝑥2 is treated in exactly the same way. density 𝑓𝑋1 (𝑥1 |𝜆) density 𝑓𝑋2 (𝑥2 |𝜆) 𝜆 𝜆 probability ≈ 𝑓𝑋1 (𝑥1 |𝜆) 𝑑𝑥1 probability ≈ 𝑓𝑋2 (𝑥2 |𝜆) 𝑑𝑥2 𝑑𝑥1 𝑑𝑥2 𝑥 𝑥1 𝑥 𝑥2 The usual relationship between density and probability for small ranges. Since the data is collected independently the joint probability is the product of the individual probabilities. Stated carefully 𝑃 (𝑋1 in range, 𝑋2 in range|𝜆) ≈ 𝑓𝑋1 (𝑥1 |𝜆) 𝑑𝑥1 ⋅ 𝑓𝑋2 (𝑥2 |𝜆) 𝑑𝑥2 Finally, using the values 𝑥1 = 2 and 𝑥2 = 3 and the formula for an exponential pdf we have 𝑃 (𝑋1 in range, 𝑋2 in range|𝜆) ≈ 𝜆e−2𝜆 𝑑𝑥1 ⋅ 𝜆e−3𝜆 𝑑𝑥2 = 𝜆2 e−5𝜆 𝑑𝑥1 𝑑𝑥2 . Now that we have a genuine probability we can look for the value of 𝜆 that maximizes it. Looking at the formula above we see that the factor 𝑑𝑥1 𝑑𝑥2 will play no role in finding the maximum. So for the MLE we drop it and simply call the density the likelihood: likelihood = 𝑓(𝑥1 , 𝑥2 |𝜆) = 𝜆2 e−5𝜆 . The value of 𝜆 that maximizes this is found just like in the example above. It is 𝜆̂ = 2/5. 5 Appendix: Properties of the MLE For the interested reader, we note several nice features of the MLE. These are quite technical and will not be on any exams. 18.05 Class 10, Maximum Likelihood Estimates , Spring 2022 8 The MLE behaves well under transformations. That is, if 𝑝̂ is the MLE for 𝑝 and 𝑔 is a one-to-one function, then 𝑔(𝑝)̂ is the MLE for 𝑔(𝑝). For example, if 𝜎̂ is the MLE for the standard deviation 𝜎 then (𝜎)̂ 2 is the MLE for the variance 𝜎2 . Furthermore, under some technical smoothness assumptions, the MLE is asymptotically unbiased and has asymptotically minimal variance. To explain these notions, note that the MLE is itself a random variable since the data is random and the MLE is computed from the data. Let 𝑥1 , 𝑥2 , … be an infinite sequence of samples from a distribution with parameter 𝑝. Let 𝑝𝑛̂ be the MLE for 𝑝 based on the data 𝑥1 , … , 𝑥𝑛 . Asymptotically unbiased means that as the amount of data grows, the mean of the MLE converges to 𝑝. In symbols: 𝐸[𝑝𝑛̂ ] → 𝑝 as 𝑛 → ∞. Of course, we would like the MLE to be close to 𝑝 with high probability, not just on average, so the smaller the variance of the MLE the better. Asymptotically minimal variance means that as the amount of data grows, the MLE has the minimal variance among all unbiased estimators of 𝑝. In symbols: for any unbiased estimator 𝑝𝑛̃ and 𝜖 > 0 we have that Var(𝑝𝑛̃ ) + 𝜖 > Var(𝑝𝑛̂ ) as 𝑛 → ∞.