Chapter 9: Some Data Analysis

9.1 Least squares curve fitting

Remove["Global`*"]

BallData={{0,1.302},{0.03333,1.411},{0.06667,1.5},{0.1,1.578},{0.1333,1.646},{0.1667,1.703},{0.2,1.745},{0.2333,1.781},{0.2667,1.807},{0.3,1.828},{0.3333,1.818},{0.3667,1.818},{0.4,1.807},{0.4333,1.776},{0.4667,1.734},{0.5,1.682},{0.5333,1.63},{0.567,1.552},{0.6,1.469},{0.6333,1.37},{0.667,1.266},{0.7,1.151},{0.733,1.026},{0.7667,0.875},{0.8,0.719},{0.8333,0.557},{0.867,0.385},{0.9,0.193},{0.9333,0.005}};

BallFig=ListPlot[BallData,PlotStyle->PointSize[0.02],Frame->True,FrameLabel->{"Time (s)","Position (m)"},PlotLabel->"Ball thrown upward"]

In experimental Physics curve fitting is an important statistical tool for analysing data and quantifying correlations between variables. The command Fit finds the parameters of a function by doing a best fit, in the least squares sense, of the function on a set of data. For example, we can fit a quadratic to the ball data like this:

BallFit=Fit[BallData,{1,x,(x^2)},x]

BallFitFig=Plot[BallFit,{x,0.,1}]

Show[BallFig,BallFitFig]

Here's another example with a linear fit and a set of artificial data

Remove["Global`*"]

Data={{0,-1},{1,2},{3,7},{4,6},{7,9},{11,14}};

DataFit=Fit[Data,{1,x},x]

DataFig=ListPlot[Data,PlotStyle->PointSize[0.02]]

FitFig=Plot[DataFit,{x,0,11}]

Show[DataFig,FitFig,AxesLabel->{"x","y"}]

As a reminder of linear algebra behind the fitting procedure, let's get the best fit parameters "by hand"

MatrixForm[Data]

t=Data[[All,1]];
MatrixForm[t]

y=Data[[All,2]];
MatrixForm[y]

The first column of our data is time t, the second is our measured quantity, y.  We have a model that says y and t are related through y= C +D t.  If this equation held for each data entry, i.e. each row of our data matrix, we could solve the equation
A x = y
where the vector y is given above, x= (({{C}, {D}})) and the matrix A is given by

A=Table[{1,t[[i]]},{i,1,Length[t]}];
MatrixForm[A]

This matrix equation represents, in this case, 6 equations and 2 unknowns, C and D.  Unless the data prefectly match our model, we can't solve such a system.  Instead, we chose to find C and D that minimize the square of the errors given by (Subscript[y, i])-C-D (Subscript[t, i]).  The result is we now have to solve the system of equations given by
(A^T)A x= (A^T) y,
where (A^T) is the transpose of A.  The system now is square, with the number of equations equal to the number of unknowns.  Formally, the solution is given by
x=((((A^T)A))^(-1))(A^T) y.
Getting Mathematica to do this, we obtain

Inverse[Transpose[A].A].Transpose[A].y//N

the same parameters we obtained from the more efficient built-in function Fit.

9.2 Data analysis: Least Square Curve Fitting with Statistics

The Fit command does not give any information on the statistics of the fit (i.e. how good the fit is). There is a more complete command called LinearModelFit that does basically the same as Fit except that statistics are also returned.

Remove["Global`*"]

BallData={{0,1.302},{0.03333,1.411},{0.06667,1.5},{0.1,1.578},{0.1333,1.646},{0.1667,1.703},{0.2,1.745},{0.2333,1.781},{0.2667,1.807},{0.3,1.828},{0.3333,1.818},{0.3667,1.818},{0.4,1.807},{0.4333,1.776},{0.4667,1.734},{0.5,1.682},{0.5333,1.63},{0.567,1.552},{0.6,1.469},{0.6333,1.37},{0.667,1.266},{0.7,1.151},{0.733,1.026},{0.7667,0.875},{0.8,0.719},{0.8333,0.557},{0.867,0.385},{0.9,0.193},{0.9333,0.005}};

lm=LinearModelFit[BallData,{1,x,(x^2)},x];

Here are some examples of getting access the information pertaining to various aspects of the fit.

lm["BestFit"]

lm["BestFitParameters"]

lm["ParameterErrors"]

lm["ParameterTable"]

lm["RSquared"]

9.3 Nonlinear Curve Fitting

In the sections above, we were able to use matrix algebra to find our fit parameters uniquely, i.e. we could solve a linear system of equations to find our parameters.  This is possible because the fitting function depended linearly on our parameters.  Let's bring up our small data set again.

Remove["Global`*"]

Data={{0,-1},{1,2},{3,7},{4,6},{7,9},{11,14}};

n1=ListPlot[Data]

We fit a function y=A+Bx to the data.  We could also have chosen to do a linear fit with a more complicated function of x, like y= A + B sin(x)+ C exp(x) by entering the following.

DataFit=Fit[Data,{1,Sin[x],Exp[x]},x]

This works because our fitting function parameters still appear linearly.  We could not use, say, y= A Sin[B ((x])^2), because the parameter B appears in a nonlinear way in the fitting function.  To do a nonlinear fit we must use different routines.

nlFit=FindFit[Data,a ((Sin[b x])^2),{{a,1},{b,0.01}},x]

Show[n1,Plot[a ((Sin[b x])^2)/.nlFit,{x,0,10}]]

Note that in calling FindFit, you are able to provide inital guesses for the parameters.  In general it is difficult to find the best fit in the nonlinear case.  Rather, routines find the best fit near the initial guess.  Different initial guesses may yield different fit parameters.  We won't go into the methods of finding parameters for nonlinear fits.

As with linear fits (or regression as it is also called), we can get information of the fit itself, the time using a function called NonlinearModelFit.

nonlinm=NonlinearModelFit[Data,a ((Sin[b x])^2),{{a,1},{b,0.01}},x]

Here are some examples of getting access the information pertaining to various aspects of the fit.

nonlinm["BestFit"]

nonlinm["BestFitParameters"]

nonlinm["ParameterErrors"]

nonlinm["ParameterTable"]

nonlinm["RSquared"]

9.4 Application of Fitting

There is a data file that you should download along with this assignment called falling_object.dat that has two columns.  The first column contains time (s) and the second contains height (m).

a) Load the data into a table and plot the height as a function of time.

Remove["Global`*"]

Directory[]

SetDirectory["~/Desktop"];

FileNames["f*"]

xt=Import["falling_object.dat","Table"];

l1=ListPlot[xt,AxesLabel->{"t(s)","y(m)"}]

b) Using the centred difference scheme, calculate and plot the velocity (do not interpolate) as a function of time.  You should see the velocity approach a terminal value.

firstD[data_,i_]:=((data[[i+1,2]]-data[[i-1,2]])/(data[[i+1,1]]-data[[i-1,1]]))

vt=Table[{xt[[i,1]],firstD[xt,i]},{i,2,Length[xt]-1}];

ListPlot[vt,AxesLabel->{"t(s)","v(m/s)"},Joined->True]

c) Calculate and plot the acceleration of the falling object directly from the height data by using the centred scheme for the second derivative.  You should see the acceleration approach zero.  The acceleration graph looks noisy. This is because the original measurements contain some uncertainty and random noise. This noise gets amplified by taking derivatives.

secondD[data_,i_]:=((data[[i+1,2]]-2data[[i,2]]+data[[i-1,2]])/(((data[[i+1,1]]-data[[i,1]]))^2))

at=Table[{xt[[i,1]],secondD[xt,i]},{i,2,Length[xt]-1}];

ListPlot[at,AxesLabel->{"t(s)","a(m/(s^2))"},Joined->True]

d)  Let's assume that the object experiences a drag force of F= - b v.  We will attempt to find b.  As usual, we'll begin with Newton's 2nd illustrious Law
F=ma
ma = -bv -mg

Aha!
If we plot the quantity ma as a function of v, the graph should be a straight line (with noise) with interecept -mg and having slope equal to -b.

Construct a table with all the v values in the first column and m a  in the second.  Values of a and v in the same row should correspond to the same time.
(The mass of the object is 0.2 kg, while acceleration due to gravity is 9.81 m/s/s.  These values should jive with our intercept value. )

m=0.2;
g=9.81;
matxt=Table[{vt[[i,2]],m at[[i,2]]},{i,1,Length[at]}];

p1=ListPlot[matxt,AxesLabel->{"-v (m/s)","ma + mg (m/(s^2))"}]

e) Make a linear fit for the data in part (d) using Fit command, and then find the value of the constant b.

linfit[x_]=Fit[matxt,{1,x},x]

f)Show both the fit line and the discrete data in single graph.

Show[p1,Plot[linfit[x],{x,-16,0},PlotStyle->{Thick,Red}]]

g) Use LinearModelFit with useful corresponding commands (ParameterTable, BestFitParameters and ParameterErrors) to obtain more information about the fits.

lm=LinearModelFit[matxt,{1,x},x]

lm["ParameterTable"]

lm["BestFitParameters"]

lm["ParameterErrors"]

The coefficient of linear air resistance is

b=0.127±0.002 kg/(s^2).

Let’s compare the intercept (-1.96 ± 0.03) value with what we expect it to be (-mg).

-m g

Golden!

h)  We're not done.  Newton's equation m a = -b v - m g can be solved analytically.  Use DSolve to find the analytical solution for the height as a function of time.  Use the values of m and g above.  Leave b, initial height and initial velocity as unknown parameters.

sol=DSolve[{m y''[t]==-b y'[t]-m g,y[0]==h,y'[0]==v0},y[t],t]//First

solution=sol[[1,2]]

j)  Use  NonlinearModelFit to find the values of b, initial height and initial velocity.

nonlinm=NonlinearModelFit[xt,solution,{{b,2},{h,10},{v0,-0.6}},t]

nonlinm["BestFitParameters"]

nonlinm["ParameterErrors"]

nonlinm["ParameterTable"]

k)  Now use the function FindFit in a similar way to NonlinearRegress in order to plot the fitting function.  Plot this fitting function over the data.  Does it seem like our model (F=-b v - m g) describes the data?

itdata=FindFit[xt,solution,{{b,2},{h,100},{v0,-0.6}},t]

p1=Plot[solution/.itdata,{t,xt[[1,1]],xt[[Length[xt],1]]},PlotStyle->{Red}]

Show[l1,p1]

9.5 Statistics

For this section and the next, you will need to download the file TUP.dat, let's say to the Desktop.   The data is from a simulation of a collection of atoms.  The columns are time, pressure and potential energy.

Remove["Global`*"]

First let's find out where we are in the directory structure.

Directory[]

Let's move ourselves to the Desktop

SetDirectory["~/Desktop"]

Let's see if TUP.dat is here

FileNames["TU*"]

Import the data into a table.  There are 3 columns: time, potential energy and pressure.

TiEnPr=Import["TUP.dat","Table"];

An abbreviated look at the data

Short[TiEnPr]

The length of the data set

L=Length[TiEnPr]

We want to make a table with two colums: time and energy.  Here's one way to do it.

EnergyTime1=Table[{TiEnPr[[i,1]],TiEnPr[[i,2]]},{i,1,L}];

or just

EnergyTime=TiEnPr[[All,1;;2]];

Let's check that the two results are the same

EnergyTime==EnergyTime1

Similarly for the pressure.

PressTime1=Table[{TiEnPr[[i,1]],TiEnPr[[i,3]]},{i,1,L}];

PressTime=TiEnPr[[All,{1,3}]];

PressTime==PressTime1

Plot the energy time series.

p1=ListPlot[EnergyTime]

Plot the pressure time series.

pp=ListPlot[PressTime]

Now we want a statistical description of the data.  We'll first the energy data into a separate list.

Energy=TiEnPr[[All,2]];

We can work out the mean energy by hand...

M=Sum[Energy[[i]],{i,1,Length[Energy]}]/Length[Energy]

Or use a built-in function

Mean[Energy]

or on the fly using the original data list

Mean[TiEnPr[[All,2]]]

Similarly for the standard deviation.

SD=Sqrt[Sum[(((Energy[[i]]-M))^2),{i,1,Length[Energy]}]/(Length[Energy]-1)]

StandardDeviation[Energy]

Minimum energy

Min[Energy]

Maximum energy

Max[Energy]

Range in energies

Max[Energy]-Min[Energy]

After the standard deviation, which is basically the width of the distribution, the next moment is the skewness... or how asymmetric the energy distribution is.

Skewness[Energy]

A plot of the distribution of energies, i.e. how many data points have an energy with a given (narrow) range.

Histogram[Energy]

We see that the energy is actually bimodal.  Le's replot the energy time series.

ListPlot[EnergyTime,Joined->True,PlotRange->{{390000,1200000},{-5.8`,-5.2`}}]

Is the pressure bimodal?

Press=TiEnPr[[All,3]];

Histogram[Press]

Not really.  Let's plot the energy as a function of pressure to see if there's a strong correlation.  First we make an appropriate table.

EnergyPress=TiEnPr[[All,{3,2}]];

ListPlot[EnergyPress,PlotRange-> {{-0.05,0.07},{-5.8,-5.2}}]

9.6 Filtering

Sometimes data can be a little noisy and it helps to smooth it out.  Smoothing out the data let's you see trends that might otherwise be obscured by the noise.

One way to do this is by using moving averages.  Mathematica's built in function MovingAverages can do the trick.

?MovingAverage

Here, we will take the average of 500 data points at a time in order to smooth out the energy and pressure time series.

smoothET=MovingAverage[EnergyTime,500];

i.e. the first entry of smoothET is the average of data points 1 to 500, while the second entry is the average of points 2 to 501 etc.

smoothPT=MovingAverage[PressTime,500];

p2=ListPlot[smoothET,PlotStyle->{Red}]

Show[p1,p2,PlotRange->{{390000,1200000},{-5.8,-5.2}}]

ListPlot[smoothPT,PlotStyle->{Red}]

Show[pp,ListPlot[smoothPT,PlotStyle->{Red}]]

A slightly different way to average data in order to discern trends is to "Decimate" the data, or break it up into a given number of block.  An average is reported for each block.  The difference with the moving average is that for a a moving average, two consecutive points in the smoothed data table are averages of almost the same underlying raw data.  In our example above, smooth_1 averaged points 1 to 500, smooth_2 averaged points 2 to 501.  With blocking, the averages are reported for separate blocks of original data.  Let's look at an example.

time=TiEnPr[[All,1]];

L=Length[Energy]

Our energy data set has 3947 elements

nblocks=10

and we want to break it up into 10 blocks.  For each block, we will calculate the average.  The loop below does this.

Do[
jstart=Round[L (i-1)/nblocks]+1;
jstop=Round[L  i /nblocks];
nj=jstop-jstart +1;
Print[jstart," ",jstop," ",nj];
DecimatePress[i]=Sum[Press[[j]],{j,jstart,jstop}]/nj;
DecimateEnergy[i]=Sum[Energy[[j]],{j,jstart,jstop}]/nj;
DecimateTime[i]=Sum[time[[j]],{j,jstart,jstop}]/nj;
,{i,1,nblocks}]

Each block has about 395 points.  Let's see what's in our blocked data.

?DecimateEnergy

Plot blocked energy versus blocked time

DecEnT=Table[{DecimateTime[i],DecimateEnergy[i]},{i,1,nblocks}]

p3=ListPlot[DecEnT,PlotStyle->{Blue},Joined->True]

Comparing moving averages and blocking...

Show[p2,p3]

9.7 Plotting

Here is an example of producing a plot more or less fit for publication.  First we generate some data.

Remove["Global`*"]

d1={1,2,2.5};
d2={2,2.3,5};

If we don't need a legend, we can just make the plot.

p1=ListPlot[{d1,d2,d2},
Frame->True,
FrameLabel->{"T","I"},
LabelStyle->{Directive[20],Black},
Joined->{True,False,True},PlotStyle->{{Thickness[0.005],Red},{Blue,PointSize[0.03]},{Dashed,Blue}},PlotRange->{{0,4},{0,6}},
AspectRatio->1/1.2,
ImageSize->{600,500}]

The are many options that can help tailor the appearance of a graph.  Entering ??ListPlot will list them, and the help section is actually helpful.  Getting plots to look just right can be tricky in Mathematica.

To include a legend, we can do the following.

p1leg=ListPlot[{d1,d2,d2},
Frame->True,
FrameLabel->{"T","I"},
LabelStyle->{Directive[25],Black},
Joined->{True,False,True},PlotStyle->{{Thickness[0.005],Red},{Blue,PointSize[0.03]},{Dashed,Blue}},PlotRange->{{0,4},{0,6}},
AspectRatio->1/1.2,
ImageSize->{600,500},
PlotLegends->Placed[LineLegend[{"Dataset 1","Dataset 2",None},LabelStyle->Directive[20,Brown]],{0.18,0.88}]]

Let's write our graphs to appropriate files.

?Export

Export["plot1.eps",p1,"eps"]

Export["plot1Legend.pdf",p1leg,"pdf"]

Directory[]

FileNames["plot*"]

9.8 Chi-squared Test

Often, one wishes to compare theory against experimental data.  One way to quantify how well the two match is the chi-squared test.  As an introduction to this topic, we will go through an example.

Remove["Global`*"]

Let's pretend we have some data.  The data is in the form of a histogram, where each column tells how many people we measured have a certain height, i.e. measure the heights of a few hundred people to the nearest inch, then count up the number of people that have the same height.

Instead of measuring a bunch of people, let's construct a data set.  We will use a normal distribution (Gaussian) with an average at 70 inches (5' 10'') and a standard deviation of 6 inches.  We will have about 300 people.  To lend some realism, we'll add some randomness to the data.

σ=6;
have=70;
Np=200;
bin=1;

?Random

height=Table[{i,Max[0,Round[((Np Exp[-((((i-have))^2)/(2 (σ^2)))])/(σ (Sqrt[2 π])))+RandomReal[NormalDistribution[0,1]]]]},{i,54,86,bin}]

So now we have a table.  The 1st column of each entry gives a height in inches, the second column gives the number of people with that height.  Let's plot our "data" .

l1=ListPlot[height,PlotRange->{{30,100},{0,Np/10}},PlotStyle->PointSize[0.02]]

It still looks roughly Gaussian.  If this were real data and we had reason to believe that heights should be normally distributed, we would fit a Gaussian to the data.

But first, we want to know how many different heights we have and how many people we've got in our sample.

Ndat=Length[height]

Ne=Sum[height[[i,2]],{i,1,Ndat}]

Now let's use the nonlinear fitting routine to fit a normal distribution to the data.  For the nonlinear case, we can provide an inital guess for our parameters.

Fitfunc[x_]:=((Ne bin)/(A Sqrt[2π]))Exp[((-(((x-B))^2))/(2 (A^2)))]
dis[x_]=Fitfunc[x]/.FindFit[height,Fitfunc[x],{{A,6},{B,70}},x]

Let's see how the fit did.

Show[l1,Plot[dis[x],{x,50,90},PlotStyle->{Red,Thick}]]

And now for the (χ^2) statistic.

(χ^2)= (∑_(i=1)^(ALL BINS))((((Observed (Subscript[Number , i])- Expected (Subscript[Number, i])))^2)/(Expected (Subscript[Number, i])))

For example, let's start by looking at the first piece of our "data"

height

This says that for a height of 57 inches, we measured

height[[4,2]]

people.  Using our fitting function (our model for the distribution), we expect

dis[57]

people.  A fraction for the expected number is okay, even necessary.  So the contribution to (χ^2) from this bin of our histogram is

((((height[[4,2]]-dis[57]))^2)/(dis[57]))

Summing up over all the bins, we get our value for (χ^2)

chi2=Sum[(((height[[i,2]]-dis[height[[i,1]]]))^2)/dis[height[[i,1]]],{i,1,Ndat}]

What does this number tell us?  Not much yet.  However, recall that the distribution of Chi-Squared values follows

Y[x2_,k_]:=(1/((2^(k/2)) Gamma[k/2]))(x2^(k/2-1))Exp[-x2/2]

where k is the number of degrees of freedom (k = the number of bins minus 1... because the number of counts in each bin add up to the total number of counts, i.e. all the bins are not independent.)

Plot[Y[x,Ndat-1],{x,0,100},AxesLabel->{"(χ^2) value","Probability density"}]

The area between, say 20 and 30, gives the probability of getting a Chi-squared value in that range (assuming the contribution to Chi-Squared from each height bin has a standard deviation of unity).

NIntegrate[Y[x,Ndat-1],{x,20,30}]

Thus, the probability of getting a chi-squared value of chi2 or higher from random variations in sampling is

NIntegrate[Y[x,Ndat-1],{x,chi2,∞}]

This integral has an analytic form.

Pr[k_,x2_]=Gamma[(k/2),(x2/2)]/Gamma[(k/2)]

where k is still the number of degrees of freedom (nbins - 1), and x2 is the obtained value of (χ^2).

Gamma is the function for the various gamma functions.

?Gamma

In out case n=Ndat, so our probability is

Pr[Ndat-1,chi2]

This says that the probability that the deviation from the model distribution arise from chance, and not some underlying physics, is very close to 1, i.e. it is very likely that our model describes the data (not too surprising since we constructed it that way).  It is almost 100 percent likely that our model is right and that the diviations we see come from sampling errors.  I.e. the probability that our model is wrong is close to zero.  So we can be confident that our model is correct.

next experiment

Now let's say we were scientists who thought that this height distribution is the same for people everywhere, and we go to a far away place and measure the heights of people with the expectation that our exact same model (including the model parameters) describes people there.

We will construct another hypothetical data set that is slightly different from the original one.  This data set represents the height distribution of people in the far away place.

σ=6;
have=70;
Np=303;
Ntall=0.1 Np;
st=2;
ht=74;
Nshort=0.1 Np;
ss=2;
hs=62;

height2=Table[{i,Max[0,Round[((Np bin Exp[-((((i-have))^2)/(2 (σ^2)))])/(σ (Sqrt[2 π])))+RandomReal[NormalDistribution[0,1]]+((Ntall bin Exp[-((((i-ht))^2)/(2 (st^2)))])/(st (Sqrt[2 π])))-((Ntall bin Exp[-((((i-hs))^2)/(2 (ss^2)))])/(ss (Sqrt[2 π])))]]},{i,54,86,bin}];
l2=ListPlot[height2,PlotRange->{{30,100},{0,25}},PlotStyle->PointSize[0.02]]

Let's perfrom the chi-squared analysis using the same fitting function as before, only rescaling because our sample size is different.

Ndat2=Length[height2]
Ne2=Sum[height2[[i,2]],{i,1,Ndat2}]
dis2[x_]=dis[x] Ne2/Ne

Show[ListPlot[height2,PlotRange->{{30,100},All},PlotStyle->PointSize[0.02]],Plot[dis2[x],{x,50,90},PlotStyle->{Red,Thick}]]

chi2=Sum[(((height2[[i,2]]-dis2[height2[[i,1]]]))^2)/dis2[height2[[i,1]]],{i,1,Ndat2}]

Pr[Ndat2,chi2]

This time we are much less confident that our model describes the height distribution.  To totally reject our model, we would seek a probability of 0.05 or less, or a chi-squared value of

xupper=x/.FindRoot[Pr[Ndat2,x]==0.05,{x,3}]

or more.

Since our probability is in between 0.05 and 0.95, we won't confidently make any conclusions about our model.