Least Squares Fitting:
Physics 338

In Chapter 8 of An Introduction to Error Analysis, by John R. Taylor, the author discusses the general theme of least-squares fitting.  This is based on the normal distribution discussion in Chapter 5.  Here, we will take the theoretical background as given, and show how the minimization of the least square difference leads to the standard results for some simple cases.

Computing the Mean

As a simple application, we will use the principle of least squares to calculate the mean of a set of numbers.

First, we will make up a data set of random numbers.  The following command creates a list of 5000 numbers selected from a normal distribution with a mean of 2.0 and a standard deviation of 0.1

rdata = RandomReal[NormalDistribution[2, 0.1], 5000];
npts = Length[rdata];
Histogram[rdata, ImageSize->Large]

Computing (χ^(2))

Next, we calculate (χ^2).   For N data points with a mean value of X, (and assuming all have the same uncertainty) the formula is

chisq1[data_, X_] := (1/(npts - 1))(∑_(i=1)^npts)(((data[[i]] - X))^2)

Of course, we don’t know X yet.  The purpose of this exercise is to find the value of X which minimizes (χ^2).  We claim that value of X is our best estimate of the mean.  Here is a plot showing (χ^2) for various guesses of X.

Plot[chisq1[rdata, X], {X, 1, 3}, AxesLabel->{"mean", "(χ^2)"}, LabelStyle -> Larger, ImageSize->Large]

FindMinimum[chisq1[rdata, X], X]

This tells us the value of X  that minimizes (χ^2), and also what that minimum value is.

Sqrt[%[[1]]] (* This will be the standard deviation *)

Mean and Standard Deviation

The usual definition for the mean results from minimizing (χ^2), and that minimum value is the same as the usual definition for the standard deviation.

Mean[rdata]

StandardDeviation[rdata]

Fitting a Straight Line

First, here is some sample data of {current, mass} pairs from the magnetic susceptibility experiment.

Data

Typing the data directly:

The data is a list of data points.  All lists are made of elements in curly braces { } separated by commas.  Each data point is itself a list, consisting of two numbers (in curly braces) separated by a comma.  To see what it looks like, display it with TableForm (or MatrixForm).

rawdata = {
{0.000,181.184},
{0.497,180.964},
{1.005,180.736},
{1.507,180.514},
{1.930,180.326},
{0.000,181.185},
{-0.504,181.409},
{-1.008,181.633},
{-1.502,181.853},
{-1.930,182.045},
{0.000,181.184}
};

npts = Length[rawdata]

Typing the data using the Classroom Assistant:

From the Palettes -> Classroom Assistant palette, look for the matrix item: (({{□, □}, {□, □}})).  (Depending on the version of Mathematica, it might be in the “Typesetting” menu, or in the “Advanced” menu.  (Use Ctrl-Enter to add a new row.)

rawdata = (({{0, 181.184}, {0.497, 180.964}, {1.005, 180.736}, {1.507, 180.514}, {1.930, 180.326}, {0, 181.185}, {-0.504, 181.409}, {-1.008, 181.633}, {-1.502, 181.853}, {-1.930, 182.045}, {0, 181.184}}));

Reading from a file:

To read the data from a text file in the same directory as your notebook, first tell Mathematica what directory to use, and then tell it   This particular file consists of two columns (current, period).  Read it in as a “Table”.  If you store your data in the same directory as the Mathematica notebook, then it is probably simplest to just set the notebook directory as the default.

SetDirectory[NotebookDirectory[]]

Typically, files might have header information (names of columns, dates, etc.)  Often they are delimited by some sort of comment symbol.  I often use ‘#’ for that purpose when I type the files in by hand.  Here is what my raw data file looks like:

FilePrint["mass_vs_current.txt"]

rawdata = Import["mass_vs_current.txt", "Table"]

The default import includes all those comment lines.  We want to select only those lines that contain the actual data.  The Selectp[ comnmand is handy for this.  In this specific case, we want only lines where there are exactly two elements, and both of them are numeric.  The following does that.  The first test, Length[#] == 2 checks that there are 2 elements.  The second test checks that the element is a vector (i.e. a list)  that consists of only numeric items.  You use ‘&&’ to mean logical “AND”.  You would use ‘||’ to mean logical OR.   The ‘&’ at the end is needed to make Mathematica set the ‘#’ to each individual item in the list.

data = Select[rawdata, Length[#] == 2 && VectorQ[#, NumericQ] &]

Mathematica also can import CSV (comma-separated-values) files.  More information is in the Import[ ] function.

Plot the data.

Be sure to give your axes meaningful labels and make them large enough to be read.  I have assigned it to a variable for later reuse.

dataplot  = ListPlot[data, AxesLabel -> {"Current (A)", "Mass (g)"}, LabelStyle->Larger, ImageSize->Large]

Theoretical Curve:  A straight-line fit.

Here we explore how to find the best-fit straight line by the method of least squares.  First, define the target function.  Use the unknown items (intercept, a0, and slope, a1) as parameters to the function.

yfit[x_, a0_, a1_] := a0 + a1 x

Calculating (χ^2)

Define a (χ^2)  function for the linear fit.  For given a0 and a1 values, compute the average of the differences squared between the data and the fit value.  Within the function, the data is in {i, m} pairs, so we pull out the current of the (i^th) value with data[[i, 1]], and the mass of the (i^th) value with data[[i, 2]].  You can use the Classroom Assistant Palette to format the sum, or you can use Mathematica’s Sum[ ] function directly.  Use whichever is easier for you to read.  The denominator has the ‘-2’ because with N data points and 2 free parameters, there are only N-2 degrees of freedom.

calculateChisq[data_, a0_, a1_] := (1/(Length[data] - 2))(∑_(i=1)^(Length[data]))(((yfit[data[[i, 1]], a0, a1]  - data[[i, 2]]))^2)

calculateChisq[data_, a0_, a1_] := Sum[(((yfit[data[[i, 1]], a0, a1]  - data[[i, 2]]))^2), {i, 1, Length[data]}] / (Length[data] - 2)

Initial Explorations

This command builds an interactive window showing the data, the current fit, and the chi squared value (as the plot title).  It draws vertical lines from each data point to the fit line.  Move the a0 and a1 sliders to minimize (χ^2).  Mathematica’s “Filling” option gets confused if the data isn’t sorted, so let’s go ahead and sort the data by the first entry in each line.

data = SortBy[data, First]

Manipulate[
fitted = Table[{data[[i, 1]], yfit[data[[i, 1]], a0, a1]}, {i, 1, Length[data]}];
Show[{ListPlot[{data,fitted },  Filling -> {1 ->  {{2}, {Red, Black}}}, PlotLegends -> {"data", "fitted"}],
Plot[yfit[x, a0, a1], {x, -2, 2}]}, PlotLabel->calculateChisq[data, a0, a1], ImageSize->Large] ,
{{a0, 181.2}, 180.7, 181.7, 0.001,  Appearance->"Open"},
{{a1, 0}, -1.00, 1.00, 0.02, Appearance->"Open"}
]

Minimizing (χ^2)

Picking the intercept a0 = 181.2 (close to the expected value), look at how (χ^2) varies with a1.

Plot[calculateChisq[data, 181.2, a1 ], {a1, -1, 1 }, AxesLabel -> {"a1", "(χ^2)"}, LabelStyle ->Larger]

Picking a1 = -0.4 (again, close to the expected value), look at how (χ^2) varies with a0.

Plot[calculateChisq[data, a0, -0.4 ], {a0, 180, 182 }, AxesLabel -> {"a0", "(χ^2)"}, LabelStyle ->Larger]

Ultimately, the “best” fit involves a two-dimensional minimization of (χ^2).  Sometimes it helps to visualize this sort of thing as a 3D plot or a density plot.

Plot3D[calculateChisq[data, a0, a1], {a0, 180.7, 181.7},{a1, -1, 0}, AxesLabel -> {"a0", "a1", "(χ^2)"}, LabelStyle->Larger, ImageSize->Large]

DensityPlot[calculateChisq[data, a0, a1], {a0, 180.7, 181.7},{a1, -1, 0}, LabelStyle->Larger, ImageSize->Large]

Mathematica can find the minimum here as well.

result = FindMinimum[calculateChisq[data, a0, a1 ], {a0, a1}]

Finally, here is the best fit curve.

Show[{
dataplot,
Plot[yfit[x, a0, a1] /. result[[2]], {x, -2, 2}, PlotStyle -> Red]
}]

Linear Model Fit

Mathematica can do this minimization automatically.  The LinearModelFit[ ] function searches the "parameter space" for the minimum value of (χ^2).    It reports a "Confidence Interval" that reflects the curvature of (χ^2) -- how much you can vary either parameter without making (χ^2) too large.

One other interesting thing to note:  Note how our (χ^2) space has a "valley" (dark purple in the DensityPlot) where you can explore what happens if you change a0 but don’t change a1, and vice-versa.  This is reflected in the LinearModelFit[ ] report.  Look at the "Correlation Matrix."  It tells you, in essence, how well correlated each parameter is with the other.  The 1's along the diagonal mean a0 is perfectly correlated with a0 (obviously) and a1 is perfectly correlated with a1.  The off-diagonal elements tell you that a0 and a1 are mostly uncorrelated -- increasing one doesn’t really affect the other.  This isn’t always true for fits.  Sometimes, a change in one parameter can be partially compensated for by a change in another.  In those cases, they are not completely independent parameters, and the Correlation Matrix elements are further from zero.

fit =LinearModelFit[data, x, x]

See the on-line help for more information on dealing with the results from LinearModelFit.  Here are some examples of things you can do with it.

fit["BestFit"]

fit["BestFitParameters"]

fit["ParameterConfidenceIntervalTable"]

fit["CorrelationMatrix"] // MatrixForm

Show[{ListPlot[data], Plot[fit[x] , {x, -2, 2}, PlotStyle -> Red, LabelStyle->Larger]}, ImageSize->Large]

fit["EstimatedVariance"]

This is the same as our (χ^(2 )).

{a0, a1} = fit["BestFitParameters"]

calculateChisq[data, a0, a1]

Interpreting the Uncertainties.

Mathematica gives you the uncertainties in the parameters, but you should also check whether the fit is reasonable.  In particular, do the differences between the data points and the fit make sense?  Is the size believable?  Are there systematic trends?

Mathematica will report the “FitResiduals”, which are the  differences between each data point and the fitted curve.  The average square of the residuals is called the “EstimatedVariance”.  (Actually you divide by N-2, since there are 2 degrees of freedom used up by the two fit parameters, a0 and a1.

fit["FitResiduals"]

Total[fit["FitResiduals"]^2]/(Length[fit["FitResiduals"]] - 2)

fit["EstimatedVariance"]

The square root of the estimated variance is thus the RMS (root mean square) error -- the “typical” amount by which the fitted line misses the data.  It has the same units as the original y values, and should be compared to the y uncertainties.

Sqrt[fit["EstimatedVariance"]]

Are differences of 0.0007 grams meaningful?   For this experiment, you can read the scale to 0.001 g, and the zero current reading in this case varied between 181.183 and 181.184.  Similarly, you probably saw fluctuations on the order of 0.001g routinely as you did the experiment.  Accordingly, a variance of 0.0007 grams seems quite plausible.

Other Fit Properties.

This is a complete list of all available properties of the fit.

fit["Properties"]

A note about the “Rsquared” value

fit["RSquared"]

This is close to 1.0.  Does that mean you have a good fit?  Not necessarily.  Consider the following parabolic data set, where each point has a “random noise” added to each point (the noise has a mean of 0 and a standard deviation of 1).

RandomReal[NormalDistribution[0, 1]]

parab = Table[{x, 1.0 + 0.5 x + 0.3 (x^2) + RandomReal[NormalDistribution[0, 1]]}, {x, 1, 20, 2}]

Let’s try fitting it with a straight line

pfit =LinearModelFit[parab, x, x]

pfit["ParameterConfidenceIntervalTable"]

Show[{ListPlot[parab], Plot[pfit[x] , {x, 0, 20}, PlotStyle -> Red, LabelStyle -> Larger]}]

pfit["RSquared"]

This is still close to 1.0, even though we’re missing a systematic trend in the data.  The (R^2) value does not distinguish between data points that scatter about a line and those that deviate systematically from the line.  You need to look at the actual graph.  Obviously, a parabola is the right fit here:

Sqrt[pfit["EstimatedVariance"]]

A typical variance of 9.4 seems unreasonable given that our data should have noise on the order of 1.

pfit2 =LinearModelFit[parab, {x, (x^2)}, x]

Show[{ListPlot[parab], Plot[{pfit[x], pfit2[x]} , {x, 0, 20}, PlotLegends -> {"Linear", "Parabolic"}]}]

Sqrt[pfit2["EstimatedVariance"]]

A typical mismatch of roughly 1 between the data and the fit makes sense, given that we added noise with a typical amplitude of about 1.