In [13]:
from mth309 import *
import matplotlib.pyplot as plt
%config InlineBackend.figure_format = 'retina'
%matplotlib inline

First cook up some data that is quadratic-ish:

In [14]:
x = random.randint(0,100,12)
x
Out[14]:
array([39, 11, 99, 50, 65, 90, 68, 34, 20, 19, 98, 51])
In [22]:
y = 100 + 200*x - 3*x*x + random.randint(-400,400,12)
y
Out[22]:
array([ 3535,  1574, -9344,  2870,   481, -6289,  -348,  3201,  3190,
        3207, -8937,  2438])
In [23]:
plt.plot(x,y,'mo');

Form the normal equations

In [24]:
X = vstack((ones_like(x),x,x*x)).T
X
Out[24]:
array([[   1,   39, 1521],
       [   1,   11,  121],
       [   1,   99, 9801],
       [   1,   50, 2500],
       [   1,   65, 4225],
       [   1,   90, 8100],
       [   1,   68, 4624],
       [   1,   34, 1156],
       [   1,   20,  400],
       [   1,   19,  361],
       [   1,   98, 9604],
       [   1,   51, 2601]])
In [25]:
XTX = dot(X.T,X)
XTX
Out[25]:
array([[       12,       644,     45014],
       [      644,     45014,   3602012],
       [    45014,   3602012, 310108358]])
In [26]:
XTy = dot(X.T,y)
XTy
Out[26]:
array([     -4422,   -1702707, -202711809])

Solve the normal equations

In [27]:
beta = linalg.solve(XTX,XTy)
beta
Out[27]:
array([158.06700616, 199.2514382 ,  -2.99099681])

Make a plot to see if the fit looks good:

In [28]:
plt.plot(x,y,'mo')
xx = linspace(min(x),max(x),100)
plt.plot(xx,beta[0]+beta[1]*xx + beta[2]*xx*xx,'k',lw=3,alpha=0.3);