Physics 218:  Fitting a Normal Distribution
to the Pendulum Data

Introduction

In this notebook, you will see how to import your data, create a histogram, and fit the normal distribution to your data.  I have used this notebook to introduce a number of powerful Mathematica functions that I hope will be useful beyond this experiment.  Please make extensive use of Mathematica’s built-in Help function to learn more about any of the functions used here.  Also feel free to add in extra steps where you print out various arrays to see what is happening.

For your final submitted work, however, please trim out all unnecessary or duplicate printouts.  Make it concise, but accurate and readable.

Importing Data

You need to tell Mathematica where to find your data file.  If you store your notebook in the same directory as your data, you can use the following trick to set Mathematica’s default directory.  Otherwise, you have to either give the full pathname or give the pathname relative to Mathematica’s default.

Note the trailing ‘;’ on each line.  If you omit them, Mathematica prints out the result of the command.  Try it and see.

SetDirectory[NotebookDirectory[]];
rawdata = Import["Pendulum-20220210-2.csv", "CSV"];

The data is in columns.  Here are the first 10 lines of the file.  The first line of the file gives the column titles.  We want the third column, but only the numeric entries.

rawdata[[Range[10]]] // TableForm

There are many ways to pick out specific array elements.  One compact way is to use the [[All]] subscript notation.  The notation rawdata[[All, 3]] means take the 3rd item in every rawdata element.  The Select[ ] command is a powerful way to select data from an array that satisfy a particular criterion.  The command loops through the array, setting ‘#’ equal to each element.  In this invocation, we only select ones which are “Numbers”.  (You have to remember to put the ‘&’ at the end inside the final square bracket.)

data = Select[rawdata[[All,3]], NumberQ[#] &];

data[[Range[10]]] (* The first 10 points *)

We can find out how many data points we have by finding the length of the array.

Length[data]

Eliminating Bad Points

ListPlot[data, PlotRange ->All, Frame -> True, LabelStyle->Larger, FrameLabel->{"cycle", "Period (s)"}, ImageSize -> Scaled[0.7]]

If there are any outlier points, they should be easy to spot here.  Most likely, someone walked past the apparatus.  More generally, it is often reasonable to compute the spread of the data and throw out spurious points that are many standard deviations away from the mean.  Don’t do this unless you have good reason to suspect those points are bad, however! 

Let’s make a list of the minimum, maximum, average, and standard deviation for our data.

{min = Min[data], max = Max[data], avg = Mean[data], σ = StandardDeviation[data]}

How far away from the mean are the min and max points?   Compare them to the standard deviation:

(max - avg)/σ

(min - avg)/σ

The maximum point is within 5 σ of the mean, so we could make a case for deleting it, but we can also just keep all the points.The minimum point is slightly farther away -- 6σ.  We could probably get rid of it, or we could keep it.  It won’t make much statistical difference either way.  Finally, in this data set, the data starts to get erratic after about 600 cycles.  The Select[ ] command below 
considers only the first 600 data points and then selects only those points within 5 σ of the mean.  (You can adjust the ‘600’ number as appropriate for a particular data set.)

Length[data]
data = Select[data[[Range[600]]], Abs[# - avg] < 5σ&];
Length[data]

Now compute our statistics.

{min = Min[data], max = Max[data], avg = Mean[data], σ= StandardDeviation[data], δT = σ/Sqrt[Length[data]]}

This is our best representation of the good data.

ListPlot[data, PlotRange -> All, Frame -> True, LabelStyle->Larger, FrameLabel->{"cycle", "Period (s)"}, ImageSize -> Scaled[0.7]]

Making a Histogram

Mathematica can make a nice histogram.  It often (but not always) chooses reasonable bins.  For this data, I have overridden the default and picked a bin size of 0.002s.

Histogram[data, {0.0002}, Frame -> True, LabelStyle->Larger, FrameLabel->{"Period (s)", "Number of Measurements"}, ImageSize -> Scaled[0.7]]

Getting a manual list of the histogram numbers is a bit more tedious.

hlist = HistogramList[data, {0.0002}];

This actually gives two separate lists, first is a list of bin boundaries, and second is the list of heights.  Interleaving them is straightforward, if verbose.  We are assigning each height to the center point of each bin.  Do that by constructing a new table.

hlist[[1]]

hlist[[2]]

hist = Table[{(hlist[[1, i]] + hlist[[1, i+1]])/2., hlist[[2, i]]}, {i, 1, Length[hlist[[1]]]-1}];

dataplot = ListPlot[hist, PlotRange -> All, Frame -> True, LabelStyle->Larger, FrameLabel->{ "Period (s)", "Number"}, ImageSize -> Scaled[0.7]]

Now fit to the (unnormalized) normal distribution.

norm[A_, x_, μ_, σ_] := A Exp[-(((x - μ))^2)/(2 (σ^2))]

The NonlinearModelFit function works best if you give it a reasonable guess for the parameters.  Fortunately, we think we know the mean and standard deviation pretty well, so just use those as the initial guesses.

gaussfit = NonlinearModelFit[hist, norm[A, x, μ, σfit], {{A, 10}, {μ, avg}, {σfit, σ}}, x]

gaussfit["BestFitParameters"] (* The sign of σfit doesn't matter since it's squared. *)

gaussfit["ParameterConfidenceIntervalTable"]

(* Compare statistics values with values from the curve fit. *)TableForm[{{avg, μ, avg - μ}, {σ, Abs[σfit], σ - Abs[σfit]}}  /. gaussfit["BestFitParameters"], TableHeadings->{{"Mean", "StandardDeviation"}, {"Statistics", "Fit", "Difference"}}]

min = Min[data]; max = Max[data]; 
Show[{dataplot,
Plot[gaussfit[x], {x, min, max}, PlotStyle -> Red, PlotRange -> All]},
 Frame -> True, LabelStyle->Larger, FrameLabel->{ "Period (s)", "Number"}, ImageSize -> Scaled[0.7]]