Project Help Notebook by Keith Stroyan

Help NoteBook Copyright 1998 by Academic Press, Inc. - All rights reserved.
Modify this Notebook to prepare your project report.

Project 32: Least Squares Fit to a Linear Function

by A. Student and B. Pupil

with technical help from Da LabMonitor and Help Notebook

Initialization

Needs["Graphics`Colors`"];

My Project Introduction

Write this last.  Begin with the basic question and state your main finding clearly, but with few technical details.  You might even modify the title of this section to make it more interesting.

Beginning Explorations or Computer Experiments

Introductory Example of Section 32.1

Given the (x,y)  data points (1,6), (2,9), (3,10), find the parameters  m  and  b  that minimized the sum of the squared errors.  (Read the print project for more details.)

Mathematica's direct solution:

Clear[xData, yData, pairData, x, y]; xData = {1, 2, 3};
yData = {6, 9, 10};
pairData = Transpose[{xData, yData}];
y := Fit[pairData, {1, x}, x];
line = Plot[y, {x, 0, 4}];
Show[line, 
        Graphics[{AbsolutePointSize[4], Red, Map[Point, pairData]}]];
y

The sumed square error function with Mathematica:

Clear[L, x, error, m, b];
L[x_] :=  m x  +  b;
error[m_, b_] = 
      Sum[(((L[xData[[i]]] - yData[[i]]))^2), {i, 1, Length[xData]}]

Partial derivatives:

D[error[m, b], m]

D[error[m, b], b]

Critical values:

Solve[{D[error[m, b], m] == 0, D[error[m, b], b] == 0}]

Section 32.2 The Closed Form Solution

Read section 32.2 and supplement the program steps below with your description.

The Coefficients a, c, d, f, g

x1^2 + x2^2 + x3^2 + ...+ xn^2 = x dot x

a = xData.xData

x1 + x2 + x3 + ...+ xn =

c = Apply[Plus, xData]
d = c

e = n

e = Length[xData]

x1 y1 + x2 y2 + ... + xn yn = x dot y

f = xData.yData

y1 + y2 + y3 + ...+ yn =

g = Apply[Plus, yData]

The Linear System A.X  = B

The matrix  A

A = {{a,d},{c,e}}

The column  B

B = {f,g}

The solution X = {m , b}

{m,b} = LinearSolve[A, B]

Print["The minimizing function is  y = ",m,"x + ",b]

32.3 The General Algebraic Solution

Solve Hints 32.5, 6, 7 and write your answers here.

Extra Section 32.4: Least Squares by FindMinimum

This section does not appear in the print version of the project.  It shows how to find the least squares line by direct minimization.

Clear[pairData, xData, yData, L, m, b, x];
pairData = {{1, 6}, {2, 9}, {3, 10}};
{xData, yData} = Transpose[pairData];
L[x_] := m x  +  b;
yDifferences = Map[L, xData] - yData;
error[m_, b_] = yDifferences . yDifferences

minValues = FindMinimum[error[m, b], {m, 1}, {b, 1}];
mMin = m /. minValues[[2, 1]];
bMin = b /. minValues[[2, 2]];
minL[x_] = L[x] /. {m -> mMin, b -> bMin}

dataPlot = 
        ListPlot[pairData, PlotStyle -> {Red, AbsolutePointSize[4]}, 
          AxesLabel -> {"x", "y"}];
minPlot = 
        Plot[minL[x], {x, 0, 4}, AxesLabel -> {"x", "y"}, 
          PlotStyle -> {Blue}];
Show[dataPlot, minPlot];

Project 32.5: More General Minimization

Following is another kind of fit minimizing the summed square errors.  This time we fit data to a logistic growth curve and minimize the error based on choices of the parameters.  (See text Exercise 21.2.6 for the general logistic formula.)

Clear[pairData, xData, yData, L, r, c, k, x];
pairData = {{0, 1}, {1, 6}, {2, 9}, {3, 10}, {4, 11}};
{xData, yData} = Transpose[pairData];
L[x_] := 
        ((k (E^(r x)))/(1 + ((k (E^(r x)))/c)));
yDifferences = Map[L, xData] - yData;
error[r_, c_, k_] = yDifferences . yDifferences

minValues = FindMinimum[error[r, c, k], {r, 2}, {c, 12}, {k, 1}];
      
rMin = r /. minValues[[2, 1]];
cMin = c /. minValues[[2, 2]];
kMin = k /. minValues[[2, 3]];
minL[x_] = L[x] /. {r -> rMin, c -> cMin, k -> kMin}

dataPlot = 
        ListPlot[pairData, PlotStyle -> {Red, AbsolutePointSize[4]}, 
          AxesLabel -> {"x", "y"}];
minPlot = 
        Plot[minL[x], {x, 0, 6}, AxesLabel -> {"x", "y"}, 
          PlotStyle -> {Blue}];
Show[dataPlot, minPlot, PlotRange -> {0, 12}];

Notice that the partial derivatives are quite messy:

D[error[r, c, k], r]

D[error[r, c, k], c]

D[error[r, c, k], k]

Exact solution of the critical equations is difficult, but approximate roots near your guess is possible:

eqs = 
      {D[error[r, c, k], r] == 0, D[error[r, c, k], c] == 0, 
        D[error[r, c, k], k] == 0}
FindRoot[eqs, {r, 2}, {c, 10}, {k, 1.5}]

{rMin, cMin, kMin}

But other "solutions" may be found...

FindRoot[eqs, {r, 1}, {c, 1}, {k, 1}]

FindRoot[eqs, {r, 19.07}, {c, 9. }, {k, 1.125}]

Let's take a critical look at these critical values.

r = 21.07;
c = 9.0;
k = 1.125;
newPlot = 
        Plot[L[x], {x, 0, 6}, AxesLabel -> {"x", "y"}, 
          PlotStyle -> {Blue}];
Show[dataPlot, newPlot, PlotRange -> All];