Least Squares Fitting:
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 some simple cases.

Clear["Global`*"]

Fitting a Straight Line

First, here is some sample data of {mass, angle} pairs from the torsional oscillator experiment.  (See the separate importing notebook for more about entering and importing data.)

Data

SetDirectory[NotebookDirectory[]]

rawdata = Import["TO-static-20250212.csv", "CSV"]

rawdata[[Range[3]]] // TableForm

data3 = Select[rawdata, VectorQ[#, NumberQ] &] (* Select only those lines for which all elements are numbers. *)

For this data, we want columns 1 and 2.

data = data3[[All, {1, 2}]]

Plot the data.

Be sure to give your axes meaningful labels

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

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 {m, θ} pairs, so we pull out the mass of the (i^th) value with data[[i, 1]], and the angle for 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) (* Fancy typeset presentation *)

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

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]

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

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}}}, PlotRange ->{{-425, 425}, {0, 5.5}}, PlotLegends -> {"data", "fitted"}],
Plot[yfit[x, a0, a1], {x, -400, 400}]}, PlotLabel-> StringForm["(χ^2) = ``", calculateChisq[data, a0, a1]], ImageSize->Large] ,
{{a0,2.5},1, 4, 0.1,  Appearance->"Open"},
{{a1, 0},-0.01, 0.01, 0.001, Appearance->"Open"}
]

Minimizing (χ^2)

Picking the intercept a0 = 3  (close to the expected value), look at how (χ^2) varies with a1.  The best value for a1 seems to be a little less that -0.004.

Plot[calculateChisq[data, 3, a1 ], {a1, -0.01,0}, AxesLabel -> {"a1", "(χ^2)"}]

Picking a1 = -0.004 (again, close to the expected value), look at how (χ^2) varies with a0.  The best value for the intercept is near 3.0 radians.

Plot[calculateChisq[data, a0, -0.004], {a0, 2, 4}, 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, 4},{a1, -0.01, 0}, AxesLabel -> {"a0", "a1", "(χ^2)"}]

DensityPlot[calculateChisq[data, a0, a1], {a0, 2, 4},{a1, -0.01, 0}]

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, -400, 400}, 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, m}, m]

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"]

You can extract the values from that table, though it’s clunky.

fit["ParameterConfidenceIntervalTableEntries"] // TableForm

For example, here is the slope and its uncertainty (the slope is the second entry in the table):

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

fit["CorrelationMatrix"] // MatrixForm

Show[{ListPlot[data], Plot[fit[m] , {m, -400, 400}, PlotStyle -> Red]}, LabelStyle -> Larger, AxesLabel -> {"Mass (g)", "Angle (radians)"}, ImageSize->Scaled[0.7]]

fit["EstimatedVariance"]

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

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

calculateChisq[data, a0, a1]

Quality of Fit

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?  Mathematica will give you the “Estimated Variance”, which is the same as our manually-calculated (χ^2).  That is, it is the average of the square difference between each point and its corresponding fit.  (Actually it divides by N-2, since there are 2 degrees of freedom used up by the two fit parameters, a0 and a1.

fit["EstimatedVariance"] (* Again, note this is the same as our χ-squared. *)

RMSE = Sqrt[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.  Are differences of 0.04 radians significant?   For this experiment, you can easily read the scale to 0.1 radian, so differences of 0.04 radians seem plausible.

Residuals

Mathematica also lists the individual deviations, known as Residuals.  These can be useful for exploring systematic trends.  You will explore this idea in a separate homework assignment.

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 (with a small amount of random noise added)

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

Let’s try fitting it with a straight line

pfit =LinearModelFit[parab, {1, x}, x]

pfit["ParameterConfidenceIntervalTable"]

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

pfit["RSquared"]

Sqrt[pfit["EstimatedVariance"]]

The (R^2) value is is 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, {1, x, (x^2)}, x]

pfit2["ParameterConfidenceIntervalTable"]

Sqrt[pfit2["EstimatedVariance"]]

Show[{ListPlot[parab], Plot[{pfit[x], pfit2[x]} , {x, 0, 20}, PlotLegends -> Automatic]}, LabelStyle->Larger, ImageSize -> Scaled[0.6]]

Summary of Linear Fit

Read or input the data and transform into appropriate form.

rawdata = Import["TO-static-20250212.csv", "CSV"];

data3 = Select[rawdata, VectorQ[#, NumberQ] &]; (* 3 columns of data *)

data = data3[[All, {1, 2}]]; (* We only want 2 columns:  mass and angle. *)

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

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

fit["BestFitParameters"]

fit["ParameterConfidenceIntervalTable"]

RMSE = Sqrt[fit["EstimatedVariance"]] (* Root Mean Square error *)

Overlay plots

finalplot = Show[{
dataPlot,
Plot[fit[m], {m, -400, 400}, PlotStyle->Red]
}, ImageSize -> Scaled[0.5]
]

Save your work

Export["TO-static-plot.pdf", finalplot]