```
In [2]:
```%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
n = 20
sigma = 1.0
xdata = np.linspace(-2, 2, n)
fdata = 3*xdata**2 + 2*xdata + 1 + np.random.randn(n)*sigma
plt.figure()
plt.plot(xdata, fdata, 'o')
plt.xlabel('x')
plt.ylabel('f')
plt.show()

```
```

This data is our training data and consists of two pairs: $(x^{(i)}, f^{(i)})$. Given this data, we want to compute a polynomial fit. In this case, we know that it's quadratic, so let's use a quadratic model of the form:

$$\hat{f}(x) = a x^2 + b x + c$$This is a simple least squares problem.

$$\text{minimize} \sum_i \left(\hat{f}(x^{(i)}) - f^{(i)} \right)^2 $$Which we can rewrite in matrix form as:

$$\text{minimize} \; || \Psi w - f ||^2 $$where

$$ \Psi = \begin{bmatrix} \text{---} & \psi(x^{(1)})^T & \text{---} \\ \text{---} & \psi(x^{(2)})^T & \text{---} \\ & \vdots \\ \text{---} & \psi(x^{(n)})^T & \text{---} \\ \end{bmatrix} $$$$ f = \begin{bmatrix} f^{(1)}\\ f^{(2)}\\ \vdots \\ f^{(n)}\\ \end{bmatrix} $$In our specific case: $w = [a, b, c]$ and $\psi = [x^2, x, 1]$. The matrix equation can be solved as a least squares solution. Notice that this simple form works in any number of dimensions (not just one). In higher dimensions we just add more terms to $\psi$.

```
In [3]:
```Psi = np.zeros((n, 3))
Psi[:, 0] = xdata**2
Psi[:, 1] = xdata
Psi[:, 2] = 1
w, residuals, rank, s = np.linalg.lstsq(Psi, fdata)
print w

```
```

```
In [4]:
```nfull = 200
xfull = np.linspace(-2, 2, nfull)
Psifull = np.zeros((nfull, 3))
Psifull[:, 0] = xfull**2
Psifull[:, 1] = xfull
Psifull[:, 2] = 1
ffull = np.dot(Psifull, w)
plt.figure()
plt.plot(xdata, fdata, 'o')
plt.plot(xfull, ffull, '--')
plt.xlabel('x')
plt.ylabel('f')
plt.show()

```
```

```
In [5]:
```def func(x):
sigma = 1.0
return (6*x-2)**2*np.sin(12*x-4) + np.random.randn(len(x))*sigma
# ---- create training data ---------
ndata = 20
xdata = np.linspace(0, 1, ndata)
fdata = func(xdata)
plt.plot(xdata, fdata, 'o')
plt.xlabel('x')
plt.ylabel('f')
plt.show()

```
```

```
In [6]:
```def getPsi(x, order):
n = len(x)
Psi = np.zeros((n, order+1))
for i in range(order+1):
Psi[:, i] = x**(order-i)
return Psi
def createpoly(x, f, order):
Psi = getPsi(x, order)
w, residuals, rank, s = np.linalg.lstsq(Psi, f)
return w
def evalpoly(x, w):
order = len(w) - 1
Psi = getPsi(x, order)
f = np.dot(Psi, w)
return f

Let's try different polynomial orders and check to see what the error is in our fit.

```
In [7]:
```ordervec = np.arange(1, 21)
error = np.zeros(20)
for idx, order in enumerate(ordervec):
# build a polynomial model from the training data
w = createpoly(xdata, fdata, order)
# test the error
fhat = evalpoly(xdata, w)
error[idx] = np.linalg.norm(fhat - fdata)
# plot error
plt.figure()
plt.plot(ordervec, error, 'o')
plt.xlabel('order of polynomial')
plt.ylabel('error')
plt.show()

```
```

```
In [8]:
```order = 20
w = createpoly(xdata, fdata, order)
nhat = 200
xhat = np.linspace(0, 1, nhat)
fhat = evalpoly(xhat, w)
plt.figure()
plt.plot(xdata, fdata, 'o')
plt.plot(xhat, fhat, '--')
plt.ylim([-10, 20])
plt.xlabel('x')
plt.ylabel('f')

```
Out[8]:
```

```
In [19]:
```div = 10 # we will divide our data into div segments
ndata = 20 # number of data points
arrlength = ndata/div # each segment should contain this much data
idxrand = np.random.permutation(n) # random index into data from 0 ... n
error = np.zeros(len(ordervec))
# iterate through polynomial orders
for i, order in enumerate(ordervec):
# iterate through divisions of data for k-holdout
for j in range(div):
# indicies of data to leave out from the random permutations
holdout = idxrand[arrlength*j:arrlength*(j+1)]
# separaet into training set and testing set
xtrain = np.delete(xdata, holdout)
ftrain = np.delete(fdata, holdout)
xtest = xdata[holdout]
ftest = fdata[holdout]
# build a polynomial model from the training data
w = createpoly(xtrain, ftrain, order)
# test the error with the validation set
fhat = evalpoly(xtest, w)
error[i] += np.linalg.norm(fhat - ftest) / div # average error across divisions
# plot error
plt.figure()
plt.plot(ordervec, error, 'o')
plt.xlabel('order of polynomial')
plt.ylabel('error')
plt.show()
# plot error
plt.figure()
plt.plot(ordervec, error, 'o')
plt.xlabel('order of polynomial')
plt.ylabel('error')
plt.ylim([0, 25])
plt.show()

```
```

```
In [20]:
```order = 4
w = createpoly(xdata, fdata, order)
nhat = 200
xhat = np.linspace(0, 1, nhat)
fhat = evalpoly(xhat, w)
plt.figure()
plt.plot(xdata, fdata, 'o')
plt.plot(xhat, fhat, '--')
plt.ylim([-10, 20])
plt.xlabel('x')
plt.ylabel('f')

```
Out[20]:
```

We can see that this results in a much predictive model than the 20th order polynomial.