Maximum Likelihood Estimators

Lecture 2 Maximum Likelihood Estimators. Matlab example. As a motivation, let us look at one Matlab example. Let us generate a random sample of size 100 from beta distribution Beta(5, 2). We will learn the definition of beta distribution later, at this point we only need to know that this isi a continuous distribution on the interval [0, 1]. This can be done by typing ’X=betarnd(5,2,100,1)’. Let us fit different distributions by using a distribution fitting tool ’dfittool’. We try to fit normal distribution and beta distribution to this sample and the results are displayed in figure 2.1. 1 3 100 samples ~ Beta(5,2) Normal fit Beta fit 0.9 100 samples ~ Beta(5,2) Normal fit Beta fit 0.8 Cumulative probability 2.5 Density 2 1.5 1 0.7 0.6 0.5 0.4 0.3 0.2 0.5 0.1 0 0.3 0.4 0.5 0.6 0.7 0.8 0.9 0 0.3 1 Data 0.4 0.5 0.6 0.7 0.8 0.9 1 Data Figure 2.1: Fitting a random sample of size 100 from Beta(5, 2). (a) Histogram of the data and p.d.f.s of fitted normal (solid line) and beta (dashed line) distributions; (b) Empirical c.d.f. and c.d.f.s of fitted normal and beta distributions. Besides the graphs, the distribution fitting tool outputs the following information: Distribution: Log likelihood: Normal 55.2571 7 Domain: Mean: Variance: Parameter mu sigma -Inf < y < Inf 0.742119 0.0195845 Estimate 0.742119 0.139945 Estimated covariance mu mu 0.000195845 sigma 6.01523e-020 Distribution: Log likelihood: Domain: Mean: Variance: Parameter a b Std. Err. 0.0139945 0.00997064 of parameter estimates: sigma 6.01523e-020 9.94136e-005 Beta 63.8445 0 < y < 1 0.741371 0.0184152 Estimate 6.97783 2.43424 Std. Err. 1.08827 0.378351 Estimated covariance of parameter estimates: a b a 1.18433 0.370094 b 0.370094 0.143149 The value ’Log likelihood’ indicates that the tool uses the maximum likelihood estimators to fit the distribution, which will be the topic of the next few lectures. Notice the ’Parameter estimates’ - given the data ’dfittool’ estimates the unknown parameters of the distribution and then graphs the p.d.f. or c.d.f. corresponding to these parameters. Since the data was generated from beta distribution, it is not surprising that beta distribution fit seems better than normal distribution fit, which is particularly clear from figure 2.1 (b), that compares how estimated c.d.f. fits the empirical c.d.f. Empirical c.d.f. is defined as n 1� Fn (x) = I(Xi � x) n i=1 where I(Xn � x) is the indicator that Xi is � x. In other words, Fn (x) is the proportion of observations below level x. One can ask several questions about this example: 1. How to estimate the unknown parameters of a distribution given the data from this distribution? 8 2. How good are these estimates, are they close to the actual ’true’ parameters? 3. Does the data come from a particular type of distribution, for example, normal or beta distribution? In the next few lectures we will study the first two questions and we will assume that we know what type of distribution the sample comes from, so we only do not know the parameters of the distribution. In the context of the above example, we would be told that the data comes from beta distribution, but the parameters (5, 2) would be unknown. Of course, in general we might not know what kind of distribution the data comes from - we will study this type of questions later when we look at the so called goodness-of-fit hypotheses tests. In particular, we will see graphs like 2.1 (b) again when we study the Kolmogorov-Smirnov goodness-of-fit test. Example. We consider a dataset of various body measurements from [1] (dataset can be dowloaded from journal’s website), including weight, height, waist girth, abdomen girth, etc. First, we use Matlab fitting tool to fit weight and waist girth of men and women (separately) with lognormal distribution, see figure 2.2 (a) and (b). Wikipedia article about normal dis­ tribution gives a reference to a 1932 book ”Problems of Relative Growth” by Julian Huxley for the explanation why the sizes of full-grown animals are approximately log-normal. One short explanation is consistency between linear and volume dimensions - if linear dimensions are lognormal and volume dimensions are proportional to cube of linear dimensions then they also are lognormal. Assumption that sizes are normal would violate this consistency, since the cube of normal is not normal. We observe, hovewer, that the fit of women’s waist with lognormal is not very accurate. Later in the class we will learn several statistical tests to decide if the data comes from a certain distribution or a family of distributions, but here is a preview of what’s to come. Chi-squared goodness-of-fit test rejects the hypothesis that the distribution of logarithms of women’s waists is normal: [h,p,stats]=chi2gof(log_women_waist) h = 1, p = 5.2297e-004 stats = chi2stat: 22.0027 df: 5 edges: [1x9 double] O: [21 44 67 60 28 18 12 10] E: [1x8 double] and so does Lilliefor’s test (adjusted Kolmogorov-Smirnov test): [h,p,stats]=lillietest(log_women_waist) h = 1, p = 0, stats = 0.0841. The same tests accept the hypotheses that other variables have lognormal distribution. Au­ thor’s in [1] suggest that we can fit women’s waist with Gamma distribution. Since Gamma 9 1 0.9 0.9 0.8 0.8 0.7 0.7 Cumulative probability Cumulative probability 1 0.6 0.5 0.4 0.3 women’s weight lognormal fit women men’s weight lognormal fit men 0.2 0.1 0 50 60 70 80 Data 90 100 0.6 0.5 0.4 0.3 women’s waist lognormal fit women men’s waist lognormal fit men 0.2 0.1 0 110 60 70 80 90 100 110 Data 1 0.9 Cumulative probability 0.8 0.7 0.6 0.5 0.4 0.3 women’s waist (shifted) Gamma fit normal fit 0.2 0.1 0 5 10 15 20 25 30 35 40 45 Data Figure 2.2: Fitting weight (upper left) and waist girth (upper right) with lognormal distribution. Lower left: fitting women’s waist with shifted Gamma and normal distributions. does not have a translation (shift) parameter, when we fit Gamma distribution we can either add to it a shift parameter or instead shift all data to start at zero. In figure 2.2 (c) we fit Gamma and, for the sake of illustration, normal distribution, to women’s waist sample. As we can see, Gamma fits the data better than lognormal and much better than normal. To find the parameters of fitted Gamma distribution we use Matlab ’gamfit’ function: param=gamfit(women_waist_shift) param = 2.8700 4.4960. Chi-squared goodness-of-fit test for a specific (fitted) Gamma distribution: [h,p,stats]=chi2gof(women_waist_shift,’cdf’,@(z)gamcdf(z,param(1),param(2))) 10 h = 0, p = 0.9289, stats = chi2stat: 2.4763, df: 7 accepts the hypothesis that the sample has Gamma distribution �(2.87, 4.496). This test is not ’accurate’ in some sense, which will be explained later. One can also check that Gamma distribution fits well other variables - men’s waist girth, weight of men and weight of women. Let us consider a family of distributions P� indexed by a parameter (which could be a vector of parameters) ϕ that belongs to a set �. For example, we could consider a family of normal distributions N (�, α 2 ) in which case the parameter would be ϕ = (�, α 2 ) - the mean and variance of the distribution. Let f (X |ϕ) be either a probability function (in case of discrete distribution) or a probability density function (continuous case) of the distribution P� . Suppose we are given an i.i.d. sample X1 , . . . , Xn with unknown distribution P� from this family, i.e. parameter ϕ is unknown. A likelihood function is defined by �(ϕ) = f (X1 |ϕ) × . . . × f (Xn |ϕ). We think of the sample X1 , . . . , Xn as given numbers and we think of � as a function of the parameter ϕ only. The likelihood function has a clear interpretation. For example, if our distributions are discrete then the probability function f (x|ϕ) = P� (X = x) is the probability to observe a point x and the likelihood function �(ϕ) = f (X1 |ϕ) × . . . × f (Xn |ϕ) = P� (X1 ) × . . . × P� (Xn ) = P� (X1 , . . . , Xn ) is the probability to observe the sample X1 , . . . , Xn when the parameters of the distribution are equal to ϕ. In the continuous case the likelihood function �(ϕ) is the probability density function of the vector (X1 , . . . , Xn ). Definition: (Maximum Likelihood Estimators.) Suppose that there exists a parameter ϕ̂ that maximizes the likelihood function �(ϕ) on the set of possible parameters �, i.e. �(ϕ̂) = max �(ϕ). ��� Then ϕ̂ is called the Maximum Likelihood Estimator (MLE). When finding the MLE it sometimes easier to maximize the log-likelihood function since �(ϕ) � maximize ≥ log �(ϕ) � maximize maximizing � is equivalent to maximizing log �. Log-likelihood function can be written as log �(ϕ) = n � i=1 log f (Xi |ϕ). Let us give several examples of computing the MLE. 11 Example 1. Bernoulli distribution B(p). X = {0, 1}, P(X = 1) = p, P(X = 0) = 1 − p, p ∞ [0, 1]. Probability function in this case is given by � p, x=1 f (x|p) = = px (1 − p)1−x . 1 − p, x = 0 Likelihood function is �(p) = f (X1 |p)f (X2 |p) . . . f (Xn|p) = p# of 1’s (1 − p)# of 0’s = pX1 +...+Xn (1 − p)n−(X1 +...+Xn ) and the log-likelihood function is log �(p) = (X1 + . . . + Xn ) log p + (n − (X1 + . . . + Xn )) log(1 − p). To maximize this over p ∞ [0, 1] let us find the critical point (log �(p))� = 0, 1 1 (X1 + . . . + Xn ) − (n − (X1 + . . . + Xn )) = 0. p 1−p Solving this for p gives, X1 + . . . + X n ¯ =X n ¯ in the sample is the MLEstimator of the and, therefore, the proportion of successes p̂ = X unknown true probability of success, which is a very natural and intuitive estimator. For example, by law of large numbers, we know that p= ¯ � EX1 = p X in probability (we will recall this definition in the next lecture), which means that our estimate will approximate the unknown parameter p well when we get more and more data. Remark. In each example, once we compute the estimate of parameters, we can try to prove directly, using the explicit form of the estimate, that it approximates well the unkown parameters, as we did in Example 1. However, in the next lecture we will describe in a general setting that MLE has ’good properties’. Example 2. Normal distribution N(�, α 2 ). The p.d.f. of normal distribution is f (X|(�, α 2 )) = ≤ (X−�)2 1 e− 2�2 . 2�α and, therefore, likelihood function is �(�, α 2 ) = n � i=1 ≤ 12 (Xi −�)2 1 e− 2�2 . 2�α and log-likelihood function is n � � 1 (Xi − �)2 � log ≤ − log α − 2α 2 2� i=1 n 1 1 � = n log ≤ − n log α − 2 (Xi − �)2 . 2α i=1 2� 2 log �(�, α ) = 2 We want to maximize the log-likelihood with � respect2 to −→ < � < → and α > 0. First, obviously, for any α we need to minimize (Xi − �) over �. The critical point condition is n n � d � (Xi − �)2 = −2 (Xi − �) = 0 d� i=1 i=1 ¯ We can plug this estimate in the log-likelihood and solving this for � we get that � ˆ = X. and it remains to maximize n 1 1 � n log ≤ − n log α − 2 (Xi − X̄)2 2α i=1 2� over α. The critical point condition reads, − n 1 � + 3 (Xi − X̄)2 = 0 α α and solving this for α we obtain that the MLE of α 2 is n 1� α̂ = (Xi − X̄)2 . n i=1 2 The normal distribution fit in figure 2.1 corresponds to these parameters (�, ˆ α̂ 2 ). Exercise. Generate a normal sample in Matlab and fit it with a normal distribution using ’dfittool’. Then plot a p.d.f. or c.d.f. corresponding to MLE above and compare this with ’dfittool’. Let us give one more example of MLE. Uniform distribution U [0, ϕ] on the interval [0, ϕ]. This distribution has p.d.f. � 1 , 0 � x � ϕ, � f (x|ϕ) = 0, otherwise. The likelihood function �(ϕ) = n � i=1 f (Xi |ϕ) = = 1 I(X1 , . . . , Xn ∞ [0, ϕ]) ϕn 1 I(max(X1 , . . . , Xn ) � ϕ). ϕn 13 Here the indicator function I(A) equals to 1 if event A happens and 0 otherwise. What the indicator above means is that the likelihood will be equal to 0 if at least one of the factors is 0 and this will happen if at least one observation Xi will fall outside of the ’allowed’ interval [0, ϕ]. Another way to say it is that the maximum among observations will exceed ϕ, i.e. �(ϕ) = 0 if ϕ < max(X1 , . . . , Xn ), and �(ϕ) = 1 if ϕ � max(X1 , . . . , Xn ). ϕn Therefore, looking at the figure 2.3 we see that ϕ̂ = max(X1 , . . . , Xn ) is the MLE. 2 1.8 1.6 �(ϕ) 1.4 1.2 1 0.8 0.6 0.4 PSfrag replacements 0.2 0 0 0.5 1 1.5 2 2.5 max(X1 , . . . , Xn ) 3 ϕ Figure 2.3: MLE for the uniform distribution. Sometimes it is not so easy to find the maximum of the likelihood function as in the examples above and one might have to do it numerically. Also, MLE does not always exist. Here is an example: let us consider uniform distribution U[0, ϕ) and define the density by � 1 , 0 � x < ϕ, � f (x|ϕ) = 0, otherwise. The difference is that we ’excluded’ the point ϕ by setting f (ϕ |ϕ) = 0. Then the likelihood function is �(ϕ) = n � i=1 f (Xi |ϕ) = 1 I(max(X1 , . . . , Xn ) < ϕ) ϕn 14 and the maximum at the point ϕ̂ = max(X1 , . . . , Xn ) is not achieved. Of course, this is an artificial example that shows that sometimes one needs to be careful. References: [1] Grete Heinz, Louis J. Peterson, Roger W. Johnson, Carter J. Kerk, (2003) “Exploring Relationships in Body Dimensions“. Journal of Statistics Education, Volume 11, Number 2. 15