c15-7 (779591), страница 2
Текст из файла (страница 2)
Further reproduction, or any copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).Notice that the ψ function occurs as a weighting function in the generalizednormal equations (15.7.5). For normally distributed errors, equation (15.7.6) saysthat 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 pointsget the same relative weight, with only the sign information used.
Finally, whenthe tails are even larger, (15.7.10) says the ψ increases with deviation, then startsdecreasing, so that very deviant points — the true outliers — are not counted at allin the estimation of the parameters.This general idea, that the weight given individual points should first increasewith deviation, then decrease, motivates some additional prescriptions for ψ whichdo not especially correspond to standard, textbook probability distributions. Twoexamples areAndrew’s sinesin(z/c)|z| < cπψ(z) =(15.7.11)0|z| > cπ15.7 Robust Estimation703Fitting a Line by Minimizing Absolute DeviationOccasionally there is a special case that happens to be much easier than issuggested by the general strategy outlined above.
The case of equations (15.7.7)–(15.7.8), when the model is a simple straight liney(x; a, b) = a + bx(15.7.13)and where the weights σi are all equal, happens to be such a case. The problem isprecisely the robust version of the problem posed in equation (15.2.1) above, namelyfit a straight line through a set of data points. The merit function to be minimized isNX|yi − a − bxi |(15.7.14)i=1rather than the χ2 given by equation (15.2.2).The key simplification is based on the following fact: The median cM of a setof numbers ci is also that value which minimizes the sum of the absolute deviationsX|ci − cM |i(Proof: Differentiate the above expression with respect to cM and set it to zero.)It follows that, for fixed b, the value of a that minimizes (15.7.14) isa = median {yi − bxi }(15.7.15)Equation (15.7.5) for the parameter b is0=NXxi sgn(yi − a − bxi )(15.7.16)i=1(where sgn(0) is to be interpreted as zero). If we replace a in this equation by theimplied function a(b) of (15.7.15), then we are left with an equation in a singlevariable 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 thediscontinuities in equation 15.7.16.)Here is a routine that does all this.
It calls select (§8.5) to find the median.The bracketing and bisection are built in to the routine, as is the χ2 solution thatgenerates the initial guesses for a and b. Notice that the evaluation of the right-handside of (15.7.16) occurs in the function rofunc, with communication via global(top-level) variables.Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-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 copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).algorithm exemplified in amoeba §10.4 or amebsa in §10.9. Those algorithms makeno assumptions about continuity; they just ooze downhill and will work for virtuallyany sane choice of the function ρ.It is very much to your (financial) advantage to find good starting values,however.
Often this is done by first fitting the model by the standard χ2 (nonrobust)techniques, e.g., as described in §15.4 or §15.5. The fitted parameters thus obtainedare then used as starting values in amoeba, now using the robust choice of ρ andminimizing the expression (15.7.3).704Chapter 15.Modeling of Data#include <math.h>#include "nrutil.h"int ndatat;float *xt,*yt,aa,abdevt;ndatat=ndata;xt=x;yt=y;for (j=1;j<=ndata;j++) {As a first guess for a and b, we will find the leastsx += x[j];squares fitting line.sy += y[j];sxy += x[j]*y[j];sxx += x[j]*x[j];}del=ndata*sxx-sx*sx;aa=(sxx*sy-sx*sxy)/del;Least-squares solutions.bb=(ndata*sxy-sx*sy)/del;for (j=1;j<=ndata;j++)chisq += (temp=y[j]-(aa+bb*x[j]),temp*temp);sigb=sqrt(chisq/del);The standard deviation will give some idea of howb1=bb;big an iteration step to take.f1=rofunc(b1);b2=bb+SIGN(3.0*sigb,f1);Guess bracket as 3-σ away, in the downhill direction known from f1.f2=rofunc(b2);if (b2 == b1) {*a=aa;*b=bb;*abdev=abdevt/ndata;return;}while (f1*f2 > 0.0) {Bracketing.bb=b2+1.6*(b2-b1);b1=b2;f1=f2;b2=bb;f2=rofunc(b2);}sigb=0.01*sigb;Refine until error a negligible number of standardwhile (fabs(b2-b1) > sigb) {deviations.bb=b1+0.5*(b2-b1);Bisection.if (bb == b1 || bb == b2) break;f=rofunc(bb);if (f*f1 >= 0.0) {f1=f;b1=bb;} else {f2=f;b2=bb;}}*a=aa;*b=bb;Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-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 copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).void medfit(float x[], float y[], int ndata, float *a, float *b, float *abdev)Fits y = a + bx by the criterion of least absolute deviations. The arrays x[1..ndata] andy[1..ndata] are the input experimental points.
The fitted parameters a and b are output,along with abdev, which is the mean absolute deviation (in y) of the experimental points fromthe fitted line. This routine uses the routine rofunc, with communication via global variables.{float rofunc(float b);int j;float bb,b1,b2,del,f,f1,f2,sigb,temp;float sx=0.0,sy=0.0,sxy=0.0,sxx=0.0,chisq=0.0;15.7 Robust Estimation705*abdev=abdevt/ndata;}extern int ndatat;Defined in medfit.extern float *xt,*yt,aa,abdevt;float rofunc(float b)Evaluates the right-hand side of equation (15.7.16) for a given value of b. Communication withthe routine medfit is through global variables.{float select(unsigned long k, unsigned long n, float arr[]);int j;float *arr,d,sum=0.0;arr=vector(1,ndatat);for (j=1;j<=ndatat;j++) arr[j]=yt[j]-b*xt[j];if (ndatat & 1) {aa=select((ndatat+1)>>1,ndatat,arr);}else {j=ndatat >> 1;aa=0.5*(select(j,ndatat,arr)+select(j+1,ndatat,arr));}abdevt=0.0;for (j=1;j<=ndatat;j++) {d=yt[j]-(b*xt[j]+aa);abdevt += fabs(d);if (yt[j] != 0.0) d /= fabs(yt[j]);if (fabs(d) > EPS) sum += (d >= 0.0 ? xt[j] : -xt[j]);}free_vector(arr,1,ndatat);return sum;}Other Robust TechniquesSometimes you may have a priori knowledge about the probable values andprobable uncertainties of some parameters that you are trying to estimate from a dataset.
In such cases you may want to perform a fit that takes this advance informationproperly into account, neither completely freezing a parameter at a predeterminedvalue (as in lfit §15.4) nor completely leaving it to be determined by the data set.The formalism for doing this is called “use of a priori covariances.”A related problem occurs in signal processing and control theory, where it issometimes desired to “track” (i.e., maintain an estimate of) a time-varying signal inthe presence of noise. If the signal is known to be characterized by some numberof parameters that vary only slowly, then the formalism of Kalman filtering tellshow the incoming, raw measurements of the signal should be processed to producebest parameter estimates as a function of time.
For example, if the signal is afrequency-modulated sine wave, then the slowly varying parameter might be theinstantaneous frequency. The Kalman filter for this case is called a phase-lockedloop and is implemented in the circuitry of good radio receivers [3,4] .Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-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 copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited.
To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).#include <math.h>#include "nrutil.h"#define EPS 1.0e-7Modeling of DataChapter 15.706CITED 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: AcademicPress). [2]Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-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 copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).Jazwinski, A.
H. 1970, Stochastic Processes and Filtering Theory (New York: AcademicPress). [4]Bryson, A. E., and Ho, Y.C. 1969, Applied Optimal Control (Waltham, MA: Ginn). [3].
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
