Non-linear least squares fitting tool

INSTRUCTIONS:  Each of the cells below can be evaluated one at a time.  To repeat or redo a step, just reevaluate the cell.  If you are unfamiliar with Mathematica, a notebook cell (delimited by a square bracket at the right) is evaluated by selecting the cell and either applying the Shift-Enter key combination or the Enter(numeric) key on the numeric keypad at the far right-hand side of your keyboard.  To select a cell, either place your cursor within the input region of the cell or click on the cell bracket at the right with your mouse (usually the second bracket in from the right edge of the frame).  When asked if you want to evaluate the initialization cells, choose "Yes".  This script REQUIRES AT LEAST MATHEMATICA 7 to run.

Initialization Cells (don't modify these)

Avoid editing the content of this cell unless you know what you are doing.  Hide it if you don't want to see it.

(* Clear any past stuff out of memory to avoid carryover *)
Clear["`*"] 

(* Set the directory to find the data file to that in which the notebook file is located *)
(* If the notebook was specified as an argument to the Mathematica command and we were already in the directory where the notebook file is found we don't want to do this *)
TestFileName="FileName"/.NotebookInformation[EvaluationNotebook[]];
If[Length[TestFileName]>0,
SetDirectory[ToFileName[Extract[TestFileName,{1},FrontEnd`FileName]]];]
Off[General::spell1];

(* Define the fitting function - if the data in incoming includes 3 values/record use the third value as the weights for fitting.  If the 3rd column is the error in the second column the weight is the inverse of the square of that value. *)
nicefit[incoming_, func_,params_,col3errors_,vars_]:= Module[{datax,datay,dataerr,nlm, fitresults,fiterrors,haveerrors,ptable},
haveerrors = Length[Transpose[incoming]] == 3;
If[haveerrors,
{datax,datay,dataerr}=Transpose[incoming];
If[col3errors,
newweights = Quiet[(1/(dataerr*dataerr))/.ComplexInfinity->10^(-20)];
newweights=newweights/Max[newweights];
nlm=NonlinearModelFit[Transpose[{datax,datay}],func,params,vars,Weights->Quiet[(1/(dataerr*dataerr))/.ComplexInfinity->10^(-20)],VarianceEstimatorFunction->Automatic],nlm = NonlinearModelFit[Transpose[{datax,datay}],func,params,vars,Weights->dataerr]],
{datax,datay} = Transpose[incoming];
nlm = NonlinearModelFit[Transpose[{datax,datay}],func,params,vars]
];
(*fitresults = nlm["BestFitParameters"];*)(*fiterrors = nlm["ParameterErrors"];*)
(*Print[Flatten[{fitresults}]];*)(*,Table[σ[params[[i,1]]]->fiterrors[[i]],{i,Length[params]}]}]];*)
(*ptable=nlm["ParameterTable"];
Print[Flatten[ptable]];
ptable=nlm["ParameterConfidenceIntervalTable",ConfidenceLevel->0.95];
Print["95% Confidence Interval:"];
Print[Flatten[ptable]];
ptable=nlm["ParameterConfidenceIntervalTable",ConfidenceLevel->0.63];
Print["63% Confidence Interval:"];
Print[Flatten[ptable]];*)
{nlm["BestFit"],nlm["AdjustedRSquared"],nlm}
]

(* Draw a nice plot with an overlaid function.  Show error bars if there are 3 values/record and showerror is True *)
fitresplot[incoming_,func_,rsquared_,showequation_,Xlabel_,Ylabel_,fontpoints_,logx_:False,logy_:False, showerror_:True]:= Module[{datax,datay,dataerr,newx,newy,newerr,dx,errorbar,haveerrors,plotstuff,standopts,myplot,oplot,cplot},
haveerrors = Length[Transpose[incoming]] == 3;
If[haveerrors,
(* have errors in y-values *)
{datax,datay,dataerr} = Transpose[Select[incoming,If[logx, #[[1]]>0, True] && If[logy,#[[2]]>0,True]&]];
If[showerror,
errorbar[x_,y_,dy_,dx_]:={Line[{{x,y-dy},{x,y+dy}}],(*Line[{{x-dx,y-dy},{x+dx,y-dy}}],Line[{{x-dx,y+dy},{x+dx,y+dy}}],*)Point[{x,y}]},
errorbar[x_,y_,dy_,dx_] :={Point[{x,y}]}
];
newx = If[logx,Log/@ datax,datax];
newerr = If[logy,dataerr/Abs[datay],dataerr];
newy = If[logy,Log/@ datay,datay];
dx = (Max[newx]-Min[newx])/Length[newx]/4;
plotstuff = Map[errorbar[Sequence@@#,dx]&,Transpose[{newx,newy,newerr}]],

(* don't have error in y-values *)
{datax,datay} = Transpose[Select[incoming,If[logx, #[[1]]>0, True] && If[logy,#[[2]]>0,True]&]];
errorbar[x_,y_,dx_]:={Point[{x,y}]};
newx = If[logx,Log/@ datax,datax];
newy = If[logy,Log/@ datay,datay];
dx = (Max[newx]-Min[newx])/Length[newx]/4;
plotstuff = Map[errorbar[Sequence@@#,dx]&,Transpose[{newx,newy}]]
];
myplot := Which[Not[logx]&&Not[logy],Plot,logx&&Not[logy],LogLinearPlot,Not[logx]&&logy,LogPlot,logx&&logy,LogLogPlot];
oplot = Show[ Graphics[{PointSize[0.015],Thickness[0.002],plotstuff}]];
cplot = myplot[func,{x,datax[[1]],datax[[-1]]},PlotRange->All,AxesLabel->{Xlabel,Ylabel},PlotStyle->{Green}];
If[haveerrors,
If[showequation,
standopts = {BaseStyle->{FontSize->fontpoints,FontFamily->"Helvetica"},Axes->True,AspectRatio->1/GoldenRatio,PlotRange->All,AxesOrigin->{Min[newx],Min[newy-newerr]},PlotLabel->StringJoin[ToString[func,StandardForm],rsquared],AxesLabel->{Xlabel,Ylabel}},
standopts = {BaseStyle->{FontSize->fontpoints,FontFamily->"Helvetica"},Axes->True,AspectRatio->1/GoldenRatio,PlotRange->All,AxesOrigin->{Min[newx],Min[newy-newerr]},AxesLabel->{Xlabel,Ylabel}}
];,
If[showequation,
standopts = {BaseStyle->{FontSize->fontpoints,FontFamily->"Helvetica"},Axes->True,AspectRatio->1/GoldenRatio,PlotRange->All,AxesOrigin->{Min[newx],Min[newy]},PlotLabel->StringJoin[ToString[func,StandardForm],rsquared],AxesLabel->{Xlabel,Ylabel}},
standopts = {BaseStyle->{FontSize->fontpoints,FontFamily->"Helvetica"},Axes->True,AspectRatio->1/GoldenRatio,PlotRange->All,AxesOrigin->{Min[newx],Min[newy]},AxesLabel->{Xlabel,Ylabel}}
];
];
Show[oplot,cplot,standopts]
]

(* If we only have two points we will do a manual fit of a straight line to those points *)
(* Parameter uncertainty determined by standard error propagation formula *)
twopointfit[incoming_, col3errors_]:= Module[{datax,datay,dataerr,haveerrors,sol,y1,y2,x1,x2,a,b,σy1,σy2,σa,σb,values},
sol=Solve[{y1==a+b*x1,y2==a+b*x2},{a,b}];
a=a/.sol[[1]];
b=b/.sol[[1]];
σa=Sqrt[(D[a,y1]*σy1)^2+(D[a,y2]*σy2)^2];
σb=Sqrt[(D[b,y1]*σy1)^2+(D[b,y2]*σy2)^2];
haveerrors = Length[Transpose[incoming]] == 3;
If[haveerrors,
{datax,datay,dataerr}=Transpose[incoming];
values={x1->datax[[1]],x2->datax[[2]],y1->datay[[1]],y2->datay[[2]],σy1->dataerr[[1]],σy2->dataerr[[2]]};
a=a/.values;
b=b/.values;
If[col3errors,
σa=σa/.values;
σb=σb/.values;,
σa=0;
σb=0;];
,
{datax,datay} = Transpose[incoming];
values={x1->datax[[1]],x2->datax[[2]],y1->datay[[1]],y2->datay[[2]]};
a=a/.values;
b=b/.values;
σa=0;
σb=0;
];
{a+b*x,a,b,σa,σb}
]

Input the experimental data.

Modify the entry below to set datafilename to the name of the file containing your data, which must be located in the same directory as NLSQfit240.nb.  The data file should be an ascii text file with with either two or three columns: either (x y) or (x y yerror).  The columns must be separated by white space (spaces or tabs) - comma-separated values can not be used.  If you have a single column of data you will need to define the values of xinit and xstep below to generate the array of x values needed for the fit.

If your columns aren’t in the given order you will need to edit the column definitions below (xcol, ycol, errorcol).

datafilename = " -- put data file name here --";

incoming = ReadList[datafilename,Number,RecordLists->True];
noxval = Length[Transpose[incoming]] == 1;
(* If you have single column data (y value only) you have to enter an initial value and a step size for the x array to be generated! *)
(* The default will just generate point numbers starting at 0, i.e., an array of integers from 0 to N-1 *)
(* The data have to be reread because we don't want a single column to be sorted *)
If[noxval,
xinit = 0;
xstep = 1;
xend = xinit+xstep*(Length[incoming]-1);
newx = Table[xx,{xx,xinit,xend,xstep}];
incoming = Transpose[{newx,Transpose[incoming][[1]]}];
, (* else, xvals are specified in file *)
(* change these values if your data are not in the standard columns *)
(* set errorcol to 0 if you don't want to use any uncertainty data *)
xcol=1;
ycol=2;
errorcol=3;
If[Length[Transpose[incoming]]< 3||errorcol==0,
incoming=Transpose[{Transpose[incoming][[xcol]],Transpose[incoming][[ycol]]}]//Sort
, (* else *)
incoming=Transpose[{Transpose[incoming][[xcol]],Transpose[incoming][[ycol]],Transpose[incoming][[errorcol]]}]//Sort
];
];

Define the theoretical model curve and plot the initial guess.

Modify the theory function as desired, but keep x as your independent variable.  Modify params to reflect your fitting-parameter names.  Place an initial guess for each fitting parameter in the curly brackets with each parameter name (eg. "{param,value}").  In general a nonlinear fit will fail unless your initial guess is fairly reasonable.  You may need to adjust the initial values to get a curve that is in the neighborhood of your data.

Experimental data points include error bars if yerror is specified in the input data file.  They can be supressed by setting showerror equal to False.  Set logx or logy equal to True to get a log scale on the repsective horizontal or vertical axis.  If column 3 in the file is the estimated errors in the column 2 values, setting Col3AreErrors to True will correctly treat those values as weights in the fit.

logx=False;logy=False;showerror=True;
Col3AreErrors=True;
(* If you want labels on the axes of your graph you can enter them in these two variables *)
XLabel=" ";
YLabel=" ";
If[Length[incoming]>2,
theory=a+b*x^c;
  params={{a,2},{b,-1},{c,0.5}};
  fitresplot[incoming,theory/.((#[[1]]->#[[2]]&)/@params)," ",True,XLabel,YLabel,12,logx,logy,showerror],
(* This is for a straight line fit to two points *)
{theory,a,b,σa,σb}=twopointfit[incoming,Col3AreErrors];
]

Perform Non-Linear Least Squares fit and plot it together with the data

(* Enter the value False on the next line if you do not want the fitting equation shown on the plot - if you are putting the plot into a document you generally don't want the equation shown since it will be in the caption and/or the body of the text *)
ShowEquation=True;
(* Enter the desired font size in FontPoints.  On a plot for publication the typical final size is a width of 3 3/8 inches (the width of a single column).  At the final width the text must be at least 8 pt size to be considered legible. *)
FontPoints=12;
If[Length[incoming]>2,
fitinformation = nicefit[incoming,theory,params,Col3AreErrors,{x}];
fittedfunction=fitinformation[[1]];
nlm=fitinformation[[3]];
ptable=nlm["ParameterTable"];
Print[Flatten[ptable]];
ptable1=nlm["ParameterConfidenceIntervalTable",ConfidenceLevel->0.6827];
ptable2=nlm["ParameterConfidenceIntervalTable",ConfidenceLevel->0.9545];
Print["68.27% Confidence Interval (1σ):\n",Flatten[ptable1],"\n95.45% Confidence Interval (2σ):\n",Flatten[ptable2]];
rsquared=StringJoin[", (R^2)=",ToString[fitinformation[[2]]]];,
(* This is for a straight line fit to two points *)
fittedfunction=theory;
Print[Grid[{{" ","Estimate","Standard Error"},{"a",a,σa},{"b",b,σb}}]];
rsquared=" ";
]
fitresplot[incoming,fittedfunction,rsquared,ShowEquation,XLabel,YLabel,FontPoints,logx,logy,showerror]

Troubleshooting

Students commonly experience the following pitfalls when using NLSQfit240.nb for the first time.
- used a modified copy of NLSQfit240.nb leftover from previous students instead of starting with a fresh copy 
- did not copy NLSQfit240.nb into the same folder as their data
- tried to proceed without first evaluating the initialization cell
- formatted the data file incorrectly (columns must be separated by white space with one set of values/line)
- didn't use square brackets when calling square roots (e.g. Sqrt[x]) and other pre-defined functions in Mathematica
- attempted to use a log-scale plot for data with negative values.