Linear Regression Didn’t Give the Best Curve

Take a known curve:

f[x_]=0.05*Exp[1.45*x]

Plot[f[x],{x,0,5},AxesOrigin->{0,0}]

Generate data that is “close” to this curve.  The data:

xd={0.0,.5,1.0,1.5,2.0,2.5,3.0,3.5,4.0,4.5}
ydt=f[xd]

yd={0.0,0.105,0.215, 0.44, 0.91, 1.85, 3.85, 8.0, 16.5, 34.0}

Note that the y-values (yd) are rounded versions of the original function.

Graph of the data points:

dpts=Table[{xd[[i]],yd[[i]]},{i,1,10}]

dataplot=ListPlot[dpts,PlotStyle->Directive[Blue,PointSize[Medium]],Frame->True,FrameLabel->{"x","y"},Axes->None,PlotRange->{{-.25,5},{-5,40}},ImageSize->240]

Best square-error parabola

parabola=Fit[dpts,{x,x^2},x]

gp[x_]=parabola

parplot=Plot[parabola,{x,0,4.5},PlotStyle->Black]

Show[parplot,dataplot,Frame->True,Axes->None,FrameLabel->{"x","y"},Axes->None,PlotRange->{{-.25,5},{-5,40}},ImageSize->240]

errparabola=Sum[(yd[[i]]-gp[xd[[i]]])^2,{i,1,10}]

Conclusion:  Not a bad fit for using a quadratic model.  Note that we assume the model curve goes through the origin.  Square error of 102.428 is large compared to the indroduced “error” in the data.

Best fit to y=(Ax^b) using Log-Log data.  Note that the first data point is dropped.

lldatapts=Table[{Log[xd[[i]]],Log[yd[[i]]]},{i,2,10}]

lldataplot=ListPlot[lldatapts,PlotStyle->Directive[Red,PointSize[Large]],Frame->True,FrameLabel->{"ln(x)","ln(y)"},Axes->None,PlotRange->{{-.8,1.6},{-2.5,4.5}},ImageSize->240]

llfit=Fit[lldatapts,{1,x},x]

llfitplot=Plot[llfit,{x,-0.8,1.6},PlotStyle->Directive[Black,Dashed]]

Show[lldataplot,llfitplot,ImageSize->240]

It is easy to see that the data does not really fit this type of model, but let’s convert back to x vs y to see how bad this common technique really was.

loga=llfit/.{x->0}
a=Exp[loga]
b=Coefficient[llfit,x]

gr[x_]=a*x^b

ratplot=Plot[gr[x],{x,0,4.5},PlotStyle->Directive[Black,Dashed]]

Show[dataplot,ratplot,ImageSize->240]

errrat=Sum[(yd[[i]]-gr[xd[[i]]])^2,{i,1,10}]

Here we get a “reasonable” fit for the left half of the graph, but not the right.  Also notice that the square error is four times worst than the quadratic model.  Note that this is the result given by some modern calculaters using their power-regression function.

Best fit to y=(Ax^b) using multivariable calculus (Levenberg-Marquardt modified gradient method).

sols=FindFit[dpts,A*x^B,{A,B},x]

Note that the FindFit command, which uses a robust and well known Levenberg-Marquardt algorithm, still has problems with the (0, 0) data point.  Since the model passes through the point (0, 0), removing that point does not introduce additional error.  Removing that data point gives.

dpts2to10=Table[dpts[[i]],{i,2,10}]

sols=FindFit[dpts2to10,A*x^B,{A,B},x]

gr2[x_]=A*x^B/.sols

ratplot2=Plot[gr2[x],{x,0,4.8},PlotStyle->Black]

Show[dataplot,ratplot2,ImageSize->240]

Show[dataplot,ratplot,ratplot2,ImageSize->240]

errrat2=Sum[(yd[[i]]-gr2[xd[[i]]])^2,{i,1,10}]

The less commonly used multivariable calculus method does a much better job of fitting this model.  This is the actual least square error for the power function model.

Best fit to y=S (e^(T x)) using semi-log data.  Note that the first data point is dropped as per common convention.  In this case, the model does not pass through (0, 0), so additional error is introduced.

sldatapts=Table[{xd[[i]],Log[yd[[i]]]},{i,2,10}]

sldataplot=ListPlot[sldatapts,PlotStyle->Directive[Green,PointSize[Large]],AxesOrigin->{0,0},Axes->False,Frame->True,FrameLabel->{"x","ln(y)"},ImageSize->240]

slfit=Fit[sldatapts,{1,x},x]

slfitplot=Plot[slfit,{x,0,5},PlotStyle->Directive[Black,Dashed]]

Show[sldataplot,slfitplot,ImageSize->240]

This time the data does seem to fit the model.

logs=slfit/.{x->0}
s=Exp[logs]
t=Coefficient[slfit,x]

ge[x_]=s*Exp[t*x]

expplot=Plot[ge[x],{x,0,4.9},PlotStyle->Directive[Black,Dotted]]

Show[dataplot,expplot,ImageSize->240]

errexp=Sum[(yd[[i]]-ge[xd[[i]]])^2,{i,1,10}]

Here the square error appears to be good.  This is the exponential regression given by some modern calculators.  However, the multivariable calculus method will turn out to be much better.

Best fit to y=s (e^tx) using calculus.

solse=FindFit[dpts,S*Exp[T*x],{S,T},x]

ge2[x_]=S*Exp[T*x]/.solse

expplot2=Plot[ge2[x],{x,0,4.9},PlotStyle->Directive[Black,Dotted]]

Show[dataplot,expplot2,ImageSize->240]

There is not visual difference between this fit and the one using semi-log linearization.  But looking at the square error:

errexp2=Sum[(yd[[i]]-ge2[xd[[i]]])^2,{i,1,10}]

Note that this gives a square error eight times smaller then when we linearize the data.  Again the calculus method gives the least square error.