678 Chapter 15.Modeling of Data an extra array of size N x M to store the whole design matrix.This storage is overwritten by the matrix U.Storage is also required for the M x M matrix V,but this is instead of the same-sized matrix for the coefficients of the normal equations.SVD can be significantly slower than solving the normal equations; however,its great advantage,that it(theoretically)cannot fail,more than makes up for the speed disadvantage. In the routine that follows,the matrices u,v and the vector w are input as working space.The logical dimensions of the problem are ndata data points by ma basis functions (and fitted parameters).If you care only about the values a of the fitted parameters,then u,v,w contain no useful information on output.If you want probable errors for the fitted parameters,read on. 19881992 #include "nrutil.h" #define TOL 1.0e-5 Default value for single precision and vari- ables scaled to order unity void svdfit(f1oatx[],float y[l,float sig[☐，int ndata,f1oata☐，int ma, float **u,float **v,float w[],float *chisg, from NUMERICAL RECIPESI void (*funcs)(float,float []int)) Given a set of data points x[1..ndata],y[1..ndata]with individual standard deviations sig[1..ndata],use x2 minimization to determine the coefficients a[1..ma]of the fit- (Nort server 令 ting function y =>ai x afunci(x).Here we solve the fitting equations using singular value decomposition of the ndata by ma matrix,as in $2.6.Arrays u[1..ndata][1..ma], THE v[1..ma][1..ma],and w[1..ma]provide workspace on input;on output they define the Americ computer, singular value decomposition,and can be used to obtain the covariance matrix.The pro- ART gram returns values for the ma fit parameters a,and x2,chisq.The user supplies a routine funcs(x,afunc,ma)that returns the ma basis functions evaluated at x =x in the array afunc[1..ma] Programs void svbksb(float**u,float w[，f1oat**v,intm,intn,float b☐ f1oatx▣)； void svdcmp(float **a,int m,int n,float w,float **v); int j,i; to dir float wmax,tmp,thresh,sum,*b,*afunc; b=vector(1,ndata); OF SCIENTIFIC COMPUTING(ISBN afunc=vector(1,ma); 19881992 for (i=1;i<-ndata;i++){ Accumulate coefficients of the fitting ma- (*funcs)(x[i],afunc,ma); trix. tmp=1.0/sig[1]; 10-621 for (j=1;j<=ma;j++)u[i][j]=afunc[j]*tmp; b[i]=y[i]*tmp; 43108 svdcmp(u,ndata,ma,w,v); Singular value decomposition. Numerical Recipes wmax=0.0; Edit the singular values,given TOL from the for (j=1;j<=ma;j++) #define statement,between here... (outside if (w[j]wmax) wmax=w[】; thresh=TOL*wmax; North Software. 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);678 Chapter 15. Modeling of Data Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machine￾readable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). an extra array of size N × M to store the whole design matrix. This storage is overwritten by the matrix U. Storage is also required for the M × M matrix V, but this is instead of the same-sized matrix for the coefficients of the normal equations. SVD can be significantly slower than solving the normal equations; however, its great advantage, that it (theoretically) cannot fail, more than makes up for the speed disadvantage. In the routine that follows, the matrices u,v and the vector w are input as working space. The logical dimensions of the problem are ndata data points by ma basis functions (and fitted parameters). If you care only about the values a of the fitted parameters, then u,v,w contain no useful information on output. If you want probable errors for the fitted parameters, read on. #include "nrutil.h" #define TOL 1.0e-5 Default value for single precision and vari￾ables scaled to order unity. void 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 deviations sig[1..ndata], use χ2 minimization to determine the coefficients a[1..ma] of the fit￾ting function y = i ai × afunci(x). Here we solve the fitting equations using singular value 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 the singular value decomposition, and can be used to obtain the covariance matrix. The pro￾gram returns values for the ma fit parameters a, and χ2, chisq. The user supplies a routine funcs(x,afunc,ma) that returns the ma basis functions evaluated at x = x in the array afunc[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 the for (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); }