Estimation Theory: Method-of-Moments Estimation

9.07 Introduction to Statistics for Brain and Cognitive Sciences Emery N. Brown Lecture 8 Estimation Theory: Method-of-Moments Estimation I. Objectives Understand the statistical paradigm for reasoning under uncertainty. Understand the principle of estimation in the statistical paradigm. Learn how to construct method-of-moments estimates for standard probability models. Understand how we assess the uncertainty in the method-of-moments estimates in terms of confidence intervals. Having completed our section on probability theory, we will study statistics in the balance of the course. We review the logic of the statistics paradigm. II. The Statistics Paradigm Statistics is the science of making decisions under uncertainty using mathematical models derived from probability theory. Given a sample of data, x1,..., xn , a statistic is any function of the data. For example the mean, median, mode, standard deviation, variance, kurtosis, and the skewness are all statistics. A. THE STATISTICAL PARADIGM (Box, Tukey) Question Preliminary Data (Exploration Data Analysis) Models Experiment (Confirmatory Data Analysis) Model Fit Goodness-of-Fit Assessment not satisfactory Satisfactory Make an Inference Make a Decision Figure 8A. The Box-Tukey paradigm for iterative model building based on the concepts of exploratory and confirmatory data analysis. page 2: 9.07 Lecture 8: Method-of-Moments Estimation The Box-Tukey Paradigm (Figure 8A) defines an 8-step approach to experimental design statistical model building and data analysis. The iterative cycle of this approach to statistical model building, data analysis and inference is attributed to George Box. The initial steps are termed exploratory data analysis (EDA). The idea of EDA was formalized by John Tukey. In this phase the statistical model builder uses whatever knowledge is available to propose an initial model or set of models for the data to be analyzed. In this phase a set of descriptive, usually non-model based statistical tools, are used to perform a preliminary analysis of the data. This may entail analysis of data from the literature or preliminary data collected in the current laboratory. The objective of the EDA phase of the analysis is to develop a set of initial probability models that may be used to describe the stochastic (random) structure in the data. The EDA techniques we have learned thus far are data plots, histograms, five-number summaries, stem-and-leaf plots and box plots. Confirmatory data analysis (CDA) uses formal statistical model fitting and assessments of goodness-of-fit to evaluate how well a model fits the data. Once a satisfactory fit to the data is obtained as evaluated in a goodness-of-fit assessment, the scientist may go ahead to make an inference and then a decision. The first part of CDA is proposing an appropriate probability model which is why we just learned a large set of probability models and their properties. We now learn how to use these probability models and extensions of them to perform statistical analyses of data. We will develop our ability to execute the steps in this paradigm and hence, use it to conduct statistical analyses during the next seven weeks. III. Method-of-Moments Estimation A. What is Estimation? Once a probability model is specified for a given problem, the next step is estimation of the model parameters. By estimation we mean using a formal procedure to compute the model parameter from the observed data. We term any formal procedure that tells us how to compute a model parameter from a sample of data an estimator. We term the value computed from the application of the procedure to actual data the estimate. In the next series lectures we will study three methods of estimation: method-ofmoments, maximum likelihood and Bayesian estimation. Example 2.1 (continued). On the previous day of the learning experiments the animal executed 22 of the 40 trials correctly. A reasonable probability model for this problem is the binomial distribution with n = 40 and p unknown, if we assume the set of trials to be independent. We simply took p̂, the estimate of p to be 22 40 = 0.55. While this seems reasonable, in what sense is this the “best” estimate? Is this still the best estimate if the trials are not independent? Example 3.2 (continued). The MEG measurements appeared Gaussian in shape (Figure 3E). This was further supported by the Q − Q plot goodness-of-fit analysis in Figure 3H. The Gaussian distribution has two parameters µ and σ 2 . We estimated µ as x = −0.2 ⋅10−11 and σˆ 2 = 2.0 ⋅10−22 , where x = n −1 n ∑ xi i=1 page 3: 9.07 Lecture 8: Method-of-Moments Estimation σˆ 2 = n −1 n ∑ (xi − x )2 , i=1 where we recall that x and σˆ 2 are the sample mean and sample variance respectively. Again, are these the best estimates of the parameters? The plot of the estimated probability density in Figure 3E and the Q − Q plot in Figure 3H suggest that these estimates are very reasonable, i.e. highly consistent with the data. Example 8.1. Spiking Activity of Retinal Neurons (Iyengar and Liu, 1997). Retinal neurons are grown in culture under constant light and environmental conditions. The spontaneous spiking activity of the neurons is recorded. The investigators wanted to develop a statistical model that accurately describes the stochastic structure of this activity. That is, the interspike interval distribution. The gamma and inverse Gaussian distributions have been suggested as possible probability models. How do we compute the parameters of these probability models from the observed measurements? We perform first an exploratory data analysis by making the boxplot, the five-number summary and to aid intuition with the interpretation of the boxplot (Figure 8B), we include a histogram (Figure 8C). Courtesy of Springer. Used with permission. Figure 8B. Boxplot of Interspike Interval of Spiking Activity of Retinal Neurons under Constant Conditions in Slice (Iyengar and Liu, Biol. Cyber. 1997). The five-number summary for the retinal neuron spiking data is: the minimum is 2 msec, the 25th percentile is 5 msec, the median is 10 msec, the 75th percentile is 43 msec, and the maximum is 190 msec. page 4: 9.07 Lecture 8: Method-of-Moments Estimation Figure 8C. Interspike Interval Histogram of Spiking Activity of Retinal Neurons under Constant Conditions in Slice (Iyengar and Liu, Biol. Cyber. 1997). For each of these examples, how do we measure the uncertainty in these parameter estimates? B. Method-of-Moments Estimation One very straight forward, intuitively appealing approach to estimation is the method-of -moments. It is simple to apply. Recall that the theoretical moments were defined in Lecture 6 as follows Given p( x) a probability mass function we define the i th theoretical moment as µi = ∑ xij p( x j ), (8.1) j for i = 1, 2,..., and given f ( x) a probability density function, we define the i th theoretical moment as ∫ µi = xi f ( x)dx, (8.2) for i = 1, 2,...,. These are the moments about zero, or the non-central moments. Note that the first moment is simply the mean and that the variance is page 5: 9.07 Lecture 8: Method-of-Moments Estimation σ 2 = E ( x − µ )2 ∫ = ∫ (x 2 − 2 µ x + µ 2 ) f ( x)dx = ∫ x 2 f ( x)dx − µ 2 ∫ f ( x)dx = (x −µ ) 2 f ( x)dx (8.3) = µ 2 − µ12 . The same result holds for discrete random variables with the integrals changed to summations. The variance is the second central moment. It is imperative that we add at this point that these definitions only make sense if these moments exist. There are certainly cases of probability densities that are well defined for which the moments do not exist. One often quoted example is the Cauchy distribution defined as f ( x) = 1 2π (1+ x 2 ) (8.4) , for x ∈ (−∞, ∞). Given a sample of data x1,..., xn , the i th sample moment is −1 µ% i =n n ∑ xri (8.5) r=1 for i = 1, 2,3,...,. The sample moments are computed by placing mass of n −1 at each data point. Because all the numbers are finite and the sum is finite, the sample moment of each order i exists. Method-of-Moments Estimation. Assume we observe x1,..., xn , as a sample from f ( x | θ ), a probability model with a d − dimensional unknown parameter θ . Assume that the first d moments of f ( x | θ ) are finite. Equate the theoretical moments to the corresponding sample moments and estimate θ as the solution to the d − dimensional equation system µi (θ ) = µ% i (8.6) for i = 1,..., d. The estimate of θ computed from Eq. 8.6 is the method-of-moments estimate and we denote it as θˆMM . Example 8.1 (continued). To fit the gamma probability model to the retinal ganglion cell interspike model using the method-of-moments, we consider µ = E ( X ) = αβ −1 2 2 σ = E ( X − µ ) = αβ (8.7) −2 (8.8) page 6: 9.07 Lecture 8: Method-of-Moments Estimation Note that we have used the second central moment instead of the second moment about zero to simplify the computation. Because θ = (α , β ) is 2-dimensional, we compute the first sample moment and the second sample central moment as n ∑ xr (8.9) ∑ (xr − x)2 (8.10) x = n −1 r=1 σˆ 2 = n−1 n r=1 and our 2-dimensional system of equations is x = αβ −1 2 σˆ = αβ −2 (8.11) (8.12) Solving Eqs. 8.11 and 8.12 yields the method-of-moments estimates βˆMM = x (σˆ 2 )−1 α̂ MM = x 2 (σˆ 2 )−1. (8.13) Example 2.1 (continued). To fit the binomial probability model to the rat learning data using the method-of-moments we have x = n−1 n ∑ xi and r=1 x=p (8.14) p̂MM = x (8.15) The one-dimensional equation is and p̂MM = n−1 n ∑ xr = (40)−1 22 = 0.55. (8.16) r=1 Example 3.2 (continued). To fit the Gaussian probability model to the MEG noise data using method-of-moments we have simply µ̂MM = x 2 σ̂ MM = n−1 n ∑ (xi −µ̂MM )2 (8.17) r=1 Remark 8.1. The method-of-moments estimates are simple to compute. As we showed, they may not exist. They are not unique. To see this latter point, note that we have taken the first d page 7: 9.07 Lecture 8: Method-of-Moments Estimation moments to construct the estimation equations but in principle and in practice, we could use any d theoretical and d sample moments. To see this consider x1,..., xn as a sample from a Poisson distribution with unknown parameter λ. Recall that for a Poisson distribution E( X ) = λ (8.18) Var ( X ) = λ. (8.19) and Hence, we have λMM = µ̂1 = n−1 n ∑ xr (8.20) r=1 λMM = µ̂2 − ( µ̂1 )2 = n−1 n ∑ xr2 − (n−1 r=1 n ∑ xr )2 . (8.21) r=1 Is one of these estimates “better”? We will discuss this further when we discuss maximum likelihood estimation (See Kass, Ventura and Brown, 2005). C. Uncertainty in the Method-of-Moments Estimate 1. Confidence Interval for the Mean We computed above the estimate of the probability of a correct response as p̂ = 0.55 for the learning experiment in Example 2.1. The method-of-moments procedure requires a justification. We certainly do not believe that this is the exact answer. There is obvious uncertainty and we would like to quantify that uncertainty. This leads us naturally to a discussion of the error in this estimate or alternatively a confidence interval for the true parameter. We recall from the Central Limit Theorem (Lecture 7) that if X1,..., X n is a sample from a distribution with mean µ and variance σ 2 then for n large X : N (µ , σ2 n ) (8.22) We can use this fact to construct a confidence interval for µ as a way of measuring the uncertainty in X , which for the Example 2.1 is p̂MM . Assume that σ 2 is known. Then from the 1 Central Limit Theorem, if we take Z = n 2 ( X − µ) σ then we have page 8: 9.07 Lecture 8: Method-of-Moments Estimation Pr(| Z |< 1.96) = 0.95 = Pr(−1.96 ≤ Z ≤ 1.96) 1 = Pr(−1.96 ≤ = Pr( −1.96σ n2 (X − µ) σ ≤ X −µ ≤ ≤ 1.96) 1.96σ (8.23) ) 1 n n2 1.96σ −1.96σ = Pr( 1 ≥ µ − X ≥ ) 1 2 n n2 1.96σ 1.96σ = Pr( X − 1 ≤ µ ≤ X + ) 1 n2 n2 1 2 Since 1.96 ≈ 2, we have the frequently quoted statement that a 95% confidence interval for the mean of a Gaussian distribution is X− The quantity 2σ n 1 2 ≤µ≤X+ 2σ 1 n2 (8.24) . σ is the standard error (SE) of the mean. This is simply an application of Kass’ 1 n2 2/3-95% rule applied to the Gaussian approximation from the Central Limit Theorem of the distribution of X (Lecture 3). Remark 8.2. If the random sample is from a Gaussian distribution, then the result is exact. Remark 8.3. If the random sample is from a binomial distribution (Example 2.1), then we have σ 2 p(1 − p) pˆ = x and Var ( p̂) = Var ( x ) = = . Hence, the approximate 95% confidence interval is n n 1 p̂ ± 2[ p̂(1 − p̂) 2 ] n (8.25) For the learning experiment in Example 2.1, this becomes 1 0.55(.45) 2 ] 40 0.55 ± 2× 0.078 0.55 ± 2[ Hence, the confidence interval for p is [0.394, 0.706]. Also notice that p(1 − p) ≤ (8.26) 1 4 for all p so that an often used upper bound interval is p̂ ± 2[ or 1 12 ] 4n (8.27) page 9: 9.07 Lecture 8: Method-of-Moments Estimation p̂ ± 1 1 n2 (8.28) . Example 8.2. Choosing n to set the length of the confidence interval. How many samples should be taken to insure that the length of the 95% confidence interval is less than 0.1? The 1 1 1 1 p(1 − p) 2 1 1 1 ] ≤ 4[ ] 2 = 2[ ] 2 . If we set 2[ ] 2 = 0.1, Then we n 4n n n find that n = 400 . Notice that this choice of sample size is independent of the true value of p . length of the confidence interval is 4[ Remark 8.4. If the random sample is from a Poisson distribution then we have that the confidence interval is 1 ⎛ λˆ ⎞ 2 λ̂ ± 2 ⎜⎜ ⎟⎟ . ⎝n⎠ (8.29) What are the conditions for this result to be reasonable? Remark 8.5. We can construct 90% or 99% confidence intervals, or in general intervals of any desired percentage. The chosen percentage is the confidence level. This is often confused with the significance level of a statistical test. These concepts are related but different. We will discuss this in Lecture 12. For a 90% confidence interval we have x± 1.64σ 1 n2 , (8.30) . (8.31) and for a 99% confidence interval we have x± 2.57σ 1 n2 For Gaussian and symmetric distributions we construct the confidence intervals symmetrically. In general this need not be the case. 2. The Classical Interpretation of the Confidence Interval The classical (frequentist) interpretation of the confidence interval is that stated in terms of the operating characteristics of the confidence interval. That is, if we were to estimate µ repeatedly with new data then, in the long-run (with infinitely many repetitions), about 95% of the time the 2σ interval x ± 1 would cover µ . In other words if we considered using the procedure again and n2 again, it would work in the sense that the confidence we state would correspond, more or less, to the actual probability that the confidence interval would cover the unknown values of µ . We can illustrate the long-run behavior in a simulation study. Let us return to the Example 2.1 and assume that n = 40 and suppose that the true p = 0.6. That is, the animal’s true propensity to respond correctly is 0.6 . It is easy to simulate 10,000 Binomial observations from the Binomial distribution with n = 40 and p = 0.6. We can ask, how often does the interval in Eq. 8.25 contain the true value of 0.6. The simulation takes the following steps: page 10: 9.07 Lecture 8: Method-of-Moments Estimation 1. For i =1 to 10,000 Generate X i : Binomial (n = 40, p = 0.6) Compute p̂i = X i / 40 1 Compute sei = ( p̂i (1 − p̂i ) / 40) 2 2. Sum the number of i for which p̂i − 2sei < 0.6 < p̂i + 2sei . We did this and obtained 9,454 as the sum in step 2. As a result, we conclude that to two-digit accuracy, the true value of p is covered by the interval in Eq. 8.25 95% of the time. We therefore see the sense in which the standard procedure is well calibrated for large samples. One way probability is often interpreted is in terms of bets: if we say an event has probability of 0.5, this indicates a willingness to take either side of a bet that the event will happen. If it has probability 0.95, this indicates a willingness to take either side of a bet that gives 19:1 odds the event will happen. Here, if we were to take either side of a bet that gives 19:1 odds the that the interval will cover the true value of p, in the long run (playing this game a very large number of times) we would expect to neither win nor lose money but rather to moreor-less break even. In this sense, confidence can be interpreted as making a statement about the operating characteristics of the confidence interval procedure. We will review other interpretations of the intervals. 3. Confidence Interval for the Variance of a Gaussian Distribution These results have allowed us to compute the approximate distribution of the mean of a random sample from a distribution with finite mean and variance. What is the distribution of the sample variance? While we will not answer this question generally, we will answer it for the case of a Gaussian distribution. If x1, x2 ,..., xn is a sample from a Gaussian distribution with mean µ and variance σ 2 , and we define σˆ 2 = n−1 n ∑ (xr − x)2 then r=1 nσ̂ 2 σ 2 n = ∑ (xr − x)2 σ r=1 2 : χ(2n−1) . (8.32) We establish the result in the Addendum to Lecture 8. We can use this fact to construct a confidence interval 2 2 2 Pr( χ n−1 (0.025) ≤ nσˆ2 ≤ χ n−1 (0.975) = Pr( σ = Pr( 1 2 χ n−1 (0.025) 2 nσˆ 2 χ n−1 (0.975) 2 ≥ σ2≥ nσˆ 2 ≤σ ≤ 1 2 χ n−1 (0.975) 2 nσˆ 2 χ n−1 (0.025) ) (8.33) ) Remember that the χ 2 distribution with n −1 degrees of freedom is a gamma distribution with α= (n −1) 1 and β = . 2 2 page 11: 9.07 Lecture 8: Method-of-Moments Estimation D. The Gaussian Distribution Revisited When we do not know the value of σ (as is almost always the case), we can estimate it with the sample standard deviation, usually denoted by σˆ, or s = (n −1) −1 n ∑ (xi − x )2 and use i=1 either of these in place of σ . This results in an approximate confidence interval which applies for n large. The approximate standard error of the mean x becomes se( x ) = σ̂ 1 n2 (8.34) . The approximate 95% confidence interval is X− 1.96σ̂ n 1 2 ≤µ≤X+ 1.96σ̂ (8.35) 1 n2 The approximation improves when the data are more nearly Gaussian and may not be as good in moderate sample sizes (say, 25 observations) with strongly non-Gaussian data. Hence, it is important to check the Gaussian assumption with Q − Q plots. For small samples and Gaussian data, it is possible to adjust the confidence intervals for the estimation of σ by s. The interval becomes slightly wider than the “exact” interval in Eq. 8.23. We are less certain about µ when there is added uncertainty about σ 2 . The correction uses the t − distribution. The t − distribution is a unimodal symmetric distribution which looks very much like the Gaussian. The particular t − distribution used depends on the sample size used to compute x and s 2 . The distributions are lower in the middle and have thicker tails than the Gaussian. If the sample size is n then ν = n −1 is the number of degrees of freedom, and the appropriate t − distribution to use is the t − distribution with ν degrees of freedom. Because the tails of the t − distribution are fatter, the value from the t − distribution corresponding to 1.96 on the Gaussian distribution will be larger (See Table 8.1 and Figure 8D). The 95% confidence interval based on the t − distribution if x1, x2 ,..., xn is a sample from a Gaussian distribution with mean µ and variance σ 2 is x − t0.975,n−1 s n 1 2 ≤ µ ≤ x + t0.975,n−1 s 1 n2 (8.36) As n increases t0.975 goes to 1.96. In general as n increases the t − distribution on ν degrees of freedom converges to a Gaussian distribution with mean 0 and variance 1. page 12: 9.07 Lecture 8: Method-of-Moments Estimation ν 10 15 20 25 30 0.975 2.228 2.131 2.086 2.060 2.042 40 60 120 ∞ 2.021 2.000 1.980 1.960 Table 8.1. The 0.975 quantiles of the t-distributions as a function of degrees of freedom. Figure 8D. Family of t-distributions plotted as a function of the number of degrees of freedom ν . The limiting distribution is the standard Gaussian. For sample sizes of 15 and less the t correction is important. IV. Summary Method-of-moments is a simple intuitive approach to estimating the parameters in a probability density. We will analyze further the extent to which these methods are “best” once we have studied likelihood theory (Lecture 9) and Bayesian methods (Lecture 10). We will use the bootstrap (Lecture 11) to compute estimates of uncertainty for the method-of-moments estimates. 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. Several sections of this lecture have been taken with permission from the class notes for the course Statistical Methods in Neuroscience and Psychology written by Robert Kass in the Department of Statistics at Carnegie Mellon University. page 13: 9.07 Lecture 8: Method-of-Moments Estimation Textbook Reference DeGroot MH, Schervish MJ. Probability and Statistics, 3rd edition. Boston, MA: Addison Wesley, 2002. Iyengar S, Liao Q. Modeling neural activity using the generalized inverse Gaussian distribution. Biological Cybernetics, 1997, 77: 289–295. 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. Website Reference Kass RE, Chapter 3: Estimation , 2006.