15.2 Fitting Data to a Straight Line 665 =5(1+ SSu (15.2.19) = 1 Su (15.2.20) Cov(a,6)=- Su (15.2.21) Cov(a,b) Tab (15.2.22) GaOb #include <math.h> #include "nrutil.h" void fit(float x[],float y[],int ndata,float sig[],int mwt,float *a 100 float *b,float *siga,float *sigb,float *chi2,float *q) Given a set of data points x[1..ndata],y[1..ndata]with individual standard deviations sig[1..ndata],fit them to a straight line y a bx by minimizing x2.Returned are Cambridge from NUMERICAL RECIPES IN 18881892 a,b and their respective probable uncertainties siga and sigb,the chi-square chi2,and the goodness-of-fit probability q(that the fit would have x2 this large or larger).If mwt=0 on input,then the standard deviations are assumed to be unavailable:q is returned as 1.0 and (Nor the normalization of chi2 is to unit standard deviation on all points. f America server computer, float gammg(float a,float x); to make one paper UnN电.t THE int i; 是 ART float wt,t,sxoss,sx=0.0,sy=0.0,st2=0.0,ss,sigdat; *b=0.0: if (mwt){ Accumulate sums .. ss=0.0; strictly prohibited. Programs for (i=1;i<=ndata;i++){ ..with weights t=1.0/s0R(s1g[1]): Copyright (C) ss +wt; sx +x[i]*ut: sy +y[i]*wt; nail to directcustsen else for (i=1;i<=ndata;i++){ ...or without weights. sx+=x[1]; sy+=y[1]; 2 ss=ndata; @cambridge.org 1988-1992 by Numerical Recipes OF SCIENTIFIC COMPUTING(ISBN 0-521- -431085 sxoss=sx/ss; if (mwt){ for (i=1;i<=ndata;i++){ t=(x[i]-sxoss)/sig[i]; Software. st2 +t*t; *b+=t*y[i]/s1g[1]; (outside North America) else for (is1;i<=ndata;i++){ t=x[i]-sxoss; st2 +t*t; *b+=t*y[1]; *b/=st2; Solve for a,b,a and ob *a=(sy-sx*(*b))/ss; *siga=sqrt((1.0+sx*sx/(ss*st2))/ss); *sigb=sqrt(1.0/st2);15.2 Fitting Data to a Straight Line 665 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). σ2 a = 1 S 1 + S2 x SStt (15.2.19) σ2 b = 1 Stt (15.2.20) Cov(a, b) = − Sx SStt (15.2.21) rab = Cov(a, b) σaσb (15.2.22) #include <math.h> #include "nrutil.h" void fit(float x[], float y[], int ndata, float sig[], int mwt, float *a, float *b, float *siga, float *sigb, float *chi2, float *q) Given a set of data points x[1..ndata],y[1..ndata] with individual standard deviations sig[1..ndata], fit them to a straight line y = a + bx by minimizing χ2. Returned are a,b and their respective probable uncertainties siga and sigb, the chi-square chi2, and the goodness-of-fit probability q (that the fit would have χ2 this large or larger). If mwt=0 on input, then the standard deviations are assumed to be unavailable: q is returned as 1.0 and the normalization of chi2 is to unit standard deviation on all points. { float gammq(float a, float x); int i; float wt,t,sxoss,sx=0.0,sy=0.0,st2=0.0,ss,sigdat; *b=0.0; if (mwt) { Accumulate sums ... ss=0.0; for (i=1;i<=ndata;i++) { ...with weights wt=1.0/SQR(sig[i]); ss += wt; sx += x[i]*wt; sy += y[i]*wt; } } else { for (i=1;i<=ndata;i++) { ...or without weights. sx += x[i]; sy += y[i]; } ss=ndata; } sxoss=sx/ss; if (mwt) { for (i=1;i<=ndata;i++) { t=(x[i]-sxoss)/sig[i]; st2 += t*t; *b += t*y[i]/sig[i]; } } else { for (i=1;i<=ndata;i++) { t=x[i]-sxoss; st2 += t*t; *b += t*y[i]; } } *b /= st2; Solve for a, b, σa, and σb. *a=(sy-sx*(*b))/ss; *siga=sqrt((1.0+sx*sx/(ss*st2))/ss); *sigb=sqrt(1.0/st2);