Least Squares Fitting to a Straight Line
Physics 238

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 fitting a straight line.

Fitting a Straight Line

Data

First, here is some sample data of {mass, angle} pairs from the torsional oscillator experiment.  Masses are listed as negative when they apply a negative torque.

Typing the data using the Classroom Assistant:

From the Palettes -> Classroom Assistant pallete, 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, and Ctrl-, to add a new column.)

data =(({{-400, 1.505}, {-350, 1.645}, {-300, 1.800}, {-250, 1.980}, {-200, 2.180}, {-150, 2.380}, {-100, 2.600}, {-50, 2.800}, {0, 3.005}, {50, 3.205}, {100, 3.405}, {150, 3.620}, {200, 3.830}, {250, 4.025}, {300, 4.235}, {350, 4.460}, {400, 4.645}, {0, 3.005}}));

Typing the data directly:

You can also type in 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).

data = {
{-400,1.505},{-350,1.645},{-300,1.800},{-250,1.980},{-200,2.180},{-150,2.380},{-100,2.600},{-50,2.800},
{0,3.005},
{50,3.205},{100,3.405},{150,3.620},{200,3.830},{250,4.025},{300,4.235},{350,4.460},{400,4.645},
{0, 3.005}
};

npts = Length[data]

Plot the data.

Be sure to give your axes meaningful labels.

ListPlot[data, AxesLabel -> {"Mass (g)", "Angle (radians)"}]

The default plot isn’t particularly nice, especially if you want to display it or include it in another document.  There are lots of options to tweak; the following generally work fairly well in practice.  It also shows how to assign a name to that plot, to make it easy to re-use it and layer the fit on top.

dataPlot = ListPlot[data, ImageSize -> Scaled[0.8], Frame -> True, FrameLabel -> {"Mass (g))", "Angle (radians)"}, LabelStyle -> Larger, PlotStyle -> Thick]

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 fuction, 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 Pallete 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_] := (1/(Length[data] - 2))Sum[(((yfit[data[[i, 1]], a0, a1]  - data[[i, 2]]))^2), {i, 1, Length[data]}]

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, -450, 450}]}, PlotLabel->calculateChisq[data, a0, a1], LabelStyle -> Larger, ImageSize->Scaled[0.8]] ,
{{a0, 3}, 1.5, 4.5, 0.005,  Appearance->"Open"},
{{a1, 0}, -0.01, 0.01, 0.0002, Appearance->"Open"}
]

Minimizing (χ^2)

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

Plot[calculateChisq[data, 3.0, a1 ], {a1, 0.003, 0.005}, AxesLabel -> {"a1", "(χ^2)"}]

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

Plot[calculateChisq[data, a0, 0.004 ], {a0, 2.8, 3.2 }, AxesLabel -> {"a0", "(χ^2)"}]

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, 2.8, 3.2},{a1, 0.003, 0.005}, AxesLabel -> {"a0", "a1", "(χ^2)"}]

DensityPlot[calculateChisq[data, a0, a1], {a0, 2.8, 3.2},{a1, 0.003, 0.005}]

Mathematica can find the minimum here as well.

Clear[a0, a1]
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, -450, 450}, 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, {1, 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[{dataPlot, Plot[fit[x] , {x, -450, 450}, PlotStyle -> Red]}]

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.03 radians meaningful?   For this experiment, you can read the scale to at lead ± 0.02 radians, and you can likely interpolate to 0.01 or even 0.005 radians, so variances of 0.03 radians seem larger than can be explained by simple random measurement error.    Looking carefully at the graph, it seems likely that we twisted the wire a little beyond its linear range, and are starting to see a slightly nonlinear response.

Residual Plot

Sometimes, it is helpful to look at the residuals in more detail to try to understand why the fit and data might not agree.  A residual plot can sometimes highlight trends within the data that might not be evident in the initial fit plot and that might warrant additional investigation.

residuals = Table[{data[[i, 1]], (data[[i, 2]] - fit[data[[i, 1]]])}, {i, 1, Length[data]}]

ListPlot[residuals, AxesLabel -> {"mass (g)", "Residual Angle (rad)"}, LabelStyle -> Larger, ImageSize->Large]

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:

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

Let’s try fitting it with a straight line

pfit =LinearModelFit[parab, x, x]

pfit["ParameterConfidenceIntervalTable"]

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

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:

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

Show[{ListPlot[parab], Plot[{pfit[x], pfit2[x]} , {x, 0, 20}, PlotLegends -> Automatic]}, LabelStyle->Larger, AxesLabel->{"x", "p"}]

Summary Presentation

Given the data set:

data

Fit to a straight line

fit =LinearModelFit[data, {1, x}, x]

Find the values for the slope and intercept, along with their uncertainties.

fit["BestFitParameters"]
fit["ParameterConfidenceIntervalTable"]

Look at how far off (in radians) a typical data point is from the best-fit line.  Think about this value.  Is it reasonable?

Sqrt[fit["EstimatedVariance"]]

Show the data and the fit on the same graph.

Show[{ListPlot[data], Plot[fit[x] , {x, -450, 450}, PlotStyle -> Red]}, Frame -> True, LabelStyle -> Larger, FrameLabel -> {"Mass (g)", "Angle (radians)"}]