# Linear Least Squares Curve Fitting of Straight Line Models

## Summary

Curve fitting is a useful procedure for testing theoretical models and estimating model parameters and their uncertainties. Curve fitting is carried out by choosing the set of estimated parameters that minimize an error function. The conventional error function for analytic results with calculus is the sum of the squares of residuals. The residuals are the difference between the model prediction and the experimental observation. Weighted fitting takes into account the size of the error bars where measurements with smaller error bars count more heavily towards the minimization. Weighted fits are more robust to the effects of outlying points which generally have larger error bars. The weighted mean is a least squares procedure, fit to , which is used to combine different experimental results of the same quantity. This is how CODATA and NIST calculate the accepted fundamental constants of nature from various experiments. Curve fitting is a tedious calculation and is best done with a software package like MATLAB. We illustrate how to fit the various linear models that are discussed with simulated data in MATLAB. One can either put the total data set into the curve fit or average repeated measurements first and fit those.

## The least squares error formula for curve fitting

Suppose we have a set of measurements and we assume the span some range on the x axis, but there may be repetitions. This is in fact recommended. Always repeat measurements so you can calculate their uncertainties directly from the data. We will not factor in how to use the uncertainties the best way in curve fitting yet, but we will consider it later. If we have a model we define a residual as the difference between the model and the measurement

or in general

Now we can start defining some trial error functions. We could choose to minimize the residuals directly.

The problem with this error function is that some of the errors will be positive and some will be negative. They basically cancel each other out and you can’t find the best and . If we wanted to make the signs the same then we could try absolute values.

The problem with this error function is that taking derivatives of absolute value functions is nasty. You can’t get many nice results with calculus using this error function. The next simplest modification of the error function is to take

Then we look for the least squares of the residuals. This is the convention that most people use, but it is not totally arbitrary. We will discuss in a more advanced lesson some justification for the least squares method.

## Fitting y = A + Bx

Expanding the 9 terms in for gives

Using the “S” notation we can write as

where and etc.

The way we find the minimum of is to take the partial derivatives with respect to and and set them equal to zero.

This gives two equation which can be solved for and .

This equation can be written as a matrix equation.

As long as the determinant of the first matrix on the left is not zero we can solve these linear equations with Cramer’s rule. The condition is basically that two of the be different. You can’t make any sense of the fit if all the measurements are just at one point. You need two different to define a linear function.

A common expression is .

Now we find the uncertainties in and by using the propagation of error formula.

If we have not repeated the measurements at each several times then we don’t know what the are. This is a good justification for repeating the measurements. If we used all different then we can just assume that is the same for all the points. If we use some arguments from statistics we can say

There are N points, but two degrees of freedom less because we have two fit parameters. The expression in each parenthesis is roughly equal to . So we can write

Now we can do the calculation of the uncertainty. We have

I define . It has the properties that for each . , ,

A similar calculation for gives

If we know each then it makes sense to just do the weighted calculation which we will describe later in this lesson.

## Fitting y = Bx

For a model like Ohm’s law, , we want to fit a straight line forced to pass through the origin, .

The calculus is similar.

Minimizing for using we find

The propagation of error formula gives the uncertainty of .

Again if we don’t know the uncertainties we can use

## Fitting y = A

The last model we will consider is a horizontal line.

.

So we take the first derivative and find

or

That is just the arithmetic mean for . We can find the standard deviation in the mean by the same procedure we used before.

We use for where .

So the standard deviation of the mean is given by

.

## Weighted fits and using chi squared

It is more useful to do curve fitting when you have repeated the measurements. Then you can know what the different standard deviations are. If the standard deviation is large then you don’t want to count those points as much as the ones for which the standard deviation is small. Another method of doing least squares curve fitting is to weight the residuals by . This sum is called chi squared, .

This does not change the math so much that we cannot solve for and . You just have to replace each formula with . For example, , , and so on.

There are slight differences in how you interpret chi squared as a statistic if you use each pair of data points in the equation or first average the points and calculate their standard deviation of the mean. This boils down to using or . Often you don’t have access to the individual data points if you look at someone else’s data, so you have to use the average values and the standard deviations of the mean instead of . We can then summarize the results from before with the new weighting scheme.

For .

For . We have . So we can calculate the uncertainty.

Probably the most useful weighted formula is for the weighted mean when . This allows you to average together measurements.

Taking we find

This is referred to as the weighted mean. We can find the uncertainty in the weighted mean

In the limit, that all are the same we recover the arithmetic mean and the standard deviation of the mean for the uncertainty.

This formula is quite useful. If you have a series of accurate measurements some of which are precise and others not so much you can combine them to get a more accurate and precise average based on the weighting. Most people are probably already familiar with this in terms of calculating their grades. Say 25 percent for homework, 25 percent for the midterm and 50 percent for the final. Then an average score might look like

Pretty close, let’s give them an A minus anyways.

## MATLAB functions for straight line curve fitting

### y = A + Bx

Example. Write a Matlab function that fits a straight line through a dataset.

We can easily compute all the sums in Matlab. Computing the sums by hand would be quite tedious.

function [ A sigA B sigB chi2] = lineABx( xi, yi, sigi )

%[A sigA B sigB chi2] = lineBx(xi, yi, sigi)

% Least sqaures fit to a line y = A + Bx

wi = sigi.^(-2);

Sw = sum(wi);

Swx = sum(wi.*xi);

Swxx = sum(wi.*xi.*xi);

Swy = sum(wi.*yi);

Swxy = sum(wi.*xi.*yi);

Delta = Sw*Swxx – Swx^2;

A = (Swy*Swxx – Swx*Swxy)/Delta;

B = (Sw*Swxy – Swx*Swy)/Delta;

sigA = sqrt(Swxx/Delta);

sigB = sqrt(Sw/Delta);

chi2 = sum( wi.*(yi – A- B.*xi).^2)

end

Example. Write a Matlab function that fits an unweighted straight line through a dataset.

We have

function [ A sigA B sigB] = lineABxU( xi, yi )

%[A sigA B sigB] = lineBx(xi, yi)

%Unweighted least sqaures fit to a line y = A + Bx

S = length(xi);

Sx = sum(xi);

Sxx = sum(xi.*xi);

Sy = sum(yi);

Sxy = sum(xi.*yi);

Delta = S*Sxx – Sx^2;

A = (Sy*Sxx – Sx*Sxy)/Delta;

B = (S*Sxy – Sx*Sy)/Delta;

sigma = sqrt(1/(length(xi)-2)*sum((yi – A – B.*xi).^2));

sigA = sigma*sqrt(Sxx/Delta);

sigB = sigma*sqrt(S/Delta);

end

### y = Bx

Example. Write a Matlab function to fit a straight line through the origin

We can evaluate the sums easily in Matlab.

function [ B sigB ] = lineBx( xi, yi, sigi )

%[B sigB chi2] = lineBx(xi, yi, sigi)

%Weightef least sqaured fit to a line y = Bx

wi = sigi.^(-2);

Swxy = sum(wi.*xi.*yi);

Swxx = sum(wi.*xi.*xi);

B = Swxy/Swxx;

sigB = 1/sqrt(Swxx);

chi2 = sum( wi.*(yi – B.*xi).^2)

end

Example. Write a Matlab function to fit an unweighted straight line through the origin

The code for the unweighted fit is as follows

function [ BU sigBU ] = lineBxU( xi, yi )

%[B sigB] = lineBx(xi, yi )

%Weightef least sqaured fit to a line y = Bx

Sxy = sum(xi.*yi);

Sxx = sum(xi.*xi);

BU = Sxy/Sxx;

variance = 1/length(yi)* sum( (yi – BU).^2 );

sigma = sqrt(variance);

sigBU = sigma/sqrt(Sxx);

end

### y = A

Example. Write a Matlab function to calculated the weighted mean and its uncertainty.

Implementing the theory we have

function [ A sigA chi2] = lineA( xi, yi, sigi )

%[A sigA chi2] = lineA(xi, yi, sigi)

%Least sqaures fit to a constant y = A

wi = sigi.^(-2);

Sw = sum(wi);

Swy = sum(wi.*yi);

A = Swy/Sw;

sigA = 1/sqrt(Sw);

chi2 = sum( wi.*(yi – A).^2 );

end

## MATLAB curve fitting with simulated data

Example. Illustrate a weighted linear fit to .

For the purpose of simulating, data we will assume a linear model plus random gaussian noise. The gaussian distribution noise can be generated with the command . We will

take measurements at each point for . With the model . The goal of this simulation is to compare fitting the measurements and their standard deviations with the means at each and their standard deviation of the mean. It will turn out that fitting either way gives identical estimates.

clear; figure(1); clf;

format long;

A = 10;

B = 3;

mu = 0;

sig = 3;

Ni = 5;

x = 1:1:10

allx = [];

ally = [];

allysig = [];

avgy = [];

for i = 1:length(x)

xtemp = [];

temp = [];

tempsig = [];

for j = 1:Ni

err = normrnd(mu,sig);

temp(j) = A + B*x(i) + err;

xtemp = [xtemp x(i)];

end

sigy1 = std(temp);

allx = [allx xtemp];

ally = [ally temp];

allysig = [allysig ones(1,Ni)*sigy1 ];

avgy(i) = [mean(temp)];

avgysig(i) = [sigy1/sqrt(Ni)];

end

[A1 sigA1 B1 sigB1 chi1] = lineABx(allx,ally,allysig);

[A2 sigA2 B2 sigB2 chi2] = lineABx(x, avgy, avgysig);

figure(1);

plot(allx, ally,’.k’,’markersize’,20);

hold on;

plot(x, avgy, ‘xr’,’markersize’,20);

axis( [ 0 11 0 1.1*max(ally)]);

xlabel(‘x’,’fontsize’,20);

ylabel(‘y’,’fontsize’,20);

x1 = [0 11]

y1 = [ A1 (A1 +B1*11)];

plot(x1,y1);

You can check the fit parameters A1, A2, B1, and B2 to see that they are numerically the same. Using the measurements with the standard deviations is equivalent to using the means with the standard deviations of the mean. Whichever you choose to display is up to you. In physics usually one shows the standard error of the mean because the samples are assumed to be identical and there are often a surplus of points. It might make more sense in biology to show the actual variation of the samples where they can be different and the number of measurements is less. You can modify the code by increasing to show that the parameter estimates become more accurate and precise as is increased.

Model y = 10 + 3x + normrnd(0,3), with five samples using standard deviation and measurements

Model y = 10 + 3x + normrnd(0,3), with five samples using means and standard deviation of the means. We obtain the same fit as before

Example. Show the effect of an outlying point when a linear fit is weighted or unweighted.

We will use simulated data like before but where the last point is an outlier. We will show the effect in an unweighted fit and a weighted fit which is more robust and fits the data better.

clear;

format long;

x = 1:10

y = [ 13.69 14.9 17.1 24.4 25.4 26.6 29.48 32.0 37.36 48.0 ]

sigy = [ 0.58 1.0 1.0 1.2 1.6 1.4 0.86 1.8 0.84 8.1 ]

figure(1);clf;

errorbar(x,y,sigy,’.k’,’markersize’,10);

hold on;

axis( [ 0 11 0 1.3*max(y)]);

xlabel(‘x’,’fontsize’,20);

ylabel(‘y’,’fontsize’,20);

title(‘Effect of an outlying point on unweighted fit’,’fontsize’,18);

[A1 sigA1 B1 sigB1 chi1] = lineABx(x,y, sigy);

[A2 sigA2 B2 sigB2 ] = lineABxU(x, y);

x1 = [0 11]

y1 = [ A1 (A1 +B1*11)];

plot(x1,y1,’k-‘);

x2 = [0 11]

y2 = [ A2 (A2 +B2*11)];

plot(x2,y2,’k-.’);

Model y = 10 + 3x + normrnd(0,3), with five samples, standard deviation of the mean used and an outlying point. The outlying point shifts the slope of the unweighted fit, but does not effect the weighted fit too much.

Comments are Disabled