Least Squares Fitting of Data to a Curve

1 Least Squares Fittingof data to a CurveGerald RecktenwaldPortland State UniversityDepartment of Mechanical slides are a supplement to the bookNumerical Methods withMatlab: Implementations and Applications, by Gerald W. Recktenwald,c 2000 2007, Prentice-Hall, Upper Saddle River, NJ. These slides arecopyrightc 2000 2007 Gerald W. Recktenwald. The PDF versionof these slides may be downloaded or stored or printed only fornoncommercial, educational use. The repackaging or sale of theseslides in any form, without written consent of the author, is latest version of this PDF file, along with other supplemental materialfor the book, can be found ~gerry/nmm/.Version 6, 2007page 1 Overview Fitting a line to data Geometric interpretation Residuals of the overdetermined system The normal equations Nonlinear fits via coordinate transformation Fitting arbitrary linear combinations of basis functions Mathematical formulation Solution via normal equations Solution via QR factorization Polynomial Curve fits with the built-inpolyfitfunction Multivariate fittingNMM: Least Squares Curve -Fittingpage 2 Fitting a Line to DataGivenmpairs of data :(xi,yi),i=1.

2 ,mFind the coefficients and such thatF(x)= x+ is a good fit to the dataQuestions: How do we define good fit? How do we compute and after a definition of good fit is obtained?NMM: Least Squares Curve -Fittingpage 3 Plausible FitsPlausible fits are obtained by adjusting theslope ( ) and intercept ( ). Here is agraphical representation of potential fits to aparticular set of dataWhich of the lines provides the best fit?12345612345xyNMM: Least Squares Curve -Fittingpage 4 The ResidualThe difference between the givenyivalueand the fit function evaluated atxiisri=yi F(xi)=yi ( xi+ )riis theresidualfor the data pair(xi,yi).riis the vertical distance between the knowndata and the fit : Least Squares Curve -Fittingpage 5 Minimizing the ResidualTwo criteria for choosing the best fitminimizeX|ri|orminimizeXr2iFor statistical and computational reasons choose minimization of =Pr2i =mXi=1[yi ( xi+ )]2 The best fit is obtained by the values of and that minimize.

3 NMM: Least Squares Curve -Fittingpage 6 Orthogonal Distance FitAn alternative to minimizing theresidual is to minimize the orthogonaldistance to the known as theOrthogonal Distance , , Ake Bj ork,NumericalMethods for Least Squares Problems,1996, SIAM, : Least Squares Curve -Fittingpage 7 Least Squares Fit(1)Theleast Squares fitis obtained by choosing the and so thatmXi=1r2iis a minimum. Let = r 22to simplify the and by minimizing = ( , ). The minimum requires =constant=0and =constant=0 NMM: Least Squares Curve -Fittingpage 8 Least Squares Fit(2)Carrying out the differentiation leads toSxx +Sx =Sxy(1)Sx +m =Sy(2)whereSxx=mXi=1xixiSx=mXi=1xiSxy=mX i=1xiyiSy=mXi=1yiNote:Sxx,Sx,Sxy,andSyyc an be directly computed from the given(xi,yi) , Equation (1) and (2) are two equations for the two unknowns, and.

4 NMM: Least Squares Curve -Fittingpage 9 Least Squares Fit(3)Solving equations (1) and (2) for and yields =1d(SxSy mSxy)(3) =1d(SxSxy SxxSy)(4)withd=S2x mSxx(5)NMM: Least Squares Curve -Fittingpage 10 Overdetermined System for a Line Fit(1)Now, let s rederive the equations for the fit. This will give us insight into the process orfitting arbitrary linear combinations of any two points we can write x1+ =y1 x2+ =y2or x11x21 = y1y2 But why just pick two points?NMM: Least Squares Curve -Fittingpage 11 Overdetermined System for a Line Fit(2)Writing out the x+ =yequation forallof the known points(xi,yi),i=1,..,mgives = y= :We cannot solveAc=ywith Gaussian elimination. Unless the system is consistent ( , unlessylies in the column space ofA) it is impossible to find thec=( , )Tthat exactly satisfiesallmequations.

5 The system is consistentonlyif all the data points lie along a single : Least Squares Curve -Fittingpage 12 Normal Equations for a Line FitCompute = r 22,wherer=y Ac = r 22=rTr=(y Ac)T(y Ac)=yTy (Ac)Ty yT(Ac)+cTATAc=yTy 2yTAc+ requires c= 2 ATy+2 ATAc=0or(ATA)c=ATbThis is the matrix formulation of equations (1) and (2).NMM: Least Squares Curve -Fittingpage fits a line to a set of data by solving the normal [c,R2] = linefit(x,y)% linefit Least - Squares fit of data to y = c(1)*x + c(2)%% Synopsis: c = linefit(x,y)% [c,R2] = linefit(x,y)%% Input: x,y = vectors of independent and dependent variables%% Output: c = vector of slope, c(1), and intercept, c(2) of Least sq. line fit% R2 = (optional) coefficient of determination; 0 <= R2 <= 1% R2 close to 1 indicates a strong relationship between y and xif length(y)~= length(x), error( x and y are not compatible ); endx = x(:); y = y(:); % Make sure that x and y are column vectorsA = [x ones(size(x))]; % m-by-n matrix of overdetermined systemc = (A *A)\(A *y); % Solve normal equationsif nargout>1r=y-A*c;R2 = 1 - (norm(r)/norm(y-mean(y)))^2;endNMM: Least Squares Curve -Fittingpage 14 Line Fitting ExampleStore data and perform the fit>>x=[1245];>>y=[1223];>> c = linefit(x,y)c= and plot the fit>> xfit = linspace(min(x),max(x));>> yfit = c(1)*xfit + c(2)>> plot(x,y, o ,xfit,yfit, - ).

6 data and fit functionNMM: Least Squares Curve -Fittingpage 15R2 Statistic(1)R2is a measure of how well the fit function follows the trend in the R2 : yis the value of the fit function at theknown data a line fit yi=c1xi+c2 yis the average of theyvalues y=1mXyiThen:R2=X( yi y)2X(yi y)2=1 r 22P(yi y)2 WhenR2 1the fit function follows the trend of the 0the fit is not significantly better than approximating the data by its : Least Squares Curve -Fittingpage 16 Graphical Interpretation of theR2 StatisticConsider a line fit to a data set withR2=1 r 22P(yi y)2= distances between givenydata and the Least Squares line lines show contributions to r distances between givenydata and the average of lines show contributions toP(yi y)212345612345xyNMM: Least Squares Curve -Fittingpage 17R2 Statistic: Example CalculationConsider the variation of the bulk modulusof Silicon Carbide as a function oftemperature (Cf.)

7 Example )T( C)20 500 1000 1200 1400 1500G(GP a)203 197191188186184>> [t,D,labels] = loadColData( ,6,5);>> g = D(:,1);>> [c,R2] = linefit(t,g);c= = CBulk Modulus GPaNMM: Least Squares Curve -Fittingpage 18 Fitting Transformed Non-linear Functions(1) Some nonlinear fit functionsy=F(x)can be transformed to an equation of theformv= u+ Linear Least Squares fit to a line is performed on the transformed variables. Parameters of the nonlinear fit function are obtained by transforming back to theoriginal variables. The linear Least Squares fit to the transformed equations does not yield the same fitcoefficients as a direct solution to the nonlinear Least Squares problem involving theoriginal fit :y=c1ec2x lny= x+ y=c1xc2 lny= lnx+ y=c1xec2x ln(y/x)= x+ NMM: Least Squares Curve -Fittingpage 19 Fitting Transformed Non-linear Functions(2)Considery=c1ec2x(6)Taking the logarithm of both sides yieldslny=lnc1+c2xIntroducing the variablesv=lnyb=lnc1a=c2transforms equation (6) tov=ax+bNMM: Least Squares Curve -Fittingpage 20 Fitting Transformed Non-linear Functions(3)The preceding steps are equivalent to graphically obtainingc1andc2by plotting the dataon semilog + 210 1100101xyNMM: Least Squares Curve -Fittingpage 21 Fitting Transformed Non-linear Functions(4)Considery=c1xc2.

8 Taking the logarithm of both sides yieldslny=lnc1+c2lnx(7)Introduce the transformed variablesv=lnyu=lnxb=lnc1a=c2and equation (7) can be writtenv=au+bNMM: Least Squares Curve -Fittingpage 22 Fitting Transformed Non-linear Functions(5)The preceding steps are equivalent to graphically obtainingc1andc2by plotting the dataon log-log + 210 1100101100101102xyNMM: Least Squares Curve -Fittingpage 23 Example: Fitting data toy=c1xec2xConsidery=c1xec2x. The transformationv=ln yx a=c2b=lnc1results in the linear equationv=ax+bNMM: Least Squares Curve -Fittingpage 24 Fitting Transformed Non-linear Functions(6)The preceding steps are equivalent to graphically obtainingc1andc2by plotting the dataon semilog (y/x)=c2x+ 210 1100101xyNMM: Least Squares Curve -Fittingpage uses a linearizing transformation to fity=c1xec2xto c = xexpfit(x,y)% xexpfit Least Squares fit of data to y = c(1)*x*exp(c(2)*x)%% Synopsis: c = xexpfit(x,y)%% Input: x,y = vectors of independent and dependent variable values%% Output: c = vector of coefficients of y = c(1)*x*exp(c(2)*x)z = log( ); % Natural log of element-by-element divisionc = linefit(x,z); % Fit is performed by linefitc = [exp(c(2)); c(1)].

9 % Extract parameters from transformationNMM: Least Squares Curve -Fittingpage 26 Example: Fit Synthetic DataFity=c1xec2xto synthetic data . SeedemoXexp>> % Synthetic data with noise, avoid x=0>> x0 = ;>> noise = ;>> x = linspace(x0,2,200);>> y = 5*x.*exp(-3*x);>>yn=y+noise*(rand(size(x )) );>> % Guarantee yn>0 for log(yn)>> yn = abs(yn);>> c = xexpfit(x,yn);c= = c2 = points in synthetic data setoriginalnoisy fit NMM: Least Squares Curve -Fittingpage 27 Summary of Transformations Transform(x, y) data as needed Uselinefit Transform results oflinefitback>> x = .. % original data >> y = ..>> u = .. % transform the data >> v = ..>> a = linefit(u,v)>> c = .. % transform the coefficientsNMM: Least Squares Curve -Fittingpage 28 Summary of Line Fitting (1) pairs are given:(xi,yi),i=1.

10 , The fit functiony=F(x)=c1x+c2hasn=2basis functionsf1(x)=xandf2(x)= Evaluating the fit function for each of themdata points gives an overdeterminedsystem of equationsAc=ywherec=[c1,c2]T,y=[y1,y2,.. ,ym]T,andA=2664f1(x1)f2(x1)f1(x2)f2(x2). .f1(xm)f2(xm)3775= : Least Squares Curve -Fittingpage 29 Summary of Line Fitting (2)4. The Least - Squares principle defines the best fit as the values ofc1andc2that minimize (c1,c2)= y F(x) 22= y Ac Minimizing of (c1,c2)leads to the normal equations(ATA)c=ATy,6. Solving the normal equations gives the slopec1and interceptc2of the best fit : Least Squares Curve -Fittingpage 30 Fitting Linear Combinations of Functions Definition of fit function and basis functions Formulation of the overdetermined system Solution via normal equations:fitnorm Solution via QR factorization:fitqrand\NMM: Least Squares Curve -Fittingpage 31 Fitting Linear Combinations of Functions(1)Consider the Fitting functionF(x)=cf1(x)+c2f2(x)+.