c15-3 (779587)
Текст из файла
666Chapter 15.Modeling of Data}CITED REFERENCES AND FURTHER READING:Bevington, P.R. 1969, Data Reduction and Error Analysis for the Physical Sciences (New York:McGraw-Hill), Chapter 6.15.3 Straight-Line Data with Errors in BothCoordinatesIf experimental data are subject to measurement error not only in the yi ’s, but also inthe xi ’s, then the task of fitting a straight-line modely(x) = a + bx(15.3.1)2is considerably harder. It is straightforward to write down the χ merit function for this case,χ2 (a, b) =NX(yi − a − bxi )2i=1σy2 i + b2 σx2 i(15.3.2)where σx i and σy i are, respectively, the x and y standard deviations for the ith point.
Theweighted sum of variances in the denominator of equation (15.3.2) can be understood bothas the variance in the direction of the smallest χ2 between each data point and the line withslope b, and also as the variance of the linear combination yi − a − bxi of two randomvariables xi and yi ,Var(yi − a − bxi ) = Var(yi ) + b2 Var(xi ) = σy2 i + b2 σx2 i ≡ 1/wi(15.3.3)The sum of the square of N random variables, each normalized by its variance, is thusχ2 -distributed.We want to minimize equation (15.3.2) with respect to a and b. Unfortunately, theoccurrence of b in the denominator of equation (15.3.2) makes the resulting equation forthe slope ∂χ2 /∂b = 0 nonlinear.
However, the corresponding condition for the intercept,∂χ2 /∂a = 0, is still linear and yields"#,XXa=wi (yi − bxi )wi(15.3.4)iiwhere the wi ’s are defined by equation (15.3.3). A reasonable strategy, now, is to use themachinery of Chapter 10 (e.g., the routine brent) for minimizing a general one-dimensionalfunction to minimize with respect to b, while using equation (15.3.4) at each stage to ensurethat the minimum with respect to b is also minimized with respect to 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).*chi2=0.0;Calculate χ2 .*q=1.0;if (mwt == 0) {for (i=1;i<=ndata;i++)*chi2 += SQR(y[i]-(*a)-(*b)*x[i]);sigdat=sqrt((*chi2)/(ndata-2));For unweighted data evaluate typ*siga *= sigdat;ical sig using chi2, and ad*sigb *= sigdat;just the standard deviations.} else {for (i=1;i<=ndata;i++)*chi2 += SQR((y[i]-(*a)-(*b)*x[i])/sig[i]);if (ndata>2) *q=gammq(0.5*(ndata-2),0.5*(*chi2));Equation (15.2.12).}66715.3 Straight-Line Data with Errors in Both CoordinatesbsBAr0σaa∆ χ2 = 1Figure 15.3.1.
Standard errors for the parameters a and b. The point B can be found by varying theslope b while simultaneously minimizing the intercept a. This gives the standard error σb , and also thevalue s. The standard error σa can then be found by the geometric relation σa2 = s2 + r 2 .Because of the finite error bars on the xi ’s, the minimum χ2 as a function of b willbe finite, though usually large, when b equals infinity (line of infinite slope). The angleθ ≡ arctan b is thus more suitable as a parametrization of slope than b itself. The value of χ2will then be periodic in θ with period π (not 2π!). If any data points have very small σy ’sbut moderate or large σx ’s, then it is also possible to have a maximum in χ2 near zero slope,θ ≈ 0. In that case, there can sometimes be two χ2 minima, one at positive slope and theother at negative.
Only one of these is the correct global minimum. It is therefore importantto have a good starting guess for b (or θ). Our strategy, implemented below, is to scale theyi ’s so as to have variance equal to the xi ’s, then to do a conventional (as in §15.2) linear fitwith weights derived from the (scaled) sum σy2 i + σx2 i .
This yields a good starting guess forb if the data are even plausibly related to a straight-line model.Finding the standard errors σa and σb on the parameters a and b is more complicated.We will see in §15.6 that, in appropriate circumstances, the standard errors in a and b are therespective projections onto the a and b axes of the “confidence region boundary” where χ2takes on a value one greater than its minimum, ∆χ2 = 1. In the linear case of §15.2, theseprojections follow from the Taylor series expansion1 ∂ 2 χ2∂ 2 χ2∂ 2 χ222∆χ2 ≈(∆a)+(∆b)∆a∆b(15.3.5)+2 ∂a2∂b2∂a∂bBecause of the present nonlinearity in b, however, analytic formulas for the second derivativesare quite unwieldy; more important, the lowest-order term frequently gives a poor approximation to ∆χ2 .
Our strategy is therefore to find the roots of ∆χ2 = 1 numerically, by adjustingthe value of the slope b away from the minimum. In the program below the general root finderzbrent is used. It may occur that there are no roots at all — for example, if all error bars areso large that all the data points are compatible with each other. It is important, therefore, tomake some effort at bracketing a putative root before refining it (cf. §9.1).Because a is minimized at each stage of varying b, successful numerical root-findingleads to a value of ∆a that minimizes χ2 for the value of ∆b that gives ∆χ2 = 1.
This (seeFigure 15.3.1) directly gives the tangent projection of the confidence region onto the b axis,and thus σb . It does not, however, give the tangent projection of the confidence region ontothe a axis. In the figure, we have found the point labeled B; to find σa we need to find theSample 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).σb668Chapter 15.Modeling of Datapoint A. Geometry to the rescue: To the extent that the confidence region is approximatedby an ellipse, then you can prove (see figure) that σa2 = r2 + s2 . The value of s is knownfrom having found the point B. The value of r follows from equations (15.3.2) and (15.3.3)applied at the χ2 minimum (point O in the figure), giving,Xwi(15.3.6)r2 = 1Actually, since b can go through infinity, this whole procedure makes more sense in(a, θ) space than in (a, b) space. That is in fact how the following program works. Sinceit is conventional, however, to return standard errors for a and b, not a and θ, we finallyuse the relationσb = σθ / cos2 θ(15.3.7)We caution that if b and its standard error are both large, so that the confidence region actuallyincludes infinite slope, then the standard error σb is not very meaningful.
The function chixyis normally called only by the routine fitexy. However, if you want, you can yourselfexplore the confidence region by making repeated calls to chixy (whose argument is an angleθ, not a slope b), after a single initializing call to fitexy.A final caution, repeated from §15.0, is that if the goodness-of-fit is not acceptable(returned probability is too small), the standard errors σa and σb are surely not believable.
Indire circumstances, you might try scaling all your x and y error bars by a constant factor untilthe probability is acceptable (0.5, say), to get more plausible values for σa and σb .#include <math.h>#include "nrutil.h"#define POTN 1.571000#define BIG 1.0e30#define PI 3.14159265#define ACC 1.0e-3int nn;float *xx,*yy,*sx,*sy,*ww,aa,offs;Global variables communicate withchixy.void fitexy(float x[], float y[], int ndat, float sigx[], float sigy[],float *a, float *b, float *siga, float *sigb, float *chi2, float *q)Straight-line fit to input data x[1..ndat] and y[1..ndat] with errors in both x and y, the respective standard deviations being the input quantities sigx[1..ndat] and sigy[1..ndat].Output quantities are a and b such that y = a + bx minimizes χ2 , whose value is returnedas chi2. The χ2 probability is returned as q, a small value indicating a poor fit (sometimesindicating underestimated errors).
Standard errors on a and b are returned as siga and sigb.These are not meaningful if either (i) the fit is poor, or (ii) b is so large that the data areconsistent with a vertical (infinite b) line. If siga and sigb are returned as BIG, then the dataare consistent with all values of b.{void avevar(float data[], unsigned long n, float *ave, float *var);float brent(float ax, float bx, float cx,float (*f)(float), float tol, float *xmin);float chixy(float bang);void fit(float x[], float y[], int ndata, float sig[], int mwt,float *a, float *b, float *siga, float *sigb, float *chi2, float *q);float gammq(float a, float x);void mnbrak(float *ax, float *bx, float *cx, float *fa, float *fb,float *fc, float (*func)(float));float zbrent(float (*func)(float), float x1, float x2, float tol);int j;float swap,amx,amn,varx,vary,ang[7],ch[7],scale,bmn,bmx,d1,d2,r2,dum1,dum2,dum3,dum4,dum5;xx=vector(1,ndat);yy=vector(1,ndat);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.
Характеристики
Тип файла PDF
PDF-формат наиболее широко используется для просмотра любого типа файлов на любом устройстве. В него можно сохранить документ, таблицы, презентацию, текст, чертежи, вычисления, графики и всё остальное, что можно показать на экране любого устройства. Именно его лучше всего использовать для печати.
Например, если Вам нужно распечатать чертёж из автокада, Вы сохраните чертёж на флешку, но будет ли автокад в пункте печати? А если будет, то нужная версия с нужными библиотеками? Именно для этого и нужен формат PDF - в нём точно будет показано верно вне зависимости от того, в какой программе создали PDF-файл и есть ли нужная программа для его просмотра.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
