LAB #5a: Least Squares Fitting


TABLE OF CONTENTS

Introduction

Exercise 1

Theory

Exercise 2

Example

Exercise 3



Introduction

This lab will only take a single session. It covers material from Section 2.12 of Atkinson.

This lab is concerned with approximating a set of points in the plane by the best possible straight line. This is an example of a simple unconstrained minimum problem that results in a simple 2x2 linear system of equations.

Theory

Suppose you are given a set of points in the plane, (xi,yi), for i=1,2,...,n and you want to find the best straight line, y=m*x+b, that approximates the data points. One way of determining the "best" line is to find the minimum of the usual Euclidean distance of the points from the line (measured vertically) over all possible choices of m and b. It turns out that this is not difficult, and many of you have seen it already in elementary calculus. The idea goes in the following manner.

Define a function L(m,b)=Sumi(yi-m*xi-b)2 where the function "Sumi " denotes the sum over i. The function L(m,b) is a non-negative, differentiable function, so it has a global minimum. The values of m and b at the minimum are the desired slope and intercept of the best line to fit the data. It is easy to find the global minimum of L: just take the partial derivative with respect to m and b and set them to zero. These partials are

          -2*Sumi(xi*(yi-m*xi-b)=0
          -2*Sumi(yi-m*xi-b)=0



And this can be rewritten as

           m*Sumi(xi* xi)  + b*Sumi(xi) = Sumi(xi* yi)
           m*Sumi(xi)      +  b*N       = Sumi(yi)


This, of course, is a pair of equations in the two unknowns, m and b.

Example

We will now generate a set of points and find the best straight line that fits those points.

Exercise 1: Generate a set of ten data points by setting x=1:10, and y=3*x-1. Of course, if we tried to fit these points with a straight line, the line would have slope 3 and intercept (-1), so we will jiggle the points a little by using the Matlab rand function,

y=y+10*rand(size(y))-5

The rand function returns random numbers between 0 and 1, so we multiply by 10 to get something large enough to see on a plot, and subtract 5 to get something with zero average. Plot the points.

Exercise 2: Use the Matlab sum function to compute the coefficients of the system of equations written above. The easiest way to do this is to suppose the coefficient matrix is called A, and set

A(1,2)=sum(x)

The other matrix elements A(1,1), A(2,1), and A(2,2) as well as the right hand side vector B(1) and B(2) can also be computed using simple vector operations. Solve the resulting system, yielding m and b. Are they anywhere near the values of 3 and (-1) that we used to generate the data?

Exercise 3: Plot the data points as circles (plot(x,y,'o')) and plot the line y=m*x+b, with the values for m and b that you just computed, on the same plot. Is the straight line a reasonable one, based on the data points?

Remark (not an exercise): It would be interesting to know how close you would have come to the original values of 3 and (-1) if you chose a smaller value to "jiggle" the data and how big can you make the "jiggle" before all relationship to the original values seems to be lost.

Back to the MATH2070 page.

Last revised on October 9, 2000.