f15-7

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 , defined 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 15.7 RobustEstimation 695Sample 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).narrow central peak tail of outliers least squares fit robust straight-line fit(a) (b) Figure 15.7.1. Examples where robust statistical methods are desirable: (a) A one-dimensional distributionwithatailofoutliers; statistical fluctuationsintheseoutlierscanpreventaccuratedetermination oftheposition ofthecentral peak. (b)Adistribution intwodimensions fittedtoastraight line; non-robust techniques such as least-squares fitting can have undesired sensitivity to outlying points. correlation coef ficient (14.6.1) are R-estimates in essence, if not always by formal definition. Someotherkindsofrobusttechniques,comingfromthe fieldsofoptimalcontrol andfiltering rather than from the field of mathematical statistics, are mentioned at theendofthissection. Someexampleswhererobuststatisticalmethodsaredesirable are shown in Figure 15.7.1. Estimation of Parameters by LocalM-Estimates Suppose we know that our measurement errors are not normally distributed. Then,inderivingamaximum-likelihoodformulafortheestimatedparameters aina model y(x;a), we would write instead of equation (15.1.3) P=N/productdisplay i=1{exp [−ρ(yi,y{xi;a})] ∆y} (15.7.1 ) 696 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).wherethe function ρis the negativelogarithmof the probabilitydensity. Takingthe logarithm of (15.7.1) analogously with (15.1.4), we find that we want to minimize the expression N/summationdisplay i=1ρ(yi,y{xi;a})( 15.7.2 ) Very often, it is the case that the function ρdepends not independently on its twoarguments,measured yiandpredicted y(xi),butonlyontheirdifference,atleast ifscaledbysomeweightfactors σiwhichweareabletoassigntoeachpoint. Inthis casetheM-estimateissaidtobe local,andwecanreplace(15.7.2)bytheprescription minimizeover aN/summationdisplay i=1ρ/parenleftbiggyi−y(xi;a) σi/parenrightbigg (15.7.3 ) where the function ρ(z)is a functionof a single variable z≡[yi−y(xi)]/σ i. If we now de fine the derivative of ρ(z)to be a function ψ(z), ψ(z)≡dρ(z) dz(15.7.4 ) then the generalization of (15.1.7) to the case of a general M-estimate is 0=N/summationdisplay i=11 σiψ/parenleftbiggyi−y(xi) σi/parenrightbigg/parenleftbigg∂y(xi;a) ∂a k/parenrightbigg k=1,...,M (15.7.5 ) If you compare (15.7.3) to (15.1.3), and (15.7.5) to (15.1.7), you see at once that the specialization for normally distributed errors is ρ(z)=1 2z2ψ(z)=z(normal) (15.7.6 ) If the errors are distributed as a doubleortwo-sided exponential , namely Prob {yi−y(xi)}∼ exp/parenleftbigg −/vextendsingle/vextendsingle/vextendsingle/vextendsingley i−y(xi) σi/vextendsingle/vextendsingle/vextendsingle/vextendsingle/parenrightbigg (15.7.7 ) then, by contrast, ρ(x)=|z| ψ(z)=sgn(z)(doubleexponential) (15.7.8 ) Comparing to equation (15.7.3), we see that in this case the maximum likelihood estimator is obtained by minimizing the mean absolute deviation , rather than the mean square deviation. Here the tails of the distribution, although exponentially decreasing, are asymptoticallymuch larger than any correspondingGaussian. A distribution with even more extensive —therefore sometimes even more realistic—tails is the CauchyorLorentzian distribution, Prob {y i−y(xi)}∼1 1+1 2/parenleftbiggyi−y(xi) σi/parenrightbigg2(15.7.9 ) 15.7 RobustEstimation 697Sample 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).This implies ρ(z)=l o g/parenleftbigg 1+1 2z2/parenrightbigg ψ(z)=z 1+1 2z2(Lorentzian) (15.7.10 ) Notice that the ψfunction occurs as a weighting function in the generalized normal equations (15.7.5). For normally distributed errors, equation (15.7.6) says that the more deviant the points, the greater the weight. By contrast, when tails aresomewhat more prominent, as in (15.7.7), then (15.7.8) says that all deviant points get the same relative weight, with only the sign information used. Finally, when the tails are even larger, (15.7.10) says the ψincreases with deviation, then starts decreasing , so that very deviant points —the true outliers —are not counted at all in the estimation of the parameters. This general idea, that the weight given individual points should first increase with deviation, then decrease, motivates some additional prescriptions for ψwhich do not especially correspond to standard, textbook probability distributions. Twoexamples are Andrew’s sine ψ(z)=/braceleftbigg sin(z/c) 0|z|<c π |z|>c π(15.7.11 ) Ifthemeasurementerrorshappentobenormalafterall,withstandarddeviations σ i, then it can be shown that the optimal value for the constant cisc=2.1. Tukey’s biweight ψ(z)=/braceleftbigg z(1−z2/c2)2 0|z|<c |z|>c(15.7.12 ) where the optimal value of cfor normal errors is c=6.0. NumericalCalculation of M-Estimates Tofit a model by means of an M-estimate, you first decide which M-estimate you want, that is, which matching pair ρ,ψyou want to use. We rather like (15.7.8) or (15.7.10). You then have to make an unpleasant choice between two fairly dif ficult problems. Either find the solution of the nonlinear set of Mequations (15.7.5), or else minimize the single function in Mvariables (15.7.3). Notice that the function (15.7.8) has a discontinuous ψ, and a discontinuous derivative for ρ. Such discontinuities frequently wreak havoc on both general nonlinear equation solvers and general function minimizing routines. You might now think of rejecting (15.7.8) in favor of (15.7.10), which is smoother. However,youwillfindthatthelatterchoiceisalsobadnewsformanygeneralequationsolving or minimization routines: small changes in the fitted parameters can drive ψ(z) off its peak into one or the other of its asymptotically small regimes. Therefore,differenttermsin theequationspringintooroutofaction(almostas badas analytic discontinuities). Don’t despair. If your computer budget (or, for personal computers, patience) is up to it, this is an excellent application for the downhill simplex minimization 698 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).algorithmexempli fiedin amoeba §10.4or amebsain§10.9. Thosealgorithmsmake noassumptionsaboutcontinuity;theyjustoozedownhillandwillworkforvirtuallyany sane choice of the function ρ. It is very much to your ( financial) advantage to find good starting values, however. Oftenthis is doneby firstfitting themodelbythe standard χ 2(nonrobust) techniques,e.g., as describedin §15.4or §15.5. The fitted parametersthus obtained are then used as starting values in amoeba, now using the robust choice of ρand minimizing the expression (15.7.3). Fittinga Line by MinimizingAbsolute Deviation Occasionally there is a special case that happens to be much easier than is suggested by the general strategy outlined above. The case of equations (15.7.7) – (15.7.8), when the model is a simple straight line y(x;a, b)=a+bx (15.7.13 ) and where the weights σiare all equal, happens to be such a case. The problem is preciselytherobustversionoftheproblemposedinequation(15.2.1)above,namely fit a straightlinethrougha set ofdatapoints. Themerit functionto beminimizedis N/summationdisplay i=1|yi−a−bx i| (15.7.14 ) rather than the χ2given by equation (15.2.2). The key simpli fication is based on the following fact: The median cMof a set ofnumbers ciis also thatvaluewhichminimizesthesumofthe absolutedeviations /summationdisplay i|ci−cM| (Proof: Differentiatethe aboveexpressionwith respect to cMandset it to zero.) It follows that, for fixedb, the value of athat minimizes (15.7.14)is a=median {yi−bx i} (15.7.15 ) Equation (15.7.5) for the parameter bis 0=N/summationdisplay i=1xisgn(yi−a−bx i)( 15.7.16 ) (where sgn (0)is to be interpreted as zero). If we replace ain this equation by the implied function a(b)of (15.7.15), then we are left with an equation in a single variable which can be solved by bracketing and bisection, as described in §9.1. (In fact, it is dangerous to use any fancier method of root- finding, because of the discontinuities in equation 15.7.16.) Here is a routine that does all this. It calls select(§8.5) tofind the median. The bracketing and bisection are built in to the routine, as is the χ2solution that generatestheinitial guesses for aandb. Noticethat theevaluationofthe right-hand sideof(15.7.16)occursin thefunction rofunc,withcommunicationviaa common block. To save memory, you could generate your data arrays directly into that common block, deleting them from this routine ’s calling sequence. 15.7 RobustEstimation 699Sample 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).SUBROUTINE medfit(x,y,ndata,a,b,abdev) INTEGER ndata,NMAX,ndatat PARAMETER (NMAX=1000) REAL a,abdev,b,x(ndata),y(ndata), * arr(NMAX),xt(NMAX),yt(NMAX),aa,abdevt COMMON /arrays/ xt,yt,arr,aa,abdevt,ndatat C USES rofunc Fits y=a+bxby the criterion of least absolute deviations. The arrays x(1:ndata) andy(1:ndata) are the input experimental points. The fitted parameters aandbare output, along with abdev , which is the mean absolute deviation (in y) of the experimental points from the fitted line. This routine uses the routine rofunc , with communication via a common block. INTEGER j REAL b1,b2,bb,chisq,del,f,f1,f2,sigb,sx,sxx,sxy,sy,rofunc sx=0.sy=0. sxy=0. sxx=0.do 11j=1,ndata As a first guess for aand b, we will find the least- squares fitting line. xt(j)=x(j) yt(j)=y(j) sx=sx+x(j)sy=sy+y(j) sxy=sxy+x(j)*y(j) sxx=sxx+x(j)**2 enddo 11 ndatat=ndata del=ndata*sxx-sx**2 aa=(sxx*sy-sx*sxy)/del Least-squares solutions. bb=(ndata*sxy-sx*sy)/del chisq=0. do12j=1,ndata chisq=chisq+(y(j)-(aa+bb*x(j)))**2 enddo 12 sigb=sqrt(chisq/del) The standard deviation will give some idea of how big an iteration step to take. b1=bb f1=rofunc(b1)if(sigb.gt.0.)then b2=bb+sign(3.*sigb,f1) Guess bracket as 3- σaway, in the downhill direction known from f1. f2=rofunc(b2) if(b2.eq.b1)then a=aa b=bbabdev=abdevt/ndatareturn endif 1 if(f1*f2.gt.0.)then Bracketing. bb=b2+1.6*(b2-b1) b1=b2 f1=f2b2=bbf2=rofunc(b2) goto 1 endifsigb=0.01*sigb Refine until error a negligible number of standard de- viations. 2 if(abs(b2-b1).gt.sigb)then bb=b1+0.5*(b2-b1) Bisection. if(bb.eq.b1.or.bb.eq.b2)goto 3f=rofunc(bb) if(f*f1.ge.0.)then f1=fb1=bb else f2=f b2=bb endif goto 2 endif 700 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).endif 3 a=aa b=bb abdev=abdevt/ndatareturn END FUNCTION rofunc(b) INTEGER NMAXREAL rofunc,b,EPSPARAMETER (NMAX=1000,EPS=1.e-7) C USES select Evaluates the right-hand side of equation (15.7.16) for a given value of b. Communication with the program medfit is through a common block. INTEGER j,ndata REAL aa,abdev,d,sum,arr(NMAX),x(NMAX),y(NMAX),selectCOMMON /arrays/ x,y,arr,aa,abdev,ndatado 11j=1,ndata arr(j)=y(j)-b*x(j) enddo 11 if (mod(ndata,2).eq.0) then j=ndata/2 aa=0.5*(select(j,ndata,arr)+select(j+1,ndata,arr)) else aa=select((ndata+1)/2,ndata,arr) endif sum=0.abdev=0.do 12j=1,ndata d=y(j)-(b*x(j)+aa) abdev=abdev+abs(d)if (y(j).ne.0.) d=d/abs(y(j)) if (abs(d).gt.EPS) sum=sum+x(j)*sign(1.0,d) enddo 12 rofunc=sumreturn END OtherRobust Techniques Sometimes you may have a prioriknowledge about the probable values and probable uncertainties of some parameters that you are trying to estimate from a data set. In suchcasesyou maywanttoperforma fitthattakes thisadvance informationproperly intoaccount, neither completely freezing a parameter at a predetermined value (as in lfit §15.4) nor completely leaving it to be determined by the data set. The formalism for doing this is called“use ofa prioricovariances. ” A related problem occurs in signal processing and control theory, where it is sometimes desired to “track”(i.e., maintain an estimate of) a time-varying signal in the presence of noise. Ifthesignalisknown tobecharacterized bysomenumber ofparameters thatvary onlyslowly, then the formalism of Kalman filtering tells how the incoming, raw measurements of the signal should be processed to produce best parameter estimates as a function of time. For example, ifthe signal is a frequency-modulated sine wave,then the slowly varying parameter might bethe instantaneous frequency. The Kalman filterforthis caseiscalled a phase-locked loopand is implemented in the circuitry of good radio receivers [3,4]. CITED REFERENCES AND FURTHER READING: Huber, P.J. 1981, Robust Statistics (New York: Wiley). [1] Launer, R.L., and Wilkinson, G.N. (eds.) 1979, Robustness in Statistics (New York: Academic Press). [2] Bryson, A. E., and Ho, Y.C. 1969, Applied Optimal Control (Waltham, MA: Ginn). [3] Jazwinski, A. H. 1970, Stochastic Processes and Filtering Theory (New York: Academic Press). [4]