15.4 General Linear Least Squares 675 (*funcs)(x[i],afunc,ma); ym=y[i]: if (mfit ma){ Subtract off dependences on known pieces for (j=1;j<=ma;j++) of the fitting function. if (!ia[j])ym -a[j]*afunc[j]; sig2i=1.0/SQR(sig[i]); for(j=0,1=1;1<ma;1+)( if(ia[1]){ wt=afunc[l]*sig2i; for(j+,k=0,m=1;m<=1;m+) if (ia[m])covar[j][++k]+wt*afunc[m] http://www.nr Permission is read able files beta[j门[1]+=ym*wt: 83g granted for 19881992 for (j=2;j<=mfit;j++) Fill in above the diagonal from symmetry. for (k=1;k<j;k++) 1-.200 covar [k][j]=covar[j][k]; gaussj(covar,mfit,beta,1); Matrix solution. for(j=0,1=1;1<ema;1++) if(ia[1])a[1]=beta[+j][1]; Partition solution to appropriate coefficients from NUMERICAL RECIPES IN *chisq=0.0; for (i=1;i<=ndat;i++){ Evaluate x2 of the fit. (*funcs)(x[i],afunc,ma); (North to make for (sum=0.0,j=1;j<=ma;j++)sum +a[j]*afunc[j]; *chisq +SQR((y[i]-sum)/sig[i]); America 州bMe se one paper UnN电.t THE ART covsrt(covar,ma,ia,mfit); Sort covariance matrix to true order of fitting free_vector(afunc,1,ma); coefficients. free_matrix(beta,1,ma,1,1); Programs That last call to a function covsrt is only for the purpose of spreading the covariances back into the full ma x ma covariance matrix,in the proper rows and to dir columns and with zero variances and covariances set for variables which were held frozen. The function covsrt is as follows. 1788-1982 OF SCIENTIFIC COMPUTING(ISBN v@cam #define SWAP(a,b){swap=(a);(a)=(b);(b)=swap; Numerical Recipes 10-621 void covsrt(float **covar,int ma,int ia[],int mfit) 43108 Expand in storage the covariance matrix covar,so as to take into account parameters that are being held fixed.(For the latter,return zero covariances.) int i,j,k; (outside 膜 float swap; Software. for (ismfit+1;i<=ma;i++) ying of for (j=1;j<=i;j++)covar [i][j]=covar[j][i]=0.0; k=mfit; for（j=ma;j>=1;j--)[ if(1a[j]){ for (i=1;i<=ma;i++)SWAP(covar[i][k],covar[i][j]) for (i=1;i<=ma;i++)SWAP(covar[k][i],covar[j][i]) 2215.4 General Linear Least Squares 675 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). (*funcs)(x[i],afunc,ma); ym=y[i]; if (mfit < ma) { Subtract off dependences on known pieces for (j=1;j<=ma;j++) of the fitting function. if (!ia[j]) ym -= a[j]*afunc[j]; } sig2i=1.0/SQR(sig[i]); for (j=0,l=1;l<=ma;l++) { if (ia[l]) { wt=afunc[l]*sig2i; for (j++,k=0,m=1;m<=l;m++) if (ia[m]) covar[j][++k] += wt*afunc[m]; beta[j][1] += ym*wt; } } } for (j=2;j<=mfit;j++) Fill in above the diagonal from symmetry. for (k=1;k<j;k++) covar[k][j]=covar[j][k]; gaussj(covar,mfit,beta,1); Matrix solution. for (j=0,l=1;l<=ma;l++) if (ia[l]) a[l]=beta[++j][1]; Partition solution to appropriate coefficients *chisq=0.0; a. for (i=1;i<=ndat;i++) { Evaluate χ2 of the fit. (*funcs)(x[i],afunc,ma); for (sum=0.0,j=1;j<=ma;j++) sum += a[j]*afunc[j]; *chisq += SQR((y[i]-sum)/sig[i]); } covsrt(covar,ma,ia,mfit); Sort covariance matrix to true order of fitting free_vector(afunc,1,ma); coefficients. free_matrix(beta,1,ma,1,1); } That last call to a function covsrt is only for the purpose of spreading the covariances back into the full ma × ma covariance matrix, in the proper rows and columns and with zero variances and covariances set for variables which were held frozen. The function covsrt is as follows. #define SWAP(a,b) {swap=(a);(a)=(b);(b)=swap;} void covsrt(float **covar, int ma, int ia[], int mfit) Expand in storage the covariance matrix covar, so as to take into account parameters that are being held fixed. (For the latter, return zero covariances.) { int i,j,k; float swap; for (i=mfit+1;i<=ma;i++) for (j=1;j<=i;j++) covar[i][j]=covar[j][i]=0.0; k=mfit; for (j=ma;j>=1;j--) { if (ia[j]) { for (i=1;i<=ma;i++) SWAP(covar[i][k],covar[i][j]) for (i=1;i<=ma;i++) SWAP(covar[k][i],covar[j][i]) k--; } } }