704 Chapter 15. Modeling of Data #include <math.h> #include "nrutil.h" int ndatat; float *xt,*yt,aa,abdevt; void medfit(float x[],float y[],int ndata,float *a,float *b,float *abdev) Fits y =a+ba by the criterion of least absolute deviations.The arrays x[1..ndata]and y[1..ndata]are the input experimental points.The fitted parameters a and b are output, along with abdev,which is the mean absolute deviation (in y)of the experimental points from the fitted line.This routine uses the routine rofunc,with communication via global variables. float rofunc(float b); int j; float bb,b1,b2,del,f,f1,f2,sigb,tempi f1 oat sx=0.0,sy=0.0,sxy=0.0,sxx=0.0,ch1sq=0.0; 18881992 ndatat=ndata; xt=x; 1600 (including this one) yt-yi for (j=1;j<=ndata;j++){ As a first guess for a and b,we will find the sx +x[j]; least-squares fitting line. sy +=y[j]; to any users to make from NUMERICAL RECIPES IN sxy +x[j]*y[i]: sxx +x[j]*x[j]; del=ndata*sxx-sx*sx; Cambridge University Press. THE aa=(sxx*sy-sx*sxy)/del; Least-squares solutions. (North America server computer, bb=(ndata*sxy-sx*sy)/del; one paper ART for (j=1;j<=ndata;j++) chisg +(temp=y[j]-(aa+bb*x[j]),temp*temp); 是 sigb=sqrt(chisq/del); The standard deviation will give some idea of b1-bb; how big an iteration step to take. fi=rofunc(b1): st st for their Programs if (sigb 0.0){ b2=bb+SIGN(3.0*sigb,f1); Guess bracket as 3-o away,in the downhill direction known from f1 Copyright (C) f2=rofunc(b2): 1f(b2==b1)C *a=aa; to directcustsen *b=bb; *abdev=abdevt/ndata; OF SCIENTIFIC COMPUTING(ISBN 0-521 return; wh11e(f1*f2>0.0) Bracketing. bb=b2+1.6*(b2-b1); b1=b2; v@cambridge.org personal use.Further reproduction, To order Numerical Recipes books or 1988-1992 by Numerical Recipes f1=f2: -431085 b2=bbi f2=rofunc(b2); (outside sigb=0.01*sigb; Refine until error a negligible number of standard North Software. while (fabs(b2-b1)>sigb){ deviations. bb=b1+0.5*(b2-b1); Bisection. if (bb =b1 Il bb =b2)break; f=rofunc(bb); 1f(f*f1>=0.0) f1=f; bl=bb; else f2=f; b2=bbi 2704 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). #include <math.h> #include "nrutil.h" int ndatat; float *xt,*yt,aa,abdevt; void medfit(float x[], float y[], int ndata, float *a, float *b, float *abdev) Fits y = a + bx by the criterion of least absolute deviations. The arrays x[1..ndata] and y[1..ndata] are the input experimental points. The fitted parameters a and b are output, along with abdev, which is the mean absolute deviation (in y) of the experimental points from the fitted line. This routine uses the routine rofunc, with communication via global variables. { float rofunc(float b); int j; float bb,b1,b2,del,f,f1,f2,sigb,temp; float sx=0.0,sy=0.0,sxy=0.0,sxx=0.0,chisq=0.0; ndatat=ndata; xt=x; yt=y; for (j=1;j<=ndata;j++) { As a first guess for a and b, we will find the sx += x[j]; least-squares fitting line. sy += y[j]; sxy += x[j]*y[j]; sxx += x[j]*x[j]; } del=ndata*sxx-sx*sx; aa=(sxx*sy-sx*sxy)/del; Least-squares solutions. bb=(ndata*sxy-sx*sy)/del; for (j=1;j<=ndata;j++) chisq += (temp=y[j]-(aa+bb*x[j]),temp*temp); sigb=sqrt(chisq/del); The standard deviation will give some idea of b1=bb; how big an iteration step to take. f1=rofunc(b1); if (sigb > 0.0) { b2=bb+SIGN(3.0*sigb,f1); Guess bracket as 3-σ away, in the downhill direction known from f1. f2=rofunc(b2); if (b2 == b1) { *a=aa; *b=bb; *abdev=abdevt/ndata; return; } while (f1*f2 > 0.0) { Bracketing. bb=b2+1.6*(b2-b1); b1=b2; f1=f2; b2=bb; f2=rofunc(b2); } sigb=0.01*sigb; Refine until error a negligible number of standard while (fabs(b2-b1) > sigb) { deviations. bb=b1+0.5*(b2-b1); Bisection. if (bb == b1 || bb == b2) break; f=rofunc(bb); if (f*f1 >= 0.0) { f1=f; b1=bb; } else { f2=f; b2=bb; } } }