# Practical Session 3: Ensemble Learning Techniques

*Notebook by Ekaterina Kochmar*

This practical will address the use of ensemble-based learning techniques. You will be working with the otherwise familiar settings of classification and regression tasks. In this practical, you will use the familiar datasets and will also learn how to generate artificial datasets.

## Learning objectives

In this practical you will learn about:
- simple voting classifiers using hard and soft voting strategies
- bagging and pasting ensembles
- boosting and early stopping
- popular ensemble-based learning algorithms including `RandomForests` and `AdaBoost`

In addition, you will get more practice with `sklearn's` machine learning routines, and you will learn how to implement ensembles both from scratch and using `sklearn's` implementation.

**Bibliography**: Aurelien Geron, *Hands-On Machine Learning with Scikit-Learn and TensorFlow*.

## Setup

First of all, let's import a few common modules, ensure `matplotlib` plots are inline and define the parameters for the figures. Feel free to add your own settings for the notebook here:


```python
# Common imports
import numpy as np
import os

# to make this notebook's output stable across runs
np.random.seed(42)

# To plot pretty figures
%matplotlib inline
import matplotlib
import matplotlib.pyplot as plt
plt.rcParams['axes.labelsize'] = 14
plt.rcParams['xtick.labelsize'] = 12
plt.rcParams['ytick.labelsize'] = 12
```

## Simple voting classifiers: hard and soft voting strategies

[The wisdom of the crowd](https://en.wikipedia.org/wiki/Wisdom_of_the_crowd) is the collective opinion of a group of individuals rather than that of a single expert. It has been shown that if you ask a complex question to thousands of random people and then aggregate their answers, in many cases you are likely to get an answer that is more accurate than an expert's prediction. This effect is nowadays widely used in practice and can be said to provide the basis for crowdsourcing approaches.

By analogy, if you combine decisions of multiple diverse classifiers, the resulting vote will often be better than a vote of any single predictor. Such combination of classifiers is called _ensemble_. Perhaps somewhat counterintuitively, the aggregated vote of an ensemble often achieves a higher accuracy than the best classifier in the ensemble. This works even if each predictor in the ensemble is a _weak learner_ (meaning that it works only slightly better than random guessing). The ensemble can still be a _strong learner_ achieving high accuracy, provided that there are:
- a sufficient number of weak learners, and
- their decisions are sufficiently diverse.

How is this possible? Here is a statistical analogy: suppose you have a slightly biased coin that has a $51\%$ chance of coming up heads and a $49\%$ chance of coming up tails. If you toss it $1000$ times, you will get around $510$ heads and $490$ tails, and hence a majority of heads. If you do the math or run simulations, you will find out that the probability of obtaining the majority of heads after $1000$ tosses of this coin approaches $73\%$. After $10K$ tosses, the probability of having a majority of heads is about $97\%$. This is explained by the [_law of large numbers_](https://en.wikipedia.org/wiki/Law_of_large_numbers): over a large number of tosses the ratio of heads gets closer to the probability of heads ($51\%$).


```python
import scipy, math
import scipy.special

def probability(p, n, x):
    binom = scipy.special.comb(n, x, exact=True) * math.pow(p, x) * math.pow((1-p), n-x)
    return binom

prob_maj = 0.0
for x in range(501, 1001):
    prob_maj += probability(0.51, 1000, x)
    
print(prob_maj)
```

The code below simulates $10$ independent series of coin tosses, and then plots the results. You can see in the figure below that as the number of tosses increases, all $10$ series end up close to $51\%$ and thus are consistently above $50\%$.


```python
np.random.seed(42)
heads_proba = 0.51
coin_tosses = (np.random.rand(10000, 10) < heads_proba).astype(np.int32)
cumulative_heads_ratio = np.cumsum(coin_tosses, axis=0) / np.arange(1, 10001).reshape(-1, 1)

plt.figure(figsize=(8,3.5))
plt.plot(cumulative_heads_ratio)
plt.plot([0, 10000], [0.51, 0.51], "k--", linewidth=2, label="51%")
plt.plot([0, 10000], [0.5, 0.5], "k-", label="50%")
plt.xlabel("Number of coin tosses")
plt.ylabel("Heads ratio")
plt.legend(loc="lower right")
plt.axis([0, 10000, 0.42, 0.58])
plt.show()
```

Similarly, if you combine $1000$ classifiers, each of which is individually correct only $51\%$ of the time, and then predict the majority voted class, you can get up to $73\%$ accuracy.

For this to work, the classifiers need to be independent and sufficiently diverse, i.e. they need to make independent, uncorrelated errors. In practice, this is hard to achieve because the classifiers are usually all trained on the same data.

Let's apply voting strategy to some data: this time, we'll be working with a synthetic (i.e., artificially generated) [`moons` dataset](https://scikit-learn.org/stable/modules/generated/sklearn.datasets.make_moons.html), that includes two interleaving half circles and provides a good 'toy' example for testing out classification strategies.

Let's first generate the data points and visualise them on a plot:


```python
from sklearn.datasets import make_moons
X, y = make_moons(n_samples=500, noise=0.30, random_state=42)

def plot_dataset(X, y, axes):
    plt.plot(X[:, 0][y==0], X[:, 1][y==0], "bs")
    plt.plot(X[:, 0][y==1], X[:, 1][y==1], "g^")
    plt.axis(axes)
    plt.grid(True, which='both')
    plt.xlabel(r"$x_1$", fontsize=20)
    plt.ylabel(r"$x_2$", fontsize=20, rotation=0)

plot_dataset(X, y, [-1.5, 2.5, -1, 1.5])
plt.show()
```

Next, let's split the dataset into training and test sets and train some classifiers on the training data:


```python
from sklearn.model_selection import train_test_split
X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=42)
```

The code below shows how to apply several classifiers to this data, with a couple of classifiers that we haven't used before in this course: [`RandomForestClassifier`](https://scikit-learn.org/stable/modules/generated/sklearn.ensemble.RandomForestClassifier.html) and [`Support Vector Machines`](https://scikit-learn.org/stable/modules/generated/sklearn.svm.SVC.html) (`SVC`) classifier. For the moment, don't worry about the particular details of the classifiers' implementation (if you are interested in learning more about them, check out `sklearn's` documentation; we'll discuss `RandomForestClassifier` in this practical; `Support Vector Machines` is covered in detail in the Part II [Machine Learning and Bayesian Inference](https://www.cl.cam.ac.uk/teaching/2021/MLBayInfer/) course). The point here is that we are combining *diverse* classifiers in an ensemble: in essence, `RandomForestClassifier` runs recursive partitioning on the data, while `Support Vector Machines` try to construct a linear boundary between the classes (as discussed in Lecture 3). Feel free to experiment with your own selection of classifiers from the [`sklearn's`](https://scikit-learn.org/stable/supervised_learning.html#supervised-learning) suite.


```python
from sklearn.ensemble import VotingClassifier

from sklearn.ensemble import RandomForestClassifier
from sklearn.linear_model import LogisticRegression
from sklearn.svm import SVC

log_clf = LogisticRegression(random_state=42)
rnd_clf = RandomForestClassifier(random_state=42)
svm_clf = SVC(random_state=42)

voting_clf = VotingClassifier(
    estimators=[('lr', log_clf), ('rf', rnd_clf), ('svc', svm_clf)],
    voting='hard')
voting_clf.fit(X_train, y_train)
```

Let's look into how correlated the predictions of the classifiers are: first, get the predictions on the test data, next store them as a `pandas` DataFrame, and finally apply `.corr()` function:


```python
import pandas as pd

def get_predictions(clf):
    clf.fit(X_train, y_train)
    return clf.predict(X_test)


preds = {'lr': get_predictions(log_clf), 
        'rf': get_predictions(rnd_clf), 
        'svc': get_predictions(svm_clf)}
df = pd.DataFrame(data=preds)
df[:100]
```


```python
df.corr()
```

Will this correlation in the individual classifiers' predictions be sufficient for an ensemble? Let's check this out by combining the votes with a *hard voting strategy*: the ensemble classifier will simply choose the majority class predicted by the three classifiers.

The code below prints out individual classifiers' accuracy scores along with the ensemble's accuracy score:


```python
from sklearn.metrics import accuracy_score

for clf in (log_clf, rnd_clf, svm_clf, voting_clf):
    clf.fit(X_train, y_train)
    y_pred = clf.predict(X_test)
    print(clf.__class__.__name__, accuracy_score(y_test, y_pred))
```

As you can see, the voting classifier slightly outperforms each of the three individual classifiers, including the most accurate one(s). In addition, if individual classifiers can estimate class probabilities (i.e., if they have `predict_proba()` method), you can use _soft voting_ strategy, that is, estimate highest class probability averaged over individual classifiers. In comparison to the hard voting strategy, soft voting gives more weight to highly confident votes.

`SVC` classifier doesn't estimate class probabilities by default, so you need to set the `probability` hyperparameter to `True`.


```python
log_clf = LogisticRegression(random_state=42)
rnd_clf = RandomForestClassifier(random_state=42)
svm_clf = SVC(probability=True, random_state=42)

voting_clf = VotingClassifier(
    estimators=[('lr', log_clf), ('rf', rnd_clf), ('svc', svm_clf)],
    voting='soft')
voting_clf.fit(X_train, y_train)
```

Not let's estimate the accuracy of the voting classifier in this mode:


```python
from sklearn.metrics import accuracy_score

for clf in (log_clf, rnd_clf, svm_clf, voting_clf):
    clf.fit(X_train, y_train)
    y_pred = clf.predict(X_test)
    print(clf.__class__.__name__, accuracy_score(y_test, y_pred))
```

## Bagging ensembles: bagging and pasting

One way to ensure that the classifiers are sufficiently diverse in their predicitons is to select different training algorithms. As we said before, often it is hard to ensure that the classifiers' errors are not highly correlated because the classifiers are ultimately trained on the same data.

So, another way to make sure you end up with diverse classifiers is to train them on different random subsets of the training set. There are two ways to do that:
- sampling _with replacement_ is called [_bagging_](https://link.springer.com/article/10.1023/A:1018054314350) – short for _bootstrap aggregating_
- sampling _without replacement_ is called [_pasting_](https://link.springer.com/article/10.1023/A:1007563306331)

Both bagging and pasting allow training instances to be sampled several times across multiple predictors, but only bagging allows training instances to be sampled several times for the same predictor.

When applied to a new, test instance the ensemble classifier aggregates the predictions of all predictors and estimates the statistical mode (i.e., most frequent prediction).

**Question**: Is this similar to the hard or the soft voting strategy?

Let's use `sklearn's BaggingClassifier` to train an ensemble of $500$ [`Decision Trees`](https://scikit-learn.org/stable/modules/generated/sklearn.tree.DecisionTreeClassifier.html), each trained on $100$ training instances randomly selected from the training set with replacement (for bagging, but you can change the strategy to pasting by setting `bootstrap` to `False`). In addition, both bagging and pasting strategies are good for parallelisation: predictors can all be trained in parallel on different CPU cores or even different servers. Setting `n_jobs` parameter to $-1$ tells `sklearn` to use all available CPU cores:


```python
from sklearn.ensemble import BaggingClassifier
from sklearn.tree import DecisionTreeClassifier

#use bootstrap=False for pasting
bag_clf = BaggingClassifier(
    DecisionTreeClassifier(random_state=42), n_estimators=500,
    max_samples=100, bootstrap=True, n_jobs=-1, random_state=42)
#n_jobs = use all of the available CPU cores
bag_clf.fit(X_train, y_train)
y_pred = bag_clf.predict(X_test)
```

**Question**: [`BaggingClassifier`](https://scikit-learn.org/stable/modules/generated/sklearn.ensemble.BaggingClassifier.html) automatically performs soft voting when the base classifier can estimate class probabilities (i.e., has a `predict_proba()` method). Which strategy – hard or soft voting – does `BaggingClassifier` follow with the `DecisionTreeClassifier` as the base classifier?

Let's calculate the accuracy score for the bagging classifier that combines $500$ Decision Tree estimators with a prediction of a single Decision Tree trained on the same data:


```python
print(accuracy_score(y_test, y_pred))
```


```python
tree_clf = DecisionTreeClassifier(random_state=42)
tree_clf.fit(X_train, y_train)
y_pred_tree = tree_clf.predict(X_test)
print(accuracy_score(y_test, y_pred_tree))
```

To get more insight into the results, let's also plot the decision boundary for the two classifiers:


```python
from matplotlib.colors import ListedColormap

def plot_decision_boundary(clf, X, y, axes=[-1.5, 2.5, -1, 1.5], alpha=0.5, contour=True):
    x1s = np.linspace(axes[0], axes[1], 100)
    x2s = np.linspace(axes[2], axes[3], 100)
    x1, x2 = np.meshgrid(x1s, x2s)
    X_new = np.c_[x1.ravel(), x2.ravel()]
    y_pred = clf.predict(X_new).reshape(x1.shape)
    custom_cmap = ListedColormap(['#fafab0','#9898ff','#a0faa0'])
    plt.contourf(x1, x2, y_pred, alpha=0.3, cmap=custom_cmap)
    if contour:
        custom_cmap2 = ListedColormap(['#7d7d58','#4c4c7f','#507d50'])
        plt.contour(x1, x2, y_pred, cmap=custom_cmap2, alpha=0.8)
    plt.plot(X[:, 0][y==0], X[:, 1][y==0], "yo", alpha=alpha)
    plt.plot(X[:, 0][y==1], X[:, 1][y==1], "bs", alpha=alpha)
    plt.axis(axes)
    plt.xlabel(r"$x_1$", fontsize=18)
    plt.ylabel(r"$x_2$", fontsize=18, rotation=0)
    
plt.figure(figsize=(11,4))
plt.subplot(121)
plot_decision_boundary(tree_clf, X, y)
plt.title("Decision Tree", fontsize=14)
plt.subplot(122)
plot_decision_boundary(bag_clf, X, y)
plt.title("Decision Trees with Bagging", fontsize=14)
plt.show()
```

Even though the ensemble makes a comparable number of errors on the training set as a single Decision Tree, its decision boundary is less irregular. This suggests that ensemble's predictions will likely generalise better than the single `Decision Tree`'s predictions
when applied to the new data.

### Out-of-Bag evaluation

With bagging, some of the instances get sampled several times for any given predictor, while others may never be selected at all. By default, `BaggingClassifier` samples $m$ training instances with replacement, where $m$ is the size of the training set. As $m$ grows, the ratio of the training instances that are sampled on average for each predictor approaches $1 - exp(-1)$, i.e. around $63\%$. The remaining $37\%$ of the training instances that are never used by the predictors are called *out-of-bag* (*oob*) instances.

Mathematical proof of the above: With $m$ training instances in the set, the probability of not picking an instance in a random draw is $\frac{m-1}{m}$. For sampling with replacement, the probability of not picking an instance in $m$ random draws is $(\frac{m-1}{m})^m$. In the limit of a large number of $m$, this becomes:

\begin{equation}
\lim_{m\to\infty} (1 - \frac{1}{m})^m = e^{-1} = 0.368 \approx 37\%
\end{equation}

Since a predictor never sees these instances during training, they can be used for evaluation without the need for an additional validation set or cross-validation experiments. You can evaluate the ensemble itself by averaging each predictor's performance on the oob instances.

The code below shows you how to do that: if you set `oob_score` to `True`, the bagging classifier will be evaluated on the oob instances after training, and you can output the score as `oob_score_`.


```python
bag_clf = BaggingClassifier(
    DecisionTreeClassifier(random_state=42), n_estimators=500,
    bootstrap=True, n_jobs=-1, oob_score=True, random_state=40)
bag_clf.fit(X_train, y_train)
bag_clf.oob_score_
```

This means that the bagging classifier is likely to achieve around $90\%$ accuracy on the test data, too. Let's check the results:


```python
from sklearn.metrics import accuracy_score
y_pred = bag_clf.predict(X_test)
accuracy_score(y_test, y_pred)
```

Close enough! You can also take a look into the decision process by printing out the predictions of the classifier on the training instances. Since the base estimator has a `predict_proba()` method, what you see here are the class probabilities assigned by the classifier to negative and positive classes for each instance:


```python
bag_clf.oob_decision_function_
```

## Random Forests

[`RandomForestClassifier`](https://scikit-learn.org/stable/modules/generated/sklearn.ensemble.RandomForestClassifier.html) is an ensemble of Decision Trees typically trained via the bagging method. The number of training instances (`max_samples`) is usually set to the total size of the training set. So, in fact, `RandomForestClassifier` will be roughly equivalent to the `BaggingClassifier` that takes Decision Trees as base estimators with the following parameters:


```python
bag_clf = BaggingClassifier(
    DecisionTreeClassifier(splitter="random", max_leaf_nodes=16, random_state=42),
    n_estimators=500, max_samples=1.0, bootstrap=True, n_jobs=-1, random_state=42)

bag_clf.fit(X_train, y_train)
y_pred = bag_clf.predict(X_test)
```

The key difference is that, rather than using `BaggingClassifier` and passing it a `DecisionTreeClassifier`, you can rely on the `sklearn's RandomForestClassifier` implementation, which is more convenient and optimised for Decision Trees:


```python
from sklearn.ensemble import RandomForestClassifier

rnd_clf = RandomForestClassifier(n_estimators=500, max_leaf_nodes=16, 
                                 n_jobs=-1, random_state=42)
rnd_clf.fit(X_train, y_train)

y_pred_rf = rnd_clf.predict(X_test)
```

Let's estimate the difference between the two classifiers' predictions:


```python
np.sum(y_pred == y_pred_rf) / len(y_pred)  # see to what extent predictions are identical
```

I.e., the two classifiers are almost identical in their predictions. The Random Forest algorithm introduces extra randomness when growing trees: instead of searching for the very best feature when splitting a node, it searches for the best feature among a random subset of features. As a result, the model is more diverse and in general yields better results.

Let's visualise the decision boundary for a random set of $15$ Decision Trees to get a better idea of where the diversity comes from:


```python
plt.figure(figsize=(6, 4))

for i in range(15):
    tree_clf = DecisionTreeClassifier(max_leaf_nodes=16, random_state=42 + i)
    indices_with_replacement = np.random.randint(0, len(X_train), len(X_train))
    tree_clf.fit(X[indices_with_replacement], y[indices_with_replacement])
    plot_decision_boundary(tree_clf, X, y, axes=[-1.5, 2.5, -1, 1.5], 
                           alpha=0.02, contour=False)

plt.show()
```

**Extra material**: One further way to introduce randomness into the training process is to not only select a random subset of features when splitting a node, but also use a random threshold for each feature rather than searching for the best possible thresholds as traditional Decision Trees do. Since estimating best possible threshold for each feature at every node is also one of the most time-consuming tasks, this modification to the algorithm makes it more efficient. The algorithm that exploits this idea is called *Extremely Randomised Trees ensemble*, or *Extra-Trees*, and you can use it via [`sklearn's ExtraTreeClassifier`](https://scikit-learn.org/stable/modules/generated/sklearn.tree.ExtraTreeClassifier.html) implementation.

### Feature importance

Another useful property of the Random Forests classifier is that it can help you measure the relative importance of each feature by looking at how much the tree nodes that use this particular feature reduce impurity of the nodes on average (i.e., across all trees in the forest). The results are scaled so that the sum of all feature importances is equal to $1$.

Because features used in the *iris* and *digits* datasets are highly interpretable, let's estimate feature importances on these two datasets: 


```python
from sklearn.datasets import load_iris
iris = load_iris()
rnd_clf = RandomForestClassifier(n_estimators=500, n_jobs=-1, random_state=42)
rnd_clf.fit(iris["data"], iris["target"])
for name, score in zip(iris["feature_names"], rnd_clf.feature_importances_):
    print(name, score)
```


```python
rnd_clf.feature_importances_
```

**Question**: What do these results suggest? Are they similar to your observations from the previous experiments on this dataset?

*Digits* dataset provides further interpretability of the results if we plot them. What does the figure below suggest?


```python
from sklearn import datasets
digits = datasets.load_digits()

rnd_clf = RandomForestClassifier(random_state=42)
rnd_clf.fit(digits["data"], digits["target"])
```


```python
def plot_digit(data):
    image = data.reshape(8, 8)
    plt.imshow(image, cmap = matplotlib.cm.hot,
               interpolation="nearest")
    plt.axis("off")
    
plot_digit(rnd_clf.feature_importances_)

cbar = plt.colorbar(ticks=[rnd_clf.feature_importances_.min(), rnd_clf.feature_importances_.max()])
cbar.ax.set_yticklabels(['Not important', 'Very important'])

plt.show()
```

## AdaBoost

*Boosting*, or *hypothesis boosting* is an ensemble method that can combine several weak learners into a strong learner. The general idea is to train predictors sequentially in such a way that each next predictor tries to correct its predecessor. The most popular among boosting methods are *AdaBoost* (which stands for *Adaptive Boosting*) and *Gradient Boosting*. Let's apply *AdaBoost* first.

The idea behind *AdaBoost* is as follows: in order to correct the predecessor's mistakes, the algorithm assigns higher weights to the training instances that the predecessor underfitted. I.e., the next step pays more attention to such instances and the algorithm focuses more on the hard cases.

Here is the step-by-step strategy:
- start with the first classifier (i.e., this can be a Decision Trees classifier)
- train it and use it to make predictions on the training set
- increase the relative weight of misclassified training instances
- train a second classifier with these new updated weights and make new predictions
- update weights using the new predictions
- continue until stopping criteria are satisfied.

For instance, suppose each instance's original weight $w^{(i)}$ is set to $\frac{1}{m}$, where $m$ is the number of instances. The first classifier is applied, and its error rate $r_1$ is computed on the training set using the equation below:

\begin{equation}
r_j = \frac{\sum_{\hat{y}^{(i)}_j \neq y^{(i)}} w^{(i)}}{\sum_{i=1}^m w^{(i)}}
\end{equation}

where $\hat{y}^{(i)}_j$ is the prediction of the $j$-th classifier on the $i$-th instance. The predictor's weight $\alpha_j$ is then estimated using:

\begin{equation}
\alpha_j = \eta log \frac{1-r_j}{r_j}
\end{equation}

where $\eta$ is the learning rate, a hyperparameter that defaults to $1$. The more accurate the predictor is, the higher its weight will be; if a predictor is guessing randomly, its weight will be close to $0$; and if it performs worse than random guessing it will get a high negative weight.

Next, the weights are updated and the misclassified instances are boosted as follows:

\begin{equation}
  w^{(i)}=\begin{cases}
    w^{(i)}, & \text{if $\hat y_j^{(i)} = y_j^{(i)}$}\\
    w^{(i)} exp(\alpha_j), & \text{if $\hat y_j^{(i)} \neq y_j^{(i)}$}
  \end{cases}
\end{equation}

Then all the instances weights are normalised (i.e., divided by $\sum_{i=1}^m w^{(i)}$). The new predictor is trained on the updated training instances, applied to the training set, its weight is computed, weights are updated again, and so on. The algorithm stops when either the predefined number of predictors is reached or a perfect predictor is found.

(This process should remind you about other sequential learning techniques we've discussed, e.g. Gradient Descent. The difference is, instead of tweaking a single predictor's parameters to minimise a cost function, AdaBoost adds predictors to the ensemble gradually making the prediction better.)

**Question**: We said earlier than one of the advantages of bagging and pasting strategies is that predictors can all be trained in parallel. Is the same applicable to AdaBoost?

At prediction time, AdaBoost computes the predictions of all the predictors and weighs them according to $\alpha_j$. The predicted class is then the one that receives the majority of weighted votes, i.e.:

\begin{equation}
 \hat y(x) = argmax_{k} \sum_{j=1; \hat y_j(x)=k} ^{N} \alpha_j
\end{equation}

with $N$ being the number of predictors.

`sklearn's` implementation of AdaBoost uses *Stagewise Additive Modeling using a Multiclass Exponential loss function* (or *SAMME* for short) multiclass version of the algorithm, which in a binary case simply defaults to AdaBoost. In addition, if the base estimators can estimate class probabilities (have `predict_proba()` method), `sklearn's` algorithm relies on class probabilities rather than predictions, which generally performs better. This version is called *SAMME.R*, where *R* stands for "Real".

Let's train [`sklearn's AdaBoostClassifier`](https://scikit-learn.org/stable/modules/generated/sklearn.ensemble.AdaBoostClassifier.html), for example using $200$ *Decision Stumps* (trees composed of a single decision node and two leaf nodes, `max_depth=1`):


```python
from sklearn.ensemble import AdaBoostClassifier

ada_clf = AdaBoostClassifier(
    DecisionTreeClassifier(max_depth=1), n_estimators=200,
    algorithm="SAMME.R", learning_rate=0.5, random_state=42)
ada_clf.fit(X_train, y_train)
```

And let's plot the decision boundary of the algorithm:


```python
plot_decision_boundary(ada_clf, X, y)
```

To get a bit more insight into how the decision boundaries change from one step to another and with a different learning rate, let's plot decision boundaries for $5$ consecutive predictors, in this case using `SVC` as the estimator. Note how the first classifiers usually get many instances wrong, while the following predictors are gradually getting better. The plot on the right presents the very same $5$ consecutive classifiers, but assigns half the learning rate (i.e., the misclassified instances weights are only boosted half as much at every iteration):


```python
m = len(X_train)

plt.figure(figsize=(11, 4))
for subplot, learning_rate in ((121, 1), (122, 0.5)):
    sample_weights = np.ones(m)
    plt.subplot(subplot)
    for i in range(5):
        svm_clf = SVC(kernel="rbf", C=0.05, gamma="auto", random_state=42)
        svm_clf.fit(X_train, y_train, sample_weight=sample_weights)
        y_pred = svm_clf.predict(X_train)
        sample_weights[y_pred != y_train] *= (1 + learning_rate)
        plot_decision_boundary(svm_clf, X, y, alpha=0.2)
        plt.title("learning_rate = {}".format(learning_rate), fontsize=16)
    if subplot == 121:
        plt.text(1.70,  -0.90, "1", fontsize=14)
        plt.text(-0.40, -0.35, "2", fontsize=14)
        plt.text(-0.55, -0.05, "3", fontsize=14)
        plt.text(-0.70,  0.20, "4", fontsize=14)
        plt.text(-0.85,  0.45, "5", fontsize=14)

plt.show()
```

In addition, here is the full list of parameters and attributes of the algorithm:


```python
list(m for m in dir(ada_clf) if not m.startswith("_") and m.endswith("_"))
```

E.g., you can check classification errors for each estimator in the boosted ensemble as follows:


```python
ada_clf.estimator_errors_
```

## Gradient Boosting

Gradient Boosting is another popular boosting algorithm. It also works sequentially adding new predictors to the ensemble, each one correcting the errors from its predecessor. Unlike AdaBoost, it's trying to fit each new predictor to the *residual errors* made by the previous predictor. 

Let's this time apply the boosting algorithm to a regression task. First, generate a synthetic noisy quadratic training set:


```python
np.random.seed(42)
X = np.random.rand(100, 1) - 0.5
y = 3*X[:, 0]**2 + 0.05 * np.random.randn(100)
```

Next, let's fit a single [`DecisionTreeRegressor`](https://scikit-learn.org/stable/modules/generated/sklearn.tree.DecisionTreeRegressor.html):


```python
from sklearn.tree import DecisionTreeRegressor

tree_reg1 = DecisionTreeRegressor(max_depth=2, random_state=42)
tree_reg1.fit(X, y)
```

Now train a second `DecisionTreeRegressor` on the residual errors made by the first predictor:


```python
y2 = y - tree_reg1.predict(X)
tree_reg2 = DecisionTreeRegressor(max_depth=2, random_state=42)
tree_reg2.fit(X, y2)
```

Finally, let's train a third regressor on the residual errors made by the second predictor:


```python
y3 = y2 - tree_reg2.predict(X)
tree_reg3 = DecisionTreeRegressor(max_depth=2, random_state=42)
tree_reg3.fit(X, y3)
```

That's it – you've just built an ensemble containing three trees. It can now be applied to new instances by making predictions based on adding up the predictions from all three trees:


```python
X_new = np.array([[0.8]])

y_pred = sum(tree.predict(X_new) for tree in (tree_reg1, tree_reg2, tree_reg3))

y_pred
```

Let's visualise the results by plotting the predictions from each tree on the left and the ensemble on the right. E.g., in the first row the ensemble contains only one tree, so its predictions are exactly the same as the first tree's predictions. In the second row, the second tree is trained on the residual errors from the first tree, and the ensemble's prediction is based on the sum of the predictions of the two trees. You can see how the predictions of the ensemble are gradually getting better:


```python
def plot_predictions(regressors, X, y, axes, label=None, style="r-", data_style="b.", data_label=None):
    x1 = np.linspace(axes[0], axes[1], 500)
    y_pred = sum(regressor.predict(x1.reshape(-1, 1)) for regressor in regressors)
    plt.plot(X[:, 0], y, data_style, label=data_label)
    plt.plot(x1, y_pred, style, linewidth=2, label=label)
    if label or data_label:
        plt.legend(loc="upper center", fontsize=16)
    plt.axis(axes)

plt.figure(figsize=(11,11))

plt.subplot(321)
plot_predictions([tree_reg1], X, y, axes=[-0.5, 0.5, -0.1, 0.8], 
                 label="$h_1(x_1)$", style="g-", data_label="Training set")
plt.ylabel("$y$", fontsize=16, rotation=0)
plt.title("Residuals and tree predictions", fontsize=16)

plt.subplot(322)
plot_predictions([tree_reg1], X, y, axes=[-0.5, 0.5, -0.1, 0.8], 
                 label="$h(x_1) = h_1(x_1)$", data_label="Training set")
plt.ylabel("$y$", fontsize=16, rotation=0)
plt.title("Ensemble predictions", fontsize=16)

plt.subplot(323)
plot_predictions([tree_reg2], X, y2, axes=[-0.5, 0.5, -0.5, 0.5], 
                 label="$h_2(x_1)$", style="g-", data_style="k+", data_label="Residuals")
plt.ylabel("$y - h_1(x_1)$", fontsize=16)

plt.subplot(324)
plot_predictions([tree_reg1, tree_reg2], X, y, axes=[-0.5, 0.5, -0.1, 0.8], 
                 label="$h(x_1) = h_1(x_1) + h_2(x_1)$")
plt.ylabel("$y$", fontsize=16, rotation=0)

plt.subplot(325)
plot_predictions([tree_reg3], X, y3, axes=[-0.5, 0.5, -0.5, 0.5], 
                 label="$h_3(x_1)$", style="g-", data_style="k+")
plt.ylabel("$y - h_1(x_1) - h_2(x_1)$", fontsize=16)
plt.xlabel("$x_1$", fontsize=16)

plt.subplot(326)
plot_predictions([tree_reg1, tree_reg2, tree_reg3], X, y, axes=[-0.5, 0.5, -0.1, 0.8], 
                 label="$h(x_1) = h_1(x_1) + h_2(x_1) + h_3(x_1)$")
plt.xlabel("$x_1$", fontsize=16)
plt.ylabel("$y$", fontsize=16, rotation=0)

plt.show()
```

A simpler way to train ensembles of *Gradient Boosted Regression Trees* is to use [`sklearn's GradientBoostingRegressor`](https://scikit-learn.org/stable/modules/generated/sklearn.ensemble.GradientBoostingRegressor.html) implementation. In addition to the typical hyperparameters that control the growth of the trees, it also includes hyperparameters that control the ensemble training, e.g. `n_estimators`:


```python
from sklearn.ensemble import GradientBoostingRegressor

gbrt = GradientBoostingRegressor(max_depth=2, n_estimators=3, learning_rate=1.0, random_state=42)
gbrt.fit(X, y)
```

**Question**: Another hyperparameter is `learning_rate`. What is it responsible for and what effect on the resulting ensemble does it have? 

Feel free to check `sklearn's` documentation, as well as look into the comparison of two ensembles: the one with a lower number of estimators and a higher learning rate, and a 'slow learner' with a lower learning rate but a higher number of estimators:


```python
gbrt_slow = GradientBoostingRegressor(max_depth=2, n_estimators=200, learning_rate=0.1, random_state=42)
gbrt_slow.fit(X, y)
```


```python
plt.figure(figsize=(11,4))

plt.subplot(121)
plot_predictions([gbrt], X, y, axes=[-0.5, 0.5, -0.1, 0.8], label="Ensemble predictions")
plt.title("learning_rate={}, n_estimators={}".format(gbrt.learning_rate, gbrt.n_estimators), fontsize=14)

plt.subplot(122)
plot_predictions([gbrt_slow], X, y, axes=[-0.5, 0.5, -0.1, 0.8])
plt.title("learning_rate={}, n_estimators={}".format(gbrt_slow.learning_rate, gbrt_slow.n_estimators), fontsize=14)

plt.show()
```

**Question**: What generalisation behaviour can you expect from each of these ensembles?


### Gradient Boosting with Early Stopping

The examples above demonstrate ensembles with a different number of estimators. Which one is better? How can you quantitatively estimate the optimal number of predictors?

One technique that you can use is called *early stopping* – it allows your algorithm to stop as soon as the validation error reaches a minimum. An easy way to implement this with `sklearn` is to use `staged_predict()` method as the code below demonstrates:


```python
import numpy as np
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_squared_error

X_train, X_val, y_train, y_val = train_test_split(X, y, random_state=100)

gbrt = GradientBoostingRegressor(max_depth=2, n_estimators=120, random_state=42)
gbrt.fit(X_train, y_train)

errors = [mean_squared_error(y_val, y_pred)
          for y_pred in gbrt.staged_predict(X_val)]
bst_n_estimators = np.argmin(errors)

gbrt_best = GradientBoostingRegressor(max_depth=2,n_estimators=bst_n_estimators, random_state=42)
gbrt_best.fit(X_train, y_train)
```

Let's visualise the validaton errors and the optimal number of estimators:


```python
min_error = np.min(errors)

plt.figure(figsize=(11, 4))

plt.subplot(121)
plt.plot(errors, "b.-")
plt.plot([bst_n_estimators, bst_n_estimators], [0, min_error], "k--")
plt.plot([0, 120], [min_error, min_error], "k--")
plt.plot(bst_n_estimators, min_error, "ko")
plt.text(bst_n_estimators, min_error*1.2, "Minimum", ha="center", fontsize=14)
plt.axis([0, 120, 0, 0.01])
plt.xlabel("Number of trees")
plt.title("Validation error", fontsize=14)

plt.subplot(122)
plot_predictions([gbrt_best], X, y, axes=[-0.5, 0.5, -0.1, 0.8])
plt.title("Best model (%d trees)" % bst_n_estimators, fontsize=14)

plt.show()
```

Finally, an alternative to training a large number of estimators and then looking back in order to find the optimal number of those is to allow the algorithm to stop training when the validation error does not improve over a number of consecutive iterations (e.g., $5$ in the code below). For that to work, set the `warm_start` to `True`:


```python
gbrt = GradientBoostingRegressor(max_depth=2, warm_start=True, random_state=42)

min_val_error = float("inf")
error_going_up = 0
for n_estimators in range(1, 120):
    gbrt.n_estimators = n_estimators
    gbrt.fit(X_train, y_train)
    y_pred = gbrt.predict(X_val)
    val_error = mean_squared_error(y_val, y_pred)
    if val_error < min_val_error:
        min_val_error = val_error
        error_going_up = 0
    else:
        error_going_up += 1
        if error_going_up == 5:
            break  # early stopping

print(gbrt.n_estimators)
```

# Your tasks

1. Run the code in the notebook. During the practical session, be prepared to discuss the methods and answer the questions from this notebook:

    * When applied to a new, test instance bagging classifier aggregates the predictions of all predictors and estimates the statistical mode. Is this similar to the hard or the soft voting strategy?
    * [`BaggingClassifier`](https://scikit-learn.org/stable/modules/generated/sklearn.ensemble.BaggingClassifier.html) automatically performs soft voting when the base classifier can estimate class probabilities (i.e., has a `predict_proba()` method). Which strategy – hard or soft voting – does `BaggingClassifier` follow with the `DecisionTreeClassifier` as the base classifier?
    * What do the feature importance weights obtained on the *iris* dataset suggest? Are they similar to your observations from the previous experiments on this dataset?
    * We said earlier than one of the advantages of bagging and pasting strategies is that predictors can all be trained in parallel. Is the same applicable to AdaBoost?
    * One of the hyperparameters of the boosting algorithms is learning_rate. What is it responsible for and what effect on the resulting ensemble does it have?
    * What generalisation behaviour can you expect from each of the tree regressor ensembles (using different learning rates and different number of estimators – see above)?

2. Apply ensemble techniques of your choice to one of the datasets you've worked on during the previous practicals and report your findings.

3. **Optional**: If you want more practice with these techniques, try implementing a [*stacking algorithm*](https://en.wikipedia.org/wiki/Ensemble_learning#Stacking). There is no available implementation of this approach in `sklearn` so it needs to be implemented from scratch. The idea is as follows:
    * Split a dataset of your choice into three subsets – training, validation and test.
    * Train a number of predictors on the training data and apply them to the validation data.
    * Treat the predictions of those predictors on the validation data to generate new training set: you can use the predictions of the estimators as new features (i.e., each training instance will have as many features as the number of predictors you originally used), and the validation data targets as the new training instances targets. 
    * Now train a *blender* – a classifier of your choice – on the new training data created this way. Together with the original classifiers, the blender forms a stacking ensemble.
    * Finally, apply your original classifiers to the test set, feed their predictions to the blender, and use blender's output as the prediction on the test set.


```python

```

491475db87bfc7b2ce6814dd4604453812c82006
ipynb
Jupyter Notebook
DSPNP_practical3/DSPNP_practical3.ipynb
marcus800/cl-datasci-pnp-2021
aea4a1e1aaeac895c595d67f328485157f1e2b39
Apache-2.0
2020-11-06T13:03:19.000Z
2022-02-05T17:27:23.000Z
DSPNP_practical3/DSPNP_practical3.ipynb
marcus800/cl-datasci-pnp-2021
aea4a1e1aaeac895c595d67f328485157f1e2b39
Apache-2.0
DSPNP_practical3/DSPNP_practical3.ipynb
marcus800/cl-datasci-pnp-2021
aea4a1e1aaeac895c595d67f328485157f1e2b39
Apache-2.0
2020-11-08T09:34:50.000Z
2022-02-02T21:55:08.000Z
Qwen/Qwen-72B
1. YES
2. YES