In this project, you will implement bagging for regression trees.
Before you get started, let's import a few packages that you will need.
import numpy as np
from pylab import *
from numpy.matlib import repmat
import matplotlib
import matplotlib.pyplot as plt
from scipy.io import loadmat
import time
import sys
sys.path.append('/home/codio/workspace/.modules')
from helper import *
print('You\'re running python %s' % sys.version.split(' ')[0])You're running python 3.6.8
In this project, we will work with an artificial 2D spiral dataset of size 150 (for visualization), and a high dimensional dataset ION (for a binary test classification problem).
xTrSpiral, yTrSpiral, xTeSpiral, yTeSpiral = spiraldata(150)
xTrIon, yTrIon, xTeIon, yTeIon = iondata()You will use the regression tree from a previous project. As a reminder, the following code shows you how to instantiate a decision tree:
# Create a regression tree with no restriction on its depth
# and weights for each training example to be 1.
# If you want to create a tree of max_depth k,
# then call RegressionTree(depth = k).
tree = RegressionTree(depth = np.inf)
# To fit/train the regression tree
tree.fit(xTrSpiral, yTrSpiral)
# To use the trained regression tree to predict a score for the example
score = tree.predict(xTrSpiral)
# To use the trained regression tree to make a +1/-1 prediction
pred = np.sign(tree.predict(xTrSpiral))
tr_err = np.mean((np.sign(tree.predict(xTrSpiral)) - yTrSpiral)**2)
te_err = np.mean((np.sign(tree.predict(xTeSpiral)) - yTeSpiral)**2)
print("Training error: %.4f" % np.mean(np.sign(tree.predict(xTrSpiral)) != yTrSpiral))
print("Testing error: %.4f" % np.mean(np.sign(tree.predict(xTeSpiral)) != yTeSpiral))Training error: 0.0000
Testing error: 0.0400
The following code defines a function visclassifier(), which plots the decision boundary of a classifier in 2 dimensions. Execute the following code to see what the decision boundary of your tree looks like on the spiral data set.
def visclassifier(fun, xTr, yTr):
"""
visualize decision boundary
Define the symbols and colors we'll use in the plots later
"""
yTr = np.array(yTr).flatten()
symbols = ["ko", "kx"]
marker_symbols = ['o', 'x']
mycolors = [[0.5, 0.5, 1], [1, 0.5, 0.5]]
# Get the unique values from labels array
classvals = np.unique(yTr)
plt.figure()
# Return 300 evenly spaced numbers over this interval
res = 300
xrange = np.linspace(min(xTr[:, 0]), max(xTr[:, 0]), res)
yrange = np.linspace(min(xTr[:, 1]), max(xTr[:, 1]), res)
# Repeat this matrix 300 times for both axes
pixelX = repmat(xrange, res, 1)
pixelY = repmat(yrange, res, 1).T
xTe = np.array([pixelX.flatten(), pixelY.flatten()]).T
# Test all of these points on the grid
testpreds = fun(xTe)
# reshape it back together to make our grid
Z = testpreds.reshape(res, res)
# Z[0,0] = 1 # optional: scale the colors correctly
# Fill in the contours for these predictions
plt.contourf(pixelX, pixelY, np.sign(Z), colors = mycolors)
# Create x's and o's for training set
for idx, c in enumerate(classvals):
plt.scatter(xTr[yTr == c, 0],
xTr[yTr == c, 1],
marker = marker_symbols[idx],
color = 'k')
plt.axis('tight')
plt.show() # shows figure and blocks
# Compute tree on training data
tree = RegressionTree(depth = np.inf)
tree.fit(xTrSpiral, yTrSpiral)
print("Training error: %.4f" % np.mean(np.sign(tree.predict(xTrSpiral)) != yTrSpiral))
print("Testing error: %.4f" % np.mean(np.sign(tree.predict(xTeSpiral)) != yTeSpiral))
plt_is_inline = 'inline' in matplotlib.get_backend()
if plt_is_inline:
%matplotlib notebook
visclassifier(lambda X : tree.predict(X), xTrSpiral, yTrSpiral)Training error: 0.0000
Testing error: 0.0400
<IPython.core.display.Javascript object>
CART trees are known to be high variance classifiers if trained to full depth. Equivalently, CART trees of full depth easily overfit to the training set, which prevent the trees from performing well on new unseen data. An effective way to prevent overfitting is to use Bagging (short for bootstrap aggregating). Implement the function forest, which builds a forest of m regression trees of depth = maxdepth. Each tree should be built using training data drawn by randomly sampling n examples from the training data with replacement, where n is the number of points in xTr.
We are going to keep it simple and not randomly subsample a small set of features to split on. Therefore, all trees will be constructed with a simple call to RegressionTree(depth = maxdepth). The function should output the list of trees.
Hint: You may find np.random.choice(a, b, replace = True) useful. It samples b numbers from [0, ..., a-1] with replacement.
def forest(xTr, yTr, m, maxdepth = np.inf):
"""
Creates a forest of m trees, each of depth=maxdepth.
Input:
xTr: n x d matrix of data points
yTr: n-dimensional vector of labels
m: number of trees in the forest
maxdepth: maximum depth of each tree
Output:
trees: list of decision trees of length m
"""
n, d = xTr.shape
trees = []
# YOUR CODE HERE
for x in range(m):
indices = np.random.choice(np.arange(n),n, replace = True)
xTrain = xTr[indices]
yTrain = yTr[indices]
tree = RegressionTree(depth=maxdepth)
tree.fit(xTrain,yTrain)
trees.append(tree)
return trees
raise NotImplementedError()def forest_test1():
m = 20
x = np.arange(100).reshape((100, 1))
y = np.arange(100)
trees = forest(x, y, m)
# Make sure there are m trees in the forest
return len(trees) == m
def forest_test2():
m = 20
x = np.arange(100).reshape((100, 1))
y = np.arange(100)
max_depth = 4
trees = forest(x, y, m, max_depth)
# Get the depth of all trees in the forest
depths_forest = np.array([tree.depth for tree in trees])
# Make sure that the max depth of all the trees is correct
return np.all(depths_forest == max_depth)
def forest_test3():
s = set()
def DFScollect(tree):
# Do Depth first search to all prediction in the tree
if tree.left is None and tree.right is None:
s.add(tree.prediction)
else:
DFScollect(tree.right)
DFScollect(tree.left)
m = 200
x = np.arange(100).reshape((100, 1))
y = np.arange(100)
trees = forest(x, y, m);
lens = np.zeros(m)
for i in range(m):
s.clear()
DFScollect(trees[i].root)
lens[i] = len(s)
# Check that about 63% of data is represented in each random sample
return abs(np.mean(lens) - 100 * (1 - 1 / np.exp(1))) < 2
runtest(forest_test1, 'forest_test1')
runtest(forest_test2, 'forest_test2')
runtest(forest_test3, 'forest_test3')Running Test: forest_test1 ... ✔ Passed!
Running Test: forest_test2 ... ✔ Passed!
Running Test: forest_test3 ... ✔ Passed!
# Autograder test cell worth 1 point;
# runs forest_test1# Autograder test cell worth 1 point;
# runs forest_test2# Autograder test cell worth 1 point;
# runs forest_test3Now implement the function evalforest, which should take as input a set of
Note that for bagging, we take the average over all trees weighted equally. In a later project, you will implement a different version of this function that assigns different weights to predictions of different trees.
def evalforest(trees, X):
"""
Evaluates X using trees.
Input:
trees: list of length m of RegressionTree decision trees
X: n x d matrix of data points
Output:
pred: n-dimensional vector of predictions
"""
m = len(trees)
n,d = X.shape
pred = np.zeros(n)
# YOUR CODE HERE
for tree in trees:
prediction = tree.predict(X)
pred+=prediction
pred /= m
return pred
raise NotImplementedError()def evalforest_test1():
m = 200
x = np.arange(100).reshape((100, 1))
y = np.arange(100)
trees = forest(x, y, m)
preds = evalforest(trees, x)
return preds.shape == y.shape
def evalforest_test2():
m = 200
x = np.ones(10).reshape((10, 1))
y = np.ones(10)
max_depth = 0
# Create a forest with trees depth 0.
# Since the data are all ones, each tree will return 1 as prediction.
trees = forest(x, y, m, 0)
pred = evalforest(trees, np.ones((1, 1)))[0]
return np.isclose(pred, 1) # the prediction should be equal to the sum of weights
def bagging_test1():
m = 50
xTr = np.random.rand(500, 3) - 0.5
yTr = np.sign(xTr[:, 0] * xTr[:, 1] * xTr[:, 2]) # XOR Classification
xTe = np.random.rand(50, 3) - 0.5
yTe = np.sign(xTe[:, 0] * xTe[:, 1] * xTe[:, 2])
tree = RegressionTree(depth = 5)
tree.fit(xTr, yTr)
oneacc = np.sum(np.sign(tree.predict(xTe)) == yTe)
trees = forest(xTr, yTr, m, maxdepth = 5)
multiacc = np.sum(np.sign(evalforest(trees, xTe)) == yTe)
# Check that bagging yields improvement - or doesn't get too much worse
return multiacc * 1.1 >= oneacc
def bagging_test2():
m = 50
xTr = (np.random.rand(500, 3) - 0.5) * 4
yTr = xTr[:,0] * xTr[:, 1] * xTr[:, 2] # XOR Regression
xTe = (np.random.rand(50, 3) - 0.5) * 4
yTe = xTe[:, 0] * xTe[:, 1] * xTe[:, 2]
np.random.seed(1)
tree = RegressionTree(depth = 3)
tree.fit(xTr, yTr)
oneerr = np.sum(np.sqrt((tree.predict(xTe) - yTe) ** 2))
trees = forest(xTr, yTr, m, maxdepth = 3)
multierr = np.sum(np.sqrt((evalforest(trees, xTe) - yTe) ** 2))
# Check that bagging yields improvement - or doesn't get too much worse
return multierr <= oneerr * 1.5
runtest(evalforest_test1, 'evalforest_test1')
runtest(evalforest_test2, 'evalforest_test2')
runtest(bagging_test1, 'bagging_test1')
runtest(bagging_test2, 'bagging_test2')Running Test: evalforest_test1 ... ✔ Passed!
Running Test: evalforest_test2 ... ✔ Passed!
Running Test: bagging_test1 ... ✔ Passed!
Running Test: bagging_test2 ... ✔ Passed!
# Autograder test cell worth 1 point;
# runs evalforest-test1# Autograder test cell worth 1 point;
# runs evalforest-test2# Autograder test cell worth 1 point;
# runs bagging-test1# Autograder test cell worth 1 point;
# runs bagging-test2The following script visualizes the decision boundary of an ensemble of decision tree. You might observe that the decision boundary is less rigid with an ensemble of 50 trees than with just 1 CART tree. This is to be expected as a forest is just an ensemble of CART trees and averages the predictions. Consequently, the test error should also be less with the the ensemble than with just 1 CART tree.
# Compute tree on training data
trees = forest(xTrSpiral, yTrSpiral, 50)
print("Training error: %.4f" % np.mean(np.sign(evalforest(trees, xTrSpiral)) != yTrSpiral))
print("Testing error: %.4f" % np.mean(np.sign(evalforest(trees, xTeSpiral)) != yTeSpiral))
plt_is_inline = 'inline' in matplotlib.get_backend()
if plt_is_inline:
%matplotlib notebook
# Plot results
visclassifier(lambda X : evalforest(trees, X), xTrSpiral, yTrSpiral)Training error: 0.0000 Testing error: 0.0200
The following script evaluates the test and training error of an ensemble of decision trees as we vary the number of trees.
M = 20 # max number of trees
err_trB = []
err_teB = []
alltrees = forest(xTrIon, yTrIon , M)
for i in range(M):
trees =alltrees[:i+1]
trErr = np.mean(np.sign(evalforest(trees,xTrIon)) != yTrIon)
teErr = np.mean(np.sign(evalforest(trees,xTeIon)) != yTeIon)
err_trB.append(trErr)
err_teB.append(teErr)
print("[%d]training err = %.4f\ttesting err = %.4f" % (i, trErr, teErr))
plt_is_inline = 'inline' in matplotlib.get_backend()
if plt_is_inline:
%matplotlib notebook
# Plot results
plt.figure()
line_tr, = plt.plot(range(M), err_trB, '-*', label = "Training error")
line_te, = plt.plot(range(M), err_teB, '-*', label = "Testing error")
plt.title("Ensemble of Decision Trees")
plt.legend(handles = [line_tr, line_te])
plt.xlabel("number of trees")
plt.ylabel("error")
plt.show()The next interactive demo shows a 1-dimensional curve fitted with ensembles of decision trees. We sample 100 training data points from the curve with additive noise. Note how the predicted curve becomes increasingly smooth as you add trees — adding trees should thus decrease training and testing error. The testing error may be a little higher than the training error, but there is no large gap as ensembles are quite good at avoiding overfitting.
def onclick_forest(event):
"""
Visualize forest by adding more trees
"""
global xTrain, yTrain, Q, trees
if event.key == 'shift':
Q += 10
else:
Q += 1
Q = min(Q, M)
plt.cla()
plt.xlim((0, 1))
plt.ylim((0, 1))
pTest = evalforest(trees[:Q], xTest)
pTrain = evalforest(trees[:Q], xTrain)
errTrain = np.sqrt(np.mean((pTrain - yTrain)**2))
errTest = np.sqrt(np.mean((pTest - yTest)**2))
plt.plot(xTrain[:, 0], yTrain, 'bx')
plt.plot(xTest[:, 0], pTest, 'r-')
plt.plot(xTest[:, 0], yTest, 'k-')
plt.legend(['Training data','Prediction'])
plt.title('(%i Trees) Error Tr: %2.4f, Te:%2.4f' % (Q, errTrain, errTest))
plt.show()
n = 100 # number of training points
NOISE = 0.05 # magnitude of noise
xTrain = np.array([np.linspace(0, 1, n), np.zeros(n)]).T
yTrain = 2 * np.sin(xTrain[:, 0] * 3) * (xTrain[:, 0]**2)
yTrain += np.random.randn(yTrain.size) * NOISE;
ntest = 300 # density of test points
xTest = np.array([linspace(0, 1, ntest), np.zeros(ntest)]).T
yTest = 2 * np.sin(xTest[:, 0] * 3) * (xTest[:, 0]**2)
# Hyper-parameters (feel free to play with them)
M = 100 # number of trees
depth = np.inf
trees = forest(xTrain, yTrain, M, maxdepth = depth)
Q = 0
plt_is_inline = 'inline' in matplotlib.get_backend()
if plt_is_inline:
%matplotlib notebook
# Plot interactive demo
print('Note: Rerun this cell each time you want to interact with this plot.')
print('You can either click to add one tree or shift+click to add 10 trees.')
fig = plt.figure()
cid = fig.canvas.mpl_connect('button_press_event', onclick_forest)
plt.title('Start boosting on this 1D data (click or shift+click)')
plt.plot(xTrain[:, 0], yTrain, '*')
plt.plot(xTest[:, 0], yTest, 'k-')
plt.xlim(0, 1)
plt.ylim(0, 1);Note: Rerun this cell each time you want to interact with this plot. You can either click to add one tree or shift+click to add 10 trees.
The following demo visualizes the bagged classifier. Add your own points directly on the graph with click and shift+click to see the prediction boundaries. There will be a delay between clicks as the new classifier is trained.
def onclick_classifier(event):
"""
Visualize bagged classifier by adding new points
"""
global xTrain, yTrain, w, b, M
# Create position vector for new point
pos = np.array([[event.xdata, event.ydata]])
if event.key == 'shift': # add positive point
color = 'or'
label = 1
else: # add negative point
color = 'ob'
label = -1
xTrain = np.concatenate((xTrain, pos), axis = 0)
yTrain = np.append(yTrain, label)
marker_symbols = ['o', 'x']
classvals = np.unique(yTrain)
w = np.array(w).flatten()
mycolors = [[0.5, 0.5, 1], [1, 0.5, 0.5]]
# Return 300 evenly spaced numbers over this interval
res = 300
xrange = np.linspace(0, 1, res)
yrange = np.linspace(0, 1, res)
# Repeat this matrix 300 times for both axes
pixelX = repmat(xrange, res, 1)
pixelY = repmat(yrange, res, 1).T
xTe = np.array([pixelX.flatten(), pixelY.flatten()]).T
# Get forest
trees = forest(xTrain, yTrain, M)
fun = lambda X:evalforest(trees, X)
# Test all of these points on the grid
testpreds = fun(xTe)
# Reshape it back together to make our grid
Z = testpreds.reshape(res, res)
# Z[0,0] = 1 # optional: scale the colors correctly
plt.cla()
plt.xlim((0, 1))
plt.ylim((0, 1))
# Fill in the contours for these predictions
plt.contourf(pixelX, pixelY, np.sign(Z), colors = mycolors)
for idx, c in enumerate(classvals):
plt.scatter(xTrain[yTrain == c, 0],
xTrain[yTrain == c, 1],
marker = marker_symbols[idx],
color = 'k')
plt.show()
xTrain = np.array([[5, 6]])
b = yTrIon
yTrain = np.array([1])
w = xTrIon
M = 20
plt_is_inline = 'inline' in matplotlib.get_backend()
if plt_is_inline:
%matplotlib notebook
# Plot interactive demo
print('Note: Rerun this cell each time you want to interact with this plot.')
print('Click to add a negative point or shift+click to add a positive point.')
print('You may notice a delay when adding points to the visualization.')
fig = plt.figure()
plt.xlim(0, 1)
plt.ylim(0, 1)
cid = fig.canvas.mpl_connect('button_press_event', onclick_classifier)
plt.title('Start by clicking the plot to add a negative point');Note: Rerun this cell each time you want to interact with this plot. Click to add a negative point or shift+click to add a positive point. You may notice a delay when adding points to the visualization.