More on Nonlinear Fitting
and Plotting

We now consider another test of the goodness of fit of a curve to data in the case of a non-linear function that uses what we know about the uncertainties on individual data points. We also introduce the use of the chi-square or  (χ^2) as a measure of goodness of fit and learn why it has advantages over the estimated variance we have used so far.

1. Fit Commands

This section demonstrates the commands needed to fit a  non-linear function to a data. First get the necessary packages. Enter this command once and only before using any other fitting commands.

Needs["ErrorBarPlots`"]

Now get the data from the course website and load it in. We are using the same radioactive decay data used in the notebook entitled fitting1.nb. You have to make sure you set the full path correctly.

(* on a unix machine. *)
data1 = ReadList["/home/gilfoyle/tmp/Nuke.dat",{Number,Number}];

(* on a windows machine desktop. *)
data1 = ReadList["C:\\Documents and Settings\\ggilfoyl\\Desktop\\Nuke.dat",{Number,Number}]

The command below is the same one used for non-linear fitting from the notebook NonlinearFitting70.nb. We demonstrate it here to get ourselves started.

myfit=NonlinearModelFit[data1,N0*Exp[-lambda*x],{N0,lambda},x]
myfit["EstimatedVariance"]
myfit["ParameterTable"]

2. The  (χ^2)

The   (χ^2) is a measure of goodness of fit and is  defined as

(χ^2) = (∑_(i=1)^N)(((((Subscript[y, i])- f((Subscript[x, i]))))^2)/((Subscript[σ, i])^2))

where the sum is over the data, N is the number of data points, (Subscript[σ, i])is the uncertainty on each data point, (Subscript[y, i]) are the data, and f((Subscript[x, i])) is some function representing a physics model. In order to use this method in our fitting procedure we apply a different weight to each data point. Since the data we are considering comes from radioactive decay the uncertainty on each data point to the statistical uncertainty (Sqrt[M])where M is the number of counts we have collected for a particular time bin. The NonlinearModelFit command with weights is demonstrated in the next command.

wghts = Table[N[(1/(data1[[i,2]]))],{i,1,Length[data1]}];
myfit=NonlinearModelFit[data1,N0*Exp[-lambda*xt],{N0,lambda},xt,Weights->wghts]
myfit["EstimatedVariance"]
myfit["ParameterTable"]

If our data are statistically distributed, then the reduced  (χ^2) should be equal to about one. The reduced  (χ^2) is defined as

(Subscript[(χ^2) , reduced])= (χ^2)/ν  =(1/(N-ndf)) (∑_(i=1)^N)(((((Subscript[y, i])- f((Subscript[x, i]))))^2)/((Subscript[σ, i])^2))

where N is the number of data points again and ndf is the number of degrees of freedom in the fit. The function used in the fit above had two degrees of freedom because there were two parameters that were fitted in the minimization. The quantity EstimatedVariance listed above is the reduced (χ^2) when the proper weighting is used.

3. Semi-log plots

Now make some cool, semi-log plots of the data with uncertainties and the fit. These plots are very useful when your data covers a wide range on the vertical axis like our data here. The first command creates an array which includes an element for the uncertainty on each point. Note the additional curly brackets around the first two components in the row of the array. Don't forget to load the necessary packages discussed in the Section 1 above.

data1b=Table[{{data1[[i,1]],Log[10,data1[[i,2]]]},ErrorBar[Log[10,data1[[i,2]]]-Log[10,data1[[i,2]]-(Sqrt[data1[[i,2]]])],Log[10,data1[[i,2]]+(Sqrt[data1[[i,2]]])]-Log[10,data1[[i,2]]]]},{i,1,Length[data1]}];ErrorListPlot[data1b,PlotRange->{{0,4000},{-3,4}},Frame->True,
FrameLabel->{"t(s)","Counts","",""},
BaseStyle->Large,
ImageSize->8*72,PlotStyle->RGBColor[1,0,0]]