Nonlinear Fitting

One of the fundamental tasks in experimental physics is to compare theory with data. This comparison is often done by fitting a theory-motivated curve to data. In other words, some mathematical function produced by the physical theory has unknown parameters that can be varied so the function more closely resembles the data. In this laboratory you will be given a data set and develop such a procedure to find a best fit curve to the data.

1. Fit Commands

This section demonstrates the commands needed to fit a  function to a data. First get some data (see the same place you downloaded this notebook).  The data here consist of ordered pairs of times in seconds and the number of counts in a Geiger counter. 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. Note the double back slashes. *)
data1 = ReadList["C:\\Users\\ggilfoyl\\Desktop\\Nuke.dat",{Number,Number}];

Convert the time from seconds to minutes and create a new data array. For reasons yet to be understood this has to be done or the fit to the raw data array (data1) fails.

data2 = Table[{data1[[i,1]]/60,data1[[i,2]]},{i,1,Length[data1]}];

Fit the data to an exponential function of the form

f(t) = (Subscript[N, 0])(e^(-λt))

where (Subscript[N, 0])is a scale factor and λ is the decay constant and then print out the results. See the online Mathematica help by clicking on Help in the menu bar.

myfit=NonlinearModelFit[data2,N0*Exp[-lambda*xt],{N0,lambda},xt]; (* Does the fit. *)
myfit["EstimatedVariance"] (* goodness of fit. *)
myfit["ParameterTable"] (* Prints fit parameters and uncertainties. *)

2. Plotting the Results

To display the data2 array and a curve on a single plot use the commands shown below. Notice the semicolon at the end of the line turns off plotting for the first two plots that have the data and curve. The last, combined-curve-and-data plot is displayed with the Show command.

plot1 = ListPlot[data2,
PlotStyle->{RGBColor[0,0,1],PointSize[0.015]}];
plot2 = Plot[myfit[x],{x,0,60},
PlotStyle->{Thickness[0.02],RGBColor[1,0,0]},
PlotRange->{0,5000},
Frame->True,FrameStyle->Black,
FrameLabel->{"time (min)","Counts in 10 seconds","Nuclear Decay",""},
BaseStyle->FontSize->16];
Show[{plot2,plot1}]

3. Plotting the Results with Uncertainties

We need a special software package to plot the data with uncertainties. To get that package enter the following command once and only before using the ErrorListPlot command.

Create a new data array that includes the uncertainties on each point.

data3 = Table[{data2[[i,1]],Around[data2[[i,2]],(Sqrt[data2[[i,2]]])]},{i,1,Length[data2]}];
plot3 = ListPlot[data3,PlotStyle->{RGBColor[0,0,1],PointSize[0.010]}];
Show[{plot2,plot3},ImageSize->8.0*72,PlotRange->{{10,20},{0,500}}]

4. 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/(data2[[i,2]]))],{i,1,Length[data1]}];
myfit=NonlinearModelFit[data2,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. Don't forget to load the necessary packages shown here and in Section 3 above.

plot4=ListLogPlot[data3,PlotRange->{{0,60},{1,10000}},Frame->True,
FrameLabel->{"t(s)","Counts","",""},FrameStyle->Black,
BaseStyle->Large,
ImageSize->8*72,PlotStyle->RGBColor[1,0,0]];
plot5 = LogPlot[myfit[x],{x,0,60},PlotRange->{{0,60},Automatic},Frame->True,
FrameLabel->{"t(s)","Counts","",""},FrameStyle->Black,
BaseStyle->Large,
ImageSize->8*72,PlotStyle->RGBColor[0,0,1]];
Show[plot4,plot5]