Fitting Data with Mathematica for a Nonlinear Model

In a number of experiments this semester, you will be asked to fit experimental data to a theoretical curve where a nonlinear fit is needed.  This notebook walks through an example of calculating chi-squared in such a case where LinearModelFit won’t work.

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

Clear["Global`*"]; (* Clear out any stored definitions. *)
DateString[]

Entering Data

Importing from a file

SetDirectory[NotebookDirectory[]]

FilePrint["nonlinear-modelfit-1-data.txt", 5] (* Print out the first 5 lines *)

This data file is just two columns of numbers separated by spaces.  Mathematica can import this as a simple table.

data = Import["nonlinear-modelfit-1-data.txt", "Table"];

dataplot = ListPlot[data, PlotRange->All, PlotStyle -> Red, Frame -> True, LabelStyle->Larger, FrameLabel -> {"Time (s)", "Amplitude (V)"}, ImageSize -> Scaled[0.5]]

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.  This function has 4 free parameters:  A, ω, ϕ, and xoff.

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

Clear[chisq]
chisq[data_, A_, ω_, ϕ_, xoff_] := Sum[(((data[[i, 2]] - x[A, ω, ϕ, data[[i, 1]], xoff] ))^2), {i, 1, Length[data]}] / (Length[data] - 4)

For convenience of our chisq visualization, pick Amplitude and offset values that seem plausible.  Look at varying ω and ϕ.

myxoff = 0.0044; myA= 1.5;

Manipulate[
Show[{
dataplot, 
Plot[x[myA, ω, ϕ, t, myxoff], {t, 0, 3}]
}, ImageSize -> Large, PlotLabel->chisq[data, myA, ω, ϕ, myxoff]] , 
{{ω,0.65}, 0, 6, 0.05, Appearance->"Open"},
{{ϕ, π}, 0, 2π, 0.1, Appearance->"Open"}
]

Minimizing (χ^2)

One-dimensional minimization

Picking A = 1.5 and a phase of π, look at how (χ^2) varies with ω.

myA = 1.5; myϕ = π;

Plot[chisq[data, myA, ω, myϕ, myxoff], {ω, 0,6 }, AxesLabel -> {"ω", "(χ^2)"}]

There clearly is a global minimum near 4.8, but there are also local minima near 0.7 and 2.0.  Without help, Mathematica might get trapped in a local minimum.

FindMinimum[chisq[data, myA, ω, myϕ, myxoff],ω] (* Mathematica tries starting at ω = 1 *)

FindMinimum[chisq[data, myA, ω, myϕ, myxoff],{ω, 3}] (* Tell it to start at 3 *)

This finds the second minimum.  An even better initial guess will find the global minimum.

FindMinimum[chisq[data, myA, ω, myϕ, myxoff],{ω, 4}] (* Tell it to start near 4 *)

Two-dimensional minimization

Plot3D[chisq[data, myA, ω, ϕ, myxoff], {ω, 0, 6}, {ϕ, 0, 2π}, PlotRange->All, AxesLabel->{"ω", "ϕ", "(χ^2)"},LabelStyle->Larger,  ImageSize -> Large]

Sometimes it helps to visualize this sort of thing as a 2d density plot.  The minimum near ω = 4.75 and ϕ=π is evident.

DensityPlot[chisq[data, myA, ω, ϕ, myxoff],  {ω, 0, 6}, {ϕ, 0, 2π}, PlotRange->All, FrameLabel->{"ω", "ϕ"}, LabelStyle->Larger, Mesh -> True, ColorFunction -> "Rainbow", ImageSize -> Scaled[0.6]]

Here again is a failed search.  Given our guess for amplitude and offset, find the best values for ω and ϕ:

result = FindMinimum[chisq[data, myA, ω,ϕ, myxoff], {ω, ϕ}]

Mathematica chooses initial starting values of ‘1’.  This gets caught in the local minimum on the left hand side.  However, if you give it a good starting point for at least one of the parameters,   Mathematica finds the minimum correctly.  Here we will let it find all of the parameters, but give a good starting hint.  (This gives a negative amplitude.  More on that below.)

result = FindMinimum[chisq[data, A, ω, ϕ, off], {A, {ω, 4}, ϕ,off}]

NonlinearModelFit

NonlinearModelFit[ ] performs this minimization automatically.  Though it has certain algorithms to attempt to avoid local minima, you still sometimes need to give it hints.

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

fit["BestFitParameters"]

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

Sqrt[fit["EstimatedVariance"]] (* The typical error is about 1 volt. *)

This fit is quite poor.   Give it a hint for ω, and it does much better.

fit = NonlinearModelFit[data, x[A, ω, ϕ, t, xoff], {A, { ω, 4}, ϕ,  xoff}, t];
fit["BestFitParameters"]

Sqrt[fit["EstimatedVariance"]] (* The typical error is now about 8 mV. *)

However, note the negative amplitude and near zero phase.  If all you want is the frequency, then that is fine.  However, if you want a positive amplitude, then you also have to give a hint for the phase:

fit = NonlinearModelFit[data, x[A, ω, ϕ, t, xoff], {A, { ω, 4}, {ϕ,π},  xoff}, t];
fit["BestFitParameters"]

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

Fit Results and Diagnostics

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

Best Fit Parameters

fit["BestFitParameters"]

The best fit parameters are given as a set of replacement rules.  You can assign a value to a variable with the ./ notation.  For example, to extract the phase you could do

phase = ϕ /. fit["BestFitParameters"]

Similarly, you can caluculate the period by using the replacement rule for ω:

T = (2π/ω) /. fit["BestFitParameters"]

Parameter Confidence Interval

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

fit["ParameterConfidenceIntervalTable"]

The table entries are available in the raw form, and you can extract the ones you want with the Part [[ ]] notation.  For example, the ω values are in the 2nd row, 1st and 2nd columns

fit["ParameterConfidenceIntervalTableEntries"]

{bestω, δω} = fit["ParameterConfidenceIntervalTableEntries"][[2, {1, 2}]]

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 also double 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.  

Here, the only tricky one is the -0.87 correlation between ω and the phase ϕ.  Because there are only two full oscillations in the data, you can (partially) compensate for increasing ϕ by decreasing ω.   This shows up clearly in the density plot of (χ^2), where I have fixed A and xoff to approximately the right values, but allowed ω and ϕ to vary.  The deep purple diagonal valley at the center shows how you can increase ω but decrease ϕ without changing (χ^2)much.  (You can also fiddle with the Manipulate box above to see this.)

RMS Error

The RMS error is the typical amount the fit misses a data point.  Here it is about 0.0077, which is about 0.5% of the amplitude.  Depending on the specific sensor used, this is probably a plausible uncertainty.

RMSE = Sqrt[fit["EstimatedVariance"]]

If you want to compare that  to the amplitude, and express it as a percentage, you can write:

(100 * Sqrt[fit["EstimatedVariance"]]) / (A/.  fit["BestFitParameters"])

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

fit["Properties"]

Summary of Usage

Import Data

SetDirectory[NotebookDirectory[]];

data = Import["nonlinear-modelfit-1-data.txt", "Table"];
dataplot = ListPlot[data, PlotRange->All, PlotStyle -> Red, Frame -> True, LabelStyle->Larger, FrameLabel -> {"Time (s)", "Amplitude (V)"}, ImageSize -> Scaled[0.5]];

Construct Model

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

Run NonlinearModel Fit, giving hints for parameters as needed

fit = NonlinearModelFit[data, x[A, ω, ϕ, t, xoff], 
{A, { ω, 4}, {ϕ, π}, xoff}, t]; (* Here we give initial guesses for ω and ϕ *)
fit["ParameterConfidenceIntervalTable"]

RMSE = Sqrt[fit["EstimatedVariance"]] (* Find typical error. *)

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

Extracting values (with uncertainty) from the fit

phase = ϕ /. fit["BestFitParameters"]

{bestω, δω} = fit["ParameterConfidenceIntervalTableEntries"][[2, {1, 2}]]