Allman E.S., Rhodes J.A. Mathematical Models in Biology: An Introduction
Подождите немного. Документ загружается.
326 Curve Fitting and Biological Modeling
The most common means of picking the best-fit line is the method of least
squares. The philosophy of least squares is that the line that best fits data
is the one that minimizes SSE. Geometrically, the least-squares best-fit line
is the one that minimizes the sum of squares of the vertical distances between
the data points and the fitting line – of all the lines that could possibly describe
the data trend, we consider as best the one with this geometric property.
Note that one feature of this method is that it chooses the best line by a cri-
terion using all the data points.
With this criterion, the calculation of the best-fit line ultimately turns out to
be surprisingly simple. Understanding why the calculation works as it does,
though, requires a bit more effort.
If there are only two data points, then finding the least-squares best-fit
line through them is straightforward. We know that there is a line going
exactly through any two points, and that line will have SSE = 0, the minimum
possible value.
Although the algebra to find the line through two points can be formulated
in a number of different but equivalent ways, a matrix formulation will set
the stage for later work. Suppose, for instance, the data points are (3, 2.3) and
(6, 1.7). Then, because a line has equation y = mx + b, we need to find m
and b so that
2.3 = m · 3 + b
1.7 = m · 6 + b
or
31
61
m
b
=
2.3
1.7
. (8.1)
Solving the matrix equation gives
m
b
=
31
61
−1
2.3
1.7
=
−0.2
2.9
,
so the line fitting the data is y =−0.2x + 2.9. Because we solved the matrix
equation exactly, the line goes exactly through the two data points.
Suppose now we had three data points, (3, 2.3), (6, 1.7), and (9, 1.3). The
first two are the same as above, and thus lie on the line we just found. However,
the third data point is not on that line, but rather lies above it. If we are still
trying to find a line y = mx + b to fit this data, we would like to find a
8.2. The Method of Least Squares 327
solution to
2.3 = m · 3 + b
1.7 = m · 6 + b
1.3 = m · 9 + b
or
31
61
91
m
b
=
2.3
1.7
1.3
. (8.2)
Why can’t you attempt to solve this matrix equation by finding a matrix
inverse?
Because matrix inverses can exist only for square matrices, straightforward
matrix algebra is not sufficient to solve this equation. In fact, we know there
is no solution to this matrix equation – if there were, then the three data points
would lie exactly on a line. Since we cannot hope for an exact solution to
Eq. (8.2), our aim is to instead find values for m and b that minimize SSE.
More generally, suppose we want to find the equation of a line y =
ˆ
mx +
ˆ
b
that, of all lines, best fits the data points (x
1
, y
1
), (x
2
, y
2
),..., (x
n
, y
n
). We
would like a solution, (
ˆ
m,
ˆ
b), to a system of equations:
y
1
= mx
1
+ b
y
2
= mx
2
+ b
.
.
.
y
n
= mx
n
+ b,
which can be written in matrix form as
x
1
1
x
2
1
.
.
.
.
.
.
x
n
1
m
b
=
y
1
y
2
.
.
.
y
n
. (8.3)
However, this equation is unlikely to have a solution (
ˆ
m,
ˆ
b), because the
original data points are unlikely to lie exactly on a line. Instead of solving
this exactly, we want to find the values of
ˆ
m and
ˆ
b that “almost” satisfy it, in
the precise least-squares sense.
328 Curve Fitting and Biological Modeling
Although we do not yet know
ˆ
m and
ˆ
b, consider a line y = mx + b as a
candidate for the best fit one. Let
˜
y
i
= mx
i
+ b, i = 1, 2,...,n,
denote the y-coordinates of the points on this candidate line, with x-
coordinates given by x
i
. Then, the error vector for the candidate line will
be
e(m, b) = (y
1
−
˜
y
1
, y
2
−
˜
y
2
,...,y
n
−
˜
y
n
)
= (y
1
− mx
1
− b, y
2
− mx
2
− b,...,y
n
− mx
n
− b).
The total error, using the sum of squares measure, is then
SSE(m, b) = (y
1
− mx
1
− b)
2
+ (y
2
− mx
2
− b)
2
+···+(y
n
− mx
n
− b)
2
=
n
i=1
(y
i
− mx
i
− b)
2
.
Notice that the error vector and the total error depend on the choice of m and
b for the line we consider. Our goal is to find values
ˆ
m and
ˆ
b that minimize
this number among all possible choices of m and b.
We focus our attention on
ˆ
m first. If SSE(
ˆ
m,
ˆ
b) is minimal, then for any
choice of number m, the value of SSE(m,
ˆ
b) must be equal to or larger than
it. That is,
SSE(m,
ˆ
b) ≥ SSE(
ˆ
m,
ˆ
b),
n
i=1
(y
i
− mx
i
−
ˆ
b)
2
≥
n
i=1
(y
i
−
ˆ
mx
i
−
ˆ
b)
2
.
Now consider m =
ˆ
m + for some , to bring attention to the perturbation
of m from its optimal value
ˆ
m. Substituting this expression for m into the
inequality and rearranging terms gives
n
i=1
(y
i
−
ˆ
mx
i
− x
i
−
ˆ
b)
2
≥
n
i=1
(y
i
−
ˆ
mx
i
−
ˆ
b)
2
, or
n
i=1
(y
i
−
ˆ
mx
i
−
ˆ
b − x
i
)
2
− (y
i
−
ˆ
mx
i
−
ˆ
b)
2
≥ 0.
8.2. The Method of Least Squares 329
But the individual summands simplify as
(y
i
−
ˆ
mx
i
−
ˆ
b − x
i
)
2
− (y
i
−
ˆ
mx
i
−
ˆ
b)
2
= ((y
i
−
ˆ
mx
i
−
ˆ
b) − x
i
)
2
− (y
i
−
ˆ
mx
i
−
ˆ
b)
2
= (y
i
−
ˆ
mx
i
−
ˆ
b)
2
− 2x
i
(y
i
−
ˆ
mx
i
−
ˆ
b)
+ (x
i
)
2
− (y
i
−
ˆ
mx
i
−
ˆ
b)
2
=−2x
i
(y
i
−
ˆ
mx
i
−
ˆ
b) + (x
i
)
2
.
Therefore, the inequality above is
n
i=1
(−2x
i
(y
i
−
ˆ
mx
i
−
ˆ
b) + (x
i
)
2
) ≥ 0,
or
−2
n
i=1
x
i
(y
i
−
ˆ
mx
i
−
ˆ
b)
+
2
n
i=1
x
2
i
≥ 0.
In this inequality, considering a value of sufficiently close to zero, the second
term will be of negligible size in comparison with the first, due to the
2
. Thus,
for all small values of , the first term must be nonnegative for the inequality
to be satisfied. However, since might be either positive or negative, the only
way the first term is always nonnegative is if
n
i=1
x
i
(y
i
−
ˆ
mx
i
−
ˆ
b) = 0. (8.4)
This gives us an equation
ˆ
m and
ˆ
b must satisfy to minimize SSE. Because
it is an equation in only two unknowns (Remember: all the x
i
and y
i
are data
values), it is more simply expressed as
n
i=1
x
2
i
ˆ
m +
n
i=1
x
i
ˆ
b =
n
i=1
x
i
y
i
. (8.5)
After much work, we have found one equation that relates
ˆ
m and
ˆ
b.To
find a second equation relating
ˆ
m and
ˆ
b, we reason similarly focusing on
ˆ
b.
The complete argument is left for the exercises, but it yields
n
i=1
x
i
ˆ
m + n
ˆ
b =
n
i=1
y
i
. (8.6)
We now have two equations, (8.5) and (8.6), called the normal equations,
that relate
ˆ
m and
ˆ
b. With two equations in two unknowns, we can solve for
ˆ
m and
ˆ
b and so find the least-squares line y =
ˆ
mx +
ˆ
b.
330 Curve Fitting and Biological Modeling
Before continuing with the now routine calculation of the solutions to the
normal equations, a slight detour leads to a remarkable observation about the
structure of the normal equations. We need a definition first.
Definition. If M is an m × n matrix, then the n × m matrix obtained by
interchanging the rows and columns of M is known as the transpose of M,
and is denoted by M
T
; for
M =
x
11
x
12
··· x
1n
x
21
x
22
··· x
2n
.
.
.
.
.
.
.
.
.
.
.
.
x
m1
x
m2
··· x
mn
, M
T
=
x
11
x
21
··· x
m1
x
12
x
22
··· x
m2
.
.
.
.
.
.
.
.
.
.
.
.
x
1n
x
2n
··· x
mn
.
Example. If A =
31
61
91
, then A
T
=
369
111
. Notice the first row of A
becomes the first column of A
T
, the second row of A becomes the second
column of A
T
, and so on. At the same time, the columns of A have become
the rows of A
T
.
Let’s return to the original matrix Eq. (8.3) that we would have like to have
solved to find a line through the data points. If we multiply each side of the
equation on the left by the transpose of the matrix appearing in it, we obtain
x
1
x
2
... x
n
11... 1
x
1
1
x
2
1
.
.
.
.
.
.
x
n
1
m
b
=
x
1
x
2
... x
n
11... 1
y
1
y
2
.
.
.
y
n
,
which, on multiplying the matrices, gives
x
1
2
+ x
2
2
+···+x
n
2
x
1
+ x
2
+···+x
n
x
1
+ x
2
+···+x
n
1 + 1 +···+1
m
b
=
x
1
y
1
+ x
2
y
2
+···+x
n
y
n
y
1
+ y
2
+···+y
n
,
or, more succinctly,
x
2
i
x
i
x
i
n
m
b
=
x
i
y
i
y
i
, (8.7)
where the sums range over i = 1,...,n.
8.2. The Method of Least Squares 331
Now compare Eq. (8.7) to Eqs. (8.5) and (8.6). Amazingly, these equations
are exactly the same; Eq. (8.5) is stored in the top row of Eq. (8.7), while
Eq. (8.6) is in the bottom row. This observation provides a quick way to
perform the method of least-squares fitting: To find the least-squares solution
to a matrix equation of the form of Equation (8.3), multiply each side of the
equation on the left by the transpose of the matrix and solve the resulting
system for
ˆ
m and
ˆ
b.
To apply this to finding the least-squares, best-fit line for the three data
points (3, 2.3), (6, 1.7), and (9, 1.3), from Equation (8.2), we obtain the normal
equations
369
111
31
61
91
m
b
=
369
111
2.3
1.7
1.3
, or
126 18
18 3
m
b
=
28.8
5.3
.
Now multiplying both sides of the last equation by
126 18
18 3
−1
=
1
126 · 3 − 18 · 18
3 −18
−18 126
gives
ˆ
m
ˆ
b
=
1
126 · 3 − 18 · 18
3 −18
−18 126
28.8
5.3
≈
−.1667
2.7667
.
The least-squares best-fit line for the three data points is thus
y =−.1667x + 2.7667.
Graph this line and the three data points. Does the line appear to fit the
data well?
In using the least-squares approach to fit a line to data, the most important
point is that you understand the criteria that you are using to choose the
best line – the one with the smallest SSE. Of secondary importance is the
calculation you do to actually get that line. The steps for this are:
1. Write equations you would like m and b to satisfy for all the data
points to be on the line y = mx + b, by plugging each data point
into the equation. For n data points, this gives n equations in the two
unknowns, m and b, that usually cannot be solved exactly.
332 Curve Fitting and Biological Modeling
2. Express the equations in matrix form as
A
m
b
= b.
Here, A will be a matrix and b a vector, each with numerical entries.
3. Create the normal equations by multiplying on the left of both sides of
the equation by A
T
, giving
A
T
A
m
b
= A
T
b.
4. Solve the normal equations by computing A
T
A, A
T
b, and ( A
T
A)
−1
.
The solution is
ˆ
m
ˆ
b
= (A
T
A)
−1
A
T
b.
Notice that the steps here say nothing about the real ideas behind least-
squares – that appeared only in our derivation of the normal equations. How-
ever, the steps have the nice feature that they provide a simple and straight-
forward calculation.
Most software packages and calculators will calculate a least-squares, best-
fit line (often called a regression line) at the touch of a button. Once you
understand the idea and method of calculation, these are great labor-saving
devices.
Although the matrix A was not invertible in the three data point example
above, the matrix product A
T
A was invertible. Because of the particular form
of the columns of any matrix A used in least-squares regression, the product
A
T
A is almost always invertible, ensuring that a least-squares solution can be
found using matrix algebra. In fact, provided the data has at least two points
with different x-coordinates, A
T
A will be invertible, although a proof of this
fact requires additional theory from linear algebra. Moreover, when A
T
A
is invertible, there is one and only one solution to the normal equations. This
justifies talking about the least-squares, best-fit line for a data set; one line is
genuinely better than all others in giving a smallest value for SSE.
Problems
8.2.1. Plot the three points (−1, 1), (0, 3), and (1, 4). Then, find the least-
squares, best-fit line for them, following the four steps outlined in the
text and doing all calculations without a computer. Add a graph of
the line to your plot.
8.2. The Method of Least Squares 333
Table 8.3. Population Size in Year t
t 0123 4 5
P 173 278 534 895 1553 2713
8.2.2. Find the least-squares, best-fit line to the data points (3, 120), (4, 116),
(5, 114), (6, 109), and (7, 106) by:
a. following the four steps given in the text, using a computer. The
MATLAB command A' gives the transpose of a matrix A.
b. following the first two steps and then using the MATLAB com-
mand A\b to find the least-squares solution
ˆ
x to Ax = b.
c. Using the MATLAB command polyfit. For instructions, type
help polyfit.
8.2.3. Recall from the last section that the data of Table 8.1 showed an
exponential decay that we hoped to model by an exponential formula.
Table 8.2 contains transformed data that is roughly linear.
a. Find the least-squares, best-fit line ln y =
ˆ
mt +
ˆ
b to the data in
Table 8.2.
b. Use your answer to part (a) to give an exponential curve y = ae
kt
fitting the data in Table 8.1.
Note: This approach to fitting an exponential curve, using a least-
squares, best-fit line to the transformed data, does not necessarily
give the exponential that minimizes SSE for the untransformed data.
It is, however, a standard approach to exponential curve fitting.
8.2.4. Suppose the population data in Table 8.3 is believed to be described
by the model P
t+1
= λP
t
.
a. Produce a semilog plot and explain why it justifies the choice of
the model.
b. Find the least-squares, best-fit line to the transformed data.
c. Use part (b) to find an exponential curve fitting the data.
d. Use part (c) to give a good estimate of λ for this data.
8.2.5. To produce and plot simulated data points that will be nearly on the
line y = .7x + 2.1, use the MATLAB commands
x=[1:10]', y=.7*x+2.1+.3*randn(10,1),
plot(x,y,'o').
Then A=[x,x.∧0], b=y will prepare you to perform the least-
squares, line-fitting calculation.
334 Curve Fitting and Biological Modeling
a. Enter all these commands and recover the least-squares line. Is it
exactly the line y = .7x + 2.1? Is it close? Perform this experiment
several times and summarize your results.
b. What is the effect of the number .3 in these commands? If .3isre-
placed by 3, is the line that you recover usually more or less similar
to y = .7x + 2.1? Explain why your observation is reasonable.
c. What is the effect of using fewer or more data points on the recov-
ery of the line? For instance, if the number 10 in these commands
is replaced with 3 or 30, is the line you recover usually more or
less similar to y = .7x + 2.1? Explain why your observation is
reasonable.
8.2.6. That there is exactly one straight line through any two points is well
known. However, this fact manifested itself in the fact that Equation
(8.1) was solvable and had a unique solution.
a. Explain why, through any three points in the plane, you would
expect to be able to find exactly one parabola of the form y =
ax
2
+ bx + c by expressing the equations you would need to solve
to find a, b, and c in matrix form. What is it about the matrix
equation that suggests there is probably one and only one solution?
b. Given n points, what degree polynomial y = p(x) should you
consider to be likely to find one and only one such polynomial
curve through the data points? Explain.
8.2.7. Consider the four data points (−2, 8.1), (0, 7), (10, 5.9), and (15, 5).
a. Use MATLAB to plot the four data points and the least-squares
line fitting them.
b. Calculate the mean x-coordinate
¯
x and mean y-coordinate
¯
y. Does
the least-squares line pass through (
¯
x,
¯
y)?
c. Perform several similar experiments by varying the data points.
8.2.8. Show that
ˆ
m and
ˆ
b for the least-squares, best-fit line to data must
satisfy Eq. (8.6) as follows:
a. Explain why SSE(
ˆ
m, b) ≥ SSE(
ˆ
m,
ˆ
b) for any choice of b.
b. Using b =
ˆ
b + δ, show that for all δ
n
i=1
(y
i
−
ˆ
mx
i
−
ˆ
b − δ)
2
− (y
i
−
ˆ
mx
i
−
ˆ
b)
2
≥ 0.
c. Show that for all δ
−2δ
n
i=1
(y
i
−
ˆ
mx
i
−
ˆ
b)
+ δ
2
n ≥ 0.
8.3. Polynomial Curve Fitting 335
d. Explain why part (c) shows
n
i=1
(y
i
−
ˆ
mx
i
−
ˆ
b) = 0.
e. Deduce Eq. (8.6).
8.2.9. (Calculus) The normal equations for least-squares, line fitting can
also be derived using calculus. Recall that at the minimum of a dif-
ferentiable function, the derivative must be zero.
a. Derive Eq. (8.4) by differentiating SSE(m,
ˆ
b) with respect to m
and setting the result equal to zero.
b. Derive Eq. (8.6) similarly.
8.2.10. Because the normal Eqs. (8.7) can be solved by inverting a 2 × 2
matrix, the formula for the matrix inverse leads to a formula for
ˆ
m
and
ˆ
b. Use this approach to deduce
ˆ
m =
n
x
i
y
i
−
x
i
y
i
n
x
2
i
−
x
i
2
ˆ
b =
x
2
i
y
i
−
x
i
x
i
y
i
n
x
2
i
−
x
i
2
.
8.2.11. Use the result of the last problem to show that, if
¯
x an
¯
y are the
mean x- and y-coordinates of the data points, then (
¯
x,
¯
y) lies on the
least-squares, best-fit line.
8.2.12. The least-squares solution to the equation Ac = b, for a vector of
unknowns c,isgivenby
ˆ
c = ( A
T
A)
−1
A
T
b. You might think this
formula could be simplified as
ˆ
c = (A
T
A)
−1
A
T
b
= A
−1
(A
T
)
−1
A
T
b
= A
−1
I b
= A
−1
b.
Explain the error in this reasoning.
8.3. Polynomial Curve Fitting
While least-squares fitting of straight lines to data is very commonly done,
fitting higher degree polynomials is often useful as well. If we have decided
to model some data by an nth degree polynomial y = f (t), then finding the