```
In [2]:
```%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from bias_and_variance import print_equation, format_equation
from ipywidgets import interact

We wish to evaluate the performance of a learning algorithm that takes a set of $(x,y)$ pairs as input and selects a function $g(x)$, which predicts the value of $y$ for a given $x$ (that is, for a given $(x,y)$ pair, $g(x)$ should approximate $y$).

We will evaluate the 'fit' of the functions that our learning algorithm selects using the *mean squared error* (MSE) between $g(x)$ and $y$. For a set of $n$ data points $\{ (x_1, y_1), \ldots, (x_n,y_n) \}$, the MSE associated with a function $g$ is calculated as

The set of candidate functions that we will allow our algorithm to consider is the set of $k$th-order polynomials. This means that our hypothesis space will contain all functions of the form $g(x) = p_kx^k + p_{k-1} x^{k-1} + \ldots + p_{1} x + p_0$. These functions are entirely characterized by their coefficient vector ${\bf p} = (p_k, p_{k-1}, \ldots, p_0)$, and become much more flexible as $k$ increases (think about the difference between linear ($k=1$), quadratic ($k=2$), and cubic ($k=3$) functions). If we are given a set of $(x,y)$ pairs, it is straightforward to find the $k$th-order polynomial that minimizes the MSE between $g(x)$ and the observed $y$ values.

```
In [3]:
```data = np.load("data/xy_data.npy")
# only show the first ten points, since there are a lot
data[:10]

```
Out[3]:
```

You will use this data to train and evaluate your learning algorithm.

A "learning algorithm" which finds a polynomial of given degree that minimizes the MSE on a set of $(x,y)$ coordinates is implemented in Python with the `np.polyfit`

command. You can use `p = np.polyfit(x, y, k)`

to return the coefficient vector ${\bf p}$ for the $k$th-order polynomial $g$ that best fits (i.e., has the smallest $MSE$ on) the $x$ and $y$ coordinates in `x`

and `y`

.

For example, to fit a 4th-order polynomial to the data:

```
In [4]:
```# fit the 4th order polynomial
p = np.polyfit(data[:, 0], data[:, 1], 4)
print("Vector of coefficients: " + str(p))
# display the resulting equation
print_equation(p)

```
```

`np.polyval`

. For example, if you wanted to compute $g(x)$ at $x=0$, $x=0.5$, and $x=1$, you could do:

```
In [5]:
```np.polyval(p, np.array([0, 0.5, 1]))

```
Out[5]:
```

```
In [6]:
```# first load the data
data = np.load("data/xy_data.npy")
@interact
def make_polynomial_fit_and_graph(polynomial_order=(0, 9), training_set_index=(1, 11)):
"""Finds the best-fitting polynomials for k = {0, ... , 9},
using one of eleven different training datasets.
"""
# relabel the parameters
k = polynomial_order
i = training_set_index
# pull out the x and y values
x = data[((i - 1) * 10):(i * 10), 0]
y = data[((i - 1) * 10):(i * 10), 1]
# create the figure
fig, axis = plt.subplots()
# create a range of values for x between 0 and 1
plotx = np.arange(0, 1.01, 0.01)
# find the coefficients p
p = np.polyfit(x, y, k)
# find the values of the polynomial parameterized by p and
# evaluated for the points plotx
ploty = np.polyval(p, plotx)
# plot the fitted function
axis.plot(plotx, ploty, 'b-', label="${}$".format(format_equation(p)))
# plot the original data points
axis.plot(x, y, 'ko')
# set the axis limits
axis.set_xlim(0, 1)
axis.set_ylim(0, 0.35)
# put a title on each plot
axis.set_title('Dataset #{} fitted with k = {}'.format(i, k))
# create a legend
axis.legend(loc='upper left', frameon=False)

```
```

*training* dataset and a *testing* dataset.

```
In [7]:
```def mse(k, train, test):
"""Fits a polynomial with order `k` to a training dataset, and
then returns the mean squared error (MSE) between the y-values
of the training data and the fitted polynomial, and the MSE
between the y-values of the test data and the fitted polynomial.
Your answer can be done in 6 lines of code, including the return
statement.
Parameters
----------
k : integer
The polynomial order
train : numpy array with shape (n, 2)
The training data, where the first column corresponds to the
x-values, and the second column corresponds to the y-values
test : numpy array with shape (m, 2)
The testing data, where the first column corresponds to the
x-values, and the second column corresponds to the y-values
Returns
-------
a 2-tuple consisting of the training set MSE and testing set MSE
"""
### BEGIN SOLUTION
# compute the polynomial fit
p = np.polyfit(train[:, 0], train[:, 1], k)
# compute predictions and MSE for the training data
train_prediction = np.polyval(p, train[:, 0])
train_mse = np.mean((train_prediction - train[:, 1]) ** 2)
# compute predictions and MSE for the testing data
test_prediction = np.polyval(p, test[:, 0])
test_mse = np.mean((test_prediction - test[:, 1]) ** 2)
return train_mse, test_mse
### END SOLUTION

```
In [8]:
```# load the data
data = np.load("data/xy_data.npy")
# compute the MSE
train_mse, test_mse = mse(2, data[:10], data[10:])
print("The training error is: " + str(train_mse))
print("The testing error is: " + str(test_mse))

```
```

```
In [9]:
``````
# add your own test cases here!
```

```
In [10]:
```"""Test that the `mse` function is correct."""
from numpy.testing import assert_allclose
data = np.load("data/xy_data.npy")
# use first ten, and the remaining
assert_allclose(mse(2, data[:10], data[10:]), (0.000229744955167, 0.000374298371336))
assert_allclose(mse(3, data[:10], data[10:]), (0.000169612346303, 0.000463251756094))
assert_allclose(mse(9, data[:10], data[10:]), (1.46448764925e-21, 0.337001581723), atol=1e-20)
# use half-and-half
assert_allclose(mse(2, data[:55], data[55:]), (0.00034502281024316553, 0.00037620706341530435))
assert_allclose(mse(3, data[:55], data[55:]), (0.0003378190977339938, 0.00039980736728858482))
assert_allclose(mse(9, data[:55], data[55:]), (0.00026755111091101571, 0.00061531514687572487))
# use last twenty, and the remaining
assert_allclose(mse(2, data[-20:], data[:-20]), (0.00030881029910697136, 0.00040876086505745344))
assert_allclose(mse(3, data[-20:], data[:-20]), (0.00021713262385879197, 0.00055653317636801015))
assert_allclose(mse(9, data[-20:], data[:-20]), (0.00012210662449207329, 0.00071987940235435685))
print("Success!")

```
```

```
In [11]:
```def plot_mse(axis, max_order, train, test):
"""Plot the mean squared error (MSE) for the given training and testing
data as a function of polynomial order.
* Your plot should show the MSE for 0 <= k < max_order
* There should be two lines: one black, for the training set error, and
one red, for the testing set error.
* Make sure to include labels for the x- and y- axes.
* Label the training error and testing error lines as "Training set error"
and "Testing set error", respectively. These labels will be used to
create a legend later on (and so you should NOT actually create the
legend yourself -- just label the lines).
Your answer can be done in 10 lines of code, including the return statement.
Parameters
----------
axis : matplotlib axis object
The axis on which to plot the MSE
max_order : integer
The maximum polynomial order to compute a fit for
train : numpy array with shape (n, 2)
The training data, where the first column corresponds to the
x-values, and the second column corresponds to the y-values
test : numpy array with shape (m, 2)
The testing data, where the first column corresponds to the
x-values, and the second column corresponds to the y-values
Returns
-------
numpy array with shape (max_order, 2)
The MSE for the training data (corresponding to the first column) and
for the testing data (corresponding to the second column). Each row
corresponds to a different polynomial order.
"""
### BEGIN SOLUTION
k = np.arange(0, max_order)
# compute error for all values of k
error = np.empty((max_order, 2))
for i in range(max_order):
error[i] = mse(k[i], train, test)
axis.plot(k, error[:, 0], 'k-', label="Training set error")
axis.plot(k, error[:, 1], 'r-', label="Testing set error")
axis.set_xlabel("Polynomial model order (k)")
axis.set_ylabel("Mean squared error")
return error
### END SOLUTION

`plot_mse`

function, you should be able to see the error as a function of the polynomial order for both the training set and the test set:

```
In [12]:
```# load the data
data = np.load("data/xy_data.npy")
# plot it
fig, axis = plt.subplots()
plot_mse(axis, 15, data[:20], data[20:])
axis.legend(loc='upper left')

```
Out[12]:
```

```
In [13]:
``````
# add your own test cases here!
```

```
In [14]:
```"""Is the plot_mse function correctly implemented?"""
from nose.tools import assert_equal, assert_not_equal
from numpy.testing import assert_allclose
from plotchecker import get_data
data = np.load("data/xy_data.npy")
# check that it uses the mse function
old_mse = mse
del mse
try:
fig, axis = plt.subplots()
plot_mse(axis, 9, data[:10], data[10:])
except NameError:
pass
else:
raise AssertionError("plot_mse should call mse, but it does not")
finally:
plt.close('all')
mse = old_mse
del old_mse
fig, axis = plt.subplots()
error = plot_mse(axis, 9, data[:10], data[10:])
axis.legend(loc='upper left')
# check the error
assert_equal(error.shape, (9, 2))
assert_allclose(error[0], mse(0, data[:10], data[10:]))
assert_allclose(error[4], mse(4, data[:10], data[10:]))
assert_allclose(error[8], mse(8, data[:10], data[10:]))
# check the plotted data
plotted_data = get_data(axis)
assert_equal(plotted_data.shape, (18, 2))
assert_allclose(plotted_data[:9, 0], np.arange(9))
assert_allclose(plotted_data[9:, 0], np.arange(9))
assert_allclose(plotted_data[:9, 1], error[:, 0])
assert_allclose(plotted_data[9:, 1], error[:, 1])
# check the line colors
assert axis.lines[0].get_color() in ['k', 'black', (0, 0, 0), '#000000']
assert axis.lines[1].get_color() in ['r', 'red', (1, 0, 0), '#FF0000']
# check the legend
legend_labels = [x.get_text() for x in axis.get_legend().get_texts()]
assert_equal(legend_labels, ["Training set error", "Testing set error"])
# check the axis labels
assert_not_equal(axis.get_xlabel(), "")
assert_not_equal(axis.get_ylabel(), "")
plt.close('all')
print("Success!")

```
```

`plot_mse`

function with different subsets of the data, depending on the index that is set:

```
In [15]:
```# load the data
data = np.load("data/xy_data.npy")
@interact
def visualize_mse(training_set_index=(1, 11)):
# relabel the index for convenience
i = training_set_index
# pull out the training and testing data
traindata = data[((i - 1) * 10):(i * 10)]
testdata = np.concatenate([data[:((i - 1) * 10)], data[(i * 10):]])
# plot the MSE
fig, axis = plt.subplots()
plot_mse(axis, 10, traindata, testdata)
axis.set_ylim(0, 0.01)
axis.set_title("MSE for dataset #{}".format(i))
axis.legend(loc='upper left')

```
```

0.25 points were awarded if you discussed how MSE on `traindata`

changes with the order of the polynomial model. 0.25 points were awarded if you also did so for `testdata`

. 0.25 points were awarded if you discussed whether or not you should assess overall model performance on the MSE calculated on `traindata`

. 0.25 points were awarded if you discussed why or why not this is a good idea. A sample answer is included below:

The MSE calculated on

`traindata`

continues to decrease with the polynomial model order. When calculated on`testdata`

, the MSE decreases until the model order is around 2, then increases for model orders greater than 2. Assessing your model using MSE calculated on the training data is a bad idea. This is because as the complexity of a model increases, it will be better and better able to fit the trends in the training data. Eventually the MSE will go to zero. However, what we really care about is generalization. In other words we want to know how accurately a model will predict for data it has never seen before. This information is provided by instead assessing the MSE on the testing data.

The MSE on the test set is an approximation of the *prediction error* associated with a particular learning algorithm. In class we discussed how *expected prediction error* is composed of a bias and a variance term. Both the bias and the variance depend upon the true function, $f$, that we are estimating with our learning algorithm.

- An algorithm with high
*bias*will systematically produce predictions that differ from the true function. This can happen in situations where the true function is more complex than the most complex function available to the algorithm. - An algorithm with high
*variance*will produce predictions that can vary wildly depending on the specifics of the dataset we are evaluating. This can happen in situations where the algorithm is able to return functions that are more complex than the true function.

For a more detailed explanation of the bias-variance tradeoff, please take a look at the reading by Geman et al. (1992).