f15-6

684 Chapter15. ModelingofDataSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).15.6 Confidence Limits on Estimated Model Parameters Severaltimesalreadyinthischapterwehavemadestatementsaboutthestandard errors, or uncertainties, in a set of Mestimated parameters a. We have given some formulas for computing standard deviations or variances of individual parameters (equations 15.2.9, 15.4.15, 15.4.19), as well as some formulas for covariances between pairs of parameters (equation 15.2.10; remark following equation 15.4.15;equation 15.4.20; equation 15.5.15). In this section, we want to be more explicit regarding the precise meaning of these quantitative uncertainties, and to give further information about how quantitative confidence limits on fitted parameters can be estimated. The subject can get somewhat technical, and even somewhat confusing, so we will try to makeprecise statements, even when they must be offered without proof. Figure 15.6.1 shows the conceptual scheme of an experiment that “measures” a set of parameters. There is some underlying true set of parameters a truethat are known to Mother Nature but hidden from the experimenter. These true parameters arestatistically realized,alongwithrandommeasurementerrors,asameasureddata set,whichwewillsymbolizeas D(0). Thedataset D(0)isknowntotheexperimenter. He or she fits the data to a model by χ2minimizationor some other technique, and obtainsmeasured,i.e.,fitted, valuesfortheparameters,whichweheredenote a(0). Because measurement errors have a random component, D(0)is not a unique realization of the true parameters atrue. Rather, there are infinitely many other realizations of the true parameters as “hypothetical data sets” each of which could have been the one measured, but happened not to be. Let us symbolize these byD(1),D(2),.... Each one, had it been realized, would have given a slightly different set of fitted parameters, a(1),a(2),..., respectively. These parameter sets a(i)therefore occur with some probability distribution in the M-dimensional space of all possible parametersets a. The actual measured set a(0)is one memberdrawn from this distribution. Even more interesting than the probability distribution of a(i)would be the distribution of the difference a(i)−atrue. This distribution differs from the former onebyatranslationthatputsMotherNature’struevalueattheorigin. Ifweknew this distribution, we would know everything that there is to know about the quantitative uncertainties in our experimental measurement a(0). So the name of the game is to find some way of estimating or approximating theprobabilitydistributionof a(i)−atruewithoutknowing atrueandwithouthaving available to us an infinite universe of hypothetical data sets. Monte CarloSimulationof Synthetic Data Sets Although the measured parameter set a(0)is not the true one, let us consider a fictitious world in which it wasthe true one. Since we hope that our measured parameters are not toowrong, we hope that that fictitious world is not too different from the actual world with parameters atrue. In particular, let us hope — no, let us assume— that the shape of the probability distribution a(i)−a(0)in the fictitious worldisthesame,orverynearlythesame,astheshapeoftheprobabilitydistribution 15.6ConfidenceLimitsonEstimatedModelParameters 685Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).actual data set hypothetical data set hypothetical data set hypothetical data seta3 a2 a1 fitted parameters a 0χ2 min true parameters atrueexperimental realization ...... Figure 15.6.1. A statistical universe of data sets from an underlying model. True parameters atrueare realized in a data set, from which fitted (observed) parameters a0are obtained. If the experiment were repeated many times, new data sets and new values of the fitted parameters would be obtained. a(i)−atrueintherealworld. Noticethatwearenotassumingthat a(0)andatrueare equal; they are certainly not. We are only assuming that the way in which random errors enter the experiment and data analysis does not vary rapidly as a function ofa true, so thata(0)can serve as a reasonable surrogate. Now, often, the distribution of a(i)−a(0)in thefictitious world iswithin our power to calculate (see Figure 15.6.2). If we know something about the process that generated our data, given an assumed set of parameters a(0), then we can usuallyfigure out how to simulateour own sets of “synthetic”realizations of these parametersas “syntheticdatasets. ”Theprocedureis to drawrandomnumbersfrom appropriate distributions (cf. §7.2–§7.3) so as to mimic our best understanding of theunderlyingprocessandmeasurementerrorsin ourapparatus. With suchrandomdraws,weconstructdatasetswithexactlythesamenumbersofmeasuredpoints,and preciselythesamevaluesofallcontrol(independent)variables,asouractualdataset D (0). Letuscallthesesimulateddatasets DS (1),DS (2),.... Byconstructiontheseare supposed to have exactly the same statistical relationship to a(0)as the D(i)’sh a v e toatrue. (Forthe case whereyoudon ’tknowenoughaboutwhat youare measuring to do a credible job of simulating it, see below.) Next, for each DS (j), perform exactly the same procedure for estimation of parameters, e.g., χ2minimization, as was performed on the actual data to get the parameters a(0), giving simulated measured parameters aS (1),aS (2),.... Each simulated measured parameter set yields a point aS (i)−a(0). Simulate enough data setsandenoughderivedsimulatedmeasuredparameters,andyoumapoutthedesired probability distribution in Mdimensions. In fact, the ability to do Monte Carlo simulations in this fashion has revo- 686 Chapter15. ModelingofDataSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).synthetic data set 1 synthetic data set 2 synthetic data set 3 synthetic data set 4a2 χ2 min χ2 min(s)a1 (s) a3 (s) a4 (s)Monte Carlo parameters Monte Carlo realization fittedparameters a 0actual data set Figure15.6.2. MonteCarlosimulationofanexperiment. The fittedparametersfromanactualexperiment are used assurrogates for the true parameters. Computer-generated random numbers are used tosimulatemany synthetic data sets. Each of these is analyzed to obtain its fitted parameters. The distribution of thesefitted parameters around the (known) surrogate true parameters is thus studied. lutionized many fields of modern experimental science. Not only is one able to characterize the errors of parameter estimation in a very precise way; one can also try out on the computerdifferentmethods of parameterestimation, or differentdatareduction techniques, and seek to minimize the uncertainty of the result according to any desired criteria. Offered the choice between mastery of a five-foot shelf of analyticalstatistics booksandmiddlingabilityatperformingstatistical MonteCarlo simulations, we would surely choose to have the latter skill. Quick-and-Dirty Monte Carlo: TheBootstrap Method Here is a powerful technique that can often be used when you don ’t know enough about the underlying process, or the nature of your measurement errors, to do a credible Monte Carlo simulation. Suppose that your data set consists of Nindependent and identically distributed (oriid)“data points. ”Each data point probablyconsists ofseveralnumbers,e.g.,oneormorecontrolvariables(uniformly distributed, say, in the range that you have decided to measure) and one or more associatedmeasuredvalues(eachdistributedhoweverMotherNaturechooses). “Iid” meansthatthesequentialorderofthedatapointsisnotofconsequencetotheprocess that you are using to get the fitted parameters a. For example, a χ2sum like (15.5.5) does not care in what order the points are added. Even simpler examples are the mean value of a measured quantity, or the mean of some function of the measured quantities. Thebootstrapmethod [1]usestheactualdataset DS (0),withits Ndatapoints,to generate any number of synthetic data sets DS (1),DS (2),..., also with Ndata points. The procedureis simply to draw Ndata points at a time with replacement from the 15.6ConfidenceLimitsonEstimatedModelParameters 687Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).setDS (0). Because of the replacement, you do not simply get back your original data set each time. You get sets in which a random fraction of the original points,typically ∼1/e≈37%, are replaced by duplicated original points. Now, exactly as in the previous discussion, you subject these data sets to the same estimation procedure as was performed on the actual data, giving a set of simulated measuredparameters a S (1),aS (2),.... Thesewill be distributedaround a(0)in close to the same way that a(0)is distributed around atrue. Soundslikegettingsomethingfornothing,doesn ’tit? Infact,ithastakenmore thanadecadeforthebootstrapmethodtobecomeacceptedbystatisticians. Bynow, however,enoughtheoremshavebeenprovedtorenderthebootstrapreputable(see [2] forreferences). Thebasicideabehindthebootstrapisthattheactualdataset,viewed as a probability distribution consisting of delta functions at the measured values, is inmostcasesthebest —oronly—availableestimatoroftheunderlyingprobability distribution. It takes courage, but one can often simply use thatdistribution as the basis for Monte Carlo simulations. Watch out for cases where the bootstrap ’s“iid”assumption is violated. For example,ifyouhavemademeasurementsatevenlyspacedintervalsofsomecontrol variable,thenyoucan usuallygetawaywithpretendingthattheseare “iid,”uniformly distributed over the measured range. However, some estimators of a(e.g., ones involvingFouriermethods)mightbe particularlysensitiveto all the pointson a grid beingpresent. In that case, the bootstrapis goingto givea wrongdistribution. Alsowatchoutforestimatorsthatlookatanythinglikesmall-scaleclumpinesswithinthe Ndata points, or estimators that sort the data and look at sequential differences. Obviouslythebootstrapwill fail onthese, too. (Thetheoremsjustifyingthemethodarestilltrue,butsomeoftheirtechnicalassumptionsareviolatedbytheseexamples.) For a large class of problems, however, the bootstrap does yield easy, very quick, Monte Carlo estimates of the errors in an estimated parameter set. Confidence Limits Rather than present all details of the probability distribution of errors in parameter estimation, it is common practice to summarize the distribution in theform ofconfidence limits . The full probability distribution is a function de fined on the M-dimensional space of parameters a.Aconfidence region (orconfidence interval)isjustaregionofthat M-dimensionalspace(hopefullyasmallregion)that contains a certain (hopefully large) percentage of the total probability distribution. You pointto a con fidenceregionandsay, e.g., “thereis a 99 percentchancethat the true parameter values fall within this region around the measured value. ” It is worth emphasizing that you, the experimenter, get to pick both the confidencelevel (99 percent in the above example),and the shape of the con fidence region. The only requirementis that your regiondoes include the stated percentage ofprobability. Certainpercentagesare,however,customaryinscienti ficusage: 68.3 percent (the lowest con fidence worthy of quoting), 90 percent, 95.4 percent, 99 percent,and99.73percent. Highercon fidencelevelsareconventionally “ninety-nine point nine ...nine.”As for shape, obviously you want a region that is compact and reasonably centered on your measurement a (0), since the whole purpose of a confidence limit is to inspire con fidence in that measured value. In one dimension, 688 Chapter15. ModelingofDataSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).68% confidence interval on a268% confidence interval on a168% confidence region on a1 and a2 jointly biasa(i)1 − a(0)1(s)a(i)2 − a(0)2(s) Figure 15.6.3. Con fidence intervals in 1 and 2 dimensions. The same fraction of measured points (here 68%) lies (i) between the two vertical lines, (ii) between the two horizontal lines, (iii) within the ellipse. the convention is to use a line segment centered on the measured value; in higher dimensions, ellipses or ellipsoids are most frequently used. You might suspect, correctly, that the numbers 68.3 percent, 95.4 percent, and 99.73 percent, and the use of ellipsoids, have some connection with a normaldistribution. That is true historically, but not always relevant nowadays. In general, the probability distribution of the parameters will not be normal, and the above numbers, used as levels of con fidence, are purely matters of convention. Figure 15.6.3 sketches a possible probability distribution for the case M=2. Shown are three different con fidence regions which might usefully be given, all at thesamecon fidencelevel. Thetwoverticallinesencloseaband(horizontalinterval) whichrepresentsthe68percentcon fidenceintervalforthevariable a 1withoutregard to the value of a2. Similarly the horizontal lines enclose a 68 percent con fidence interval for a2. The ellipse shows a 68 percent con fidence interval for a1anda2 jointly. Noticethattoenclosethesameprobabilityasthetwobands,theellipsemust necessarily extend outside of both of them (a point we will return to below). Constant Chi-Square Boundaries as Confidence Limits When the methodused to estimate the parameters a(0)is chi-squareminimiza- tion, as in the previoussections of this chapter, then there is a naturalchoice for the 15.6ConfidenceLimitsonEstimatedModelParameters 689Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).C B A Z′Z C′∆χ2 = 6.63 ∆χ2 = 2.71 ∆χ2 = 1.00 ∆χ2 = 2.30A′ B′ Figure 15.6.4. Con fidence region ellipses corresponding to values of chi-square larger than the fitted minimum. The solid curves, with ∆χ2=1.00,2.71,6.63project onto one-dimensional intervals AA/prime, BB/prime,CC/prime. These intervals —not the ellipses themselves —contain 68.3%, 90%, and 99% of normally distributed data. The ellipse that contains 68.3% of normally distributed data is shown dashed, and has ∆χ2=2.30. For additional numerical values, see accompanying table. shape of con fidence intervals, whose use is almost universal. For the observeddata setD(0), the value of χ2is a minimum at a(0). Call this minimum value χ2 min.I f thevector aofparametervaluesisperturbedawayfrom a(0),then χ2increases. The region within which χ2increases by no more than a set amount ∆χ2defines some M-dimensional con fidence region around a(0).I f∆χ2is set to be a large number, this will be a big region;if it is small, it will be small. Somewherein between there will be choices of ∆χ2that cause the region to contain, variously, 68 percent, 90 percent, etc. of probability distribution for a’s, as defined above. These regions are taken as the con fidence regions for the parameters a(0). Very frequently one is interested not in the full M-dimensional con fidence region,butinindividualcon fidenceregionsforsomesmallernumber νofparameters. For example, one might be interested in the con fidence interval of each parameter taken separately (the bands in Figure 15.6.3), in which case ν=1. In that case, the naturalcon fidenceregionsin the ν-dimensionalsubspace ofthe M-dimensional parameter space are the projections of the M-dimensional regions de fined byfixed ∆χ2intothe ν-dimensionalspacesofinterest. InFigure15.6.4,forthecase M=2, we show regions corresponding to several values of ∆χ2. The one-dimensional confidence interval in a2corresponding to the region bounded by ∆χ2=1lies between the lines AandA/prime. Notice that the projection of the higher-dimensional region on the lower- dimension space is used, not the intersection. The intersection would be the band between ZandZ/prime.I ti sneverused. It is shownin the figureonlyfor thepurposeof 690 Chapter15. ModelingofDataSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).makingthis cautionarypoint, that it shouldnot be confusedwith the projection. Probability Distributionof Parameters inthe NormalCase You may be wondering why we have, in this section up to now, made no connection at all with the error estimates that come out of the χ2fitting procedure, most notably the covariance matrix Cij. The reason is this: χ2minimization is a useful means for estimating parameters even if the measurement errors arenot normally distributed. While normally distributed errors are required if the χ 2 parameter estimate is to be a maximum likelihood estimator ( §15.1), one is often willing to give up that property in return for the relative convenience of the χ2 procedure. Only in extreme cases, measurement error distributions with very large “tails,”isχ2minimization abandoned in favor of more robust techniques, as will be discussed in §15.7. However,theformalcovariancematrixthatcomesoutofa χ2minimizationhas aclearquantitativeinterpretationonlyif(ortotheextentthat)themeasurementerrorsactuallyarenormallydistributed. Inthecaseof nonnormalerrors,youare “allowed” •tofit for parameters by minimizing χ 2 •touseacontourofconstant ∆χ2astheboundaryofyourcon fidenceregion •to use Monte Carlo simulation or detailed analytic calculation in deter- miningwhichcontour ∆χ2is the correct one for your desired con fidence level •to give the covariance matrix Cijas the“formal covariance matrix of thefit.” You are notallowed •to use formulas that we now give for the case of normal errors, which establish quantitative relationships among ∆χ2,Cij, and the con fidence level. Here are the key theorems that hold when (i) the measurement errors are normally distributed, and either (ii) the model is linear in its parameters or (iii) the sample size is large enough that the uncertainties in the fitted parameters ado not extendoutsidearegioninwhichthemodelcouldbereplacedbyasuitablelinearizedmodel. [Note that condition (iii) does not preclude your use of a nonlinear routine likemqrfittofindthefitted parameters.] Theorem A. χ 2 minis distributed as a chi-square distribution with N−M degrees of freedom, where Nis the number of data points and Mis the number of fittedparameters. Thisisthebasictheoremthatletsyouevaluatethegoodness-of- fit of the model, as discussed above in §15.1. We list it first to remind you that unless the goodness-of- fitis credible, the whole estimation of parameters is suspect. Theorem B. IfaS (j)is drawn from the universe of simulated data sets with actual parameters a(0), then the probability distribution of δa≡aS (j)−a(0)is the multivariate normal distribution P(δa)da1...d a M=const.×exp/parenleftbigg −1 2δa·[α]·δa/parenrightbigg da1...d a M where [α]is the curvature matrix de fined in equation (15.5.8). Theorem C. IfaS (j)is drawn from the universe of simulated data sets with actualparameters a(0),thenthequantity ∆χ2≡χ2(a(j))−χ2(a(0))isdistributedas 15.6ConfidenceLimitsonEstimatedModelParameters 691Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).achi-squaredistributionwith Mdegreesoffreedom. Herethe χ2’sareallevaluated using the fixed (actual) data set D(0). This theorem makes the connection between particular values of ∆χ2and the fraction of the probability distribution that they encloseasan M-dimensionalregion,i.e.,thecon fidencelevelofthe M-dimensional confidence region. Theorem D. Suppose that aS (j)is drawn from the universe of simulated data sets (as above), that its firstνcomponents a1,...,a νare heldfixed, and that its remaining M−νcomponents are varied so as to minimize χ2. Call this minimum value χ2 ν. Then ∆χ2 ν≡χ2 ν−χ2 minis distributed as a chi-square distribution with νdegrees of freedom. If you consult Figure 15.6.4, you will see that this theorem connectsthe projected ∆χ2regionwithacon fidencelevel. Inthe figure,apointthat is heldfixed in a2and allowedto varyin a1minimizing χ2will seek outthe ellipse whose top or bottom edge is tangent to the line of constant a2, and is therefore the line that projects it onto the smaller-dimensional space. As afirst example, let us consider the case ν=1, where we want to find the confidence interval of a single parameter, say a1. Notice that the chi-square distributionwith ν=1degreeoffreedomisthesamedistributionasthatofthesquare of a single normallydistributed quantity. Thus ∆χ2 ν<1occurs 68.3percent of the time(1- σforthenormaldistribution), ∆χ2 ν<4occurs95.4percentofthetime(2- σ for the normal distribution), ∆χ2 ν<9occurs 99.73 percent of the time (3- σfor the normal distribution), etc. In this manneryou find the ∆χ2 νthat correspondsto your desiredcon fidencelevel. (Additionalvalues aregiveninthe accompanyingtable.) Letδabe a change in the parameters whose first component is arbitrary, δa 1, but the rest of whose componentsare chosen to minimize the ∆χ2. Then Theorem D applies. The value of ∆χ2is given in general by ∆χ2=δa·[α]·δa (15.6.1 ) which follows from equation (15.5.8) applied at χ2 minwhere βk=0. Since δaby hypothesis minimizes χ2in all but its first component, the second through Mth components of the normal equations (15.5.9) continue to hold. Therefore, the solution of (15.5.9) is δa=[α]−1· c 0... 0 =[C]· c 0... 0 (15.6.2 ) where cis one arbitrary constant that we get to adjust to make (15.6.1) give the desired left-hand value. Plugging (15.6.2) into (15.6.1) and using the fact that [C] and[α]are inverse matrices of one another, we get c=δa1/C 11and ∆χ2 ν=(δa1)2/C 11 (15.6.3 ) or δa1=±/radicalbig ∆χ2ν/radicalbig C11 (15.6.4 ) At last! A relation between the con fidence interval ±δa1and the formal standarderror σ1≡√C11. Notunreasonably,we findthatthe68percentcon fidence interval is ±σ1, the 95 percent con fidence interval is ±2σ1, etc. 692 Chapter15. ModelingofDataSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).∆χ2as aFunctionofCon fidenceLevelandDegreesofFreedom ν p 123456 68.3% 1.00 2.30 3.53 4.72 5.89 7.04 90% 2.71 4.61 6.25 7.78 9.24 10.6 95.4% 4.00 6.17 8.02 9.70 11.3 12.8 99% 6.63 9.21 11.3 13.3 15.1 16.8 99.73% 9.00 11.8 14.2 16.3 18.2 20.1 99.99% 15.1 18.4 21.1 23.5 25.7 27.8 These considerations hold not just for the individual parameters ai, but also for any linear combination of them: If b≡M/summationdisplay k=1ciai=c·a (15.6.5 ) then the 68 percent con fidence interval on bis δb=±/radicalbig c·[C]·c (15.6.6 ) However,thesesimple,normal-soundingnumericalrelationshipsdo notholdin the case ν> 1[3]. In particular, ∆χ2=1is not the boundary, nor does it project onto the boundary, of a 68.3 percent con fidence region when ν> 1. If you want to calculate not con fidence intervals in one parameter, but con fidence ellipses in two parameters jointly, or ellipsoids in three, or higher, then you must follow the following prescription for implementing Theorems C and D above: •Letνbethenumberof fittedparameterswhosejointcon fidenceregionyou wishtodisplay, ν≤M. Calltheseparametersthe “parametersofinterest. ” •Letpbe the con fidence limit desired, e.g., p=0.68orp=0.95. •Find ∆(i.e., ∆χ2) such that the probability of a chi-square variable with νdegrees of freedom being less than ∆isp. For some useful values of p andν,∆is given in the table. For other values, you can use the routine gammqand a simple root- finding routine (e.g., bisection) to find∆such thatgammq (ν/2,∆/2) = 1 −p. •Take the M×Mcovariance matrix [C]=[α]−1of the chi-square fit. Copy the intersection of the νrows and columns corresponding to the parameters of interest into a ν×νmatrix denoted [Cproj]. •Invertthematrix [Cproj]. (Intheone-dimensionalcasethiswasjusttaking the reciprocal of the element C11.) •Theequationfortheellipticalboundaryofyourdesiredcon fidenceregion in the ν-dimensional subspace of interest is ∆=δa/prime·[Cproj]−1·δa/prime(15.6.7 ) where δa/primeis the ν-dimensional vector of parameters of interest. 15.6ConfidenceLimitsonEstimatedModelParameters 693Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).1 w2V(2)V(1)∆χ2 = 1a2 a1lengthlength1 w1 Figure 15.6.5. Relation of the con fidence region ellipse ∆χ2=1to quantities computed by singular value decomposition. Thevectors V(i)are unit vectors along the principal axes of the con fidence region. The semi-axes have lengths equal to the reciprocal of the singular values wi. If the axes are all scaled by some constant factor α,∆χ2is scaled by the factor α2. If you are confused at this point, you may find it helpful to compare Figure 15.6.4 and the accompanying table, considering the case M=2withν=1and ν=2. You should be able to verify the following statements: (i) The horizontal band between CandC/primecontains 99 percent of the probability distribution, so it is a confidence limit on a2alone at this level of con fidence. (ii) Ditto the band between BandB/primeat the 90 percent con fidence level. (iii) The dashed ellipse, labeledby ∆χ2=2.30,contains68.3percentoftheprobabilitydistribution,so it is a confidence region for a1anda2jointly, at this level of con fidence. Confidence LimitsfromSingular Value Decomposition Whenyouhaveobtainedyour χ2fitbysingularvaluedecomposition( §15.4),the informationaboutthe fit’sformalerrorscomespackagedinasomewhatdifferent,but generally more convenient, form. The columns of the matrix Vare an orthonormal set of Mvectors that are the principal axes of the ∆χ2=constant ellipsoids. We denote the columns as V(1)...V(M). The lengths of those axes are inversely proportionaltothecorrespondingsingularvalues w1...w M;seeFigure15.6.5. The boundaries of the ellipsoids are thus given by ∆χ2=w2 1(V(1)·δa)2+···+w2 M(V(M)·δa)2(15.6.8 ) which is the justi fication for writing equation (15.4.18) above. Keep in mind that it ismucheasier to plot an ellipsoid given a list of its vector principal axes, than given its matrix quadratic form! The formulafor the covariancematrix [C]in terms of the columns V(i)is [C]=M/summationdisplay i=11 w2 iV(i)⊗V(i) (15.6.9 ) or,in components, 694 Chapter15. ModelingofDataSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).Cjk=M/summationdisplay i=11 w2 iVjiVki (15.6.10 ) CITED REFERENCES AND FURTHER READING: Efron,B.1982, TheJackknife,theBootstrap,andOtherResamplingPlans (Philadelphia:S.I.A.M.). [1] Efron, B., and Tibshirani, R. 1986, Statistical Science vol. 1, pp. 54–77. [2] Avni, Y. 1976, Astrophysical Journal , vol. 210, pp. 642–646. [3] Lampton, M., Margon, M., andBowyer, S. 1976, Astrophysical Journal ,vol. 208, pp. 177–190. Brownlee, K.A. 1965, Statistical Theory and Methodology , 2nd ed. (New York: Wiley). Martin, B.R. 1971, Statistics for Physicists (New York: Academic Press). 15.7 Robust Estimation Theconceptof robustness hasbeenmentionedinpassingseveraltimesalready. In§14.1we notedthat the medianwas a morerobustestimator ofcentralvalue than the mean; in §14.6 it was mentionedthat rank correlationis more robust than linear correlation. The concept of outlier points as exceptions to a Gaussian model for experimental error was discussed in §15.1. The term “robust”was coined in statistics by G.E.P. Box in 1953. Various definitions of greater or lesser mathematical rigor are possible for the term, but in general, referringto a statistical estimator, it means “insensitive to small departures fromtheidealizedassumptionsforwhichtheestimatorisoptimized. ”[1,2]Theword “small”can have two different interpretations, both important: either fractionally small departures for all data points, or else fractionally large departures for a small number of data points. It is the latter interpretation, leading to the notion of outlier points, that is generally the most stressful for statistical procedures. Statisticianshavedevelopedvarioussortsofrobuststatisticalestimators. Many, if not most, can be grouped in one of three categories. M-estimates follow from maximum-likelihood arguments very much as equa- tions(15.1.5)and(15.1.7)followedfromequation(15.1.3). M-estimatesareusuallythe most relevant class for model- fitting, that is, estimation of parameters. We therefore consider these estimates in some detail below. L-estimates are“linear combinations of order statistics. ”These are most applicable to estimations of central value and central tendency, though they can occasionally be applied to some problems in estimation of parameters. Two“typical”L-estimates will give you the general idea. They are (i) the median, and (ii)Tukey’s trimean ,d efined as the weighted average of the first, second, and third quartile points in a distribution, with weights 1/4, 1/2, and 1/4, respectively. R-estimates are estimates based on rank tests. For example, the equality or inequality of two distributions can be estimated by the Wilcoxon test of computing the mean rank of one distribution in a combined sample of both distributions. The Kolmogorov-Smirnov statistic (equation 14.3.6) and the Spearman rank-order