# Polynomial fitting with the normal equations¶

In [3]:
import numpy as np
import numpy.linalg as la
import matplotlib.pyplot as pt


In this demo, we will produce data from a simple parabola as a "model" and try to recover the "unknown" parameters $\alpha$, $\beta$, and $\gamma$ using least squares.

In [6]:
alpha = 3
beta = 2
gamma = 2

def f(x):
return alpha*x**2 + beta*x + gamma

plot_grid = np.linspace(-3, 3, 100)

pt.plot(plot_grid, f(plot_grid))

Out[6]:
[<matplotlib.lines.Line2D at 0x7f6900e06748>]
In [7]:
npts = 5

np.random.seed(22)
points = np.linspace(-2, 2, npts) + np.random.randn(npts)
values = f(points) + 0.3*np.random.randn(npts)*f(points)

pt.plot(plot_grid, f(plot_grid))
pt.plot(points, values, "o")

Out[7]:
[<matplotlib.lines.Line2D at 0x7f6900de0ac8>]

Now build the Vandermonde matrix:

In [9]:
A = np.array([
np.ones(npts),
points,
points**2
]).T
print(A)

[[ 1.         -2.09194992  4.37625447]
[ 1.         -2.46335065  6.06809644]
[ 1.          1.08179168  1.17027324]
[ 1.          0.76067483  0.5786262 ]
[ 1.          1.50887086  2.27669128]]


And solve for x using the normal equations:

In [10]:
x = la.solve(A.T@A, A.T@values)
x

Out[10]:
array([ 0.13178048,  1.84869187,  3.45132033])

Lastly, pick apart x into alpha_c, beta_c, and gamma_c:

In [11]:
gamma_c, beta_c, alpha_c = x

In [14]:
print(alpha, alpha_c)
print(beta, beta_c)
print(gamma, gamma_c)

3 3.45132033448
2 1.84869187248
2 0.131780484241

In [13]:
def f_c(x):
return alpha_c*x**2 + beta_c*x + gamma_c

pt.plot(plot_grid, f(plot_grid), label="true")
pt.plot(points, values, "o", label="data")
pt.plot(plot_grid, f_c(plot_grid), label="found")
pt.legend()

Out[13]:
<matplotlib.legend.Legend at 0x7f6900cdfb70>

(Edit this cell for solution.)