c15-4 (779588), страница 3
Текст из файла (страница 3)
A plausible answer to the question “howsmall is small?” is to edit in this fashion all singular values whose ratio to thelargest singularvalue is less than N times the machine precision . (You might√argue for N , or a constant, instead of N as the multiple; that starts getting intohardware-dependent questions.)There is another reason for editing even additional singular values, ones largeenough that roundoff error is not a question.
Singular value decomposition allowsyou to identify linear combinations of variables that just happen not to contributemuch to reducing the χ2 of your data set. Editing these can sometimes reduce theprobable error on your coefficients quite significantly, while increasing the minimumχ2 only negligibly. We will learn more about identifying and treating such casesin §15.6. In the following routine, the point at which this kind of editing wouldoccur is indicated.Generally speaking, we recommend that you always use SVD techniques insteadof using the normal equations. SVD’s only significant disadvantage is that it requiresSample 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).of U with the weighted data vector (15.4.5). Though it is beyond our scope to provehere, it turns out that the standard (loosely, “probable”) errors in the fitted parametersare also linear combinations of the columns of V. In fact, equation (15.4.17) canbe written in a form displaying these errors as#"M X U(i) · b 11V(i) ±V(1) ± · · · ±V(M )(15.4.18)a=wiw1wM678Chapter 15.Modeling of Data#include "nrutil.h"#define TOL 1.0e-5void svdfit(float x[], float y[], float sig[], int ndata, float a[], int ma,float **u, float **v, float w[], float *chisq,void (*funcs)(float, float [], int))Given a set of data points x[1..ndata],y[1..ndata] with individual standard deviations2sig[1..ndata], useP χ minimization to determine the coefficients a[1..ma] of the fitting function y =ai i × afunci (x).
Here we solve the fitting equations using singularvalue decomposition of the ndata by ma matrix, as in §2.6. Arrays u[1..ndata][1..ma],v[1..ma][1..ma], and w[1..ma] provide workspace on input; on output they define thesingular value decomposition, and can be used to obtain the covariance matrix. The program returns values for the ma fit parameters a, and χ2 , chisq. The user supplies a routinefuncs(x,afunc,ma) that returns the ma basis functions evaluated at x = x in the arrayafunc[1..ma].{void svbksb(float **u, float w[], float **v, int m, int n, float b[],float x[]);void svdcmp(float **a, int m, int n, float w[], float **v);int j,i;float wmax,tmp,thresh,sum,*b,*afunc;b=vector(1,ndata);afunc=vector(1,ma);for (i=1;i<=ndata;i++) {Accumulate coefficients of the fitting ma(*funcs)(x[i],afunc,ma);trix.tmp=1.0/sig[i];for (j=1;j<=ma;j++) u[i][j]=afunc[j]*tmp;b[i]=y[i]*tmp;}svdcmp(u,ndata,ma,w,v);Singular value decomposition.wmax=0.0;Edit the singular values, given TOL from thefor (j=1;j<=ma;j++)#define statement, between here ...if (w[j] > wmax) wmax=w[j];thresh=TOL*wmax;for (j=1;j<=ma;j++)if (w[j] < thresh) w[j]=0.0;...and here.svbksb(u,w,v,ndata,ma,b,a);*chisq=0.0;Evaluate chi-square.for (i=1;i<=ndata;i++) {(*funcs)(x[i],afunc,ma);for (sum=0.0,j=1;j<=ma;j++) sum += a[j]*afunc[j];*chisq += (tmp=(y[i]-sum)/sig[i],tmp*tmp);}free_vector(afunc,1,ma);free_vector(b,1,ndata);}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 extra array of size N × M to store the whole design matrix. This storageis overwritten by the matrix U.
Storage is also required for the M × M matrixV, but this is instead of the same-sized matrix for the coefficients of the normalequations. SVD can be significantly slower than solving the normal equations;however, its great advantage, that it (theoretically) cannot fail, more than makesup for the speed disadvantage.In the routine that follows, the matrices u,v and the vector w are input asworking space. The logical dimensions of the problem are ndata data points by mabasis functions (and fitted parameters). If you care only about the values a of thefitted parameters, then u,v,w contain no useful information on output. If you wantprobable errors for the fitted parameters, read on.67915.4 General Linear Least Squares#include "nrutil.h"void svdvar(float **v, int ma, float w[], float **cvm)To evaluate the covariance matrix cvm[1..ma][1..ma] of the fit for ma parameters obtainedby svdfit, call this routine with matrices v[1..ma][1..ma], w[1..ma] as returned fromsvdfit.{int k,j,i;float sum,*wti;wti=vector(1,ma);for (i=1;i<=ma;i++) {wti[i]=0.0;if (w[i]) wti[i]=1.0/(w[i]*w[i]);}for (i=1;i<=ma;i++) {Sum contributions to covariance matrix (15.4.20).for (j=1;j<=i;j++) {for (sum=0.0,k=1;k<=ma;k++) sum += v[i][k]*v[j][k]*wti[k];cvm[j][i]=cvm[i][j]=sum;}}free_vector(wti,1,ma);}ExamplesBe aware that some apparently nonlinear problems can be expressed so thatthey are linear.
For example, an exponential model with two parameters a and b,y(x) = a exp(−bx)(15.4.21)log[y(x)] = c − bx(15.4.22)can be rewritten aswhich is linear in its parameters c and b. (Of course you must be aware that suchtransformations do not exactly take Gaussian errors into Gaussian errors.)Also watch out for “non-parameters,” as iny(x) = a exp(−bx + d)(15.4.23)Here the parameters a and d are, in fact, indistinguishable. This is a good example ofwhere the normal equations will be exactly singular, and where SVD will find a zerosingular value.
SVD will then make a “least-squares” choice for setting a balancebetween a and d (or, rather, their equivalents in the linear model derived by taking thelogarithms). However — and this is true whenever SVD gives back a zero singularvalue — you are better advised to figure out analytically where the degeneracy isamong your basis functions, and then make appropriate deletions in the basis set.Here are two examples for user-supplied routines funcs.
The first one is trivialand fits a general polynomial to a set of data: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).Feeding the matrix v and vector w output by the above program into thefollowing short routine, you easily obtain variances and covariances of the fittedparameters a.
The square roots of the variances are the standard deviations ofthe fitted parameters. The routine straightforwardly implements equation (15.4.20)above, with the convention that singular values equal to zero are recognized ashaving been edited out of the fit.680Chapter 15.Modeling of Datavoid fpoly(float x, float p[], int np)Fitting routine for a polynomial of degree np-1, with coefficients in the array p[1..np].{int j;p[1]=1.0;for (j=2;j<=np;j++) p[j]=p[j-1]*x;The second example is slightly less trivial.
It is used to fit Legendre polynomialsup to some order nl-1 through a data set.void fleg(float x, float pl[], int nl)Fitting routine for an expansion with nl Legendre polynomials pl, evaluated using the recurrencerelation as in §5.5.{int j;float twox,f2,f1,d;pl[1]=1.0;pl[2]=x;if (nl > 2) {twox=2.0*x;f2=x;d=1.0;for (j=3;j<=nl;j++) {f1=d++;f2 += twox;pl[j]=(f2*pl[j-1]-f1*pl[j-2])/d;}}}Multidimensional FitsIf you are measuring a single variable y as a function of more than one variable— say, a vector of variables x, then your basis functions will be functions of a vector,X1 (x), .
. . , XM (x). The χ2 merit function is now"#2PMNXyi − k=1 ak Xk (xi )2χ =(15.4.24)σii=1All of the preceding discussion goes through unchanged, with x replaced by x. Infact, if you are willing to tolerate a bit of programming hack, you can use the aboveprograms without any modification: In both lfit and svdfit, the only use madeof the array elements x[i] is that each element is in turn passed to the user-suppliedroutine funcs, which duly gives back the values of the basis functions at that point.If you set x[i]=i before calling lfit or svdfit, and independently provide funcswith the true vector values of your data points (e.g., in global variables), then funcscan translate from the fictitious x[i]’s to the actual data points before doing its work.CITED REFERENCES AND FURTHER READING:Bevington, P.R.
1969, Data Reduction and Error Analysis for the Physical Sciences (New York:McGraw-Hill), Chapters 8–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.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
