c15-4 (779588)
Текст из файла
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,672Chapter 15.Modeling of DataX1( )X2 ( )...XM ( )x1X1(x1)σ1X2 (x1)σ1...XM (x1)σ1x2X1(x 2)σ2X2 (x 2)σ2...XM (x 2)σ2............xNX1(xN)σN.........X2 (xN)σN...XM (xN)σNFigure 15.4.1. Design matrix for the least-squares fit of a linear combination of M basis functions to Ndata points.
The matrix elements involve the basis functions evaluated at the values of the independentvariable at which measurements are made, and the standard deviations of the measured dependent variable.The measured values of the dependent variable do not enter the design matrix.Solution by Use of the Normal EquationsThe minimum of (15.4.3) occurs where the derivative of χ2 with respect to allM parameters ak vanishes.
Specializing equation (15.1.7) to the case of the model(15.4.2), this condition yields the M equationsNMXX1 0=aj Xj (xi ) Xk (xi )k = 1, . . . , M(15.4.6)yi −σi2i=1j=1Interchanging the order of summations, we can write (15.4.6) as the matrix equationMXαkj aj = βk(15.4.7)j=1whereαkj =NXXj (xi )Xk (xi )i=1σi2or equivalently[α] = AT · A(15.4.8)an M × M matrix, andβk =NXyi Xk (xi )i=1σi2or equivalently[β] = AT · b(15.4.9)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).data pointsbasis functions15.4 General Linear Least Squares673The inverse matrix Cjk ≡ [α]−1jk is closely related to the probable (or, moreprecisely, standard) uncertainties of the estimated parameters a.
To estimate theseuncertainties, consider that"N#MMXXX yi Xk (xi )(15.4.11)[α]−1Cjkaj =jk βk =σi2i=1k=1k=1and that the variance associated with the estimate aj can be found as in (15.2.7) from2σ (aj ) =NXσi2i=1∂aj∂yi2(15.4.12)Note that αjk is independent of yi , so thatX∂aj=Cjk Xk (xi )/σi2∂yiM(15.4.13)k=1Consequently, we find that2σ (aj ) =MM XXk=1 l=1Cjk Cjl"N#X Xk (xi )Xl (xi )i=1σi2(15.4.14)The final term in brackets is just the matrix [α]. Since this is the matrix inverseof [C], (15.4.14) reduces immediately toσ 2 (aj ) = Cjj(15.4.15)In other words, the diagonal elements of [C] are the variances (squareduncertainties) of the fitted parameters a. It should not surprise you to learn that theoff-diagonal elements Cjk are the covariances between aj and ak (cf. 15.2.10); butwe shall defer discussion of these to §15.6.We will now give a routine that implements the above formulas for the generallinear least-squares problem, by the method of normal equations.
Since we wish tocompute not only the solution vector a but also the covariance matrix [C], it is mostconvenient to use Gauss-Jordan elimination (routine gaussj of §2.1) to perform thelinear algebra. The operation count, in this application, is no larger than that for LUdecomposition. If you have no need for the covariance matrix, however, you cansave a factor of 3 on the linear algebra by switching to LU decomposition, withoutSample 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).a vector of length M .The equations (15.4.6) or (15.4.7) are called the normal equations of the leastsquares problem.
They can be solved for the vector of parameters a by the standardmethods of Chapter 2, notably LU decomposition and backsubstitution, Choleksydecomposition, or Gauss-Jordan elimination. In matrix form, the normal equationscan be written as either[α] · a = [β]or asAT · A · a = AT · b(15.4.10)674Chapter 15.Modeling of Data#include "nrutil.h"void lfit(float x[], float y[], float sig[], int ndat, float a[], int ia[],int ma, float **covar, float *chisq, void (*funcs)(float, float [], int))Given a set of data points x[1..ndat], y[1..ndat] with individual standard deviationssig[1..ndat], use χ2 minimization to fitPfor some or all of the coefficients a[1..ma] ofa function that depends linearly on a, y = i ai × afunci (x).
The input array ia[1..ma]indicates by nonzero entries those components of a that should be fitted for, and by zero entriesthose components that should be held fixed at their input values. The program returns valuesfor a[1..ma], χ2 = chisq, and the covariance matrix covar[1..ma][1..ma]. (Parametersheld fixed will return zero covariances.) The user supplies a routine funcs(x,afunc,ma) thatreturns the ma basis functions evaluated at x = x in the array afunc[1..ma].{void covsrt(float **covar, int ma, int ia[], int mfit);void gaussj(float **a, int n, float **b, int m);int i,j,k,l,m,mfit=0;float ym,wt,sum,sig2i,**beta,*afunc;beta=matrix(1,ma,1,1);afunc=vector(1,ma);for (j=1;j<=ma;j++)if (ia[j]) mfit++;if (mfit == 0) nrerror("lfit: no parameters to be fitted");for (j=1;j<=mfit;j++) {Initialize the (symmetric) matrix.for (k=1;k<=mfit;k++) covar[j][k]=0.0;beta[j][1]=0.0;}for (i=1;i<=ndat;i++) {Loop over data to accumulate coefficients ofthe normal equations.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}
