Lecture Notes on Likelihood Methods

9.07 Introduction to Statistics for Brain and Cognitive Sciences Emery N. Brown Lecture 9: Likelihood Methods I. Objectives 1. Review the concept of the joint distribution of random variables. 2. Understand the concept of the likelihood function. 3. Understand the concept of maximum likelihood estimation. 4. Understand the key properties of maximum likelihood estimation. 5. Understand the concept of Fisher information. 6. Understand how to compute the uncertainty in maximum likelihood estimates. II. Joint Probability Density or Distribution of a Set of Random Variables We recall that two events E1 and E2 are independent if Pr( E1 I E2 ) = Pr( E1 ) Pr( E2 ) (9.1) Intuitively, this statement says that knowledge about E1 gives no information about E2 . In general a set of events E1 , K , En is independent if n Pr( E1 I E2 K I En ) = ∏ Pr(E ) i (9.2) i=1 It follows from Lecture 1 and Lecture 5 that if x = x1 , K , xn is a sample of independent, identically distributed observations from a pmf or a pdf f ( xi | θ ) then the joint probability density of x = x1 , K , xn is f (x | θ ) = n ∏ f ( x | θ ). i (9.3) i=1 Equation 9.3 follows directly from the definition of a pmf for a discrete random variable. Using a cdf of a continuous random variable which is differentiable, it is also easy to show. Equation 9.3 is essential for our definition of the likelihood function of a sample of independent observations. III. The Likelihood Function A. Definition Definition 9.1 (Likelihood Function). Given x = x1 , K , xn is an independent, identically distributed sample from a pdf or pmf f ( xi | θ ). Let f ( x | θ ) denote the joint pdf or pmf of the sample x = ( x1 , K , xn ) defined in Eq. 9.3. Given X = x is observed, the function of θ defined as page 2: 9.07 Lecture 9: Likelihood Methods L(θ ) = f ( x | θ ) = n ∏ f (x k |θ ) (9.4) k =1 is called the likelihood function for the parameter θ . The likelihood function or likelihood provides an objective means of assessing the “information” in a sample of data about the model parameter θ . We view the data as fixed and now study L(θ ) = f ( x | θ ) as a function of θ . For a given model it summarizes all the information in the sample about θ . Use of the likelihood will allow us to overcome some of the ambiguities we faced with the method-of-moments. B. Examples of Likelihoods Example 2.1 (continued). Binomial Likelihood Function. In our learning experiment, we observed 22 correct responses from the animal in 40 trials. We derived the pmf of the 22 observations in 40 trials as the binomial distribution ⎛ 40 ⎞ f (22 | p) = ⎜ ⎟ p 22 (1 − p)18 . ⎝ 22 ⎠ (9.5) (Notice that this distribution does not have exactly the same form as Eq. 9.3, because we allowed for all the possible ways the sample could be realized. It turns out that in our likelihood analysis, this distinction will not matter because the binomial coefficient does not alter the likelihood). Since the unknown parameter is p we have the likelihood of p is ⎛ 40 ⎞ L( p) = f (22 | p) = ⎜ ⎟ p 22 (1 − p)18 . ⎝ 22 ⎠ (9.6) In general, for the binomial pmf with n trials and probability of a correct response p, the likelihood ⎛n⎞ L( p) = f ( x | p) = ⎜ ⎟ p x (1− p) n − x . ⎝ x⎠ (9.7) page 3: 9.07 Lecture 9: Likelihood Methods Figure 9A. Binomial Likelihood Function L(p) for n = 40 and k = 5, 10, 22, 30, 35. Plotting these elementary likelihood functions is useful to gain insight into what type of information they provide. We can predict the shape of the binomial likelihood function based on discussions of the beta distribution in Lecture 3 (Eq. 3.30, Figure 3J). Example 3.2 Quantile Release Hypothesis (Poisson Likelihood). We used a Poisson pmf to model the quantile release of acetylcholine at the synapse. If we record the number of quanta released in n non-overlapping time intervals of a specified length, then the likelihood function for the Poisson parameter λ , which is the expected number of quanta released, is L (λ ) = f ( x | λ ) = n ∏ λ xi e−λ i=1 xi ! n ∑ xi = (9.8) λ i=1 e−nλ n ∏x ! i i=1 We see in this example, a very important feature of the likelihood. Namely, that it has reduced n or summarized the data from all n observations to simply Sn = ∑x , i the sum of the n i=1 observations. For a large class of probability models, this is a typical feature. page 4: 9.07 Lecture 9: Likelihood Methods Figure 9B Poisson Likelihood Functions L(λ ) . The shape of the Poisson likelihood could be predicted based on results in Lecture 3. Note that the essential features of this likelihood are described by L(λ )∝ λ Sn e−nλ (9.9) which is proportional to a gamma probability density (Eq. 3.26) with α = (Sn −1) and β = n (Figure 3I). Example 3.2 (continued) Magnetoencephalogram Noise Data (Gaussian Likelihood Function). If we view the MEG data as a random sample from a Gaussian probability model then in this problem there are two unknown parameters: µ and σ 2 . page 5: 9.07 Lecture 9: Likelihood Methods L( µ , σ 2 ) = f ( x | µ , σ 2 ) n = ∏ f ( x | µ ,σ i 2 ) i=1 1 2 ⎪⎧ (xi − µ ) ⎪⎫ ⎛ 1 ⎞2 = exp − ⎨ ⎬ ⎜ ⎟ 2 2σ 2 ⎭⎪ ⎪⎩ i=1 ⎝ 2πσ ⎠ ⎧ n ⎫ (xi − µ ) 2 ⎪ ⎪ n ⎪ ⎪ ⎛ 1 ⎞2 =⎜ exp ⎨− i=1 ⎬ 2 ⎟ 2 2σ ⎝ 2πσ ⎠ ⎪ ⎪ ⎪ ⎪ ⎩ ⎭ n ∏ (9.10) ∑ If we expand the sum in the exponent, we see that the data have been compressed into n two distinct terms: ∑ i=1 n xi and ∑x i 2 . That is, we can rewrite Eq. 9.10 as i=1 ⎪⎧ 1 L(µ , σ 2 ) = exp ⎨− 2 ⎩⎪ 2σ n ∑x i i=1 2 + µ σ2 n n⎡ ∑ x − 2 ⎢⎣log(2πσ i i=1 2 )+ µ 2 ⎤ ⎪⎫ ⎥⎬ σ 2 ⎦ ⎪⎭ A crucial feature is that the coefficients on these two statistics involve the unknown parameters. This data reduction is an important consequence of the likelihood analysis. These statistics are called sufficient statistics in that it is possible to show that they contain all the information in the data about the parameters µ and σ 2 . We state this concept formally. Definition 9.2 A statistic T ( X ) is sufficient for a parameter θ if the probability density of the data conditional on the statistic is independent of θ . That is, if f ( X1, ..., X n | T ( X ),θ ) = h( X ) , then T ( X ) is sufficient. If there are two parameters, then the smallest number of sufficient statistics is two. Similarly, in the case of the Poisson likelihood above, the number of sufficient statistics was one consistent with the fact that there is only one parameter to estimate. There are very precise theorems that tell us how to identify the sufficient statistics for parameters in a probability model. We will not discuss them here. Details can be found in Pawitan (2001), DeGroot and Schervish (2002) and in Rice (2007). Likelihood analyses in general are based on the sufficient statistics. When method-ofmoments and likelihood analyses differ, it is usually because the sample moments are not the sufficient statistics for the parameters to be estimated. We will illustrate this point below when we define maximum likelihood estimation. C. Maximum Likelihood Estimation Ideally, we would like to analyze the likelihood function for every problem. This gets to be especially challenging when the number of parameters is large. Therefore, we would like to derive a summary that gives us a sense about the shape of the likelihood. This is what the maximum likelihood estimate and the Fisher information provide. We describe them in turn. page 6: 9.07 Lecture 9: Likelihood Methods Definition 9.3. Maximum Likelihood Estimation. For each sample point x = ( x1 , ..., xn ) let θˆ( x) be a parameter value of which L(θ ) = L(θ | x) attains a maximum as a function of θ for fixed x. θˆ( x) is a maximum likelihood (ML) estimator (MLE) of the parameter θ . In particular, if L(θ ) is differentiable, we can consider on ∂L(θ ) = 0 and check the conditions ∂θ ∂ 2 L(θ ) to be sure that the estimate defines a maximum. Remember that this would mean ∂θ 2 that for a one-dimensional problem, verifying that the second derivative was negative and for the d − dimensional problem verifying the condition that the determinant of the Hessian is negative definite. In likelihood analyses, it is usually easier to work with log L(θ ) instead of L(θ ). ∂ log L(θ ) is called the score function. Hence, computing the MLE can often be ∂θ formulated as finding θˆ which solves the score equation The function ∂ log L(θ ) = 0. ∂θ (9.11) Two challenging problems for finding maximum likelihood estimates of parameters in complex high dimensional models is finding the solution to Eq. 9.11 numerically and establishing that this solution achieves a global maximum. Example 9.1. Suppose that x1 , K , xn are independent, identically distributed observations from a gamma distribution with parameters α and β . If α is known then we have the likelihood function is L(α , β ) = n ∏ k =1 β α α −1 − β xk xk e Γ(α ) (9.12) β nα = exp{(α −1) log xk − β xk } Γ(α ) n k =1 k =1 n n ∑ ∑ and the log likelihood function is log L(α , β ) = −nlog Γ(α ) + nα log β + (α −1) n n ∑log( x ) − β ∑ x k k =1 k (9.13) k =1 Differentiating with respect to β n ∂ log L(α , β ) nα = − xk ∂β β k =1 ∑ and solving yields (9.14) page 7: 9.07 Lecture 9: Likelihood Methods 0= α −x β α β̂ = x (9.15) (9.16) . Notice also that if α is known, then the second derivative of the log likelihood evaluated at β̂ is ∂ 2 log L(α , β ) ∂β 2 =− β̂ nα nx =− 2 ˆ α β This is clearly negative since α is positive and x is positive. Hence, the log likelihood (and the likelihood) is maximized at β̂ . If α = 1 we have β̂ = x −1 is the maximum likelihood estimate for an exponential model. If α is unknown, then there is no closed form solution for either α or β . The estimates must then be found numerically. To see this we differentiate Eq. 9.13 with respect to α to obtain the log likelihood for α n ∂ log L(α , β ) Γ′(α ) = −n + n log β + log( xk ) ∂α Γ(α ) k =1 ∑ (9.17) Substituting for β from Eq. 9.16 we obtain n ∂ log L(α ) Γ′(α ) α = −n + n log + log( xk ) ∂α Γ(α ) x k =1 ∑ (9.18) Setting the left hand side of Eq. 9.18 equal to zero and solving for the MLE of α requires a numerical procedure such as Newton’s method. Plausible starting values for computing the ML estimates can be obtained from the method-of-moments estimates. We can see in this case that the method-of-moments estimates and the maximum likelihood estimates are different. The distinction between the maximum likelihood and the method-of-moments estimates comes about because for the gamma distribution, the sufficient statistics for α and β are n ∑ log( x ) k k =1 n and ∑ x . This is suggested by how the likelihood summarizes the data sample as indicated by k k =1 the terms in the log likelihood function in Eq. 9.13. Given the definition of the sufficient statistics, this shows that the simple method-of-moments estimates are not “efficient” (i.e., they must have some loss of information). We will make this concept more precise below. Example 3.2 (continued) Let x1 , K , xn be a random sample of MEG background noise believed to obey a Gaussian probability density with unknown µ and σ 2 . It is easy to show that x = n−1 ∑ n x k =1 k is the maximum likelihood (ML) estimate of µ and σˆ 2 = n−1 ∑ n k =1 (xk − x )2 is the maximum likelihood estimate of σ 2 . This is straightforward to show by differentiating the page 8: 9.07 Lecture 9: Likelihood Methods Gaussian log likelihood, equating the gradient to zero and solving for µ and σ 2 . Here the method-of-moments and the maximum likelihood estimates are the same. In this case, the sufficient statistics are the first two sample moments. Example 9.2. Inverse Gaussian Distribution. If x1 , K , xn is a random sample from an inverse Gaussian distribution with parameters µ and λ , then the likelihood is 1 ⎛ λ ⎞2 ⎧⎪ λ (xk − µ ) 2 ⎫⎪ exp L( µ , λ ) = ⎜⎜ ⎟ ⎨− ⎬ 3⎟ 2xk µ 2 ⎭⎪ ⎩⎪ k =1 ⎝ 2π xk ⎠ n ∏ 1 ⎧⎪ n λ (xk − µ ) 2 ⎫⎪ ⎛ λ ⎞2 = exp ⎜⎜ ⎟ ⎨− ⎬ 3⎟ 2 ⎪⎭ ⎩⎪ k =1 2xk µ k =1 ⎝ 2π xk ⎠ n ∑ ∏ ⎪⎧ 1 ⎡ λ = exp ⎨− ⎢ 2 ⎪⎩ 2 ⎢⎣ µ n ∑x k k =1 +λ n ∑x −1 k k =1 (9.19) ⎡ ⎛ λ ⎞ λ⎤ − n ⎢log ⎜ ⎟+2 ⎥+ 2 π µ⎦ ⎠ ⎣ ⎝ n ⎤ ⎪⎫ ∑ log x ⎥⎥⎦ ⎬⎪. k =1 3 k ⎭ Recall that for this probability model the mean is µ and the variance is µ 3 / λ. What are the maximum likelihood estimates of µ and λ ? Are they the same as the method-of-moments estimate? The maximum likelihood estimates are µ̂ML = n−1 n ∑x (9.20) i i=1 λ̂ML −1 = n−1 n 1 i=1 i ∑ ( x − µ̂ 1 ) (9.21) ML whereas the method-of-moments estimates are n ∑x i (9.22) λˆMM = x 3 / σ̂ 2 . (9.23) µ̂MM = n−1 i=1 Can you venture a guess as to what are the sufficient statistics for µ and λ ? (Remember the maximum likelihood estimates are always functions of the sufficient statistics. If you can derive the last line in Eq. 9.19, then the extra credit problem on Homework Assignment 7 is straightforward). Invariance Property of the Maximum Likelihood Estimator. Suppose that a distribution is indexed by a parameter θ . Let T (θ ) be a function of θ , then if θˆML is the maximum likelihood estimator of θ then T(θˆML ) is the maximum likelihood estimate of T (θ ). For example, if θ = λ , the scale parameter of the inverse Gaussian distribution the MLE λˆML is given in Eq. 9.21. If we are interested in T(λ ) = cos(λ ) then the MLE of T (λ ) is T (λˆ ). This is one of the most important properties of maximum likelihood estimates. ML page 9: 9.07 Lecture 9: Likelihood Methods D. Fisher Information (Pawitan, 2001). The MLE provides a point (or single number summary) estimate. “In general a single number is not sufficient to represent a likelihood function. If the log likelihood is well approximated by a quadratic function, then we need at least two quantities to represent it: the location of its maximum the (MLE) and the curvature at the maximum.” The curvature will be the Fisher information. When this quadratic approximation is valid we call the likelihood function regular. Fortunately, when sample sizes are sufficiently large, the likelihood function generally does become regular. To restate this critical requirement, regular problems are those in which we can approximate the log likelihood function around the MLE by a quadratic function. To define the Fisher information we first recall how we compute the Taylor series expansion of a function. The Taylor series expansion of a function f ( x) about a point x0 is defined as f ( x) = f ( x0 ) + f '( x0 )( x − x0 ) + f ''( x0 )( x − x0 )2 + .... 2! (9.24) Let us assume that our log likelihood function is differentiable and let us expand the log likelihood function in a Taylor series about the ML estimate θˆML . This gives log L(θˆML ) ''(θ − θˆML )2 log L(θ ) = log L(θˆML ) + log L(θˆML ) '(θ − θˆML ) + + ... 2! (9.25) Because log L(θˆML )' = 0 since θˆML is the ML estimate of θ , we obtain log L(θ ) ≈ log L(θˆML ) + log L(θˆML ) ''(θ − θˆML )2 . 2! (9.26) Therefore, we can approximate L(θ ) for θ close to θˆ as L(θ ) ≈ L(θˆML ) exp{− 12 I (θˆML )(θ − θˆML )2 } (9.27) I (θˆML ) = −log L(θˆML )". (9.28) where is the observed Fisher information. It is a number because it is evaluated at the MLE. I (θˆML ) measures the curvature of the likelihood around the MLE. Eq. 9.27 implies that near the MLE the likelihood can be approximated as a Gaussian distribution. We can gain intuition about the meaning of the Fisher information by considering the following Gaussian example. Example 3.2 (continued). If we assume σ 2 is known, then ignoring irrelevant constants we get log L( µ ) = − and 1 n ∑ (xi −µ )2 2σ 2 i=1 (9.29) page 10: 9.07 Lecture 9: Likelihood Methods n ∂ log L( µ ) 1 = (xi − µ ) ∂µ σ 2 i=1 ∑ ∂ 2 log L( µ ) ∂µ 2 =− n (9.30) (9.31) σ2 n . Hence, V ar( µˆ ) = n−1σ 2 = I ( µˆ )−1. Thus a key point is: A large amount of σ2 information implies smaller variance. Thus I ( µˆ ) = Example 2.1 (continued). For the binomial model of the learning experiment, the log likelihood and score functions (ignoring irrelevant constants) are log L( p) = k log p + (n − k )log(1 − p) (9.32) ∂ log L( p) k n − k = − ∂p p 1− p (9.33) The MLE is p̂ = k n (9.34) and I ( p) = − ∂ 2 log L( p) ∂p 2 = k p 2 + n−k (1− p)2 (9.35) so that at the MLE the observed Fisher information is I ( pˆ ) = n . p̂(1 − p̂) (9.36) See the Addendum 1 to Lecture 9 to understand the details of the last step. Remark 9.1. In principle, we can judge the quadratic approximation to the log likelihood by plotting the true log likelihood with the quadratic approximation superimposed. Remark 9.2. For the Gaussian distribution from Example 3.2, the quadratic approximation is exact. A practical rule for nearly all likelihood applications: a reasonably regular likelihood means θˆML is approximately Gaussian. Remark 9.3. The result in Eq. 9.27 suggests that an approximate 95% confidence interval for θ may be constructed as θˆML ± 2se(θˆML ) (9.37) page 11: 9.07 Lecture 9: Likelihood Methods where se(θˆML ) = [ I (θˆML )]− 2 . 1 Remark 9.4. Reporting θˆML and I (θˆML ) is a practical alternative to plotting the likelihood function in every problem. The quantity we derived above was the observed Fisher information. It is an estimate of the Fisher information which we define below. Definition 9.4. If x = x1,..., xn is a sample with joint probability density f ( x | θ ), then the Fisher information in x for the parameter θ is ∂ log f ( x | θ ) 2 ] ∂θ ∂ log f ( x | θ ) 2 = [ ] f ( x | θ ) dx. ∂θ I (θ ) = E[ ∫ (9.38) Equivalently, I (θ ) = −E[ =− ∫ ∂ 2 log f ( x | θ ) ∂θ 2 ∂ 2 log f ( x | θ ) ∂θ 2 ] (9.39) f ( x | θ )dx. Now that we have defined the Fisher information we can use it to make an important general statement about the performance of an estimator. Theorem 9.1 (Cramer-Rao Lower Bound). Suppose x = x1 , K , xn is a sample from f ( x | θ ), and T ( x) is an estimator of θ and Eθ [(T ( x))] is a differentiable function of θ . Suppose also that d df ( x | θ ) h( x) f ( x | θ )dx = h( x) dx, dθ dθ ∫ ∫ (9.40) for all h( x) with Eθ h( x) < ∞. Then dEθ (T ( x)) 2 ) dθ V ar(T ( x)) ≥ . ∂ log f ( x | θ ) 2 ) Eθ ( ∂θ ( (9.41) Equation 9.41 is the Cramer-Rao Lower Bound (CRLB) for the estimator T ( x). For a proof see Casella and Berger (1990). CRLB gives the lowest bound on the variance of an estimator. If the estimate is unbiased, then the numerator is 1 and the denominator is the Fisher information. If θ is a d ×1 vector then the Fisher information is a d × d matrix given by page 12: 9.07 Lecture 9: Likelihood Methods I (θ ) = Eθ [( ∂ log f ( x | θ ) T ∂ log f ( x | θ ) ∂ 2 log f ( x | θ ) ] = −Eθ [ ]. ∂θ ∂θ ∂θ∂θ T (9.42) We will make extensive use of the Fisher information to derive confidence intervals for our estimates and to make inferences in our problems. E. Criteria for Evaluating an Estimator If T ( x) is an estimator of θ , then there are several criteria that can be used to evaluate its performance. Four commonly used criteria to evaluate the performance of estimators are: 1. Mean-Squared Error (MSE) Eθ [T ( x) − θ ]2 The smaller the MSE the better the estimator. 2. Unbiasedness Eθ (T ( x)) = θ The bias of an estimator is bT (θ ) = E(T ( x)) −θ . An estimator is unbiased if its expected value is the parameter it estimates. 3. Consistency p T ( x) ⎯⎯ →θ as n → ∞ Consistency means that as the sample size increases the estimator converges in mean-squared error, probability or almost surely to the true value of the parameter. 4. Efficiency Achieves a minimum variance (Cramer-Rao Lower Bound). The efficiency of an estimator is usually characterized in terms of its variance relative to another estimator. The Cramer-Rao Lower bound gives a characterization of how small the variance of an estimator can be. F. Summary of Key Properties of Maximum Likelihood Estimation 1. Maximum likelihood estimates are generally biased. 2. Maximum likelihood estimates are consistent, hence, they are asymptotically unbiased. 3. Maximum likelihood estimates are asymptotically efficient. 4. The variance of the maximum likelihood estimate may be approximated by the reciprocal (inverse) of the Fisher information. 5. The maximum likelihood estimate θˆML is asymptotically Gaussian with mean θ and variance I(θ )−1. page 13: 9.07 Lecture 9: Likelihood Methods 6. If θˆML is the maximum likelihood estimate of θ then h(θˆML ) is the maximum likelihood estimate of h(θ ). 7. Given that θˆML : N(θ , I (θ )−1 ), by the Invariance Principle for the MLE and a Taylor Series expansion we have h(θˆML ) ≈ N (h(θ ), h '(θ )2 [ I (θ )]−1 ) (9.43) and an approximate 95% confidence interval for h(θ ) is 1 h(θˆML ) ±1.96h '(θˆML )[ I (θˆML )]− 2 . (9.44) See Addendum 2 to Lecture 9 for an explanation and illustration of this result. III. Summary The likelihood function and maximum likelihood estimation are two of the most important theoretical and practical concepts in statistics. Acknowledgments I am grateful to Uri Eden for making the figures, to Jim Mutch for careful proofreading and comments and to Julie Scott for technical assistance. Textbook References Casella, G, Berger RL. Statistical Inference. Belmont, CA: Duxbury, 1990. DeGroot MH, Schervish MJ. Probability and Statistics, 3rd edition. Boston, MA: Addison Wesley, 2002. Pawitan Y. In All Likelihood: Statistical Modeling and Inference Using Likelihood. London: Oxford, 2001. Rice JA. Mathematical Statistics and Data Analysis, 3rd edition. Boston, MA, 2007. Literature Reference Kass RE, Ventura V, Brown EN. Statistical issues in the analysis of neuronal data. Journal of Neurophysiology 2005, 94: 8-25.