Fitting Data with Mathematica

In a number of experiments this semester, you will be asked to fit experimental data to a theoretical curve.  This notebook walks through an example of doing that in Mathematica.

For this notebook, we will look at fitting a simple oscillating function to some data.

Entering Data

There are several ways to enter data in Mathematica.  If you only have a small number of data points, you can just enter them directly.  You can use the Palettes -> Classroom Assistant -> Typesetting -> Matrix entry, or you can just enter in the data by hand.

Entering by hand

If you enter it by hand, note that each data point is a list, or pair:  { t, x } .  As always in Mathematica, lists are enclosed in curly braces.  The complete data set is just a list of those data pairs:  {  { t1, x1 }, { t2, x2 } }, etc.  See below for an example data set.  You can also enter data tables with the Matrix tool from the Classroom pallette.

Importing from a file

First, tell Mathematica where the data can be found.  You will have to change this for your own set-up.

SetDirectory["E:/218/2011/lab-svn/labs/torsional"]

If your Mathematica notebook and your data file are in the same directory, this next trick sometimes works.  However, you have to save the notebook first (even if it’s just an empty notebook).

SetDirectory[NotebookDirectory[]]

Importing from a plain text file.

If your data is a simple text file with two columns of numbers, you may import it as a "Table".

data = Import["curvefit-data.txt", "Table"]

Importing from a CSV file.

If you have used LoggerPro and exported the data to a CSV file (Comma-separated-values) then Mathematica can easily import the data.

data = Import["curvefit-data.csv"]

LoggerPro included titles for the two columns of data ("Time" and "Potential", in this case).  You need to remove them from the data before you can use it.  There are many, many ways to do this in Mathematica.  In the on-line help, look for the tutorial on "Getting Pieces of Lists."  (Of course you can also simply edit the original text file as well.)

The most direct is with the Drop[ ] command.  The following drops the first element from the 'data' array and assigns it back to the data array.

data = Drop[data, 1]

Another more general way to select parts of a list is to use the Select[ ] function.  This allows you to select elements from a list that meet some criterion.  See the tutorial on Finding Expressions that Match a Pattern" for more on the Select[ ] function.  The syntax is likely unfamiliar, but the following looks at all the elements of data, and selects only those for which all the elements are numbers. 

The Select function loops through our data and assigns the symbol ‘#’ to each element in turn.   It then selects only the elements for which our test is True.  Since each element is a pair, or list of items, we need to look at both of them.  The function VectorQ does that.  It applies the “NumberQ” test to each of the elements of “#” and returns True only if all of them are numbers. (Be sure to include the '&'.  Look up the online help for Function[ ] to learn more about what it is doing.)

data = Import["curvefit-data.csv"]
data = Select[data, VectorQ[#, NumberQ] &]

Plotting the data

This function plots the data, but also saves a copy of the plot in the variable 'dataplot' for later use.  I find the default labeling to be too small; this command makes it a bit larger.

dataplot = ListPlot[data, PlotRange->All, AxesLabel -> {"t", "x"}, LabelStyle -> Larger, PlotStyle -> Red]

Fitting a Model to the Data

For this exercise, we will try to fit this to a cosine function, allowing for the possibility that the average value might not be zero, due to some offset in the sensor.

x[A_, ω_, ϕ_, t_, xoff_] := A Cos[ω t + ϕ] + xoff

The basic function to use is NonlinearModelFit[ ].  Consult the online help for more details.  This notebook will just illustrate the basic usage.

fit = NonlinearModelFit[data, x[A, ω, ϕ, t, xoff], {A,  ω, ϕ, xoff}, t]

It is also easy to generate a plot from the fit.  This pair of lines generates the plot, gives it a name, and then combines it with the data plot.  The Show[] command shows them both together.  It is easy to make the graph larger -- the ImageSize -> Scaled[0.50] makes the image equal to 1/2 the screen width.  Adjust as you see fit.

fitplot = Plot[fit[t], {t, 0, 3}, PlotRange->All];
Show[{dataplot, fitplot}, ImageSize -> Scaled[0.50]]

This fit is quite poor.  Mathematica performs the fit by starting with an initial guess, and then seeing how the fit improves as it changes the guess.  (LoggerPro works in a similar way.)  You can give it a hand by suggesting a better initial guess for any or all of the parameters.  For example, looking at the graph, it is clear that the period is a little over 1.0 seconds.  Specifying an initial guess for  ω of 2π/1.0 is likely to be closer to the desired result.

fit = NonlinearModelFit[data, x[A, ω, ϕ, t, xoff], {A, { ω, 2π/1.0}, ϕ, xoff}, t]

fitplot = Plot[fit[t], {t, 0, 3}, PlotRange->All];
Show[{dataplot, fitplot}, ImageSize -> Scaled[0.50]]

Fit Diagnostics

NonlinearModelFit returns a variety of measures for how "good" the fit is.    The online help has more details, but three particularly useful ones are shown here.

Best Fit Parameters

fit["BestFitParameters"]

Parameter Confidence Interval

This item gives both the uncertainty and the 95% confidence-level interval.

fit["ParameterConfidenceIntervalTable"]

Correlation Matrix

The "CorrelationMatrix" tells you about the correlations among the fit parameters.  Low values mean the parameters are poorly correlated.  This is good -- it means they are mostly independent.  On the other hand, high values mean the parameters are highly correlated.  This is often a sign of trouble in the proposed fit.  For example, if you are trying to fit the parameters A and B in the function
		f[x] = A(x/B)
then the fitted parameters for A and B will be highly correlated.  You could double A and B, and still get the same answer.

fit["CorrelationMatrix"] // MatrixForm

For example, the 3rd element on the first row (-0.000416606) means that changes in the first parameter (A) are quite weakly correlated with changes in the third parameter (ϕ).  The 1's along the diagonals mean each parameter is perfectly correlated with itself, as expected.   Checking the correlation matrix is a good way to check whether your proposed fitting function really has the independent parameters you thought it did.

The resulting fit object has many other “Properties” you can query.  Here is a list of all of them:

fit["Properties"]