# create a vector of values for problem 2 on Project 1
P = vector(RDF, 11); P[0] = 25.7; P[1] = 18.5; P[2] = 14.3; P[3] = 12.2; P[4] = 11.1; P[5] = 9.5; P[6] = 9.1; P[7] = 8.5; P[8] = 8.1; P[9] = 7.5; P[10] = 7.4;
|
|
# prepare lists of x and y coordinates of n points
n=10;
x = [P[i] for i in range(n)];
y = [P[i+1] for i in range(n)];
(x,y)
([25.7, 18.5, 14.3, 12.2, 11.1, 9.5, 9.1, 8.5, 8.1, 7.5], [18.5, 14.3, 12.2, 11.1, 9.5, 9.1, 8.5, 8.1, 7.5, 7.4]) ([25.7, 18.5, 14.3, 12.2, 11.1, 9.5, 9.1, 8.5, 8.1, 7.5], [18.5, 14.3, 12.2, 11.1, 9.5, 9.1, 8.5, 8.1, 7.5, 7.4]) |
# prepare lists of products xx, xy and and list of 1's
xx = [x[i]*x[i] for i in range(n)];
xy = [P[i]*P[i+1] for i in range(n)];
one= [1 for i in range(n)];
(xx, xy, one)
([660.49, 342.25, 204.49, 148.84, 123.21, 90.25, 82.81, 72.25, 65.61, 56.25], [475.45, 264.55, 174.46, 135.42, 105.45, 86.45, 77.35, 68.85, 60.75, 55.5], [1, 1, 1, 1, 1, 1, 1, 1, 1, 1]) ([660.49, 342.25, 204.49, 148.84, 123.21, 90.25, 82.81, 72.25, 65.61, 56.25], [475.45, 264.55, 174.46, 135.42, 105.45, 86.45, 77.35, 68.85, 60.75, 55.5], [1, 1, 1, 1, 1, 1, 1, 1, 1, 1]) |
#create a matrix which encodes the normal equation for least squares fitting of the points (x,y)
A =matrix(RDF, 2,3, [[sum(xx), sum(x), sum(xy)], [sum(x), sum(one), sum(y)]]);
A
[1846.45 124.5 1504.23] [ 124.5 10.0 106.2] [1846.45 124.5 1504.23] [ 124.5 10.0 106.2] |
# solve the system of equations by Gaussian elimination
# the value for the variable associated to the ith column is the last entry in the ith row
# here, for example, we'll get m/b in the last entry of the first row/seconf row
# this function works even if there are no solutions or infinitely many but further interpretation is needed
AE = A.echelon_form();
AE
[ 1.0 0.0 0.614118242388] [ 0.0 1.0 2.97422788226] [ 1.0 0.0 0.614118242388] [ 0.0 1.0 2.97422788226] |
# WARNING: indexing starts at ZERO in SAGE so to get m and b out of AE we need to say
m= AE[0,2];# not [1,3] !!
b= AE[1,2];
(m,b)
(0.614118242388, 2.97422788226) (0.614118242388, 2.97422788226) |
Q = [[x[i],y[i]] for i in range(n)]
|
|
# See http://www.sagemath.org/doc/reference/sage/plot/scatter_plot.html
# create a "scatter_plot" object from our list
pplot= scatter_plot(Q);
|
|
# Another trick: I need to define a "variable" as a var to get SAGE to treat it as one
z = var('z')
#now I can define the line I guessed when demonstrating plotting
guessline= plot((4/7)*(z-8) + 8, (z,0,27),rgbcolor=hue(0.4))
# also the least squares which is make longer so I can tell which it is
leastsquaresline= plot(m*z+b, (z,-3,30),rgbcolor=hue(0.8))
# finally plot the points and both lines
show(pplot+guessline+leastsquaresline)
|
|
# conclusion: my guess was not bad, and running the line through (8,8) was
# pretty much right on, but the slope was a bit too small to give optimal fit
|
|