c15-3 (779587), страница 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).i15.3 Straight-Line Data with Errors in Both Coordinates669}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).sx=vector(1,ndat);sy=vector(1,ndat);ww=vector(1,ndat);avevar(x,ndat,&dum1,&varx);Find the x and y variances, and scaleavevar(y,ndat,&dum1,&vary);the data into the global variablesscale=sqrt(varx/vary);for communication with the funcnn=ndat;tion chixy.for (j=1;j<=ndat;j++) {xx[j]=x[j];yy[j]=y[j]*scale;sx[j]=sigx[j];sy[j]=sigy[j]*scale;ww[j]=sqrt(SQR(sx[j])+SQR(sy[j]));Use both x and y weights in first}trial fit.fit(xx,yy,nn,ww,1,&dum1,b,&dum2,&dum3,&dum4,&dum5);Trial fit for b.offs=ang[1]=0.0;Construct several angles for referang[2]=atan(*b);ence points, and make b an anang[4]=0.0;gle.ang[5]=ang[2];ang[6]=POTN;for (j=4;j<=6;j++) ch[j]=chixy(ang[j]);mnbrak(&ang[1],&ang[2],&ang[3],&ch[1],&ch[2],&ch[3],chixy);Bracket the χ2 minimum and then locate it with brent.*chi2=brent(ang[1],ang[2],ang[3],chixy,ACC,b);*chi2=chixy(*b);*a=aa;*q=gammq(0.5*(nn-2),*chi2*0.5);Compute χ2 probability.for (r2=0.0,j=1;j<=nn;j++) r2 += ww[j];Save the inverse sum of weights atr2=1.0/r2;the minimum.bmx=BIG;Now, find standard errors for b aspoints where ∆χ2 = 1.bmn=BIG;offs=(*chi2)+1.0;for (j=1;j<=6;j++) {Go through saved values to bracketif (ch[j] > offs) {the desired roots.
Note periodd1=fabs(ang[j]-(*b));icity in slope angles.while (d1 >= PI) d1 -= PI;d2=PI-d1;if (ang[j] < *b) {swap=d1;d1=d2;d2=swap;}if (d1 < bmx) bmx=d1;if (d2 < bmn) bmn=d2;}}if (bmx < BIG) {Call zbrent to find the roots.bmx=zbrent(chixy,*b,*b+bmx,ACC)-(*b);amx=aa-(*a);bmn=zbrent(chixy,*b,*b-bmn,ACC)-(*b);amn=aa-(*a);*sigb=sqrt(0.5*(bmx*bmx+bmn*bmn))/(scale*SQR(cos(*b)));*siga=sqrt(0.5*(amx*amx+amn*amn)+r2)/scale;Error in a has additional piece} else (*sigb)=(*siga)=BIG;r2.*a /= scale;Unscale the answers.*b=tan(*b)/scale;free_vector(ww,1,ndat);free_vector(sy,1,ndat);free_vector(sx,1,ndat);free_vector(yy,1,ndat);free_vector(xx,1,ndat);670Chapter 15.Modeling of Data#include <math.h>#include "nrutil.h"#define BIG 1.0e30extern int nn;extern float *xx,*yy,*sx,*sy,*ww,aa,offs;b=tan(bang);for (j=1;j<=nn;j++) {ww[j] = SQR(b*sx[j])+SQR(sy[j]);sumw += (ww[j] = (ww[j] < 1.0/BIG ? BIG : 1.0/ww[j]));avex += ww[j]*xx[j];avey += ww[j]*yy[j];}avex /= sumw;avey /= sumw;aa=avey-b*avex;for (ans = -offs,j=1;j<=nn;j++)ans += ww[j]*SQR(yy[j]-aa-b*xx[j]);return ans;}Be aware that the literature on the seemingly straightforward subject of this sectionis generally confusing and sometimes plain wrong.
Deming’s [1] early treatment is sound,but its reliance on Taylor expansions gives inaccurate error estimates. References [2-4] arereliable, more recent, general treatments with critiques of earlier work. York [5] and Reed [6]usefully discuss the simple case of a straight line as treated here, but the latter paper hassome errors, corrected in [7].
All this commotion has attracted the Bayesians [8-10] , whohave still different points of view.CITED REFERENCES AND FURTHER READING:Deming, W.E. 1943, Statistical Adjustment of Data (New York: Wiley), reprinted 1964 (New York:Dover). [1]Jefferys, W.H. 1980, Astronomical Journal, vol. 85, pp. 177–181; see also vol.
95, p. 1299(1988). [2]Jefferys, W.H. 1981, Astronomical Journal, vol. 86, pp. 149–155; see also vol. 95, p. 1300(1988). [3]Lybanon, M. 1984, American Journal of Physics, vol. 52, pp. 22–26. [4]York, D. 1966, Canadian Journal of Physics, vol. 44, pp. 1079–1086. [5]Reed, B.C. 1989, American Journal of Physics, vol. 57, pp. 642–646; see also vol. 58, p. 189,and vol. 58, p. 1209. [6]Reed, B.C. 1992, American Journal of Physics, vol. 60, pp. 59–62. [7]Zellner, A. 1971, An Introduction to Bayesian Inference in Econometrics (New York: Wiley);reprinted 1987 (Malabar, FL: R. E. Krieger Pub. Co.). [8]Gull, S.F. 1989, in Maximum Entropy and Bayesian Methods, J.
Skilling, ed. (Boston: Kluwer). [9]Jaynes, E.T. 1991, in Maximum-Entropy and Bayesian Methods, Proc. 10th Int. Workshop,W.T. Grandy, Jr., and L.H. Schick, eds. (Boston: Kluwer). [10]Macdonald, J.R., and Thompson, W.J. 1992, American Journal of Physics, vol. 60, pp. 66–73.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).float chixy(float bang)Captive function of fitexy, returns the value of (χ2 − offs) for the slope b=tan(bang).Scaled data and offs are communicated via the global variables.{int j;float ans,avex=0.0,avey=0.0,sumw=0.0,b;15.4 General Linear Least Squares67115.4 General Linear Least Squaresy(x) = a1 + a2 x + a3 x2 + · · · + aM xM −1(15.4.1)is a polynomial of degree M − 1. Or, the functions could be sines and cosines, inwhich case their general linear combination is a harmonic series.The general form of this kind of model isy(x) =MXak Xk (x)(15.4.2)k=1where X1 (x), .
. . , XM (x) are arbitrary fixed functions of x, called the basisfunctions.Note that the functions Xk (x) can be wildly nonlinear functions of x. In thisdiscussion “linear” refers only to the model’s dependence on its parameters ak .For these linear models we generalize the discussion of the previous sectionby defining a merit function"#2PMNXyi − k=1 ak Xk (xi )2(15.4.3)χ =σii=1As before, σi is the measurement error (standard deviation) of the ith data point,presumed to be known. If the measurement errors are not known, they may all (asdiscussed at the end of §15.1) be set to the constant value σ = 1.Once again, we will pick as best parameters those that minimize χ2 . There areseveral different techniques available for finding this minimum.
Two are particularlyuseful, and we will discuss both in this section. To introduce them and elucidatetheir relationship, we need some notation.Let A be a matrix whose N × M components are constructed from the Mbasis functions evaluated at the N abscissas xi , and from the N measurement errorsσi , by the prescriptionXj (xi )(15.4.4)Aij =σiThe matrix A is called the design matrix of the fitting problem.
Notice that in generalA has more rows than columns, N ≥M , since there must be more data points thanmodel parameters to be solved for. (You can fit a straight line to two points, but not avery meaningful quintic!) The design matrix is shown schematically in Figure 15.4.1.Also define a vector b of length N byyi(15.4.5)bi =σiand denote the M vector whose components are the parameters to be fitted,a1 , .
. . , aM , by a.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).An immediate generalization of §15.2 is to fit a set of data points (xi , yi ) to amodel that is not just a linear combination of 1 and x (namely a + bx), but rather alinear combination of any M specified functions of x. For example, the functionscould be 1, x, x2, . .
. , xM −1, in which case their general linear combination,.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
