18 In-Class Assignment: Least Squares Fit (LSF)
Contents
18 In-Class Assignment: Least Squares Fit (LSF)#
0. Final Exam Review Survey#
The two class periods during the week of 4/24-4/28 will be used to review for the final exam. There will not be any PCAs for next week, and there are no graded ICAs next week either, but instructors will likely have some sort of review activity for you to work on in class. To help instructors decide what material to focus on during class next week, please fill out the following survey by this Friday 4/21 https://forms.office.com/r/XfSSM3mQnt.
1. LSF Pre-class Review#
19–LSF_pre-class-assignment.ipynb
2. Finding the best solution in an overdetermined system#
Let
Note we substitute
In other words,
We also know that the dot product is zero if two vectors are orthogonal. So we have
$
The columns of
That is:
$
The above equation is called the least squares solution to the original equation
The matrix
Example: A researcher has conducted experiments of a particular Hormone dosage in a population of rats. The table shows the number of fatalities at each dosage level tested. Determine the least squares line and use it to predict the number of rat fatalities at hormone dosage of 22.
Hormone level |
20 |
25 |
30 |
35 |
40 |
45 |
50 |
---|---|---|---|---|---|---|---|
Fatalities |
101 |
115 |
92 |
64 |
60 |
50 |
49 |
%matplotlib inline
import matplotlib.pylab as plt
import numpy as np
import sympy as sym
import time
from answercheck import checkanswer
sym.init_printing(use_unicode=True)
H = [20,25,30,35,40,45,50]
f = [101,115, 92,64,60,50,49]
plt.scatter(H,f)
plt.xlabel('Hormone Level')
plt.ylabel('Fatalities')
We want to determine a linear model that is expressed by the following equation
to approximate the connection between Hormone dosage (
✅DO THIS: Rewrite the system of equations in matrix form numpy
matrices A
and y
using the data from above:
# put your code here
checkanswer.matrix(A,'50278283443c1d502428e602f0d363db');
checkanswer.matrix(y, '37e2ed57a1516fb4a17eb2a3e9e99d2d');
✅ QUESTION: Calculate the square matrix AtA
) and the modified right hand side vector as Aty
):
#Put your code here
# AtA=
# Aty=
checkanswer.matrix(AtA,'3f579409e55e9cc409bc727c35e09c74');
checkanswer.matrix(Aty,'1dec77018852144a8b90600161c8bff6');
✅QUESTION: Find the least squares solution by solving
# Put your code here
# x=
checkanswer.matrix(x,'3213b988975659ed46496a07542eff33');
✅QUESTION: Given the solution above, define the two scalars y-intercept a
and slope b
.
#put your code here
#a=
#b=
checkanswer.float(a,'87e53cb122536f53f434cdbccf0aca94');
checkanswer.float(b,'f7dfa5a5ed0f65f0e0a87f29efd9cd74');
The following code will plot the original data and the linear model whose coefficients we found by performing a least squares fit.
H = [20,25,30,35,40,45,50]
f = [101,115, 92,64,60,50,49]
plt.scatter(H,f)
H2 = np.linspace(np.min(H), np.max(H))
f2 = a + b*H2
plt.plot(H2, f2)
✅QUESTION: Now, let’s fit a quadratic model instead of a linear model to our data, i.e. $
To do this, the vector of unknowns is now
Again, put this system of equations into matrix form
# put your code here
checkanswer.float(a,'8acc26f9f9da4edfa44529a811ebb8a9')
checkanswer.float(b,'b1e317387458d5471bd7cec94e4959fa')
checkanswer.float(c,'451fd8dc5a793f8483160363b42fa8cd')
✅QUESTION: Now, plot the original data along with the curve for the quadratic model
# put your code here
✅QUESTION: Now, see if you can generalize your code above to fit a degree-
# put your code here
d = 3 # Set d to be 3, 4, 5, or 6
# Form a len(H)x(d+1) matrix A whose entries are A[i,j] = H[i]**j
# Solve Ax = y using least squares
# Plot the original data points along with the curve for the polynomial model
✅QUESTION: The interactive function below allows you to adjust the degree of the polynomial model as well as the limits of the x_min
and x_max
parameters. Do you think that an eighth-order polynomial is a good model for this dataset? Why or why not?
from ipywidgets import interact, fixed
import ipywidgets as widgets
@interact(x=fixed(H), y=fixed(f), degree=widgets.IntSlider(min=1, max=8, step=1, value=8), x_min=widgets.IntSlider(min=min(H)-10, max=min(H), step=1, value=min(H)), x_max=widgets.IntSlider(min=max(H), max=max(H)+10, step=1, value=max(H)))
def graphPolyN(x, y, x_min, x_max, degree):
p = np.polyfit(x, y, degree)
f = np.poly1d(p)
x_pred = np.linspace(x_min, x_max, 1000)
y_pred = f(x_pred)
plt.scatter(x, y, color="red")
plt.plot(x_pred, y_pred)
Put your answer to the above question here
3. Pseudoinverse#
In this class we often talk about solving problems of the form:
Currently we have determined that this problem becomes very nice when the
Since
Now, let us consider a a more general problem where the
Assuming we can find the
How do we know there is a Psudoinverse#
Assuming the general case of a
The rowspace of
is in with dimensionThe columnspace of
is in also with dimension .The nullspace of
is in with dimensionThe nullspace of
is in with dimension .
Because the rowspace of
For any
in the rowspace, we have that is one point in the columnspace. If is another point in the row space different from , we have (The mapping is one-to-one).For any
in the columnspace, we can find in the rowspace such that (The mapping is onto).
The above is not really a proof but hopefully there is sufficient information to convince yourself that this is true.
How to compute pseudoinverse#
We want to find the
Let’s apply SVD on
: $ U m\times m V^\top n\times n \Sigma m\times n U_1 m\times r U_2 m\times (m-r) \Sigma_1 r\times r V_1^\top r\times n V_2^\top (n-r)\times n$.The columnspace of
is the columnspace of , and columnspace of is the nullspace of .The rowspace of
is the rowspace of , and rowspae of is the nullspace of .
If
is in the rowspace of , we have that . We have .If we define a matrix
, we have that . That is is is in the rowspace of .
The matrix
is the pseudoinverse of matrix . $ $
Example 1: Let $
A = np.matrix([[1,0,1],[1,2,3]])
✅TODO: Calculate the pseudoinverse numpy.linalg
function pinv
:
# put your code here
# Apinv=
checkanswer.matrix(Apinv,'7724dc0d898980cb5a80dac73bfa7aa7');
✅DO THIS: Compute
#put your code here
✅QUESTION: If
Put your answer to the above question here
✅QUESTION: If
Put your answer to the above question here
Left inverse is pseudoinverse#
We can compute the left inverse of
In this case, the SVD of
The left inverse of
Example 2: Let $
A = np.matrix([[1,1],[0,2],[1,3]])
A
✅DO THIS: Calculate the pseudoinverse
#Put your answer here
# pInv_A =
checkanswer.matrix(pInv_A,'fda9d63c9e6cdc7bdb2acc57c4bdf490');
✅DO THIS: Calculate the left inverse of
#Put your anaswer here
# leftInv_A =
checkanswer.matrix(leftInv_A,'fda9d63c9e6cdc7bdb2acc57c4bdf490');
Written by Dr. Dirk Colbry with interactive code David Yonkers, Michigan State University. Some modifications by Dr. Santhosh Karnik, Michigan State University
This work is licensed under a Creative Commons Attribution-NonCommercial 4.0 International License.