Introduction

One way to think about the interpolating polynomial discussed before is that it was a function of the form

(Subscript[p, n-1])(x)=(Subscript[a, 0])+(Subscript[a, 1])x+(Subscript[a, 2])(x^2)+...+(Subscript[a, n-1])(x^(n-1))=(∑_(j=0)^(n-1))(Subscript[a, j])(x^j)

with the (Subscript[a, j]) chosen so that the graph of the polynomial is forced to pass through the n data points.  The order of the interpolating polynomial is tied to the number of data points, but as the order of the polynomial increases so do the number of spurious oscillations it may inflict on the approximating function g.  Is it possible to choose a lower-order polynomial that is, in some sense, a good representation of the data?  This is where linear regression, or minimizing least square error, comes in.

To start, let's look at a more general problem.  Suppose we wish to approximate a function f  by a function g of the form

g(x)=(∑_(k=1)^m)(Subscript[a, k])(Subscript[∅, k])(x)

where the functions (Subscript[∅, k])(x) are known and we must choose the coefficients (Subscript[a, k]) to satisfy some criteria.  In the case of interpolating polynomials, (Subscript[∅, k])(x)=(x^(k-1)) and the criteria is that g and f must agree at m data points.  Another example would be approximation of f  by trigonometric polynomials, i.e.

g(x)=(∑_(k=0)^n)((Subscript[a, k])sin(k x)+(Subscript[b, k])cos(k x))

where the criteria is to minimize the maximum error between f and g on the interval [0, 2π].  That is, minimize the quantity

E(a) = max( | f(x) - g(x) | )

where the maximun is taken over the the interval [0, 2π] and a is a vector with components of the unknown coefficients.

The "least squares" idea is to minimize the "distance" between the approximating function, g, and the data.  In this case the method of measuring distance, or error, is to use the function

E(a)=(∑_(l=1)^n)(((f((Subscript[x, l]))-g((Subscript[x, l]))))^2)

where again a is really a vector with components of the unkown coefficients of the approximating function g.  Our goal is to make this error function as small as possible through the choice of the coefficients.  To accomplish this minimization we use a result from multivariable calculus:  At a local extreme value of E the following condition must hold

(∂/(∂(Subscript[a, j])))E=0

for j = 1, 2, ..., m.  The solution of these equations provides the coefficients (Subscript[a, k]) that minimize the error function E and give us the approximate representation we want.  It can be shown that the critical point thus obtained does indeed occur at a minimum.

As an example of how this works, consider the problem of finding a line,

g(x) = a x + b

that best represents a set of data, ((Subscript[x, i]),(Subscript[y, i])) for i = 1, 2, ..., n.  In this case, the unknown coefficients are a and b, and the square error function is

E(a,b)=(∑_(i=1)^n)((((Subscript[y, i])-g((Subscript[x, i]))))^2)=(∑_(i=1)^n)((Subscript[y, i])-a (Subscript[x, i])-b)

The partial derivatives to be set to zero are

(∂/(∂a))E=-2((∑_(i=1)^n)(Subscript[x, i])((Subscript[y, i])-a (Subscript[x, i])-b))=2(a((∑_(i=1)^n)((Subscript[x, i])^2))+b((∑_(i=1)^n)(Subscript[x, i]))-(∑_(i=1)^n)(Subscript[x, i])(Subscript[y, i]))

and

(∂/(∂b))E=-2((∑_(i=1)^n)((Subscript[y, i])-a (Subscript[x, i])-b))=2(a((∑_(i=1)^n)(Subscript[x, i]))+n b-(∑_(i=1)^n)(Subscript[y, i])).

Setting these two partial derivatives to zero gives a system of linear equations in the unknown constats a and b.  Solving the system for a and b gives,

a=(1/D)(n(∑_(i=1)^n)(Subscript[x, i])(Subscript[y, i])-((∑_(i=1)^n)(Subscript[x, i]))((∑_(i=1)^n)(Subscript[y, i])))

and then

b=(1/n)(((∑_(i=1)^n)(Subscript[y, i]))-a((∑_(i=1)^n)(Subscript[x, i])))

where

D=n(∑_(i=1)^n)((Subscript[x, i])^2)-((∑_(i=1)^n)(Subscript[x, i]))((∑_(i=1)^n)(Subscript[y, i])).

In general, fitting data with m functions generates a linear system of m equations for the coefficients in the definition of g(x) above.

How Mathematica Does it

First we need to discuss some commands that will allow us to bring data into Mathematica from external data files.  Many times data collection involves electronic recorders or manual data entry into files.  These files are stored on a computer or network usually using a directory structure.  So the first opperation needed is to tell Mathematica where to look for the file.  The Directory command gives the current location that Mathematica will search for a file.  The SetDirectory command will change that location.  In the following commands we show what the default "current directory" is when Mathematica starts, then sets the "current directory" to one where a data file resides, and finally we use the Import command to enter a data file into Mathematica.

olddir=Directory[]
newdir=SetDirectory["T:\\white\\classes\\Computer Aided Math\\ClassNotes\\Mathematica 7"]
data=Import["regdata","Table"]

This data looks like a list of ordered pairs.  To see how many data points there are we use the Dimensions command:

dim=Dimensions[data]

So we have 50 data points.  Let's graph this data to see what we have.

dataplot=ListPlot[data]

Just as a fun exercise, let's take a look at the interpolating polynomial for this data.

interpy=InterpolatingPolynomial[data,x];

We suppressed the output because it is a polynomial of degree 49.  Graphing this polynomial over the same range of values gives,

interpplot=Plot[interpy,{x,0,26},PlotRange->{0,5},PlotStyle->Red];

Show[interpplot,dataplot]

Well, the polynomial does go through all the data, but it sure doesn't look like the function the data comes from.  At a quick glance, the data kind of looks like a quadratic.  If we wanted to find the best fit quadratic function (in a least-square sense) for this data we could use the Mathematica Fit command.

quadfity=Fit[data,{1,x,x^2},x]

Here the independent variable is x and the Fit command returns the "best" linear combination of the basis functions given as a list in it's second arguement.  By "best" we mean that the square-error function given in the introduction is minimized.  Now let's plot the quadratic with the data.

quadplot=Plot[quadfity,{x,0,26},PlotRange->{0,5},PlotStyle->Green];
Show[quadplot,dataplot]

Looking at both the data and the best fit quadratic reveals that the data may not be quadratic in nature (the data and the quadratic don't look symmetric).  Let's try a cubic equation:

cubefity=Fit[data,{1,x,x^2,x^3},x]

cubeplot=Plot[cubefity,{x,0,26},PlotRange->{0,5},PlotStyle->Blue];
Show[quadplot,cubeplot,dataplot]

Hopefully, it is obvious that the best fit cubic does a better job of approximating the data then the best fit quadratic.  This should not be surprising, because there are more degrees of freedom (more parameters) in the cubic than there are in the quadratic.  We could see this by calculating the square error in each case as follows.

quaderr=Sum[(data[[k,2]]-quadfity/.{x->data[[k,1]]})^2,{k,1,dim[[1]]}]

cubeerr=Sum[(data[[k,2]]-cubefity/.{x->data[[k,1]]})^2,{k,1,dim[[1]]}]

In actuality, the data came from the function

f(x)=((Exp[-x/20] x (26-x))/25)+((RandReal[])/2)

where the RandReal command generates a random number between 0 and 1.  Often it is the case that real world data has some random noise in it.

The idea behind examples in linear regression is to graph the data and try to get a feeling for how the data is behaving.  Does the data look linear, quadratic, etc.  Then choose a set of basis function which might be powers of x, exponential, or trigonometric functions.  Finally look at the best fit solution and decides if it truly describes the data.  If not, then try a different combinations of functions.

Exercise 8.

For problems 1-3 use the given data and the fit command to find the best linear and cubic polynomial fits.  Compare your results to the interpolating polynomials in the homework on lists.

1.

data1={{-1.0,1.0},{-0.75,1.0},{-0.5,1.0},{0.5,-1.0},{0.75,-1.0},{1.0,-1.0}};

2.

data2={{-1.0,1.0},{-0.75,1.0},{-0.1,1.0},{0.1,-1.0},{0.75,-1.0},{1.0,-1.0}};

3.

data3={{-1.0,1.0},{-0.75,1.0},{-0.01,1.0},{0.01,-1.0},{0.75,-1.0},{1.0,-1.0}};

4.

Import the data from the file "regressiondata.txt".  Using techniques similar to that above, graph the data then determine a good set of functions to fit the data.  Use the fit command to find the function and plot the data and function on the same graph.  Use another set of functions and repeat the problem.  You may want to include exponential and/or trigonometric functions along with power functions in one of your attempts.