<jupyter_start><jupyter_text>Ensemble learning: when many are better that the oneIn this notebook, we will go in depth into algorithms which combine severalsimple learners (e.g. decision tree, linear model, etc.). We willsee that combining simple learners will result in a more powerful and robustlearner.We will focus on two families of ensemble methods:* ensemble using bootstrap (e.g. bagging and random-forest);* ensemble using boosting (e.g. adaptive boosting and gradient-boosting  decision tree). Benefit of ensemble method at a first glanceIn this section, we will give a quick demonstration on the power of combiningseveral learners instead of fine-tuning a single learner.We will start by loading the "California Housing" dataset.<jupyter_code>from sklearn.datasets import fetch_california_housing

california_housing = fetch_california_housing(as_frame=True)
df = california_housing.frame
X, y = california_housing.data, california_housing.target<jupyter_output><empty_output><jupyter_text>In this dataset, we want to predict the median house value in some districtin California based on demographic and geographic data.<jupyter_code>df.head()<jupyter_output><empty_output><jupyter_text>We start by learning a single decision tree regressor. As we previouslypresented in the "tree in depth" notebook, this learner needs to be tuned toovercome over- or under-fitting. Indeed, the default parameters will notnecessarily lead to an optimal decision tree. Instead of using the defaultvalue, we should search via cross-validation the optimal value of theimportant parameters such as `max_depth`, `min_samples_split`, or`min_samples_leaf`.We recall that we need to tune these parameters, as decision treestend to overfit the training data if we grow deep trees, but there are norules on how to limit the parameters. Thus, not making a search could lead usto have an underfitted model.First, let's keep a set of data to test our final model.<jupyter_code>from sklearn.model_selection import train_test_split

X_train, X_test, y_train, y_test = train_test_split(
    X, y, random_state=0,
)<jupyter_output><empty_output><jupyter_text>We will first make a grid-search to fine-tune the parameters that wementioned earlier.<jupyter_code>from sklearn.model_selection import GridSearchCV
from sklearn.tree import DecisionTreeRegressor

param_grid = {
    "max_depth": [3, 5, 8, None],
    "min_samples_split": [2, 10, 30, 50],
    "min_samples_leaf": [0.01, 0.05, 0.1, 1],
}
cv = 3
tree = GridSearchCV(
    DecisionTreeRegressor(random_state=0),
    param_grid=param_grid,
    cv=cv,
    n_jobs=-1,
)

tree.fit(X_train, y_train)<jupyter_output><empty_output><jupyter_text>We can create a dataframe storing the important information collected duringthe tuning of the parameters and investigate the results.<jupyter_code>import pandas as pd

cv_results = pd.DataFrame(tree.cv_results_)
interesting_columns = [
    "param_max_depth",
    "param_min_samples_split",
    "param_min_samples_leaf",
    "mean_test_score",
    "rank_test_score",
    "mean_fit_time",
]
cv_results = cv_results[interesting_columns].sort_values(by="rank_test_score")
cv_results<jupyter_output><empty_output><jupyter_text>From theses results, we can see that the best parameters is the combinationwhere the depth of the tree is not limited and the minimum number of samplesto create a leaf is also equal to 1 (the default values) and theminimum number of samples to make a split of 50 (much higher than the defaultvalue.It is interesting to look at the total amount of time it took to fit allthese different models. In addition, we can check the performance of theoptimal decision tree on the left-out testing data.<jupyter_code>total_fitting_time = (cv_results["mean_fit_time"] * cv).sum()
print(
    f"Required training time of the GridSearchCV: "
    f"{total_fitting_time:.2f} seconds"
)
print(
    f"Best R2 score of a single tree: {tree.best_score_:.3f}"
)<jupyter_output>Required training time of the GridSearchCV: 10.39 seconds
Best R2 score of a single tree: 0.691<jupyter_text>Hence, we have a model that has an $R^2$ score below 0.7. The amount of timeto find the best learner depends on the number of folds used during thecross-validation in the grid-search multiplied by the number of parameters.Therefore, the computational cost is quite high.Now we will use an ensemble method called bagging. More details about thismethod will be discussed in the next section. In short, this method will usea base regressor (i.e. decision tree regressors) and will train several ofthem on a slightly modified version of the training set. Then, thepredictions of all these base regressors will be combined by averaging.Here, we will use 50 decision trees and check the fitting time as well asthe performance on the left-out testing data. It is important to note thatwe are not going to tune any parameter of the decision tree.<jupyter_code>from time import time
from sklearn.ensemble import BaggingRegressor

base_estimator = DecisionTreeRegressor(random_state=0)
bagging_regressor = BaggingRegressor(
    base_estimator=base_estimator, n_estimators=50, random_state=0,
)

start_fitting_time = time()
bagging_regressor.fit(X_train, y_train)
elapsed_fitting_time = time() - start_fitting_time

print(f"Elapsed fitting time: {elapsed_fitting_time:.2f} seconds")
print(f"R2 score: {bagging_regressor.score(X_test, y_test):.3f}")<jupyter_output>Elapsed fitting time: 4.92 seconds
R2 score: 0.792<jupyter_text>We can see that the computation time is much shorter for training the fullensemble than for the parameter search of a single tree. In addition, thescore is significantly improved with a $R^2$ close to 0.8. Furthermore, notethat this result is obtained before any parameter tuning. This shows themotivation behind the use of an ensemble learner: it gives a relatively goodbaseline with decent performance without any parameter tuning.Now, we will discuss in detail two ensemble families: bagging andboosting. BaggingBagging stands for Bootstrap AGGregatING. It uses bootstrap samplesto learn several models. At predict time, the predictions of each learnerare aggregated to give the final predictions.Let's define a simple dataset (which we have used before in a previousnotebook).<jupyter_code>import numpy as np
import matplotlib.pyplot as plt

rng = np.random.RandomState(0)


def generate_data(n_samples=50, sorted=False):
    x_max, x_min = 1.4, -1.4
    len_x = x_max - x_min
    x = rng.rand(n_samples) * len_x - len_x / 2
    noise = rng.randn(n_samples) * 0.3
    y = x ** 3 - 0.5 * x ** 2 + noise
    if sorted:
        sorted_idx = np.argsort(x)
        x, y = x[sorted_idx], y[sorted_idx]
    return x, y


x, y = generate_data(n_samples=50)

plt.scatter(x, y,  color='k', s=9)
plt.xlabel("Feature")
_ = plt.ylabel("Target")<jupyter_output><empty_output><jupyter_text>The link between our feature and the target to predict is non-linear.However, a decision tree is capable of fitting such data<jupyter_code>tree = DecisionTreeRegressor(max_depth=3, random_state=0)
tree.fit(x.reshape(-1, 1), y)

grid = np.linspace(np.min(x), np.max(x), num=300)
y_pred = tree.predict(grid.reshape(-1, 1))

plt.scatter(x, y, color="k", s=9)
plt.plot(grid, y_pred, label="Tree fitting")
plt.xlabel("Feature")
plt.ylabel("Target")
plt.legend()<jupyter_output><empty_output><jupyter_text>Let's see how we can use bootstraping to learn several trees. Bootstrap sampleA bootstrap sample corresponds to a resampling, with replacement, of theoriginal dataset, a sample that is the same size as theoriginal dataset. Thus, the bootstrap sample will contain somedata points several times while some of the original data points willnot be present.We will create a function that given `x` and `y` will return a bootstrapsample `x_bootstrap` and `y_bootstrap`.<jupyter_code>def bootstrap_sample(x, y):
    bootstrap_indices = rng.choice(
        np.arange(y.shape[0]), size=y.shape[0], replace=True,
    )
    x_bootstrap_sample = x[bootstrap_indices]
    y_bootstrap_sample = y[bootstrap_indices]
    return x_bootstrap_sample, y_bootstrap_sample<jupyter_output><empty_output><jupyter_text>We will generate 3 bootstrap samples and qualitatively check the differencewith the original dataset.<jupyter_code>n_bootstrap = 3
_, axs = plt.subplots(
    ncols=n_bootstrap, figsize=(16, 6), sharex=True, sharey=True
)

for idx, (ax, _) in enumerate(zip(axs, range(n_bootstrap))):
    x_bootstrap_sample, y_bootstrap_sample = bootstrap_sample(x, y)
    ax.scatter(
        x_bootstrap_sample, y_bootstrap_sample,
    )
    ax.set_title(f"Bootstrap sample #{idx}")
    ax.set_ylabel("Target")
    ax.set_xlabel("Feature")<jupyter_output><empty_output><jupyter_text>We observe that the 3 generated bootstrap samples are all different. Toconfirm this intuition, we can check the number of unique samples in thebootstrap samples.<jupyter_code># we need to generate a larger set to have a good estimate
x_huge, y_huge = generate_data(n_samples=10000)
x_bootstrap_sample, y_bootstrap_sample = bootstrap_sample(x_huge, y_huge)

print(
    f"Percentage of samples present in the original dataset: "
    f"{np.unique(x_bootstrap_sample).size / x_bootstrap_sample.size * 100:.1f}"
    f"%"
)<jupyter_output>Percentage of samples present in the original dataset: 63.2%<jupyter_text>Theoretically, 63.2% of the original data points of the original dataset willbe present in the bootstrap sample. The other 36.8% are just repeatedsamples.So, we are able to generate many datasets, all slightly different. Now, wecan fit a decision tree to each of these datasets and each decisiontree shall be slightly different as well.<jupyter_code>_, axs = plt.subplots(
    ncols=n_bootstrap, figsize=(16, 6), sharex=True, sharey=True,
)

forest = []
for idx, (ax, _) in enumerate(zip(axs, range(n_bootstrap))):
    x_bootstrap_sample, y_bootstrap_sample = bootstrap_sample(x, y)
    ax.scatter(x_bootstrap_sample, y_bootstrap_sample)

    forest.append(
        DecisionTreeRegressor(max_depth=3, random_state=0).fit(
            x_bootstrap_sample.reshape(-1, 1), y_bootstrap_sample
        )
    )

    grid = np.linspace(np.min(x), np.max(x), num=300)
    y_pred = forest[-1].predict(grid.reshape(-1, 1))
    ax.plot(grid, y_pred, linewidth=3, label="Fitted tree")

    ax.legend()
    ax.set_ylabel("Target")
    ax.set_xlabel("Features")
    ax.set_title(f"Bootstrap #{idx}")<jupyter_output><empty_output><jupyter_text>We can plot these decision functions on the same plot to see the difference.<jupyter_code>_, ax = plt.subplots()
ax.scatter(x, y, color="k", alpha=0.4)
y_pred_forest = []
for tree_idx, tree in enumerate(forest):
    y_pred = tree.predict(grid.reshape(-1, 1))
    ax.plot(
        grid,
        y_pred,
        "--",
        label=f"Tree #{tree_idx} predictions",
        linewidth=3,
        alpha=0.8,
    )
    y_pred_forest.append(y_pred)

plt.xlabel("Feature")
plt.ylabel("Target")
_ = plt.legend()<jupyter_output><empty_output><jupyter_text>AggregatingOnce our trees are fitted and we are able to get predictions for each ofthem, we also need to combine them. In regression, the most straightforwardapproach is to average the different predictions from all learners. We canplot the averaged predictions from the previous example.<jupyter_code>_, ax = plt.subplots()
ax.scatter(x, y, color="k", alpha=0.4)
y_pred_forest = []
for tree_idx, tree in enumerate(forest):
    y_pred = tree.predict(grid.reshape(-1, 1))
    ax.plot(
        grid,
        y_pred,
        "--",
        label=f"Tree #{tree_idx} predictions",
        linewidth=3,
        alpha=0.5,
    )
    y_pred_forest.append(y_pred)

# Averaging the predictions
y_pred_forest = np.mean(y_pred_forest, axis=0)
ax.plot(
    grid,
    y_pred_forest,
    "-",
    label="Averaged predictions",
    linewidth=3,
    alpha=0.8,
)

plt.xlabel("Feature")
plt.ylabel("Target")
_ = plt.legend()<jupyter_output><empty_output><jupyter_text>The unbroken red line shows the averaged predictions, which would be thefinal preditions given by our 'bag' of decision tree regressors. Random forestA popular machine-learning algorithm is the random forest. A Random forestis a modification of the bagging algorithm. In bagging, any classifier orregressor can be used. In a random forest, the base classifier or regressormust be a decision tree. In our previous example, we used a decisiontree but we could have used a linear model as the regressor for ourbagging algorithm.In addition, random forest is different from bagging when used withclassifiers: when searching for the best split, only a subset of the originalfeatures are used. By default, this subset of feature is equal to the squareroot of the total number of features. In regression, the total number ofavailable features will be used.We will illustrate the usage of a random forest and compare it with thebagging regressor on the "California housing" dataset.<jupyter_code>X, y = california_housing.data, california_housing.target
X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=0)

from sklearn.ensemble import RandomForestRegressor

random_forest = RandomForestRegressor(
    n_estimators=100, random_state=0, n_jobs=-1
)
bagging = BaggingRegressor(
    base_estimator=DecisionTreeRegressor(random_state=0),
    n_estimators=100,
    n_jobs=-1,
)

random_forest.fit(X_train, y_train)
bagging.fit(X_train, y_train)

print(
    f"Performance of random forest: {random_forest.score(X_test, y_test):.3f}"
)
print(f"Performance of bagging: {bagging.score(X_test, y_test):.3f}")<jupyter_output>Performance of random forest: 0.794
Performance of bagging: 0.792<jupyter_text>Notice that we don't need to provide a `base_estimator` parameter to`RandomForestRegressor`, it is always a tree classifier. Also note that thescores are almost identical. This is because ourproblem is a regression problem and therefore, the number of features usedin random forest and bagging is the same.For classification problems, we would need to pass a tree model instancewith the parameter `max_features="sqrt"` to `BaggingRegressor` if we wantedit to have the same behaviour as the random forest classifier. Classifiers detailsUp to now, we have only focused on regression problems. There is a littledifference between regression and classification.First, the `base_estimator` should be chosen in line with the problem thatis solved: use a classifier with a classification problem and a regressorwith a regression problem.Then, the aggregation method is different in regression and classification:the averaged predictions is computed in regression while the majority class(weighted by the probabilities) is predicted in classification. SummaryWe saw in this section two algorithms that use bootstrap samples to createan ensemble of classifiers or regressors. These algorithms train severallearners on different bootstrap samples. The predictions are thenaggregated. This operation can be done in a very efficient manner since thetraining of each learner can be done in parallel. BoostingWe recall that bagging builds an ensemble in a parallel manner: each learneris trained independently from each other. The idea behind boosting isdifferent. The ensemble is a sequence of learners where the`Nth` learner requires all previous learners, from 1 to `N-1`.Intuitively, bagging adds learners to the ensemble to correct themistakes of the previous learners. We will start with analgorithm named Adaptive Boosting (AdaBoost) to get some intuition about themain ideas behind boosting. Adaptive Boosting (AdaBoost)We will first focus on AdaBoost, which we will use for a classificationproblem.We will load the "penguin" dataset used in the "tree in depth" notebook.We will predict penguin species from the features culmen length and depth.<jupyter_code>data = pd.read_csv("../datasets/penguins.csv")

# select the features of interest
culmen_columns = ["Culmen Length (mm)", "Culmen Depth (mm)"]
target_column = "Species"

data = data[culmen_columns + [target_column]]
data[target_column] = data[target_column].str.split().str[0]
data = data.dropna()

X, y = data[culmen_columns], data[target_column]<jupyter_output><empty_output><jupyter_text>In addition, we are also using on the function used the previous"tree in depth" notebookto plot the decision function of the tree.<jupyter_code>import seaborn as sns


def plot_decision_function(X, y, clf, sample_weight=None, fit=True, ax=None):
    """Plot the boundary of the decision function of a classifier."""
    from sklearn.preprocessing import LabelEncoder

    if fit:
        clf.fit(X, y, sample_weight=sample_weight)

    # create a grid to evaluate all possible samples
    plot_step = 0.02
    feature_0_min, feature_0_max = (
        X.iloc[:, 0].min() - 1,
        X.iloc[:, 0].max() + 1,
    )
    feature_1_min, feature_1_max = (
        X.iloc[:, 1].min() - 1,
        X.iloc[:, 1].max() + 1,
    )
    xx, yy = np.meshgrid(
        np.arange(feature_0_min, feature_0_max, plot_step),
        np.arange(feature_1_min, feature_1_max, plot_step),
    )

    # compute the associated prediction
    Z = clf.predict(np.c_[xx.ravel(), yy.ravel()])
    Z = LabelEncoder().fit_transform(Z)
    Z = Z.reshape(xx.shape)

    # make the plot of the boundary and the data samples
    if ax is None:
        _, ax = plt.subplots()
    ax.contourf(xx, yy, Z, alpha=0.4)
    sns.scatterplot(
        data=pd.concat([X, y], axis=1),
        x=X.columns[0],
        y=X.columns[1],
        hue=y.name,
        ax=ax,
    )<jupyter_output><empty_output><jupyter_text>We will purposely train a shallow decision tree. Since the tree is shallow,it is unlikely to overfit and some of the training examples will even bemisclassified on the training set.<jupyter_code>from sklearn.tree import DecisionTreeClassifier

tree = DecisionTreeClassifier(max_depth=2, random_state=0)

_, ax = plt.subplots()
plot_decision_function(X, y, tree, ax=ax)

# find the misclassified samples
y_pred = tree.predict(X)
misclassified_samples_idx = np.flatnonzero(y != y_pred)

ax.plot(
    X.iloc[misclassified_samples_idx, 0],
    X.iloc[misclassified_samples_idx, 1],
    "+k",
    label="Misclassified samples",
)
ax.legend()<jupyter_output><empty_output><jupyter_text>We observe that several samples have been misclassified by theclassifier.We mentioned that boosting relies on creating a new classifier which tries tocorrect these misclassifications. In scikit-learn, learners support aparameter `sample_weight` which forces the learner to pay more attention tosamples with higher weights, during the training.This parameter is set when calling`classifier.fit(X, y, sample_weight=weights)`.We will use this trick to create a new classifier by 'discarding' allcorrectly classified samples and only considering the misclassified samples.Thus, mosclassified samples will be assigned a weight of 1 while wellclassified samples will assigned to a weight of 0.<jupyter_code>sample_weight = np.zeros_like(y, dtype=int)
sample_weight[misclassified_samples_idx] = 1

_, ax = plt.subplots()
plot_decision_function(X, y, tree, sample_weight=sample_weight, ax=ax)

ax.plot(
    X.iloc[misclassified_samples_idx, 0],
    X.iloc[misclassified_samples_idx, 1],
    "+k",
    label="Previous misclassified samples",
)
ax.legend()<jupyter_output><empty_output><jupyter_text>We see that the decision function drastically changed. Qualitatively,we see that the previously misclassified samples are now correctlyclassified.<jupyter_code>y_pred = tree.predict(X)
newly_misclassified_samples_idx = np.flatnonzero(y != y_pred)
remaining_misclassified_samples_idx = np.intersect1d(
    misclassified_samples_idx, newly_misclassified_samples_idx
)

print(
    f"Number of samples previously misclassified and still misclassified: "
    f"{len(remaining_misclassified_samples_idx)}"
)<jupyter_output>Number of samples previously misclassified and still misclassified: 0<jupyter_text>However, we are making mistakes on previously well classified samples. Thus,we get the intuition that we should weight the predictions of each classifierdifferently, most probably by using the number of mistakes each classifieris making.So we could use the classification error to combine both trees.ensemble_weight = [    (y.shape[0] - len(misclassified_samples_idx)) / y.shape[0],    (y.shape[0] - len(newly_misclassified_samples_idx)) / y.shape[0],]ensemble_weight The first classifier was 94% accurate and the second one 69% accurate.Therefore, when predicting a class, we should trust the first classifierslightly more than the second one. We could use these accuracy values toweight the predictions of each learner.To summarize, boosting learns several classifiers, each of which willfocus more or less on specific samples of the dataset. Boosting is thusdifferent from bagging: here we never resample our dataset, we just assigndifferent weights to the original dataset.Boosting requires some strategy to combine the learners together:* one needs to define a way to compute the weights to be assigned  to samples;* one needs to assign a weight to each learner when making predictions.Indeed, we defined a really simple scheme to assign sample weights andlearner weights. However, there are statistical theory for how thesethese sample and learner weights can be optimally calculated.FIXME: I think we should add a reference to ESL here.We will use the AdaBoost classifier implemented in scikit-learn andlook at the underlying decision tree classifiers trained.<jupyter_code>from sklearn.ensemble import AdaBoostClassifier

adaboost = AdaBoostClassifier(
    base_estimator=DecisionTreeClassifier(max_depth=3, random_state=0),
    n_estimators=3,
    algorithm="SAMME",
    random_state=0,
)
adaboost.fit(X, y)

_, axs = plt.subplots(ncols=3, figsize=(16, 6))

for ax, tree in zip(axs, adaboost.estimators_):
    plot_decision_function(X, y, tree, fit=False, ax=ax)

print(f"Weight of each classifier: {adaboost.estimator_weights_}")
print(f"Error of each classifier: {adaboost.estimator_errors_}")<jupyter_output>Weight of each classifier: [3.58351894 3.46901998 3.03303773]
Error of each classifier: [0.05263158 0.05864198 0.08787269]<jupyter_text>We see that AdaBoost has learnt three different classifiers each of whichfocuses on different samples. Looking at the weights of each learner, we seethat the ensemble gives the highest weight to the first classifier. Thisindeed makes sense when we look at the errors of each classifier. The firstclassifier also has the highest classification performance.While AdaBoost is a nice algorithm to demonsrate the internal machinery ofboostingalgorithms, it is not the most efficient machine-learning algorithm.The most efficient algorithm based on boosting is the gradient-boostingdecision tree (GBDT) algorithm which we will discuss now. Gradient-boosting decision tree (GBDT)Gradient-boosting differs from AdaBoost due to the following reason: insteadof assigning weights to specific samples, GBDT will fit a decisiontree on the residuals (hence the name "gradient") of the previous tree.Therefore, each new added tree in the ensemble predicts the error made by theprevious learner instead of predicting the target directly.In this section, we will provide some intuition about the way learnersare combined to give the final prediction. In this regard, let's go backto our regression problem which is more intuitive for demonstrating theunderlying machinery.<jupyter_code>x, y = generate_data(sorted=True)

plt.scatter(x, y, color="k", s=9)
plt.xlabel("Feature")
_ = plt.ylabel("Target")<jupyter_output><empty_output><jupyter_text>As we previously discussed, boosting will be based on assembling a sequenceof learners. We will start by creating a decision tree regressor. We will fixthe depth of the tree so that the resulting learner will underfit the data.<jupyter_code>tree = DecisionTreeRegressor(max_depth=3, random_state=0)
tree.fit(x.reshape(-1, 1), y)

grid = np.linspace(np.min(x), np.max(x), num=200)
y_pred_grid_raw = tree.predict(grid.reshape(-1, 1))

plt.scatter(x, y, color="k", s=9)
plt.xlabel("Feature")
plt.ylabel("Target")
line_predictions = plt.plot(grid, y_pred_grid_raw, "--")

y_pred_raw = tree.predict(x.reshape(-1, 1))
for idx in range(len(y)):
    lines_residuals = plt.plot(
        [x[idx], x[idx]], [y[idx], y_pred_raw[idx]], color="red",
    )

_ = plt.legend(
    [line_predictions[0], lines_residuals[0]], ["Fitted tree", "Residuals"]
)<jupyter_output><empty_output><jupyter_text>Since the tree underfits the data, its accuracy is farfrom perfect on the training data. We can observe this in the figure bylooking at the difference between the predictions and the ground-truth data.We represent these errors, called "Residuals", by unbroken red lines.Indeed, our initial tree was not expressive enough to handle the complexityof the data, as shown by the residuals.In a gradient-boosting algorithm, the idea is to create a second treewhich, given the same data `x`, will try to predict the residuals instead ofthe target `y`. We would therefore have a tree that is able to predict theerrors made by the initial tree.Let's train such a tree.<jupyter_code>residuals = y - y_pred_raw

tree_residuals = DecisionTreeRegressor(max_depth=5, random_state=0)
tree_residuals.fit(x.reshape(-1, 1), residuals)

y_pred_grid_residuals = tree_residuals.predict(grid.reshape(-1, 1))

plt.scatter(x, residuals, color="k", s=9)
plt.xlabel("Feature")
plt.ylabel("Residuals")
line_predictions = plt.plot(grid, y_pred_grid_residuals, "--")

y_pred_residuals = tree_residuals.predict(x.reshape(-1, 1))
for idx in range(len(y)):
    lines_residuals = plt.plot(
        [x[idx], x[idx]], [residuals[idx], y_pred_residuals[idx]], color="red",
    )

_ = plt.legend(
    [line_predictions[0], lines_residuals[0]], ["Fitted tree", "Residuals"]
)<jupyter_output><empty_output><jupyter_text>We see that this new tree only manages to fit some of the residuals.We will focus on the last sample in `x` andexplain how the predictions of both trees are combined.<jupyter_code>_, axs = plt.subplots(ncols=2, figsize=(12, 6), sharex=True)

axs[0].scatter(x, y, color="k", s=9)
axs[0].set_xlabel("Feature")
axs[0].set_ylabel("Target")
axs[0].plot(grid, y_pred_grid_raw, "--")

axs[1].scatter(x, residuals, color="k", s=9)
axs[1].set_xlabel("Feature")
axs[1].set_ylabel("Residuals")
plt.plot(grid, y_pred_grid_residuals, "--")

for idx in range(len(y)):
    axs[0].plot(
        [x[idx], x[idx]], [y[idx], y_pred_raw[idx]], color="red",
    )
    axs[1].plot(
        [x[idx], x[idx]], [residuals[idx], y_pred_residuals[idx]], color="red",
    )

axs[0].scatter(x[-1], y[-1], color="tab:orange", label="Sample of interest")
axs[1].scatter(
    x[-1], residuals[-1], color="tab:orange", label="Sample of interest"
)

axs[0].set_xlim([1.1, 1.4])
axs[1].set_xlim([1.1, 1.4])
axs[0].set_ylim([0, 2])
axs[1].set_ylim([-0.7, 0.7])
axs[0].legend()
axs[1].legend()<jupyter_output><empty_output><jupyter_text>For our sample of interest, our initial tree is making an error (smallresidual). Whenfitting the second tree, the residual in this case is perfectly fitted andpredicted. We will quantitatively check this prediction using the fittedtree. First, let's check the prediction of the initial tree and compare itwith the true value.<jupyter_code>x_max = x[-1]
y_true = y[-1]

print(f"True value to predict for f(x={x_max:.3f}) = {y_true:.3f}")

y_pred_first_tree = tree.predict([[x_max]])[0]
print(
    f"Prediction of the first decision tree for x={x_max:.3f}: "
    f"y={y_pred_first_tree:.3f}"
)
print(f"Error of the tree: {y_true - y_pred_first_tree:.3f}")<jupyter_output>True value to predict for f(x=1.384) = 1.676
Prediction of the first decision tree for x=1.384: y=1.555
Error of the tree: 0.121<jupyter_text>As we visually observed, we have a small error. Now, we can use the secondtree to try to predict this residual.<jupyter_code>print(
    f"Prediction of the residual for x={x_max:.3f}: "
    f"{tree_residuals.predict([[x_max]])[0]:.3f}"
)<jupyter_output>Prediction of the residual for x=1.384: 0.121<jupyter_text>We see that our second tree is capable of prediting the exact residual(error) of our first tree. Therefore, we can predict the value of`x` by summing the prediction of the all trees in the ensemble.<jupyter_code>y_pred_first_and_second_tree = (
    y_pred_first_tree + tree_residuals.predict([[x_max]])[0]
)
print(
    f"Prediction of the first and second decision trees combined for "
    f"x={x_max:.3f}: y={y_pred_first_and_second_tree:.3f}"
)
print(f"Error of the tree: {y_true - y_pred_first_and_second_tree:.3f}")<jupyter_output>Prediction of the first and second decision trees combined for x=1.384: y=1.676
Error of the tree: 0.000<jupyter_text>We chose a sample for which only two trees were enough to make the perfectprediction. However, we saw in the previous plot that two trees were notenough to correct the residuals of all samples. Therefore, one needs toadd several trees to the ensemble to successfully correct the error.We will compare the performance of random-forest and gradient boosting onthe California housing dataset.<jupyter_code>from sklearn.ensemble import GradientBoostingRegressor

X, y = california_housing.data, california_housing.target
X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=0,)

gradient_boosting = GradientBoostingRegressor(n_estimators=200)
start_time = time()
gradient_boosting.fit(X_train, y_train)
fit_time_gradient_boosting = time() - start_time

random_forest = RandomForestRegressor(n_estimators=200, n_jobs=-1)
start_time = time()
random_forest.fit(X_train, y_train)
fit_time_random_forest = time() - start_time

print(
    f"The performance of gradient-boosting are:"
    f"{gradient_boosting.score(X_test, y_test):.3f}"
)
print(f"Fitting time took: {fit_time_gradient_boosting:.2f} seconds")

print(
    f"The performance of random-forest are:"
    f"{random_forest.score(X_test, y_test):.3f}"
)
print(f"Fitting time took: {fit_time_random_forest:.2f} seconds")<jupyter_output>The performance of gradient-boosting are:0.804
Fitting time took: 7.04 seconds
The performance of random-forest are:0.796
Fitting time took: 7.85 seconds<jupyter_text>In term of computation performance, the forest can be parallelized and willbenefit from the having multiple CPUs. In terms of scoring performance, bothalgorithms lead to very close results. Parameter consideration with random forest and gradient-boostingIn the previous section, we did not discuss the parameters of random forestand gradient-boosting. However, there are a couple of things to keep in mindwhen setting these parameters. Random forestThe main parameter to tune with random forest is the `n_estimators`parameter. In general, the more trees in the forest, the better theperformance will be. However, it will slow down the fitting and predictiontime. So one has to balance compute time and performance when setting thenumber of estimators when putting such learner in production.The `max_depth` parameter could also be tuned. Sometimes, there is no needto have fully grown trees. However, be aware that with random forest, treesare generally deep since we are seeking to overfit the learners on thebootstrap samples because this will be mitigated by combining them.Assembling underfitted trees (i.e. shallow trees) might also lead to anunderfitted forest.<jupyter_code>param_grid = {
    "n_estimators": [10, 20, 30],
    "max_depth": [3, 5, None],
}
grid_search = GridSearchCV(
    RandomForestRegressor(n_jobs=-1), param_grid=param_grid, n_jobs=-1,
)
grid_search.fit(X_train, y_train)

columns = ["params", "mean_test_score", "rank_test_score"]
cv_results = pd.DataFrame(grid_search.cv_results_)
cv_results[columns].sort_values(by="rank_test_score")<jupyter_output><empty_output><jupyter_text>We can observe that in our grid-search, the largest `max_depth` together withlargest `n_estimators` led to the best performance. Gradient-boosting decision treeFor gradient-boosting, parameter tuning is a combination of severalparameters instead of setting one after the other each parameter. Theimportant parameters are `n_estimators`, `max_depth`, and `learning_rate`.Let's first discuss the `max_depth` parameter. We saw in the section ongradient-boosting that the algorithm fits the error of the previous treein the ensemble. Thus, fitting fully grown trees will be detrimental. Indeed,the first tree of the ensemble would perfectly fit (overfit) the data andthus no subsequent tree would be required, since there would be no residuals.Therefore, the tree used in gradient-boosting should have a low depth,typically between 3 to 8 levels.With this consideration in mind, the deeper the trees, the faster theresiduals will be corrected and less learners are required. So `n_estimators`should be increased if `max_depth` is lower.Finally, we have overlooked the impact of the `learning_rate` parameter uptill now. Whenfitting the residuals one could choose if the tree should try to correct allpossible errors or only a fraction of them. The learning-rate allows you tocontrol this behaviour. A small learning-rate value would only correct theresiduals of very few samples. If a large learning-rate is set (e.g., 1),we would fit the residuals of all samples. So, with a very lowlearning-rate, we will need more estimators to correct the overall error.However, a too large learning-rate tends to obtain an overfitted ensemble,similar to having a too large tree depth.<jupyter_code>param_grid = {
    "n_estimators": [10, 30, 50],
    "max_depth": [3, 5, None],
    "learning_rate": [0.1, 1],
}
grid_search = GridSearchCV(
    GradientBoostingRegressor(), param_grid=param_grid, n_jobs=-1,
)
grid_search.fit(X_train, y_train)

cv_results = pd.DataFrame(grid_search.cv_results_)
cv_results[columns].sort_values(by="rank_test_score")<jupyter_output><empty_output><jupyter_text>Accelerating gradient-boostingWe previously mentioned that random-forest is an efficient algorithm sinceeach tree of the ensemble can be fitted at the same time independently.Therefore, the algorithm scales efficiently with both the number of CPUs andthe number of samples.In gradient-boosting, the algorithm is a sequential algorithm. It requiresthe `N-1` trees to have been fit to be able to fit the tree at stage `N`.Therefore, the algorithm isquite computationally expensive. The most expensive part in this algorithm isthe search for the best split in the tree which is a brute-forceapproach: all possible split are evaluated and the best one is picked. Weexplained this process in the notebook "tree in depth", whichyou can refer to.To accelerate the gradient-boosting algorithm, one could reduce the number ofsplits to be evaluated. As a consequence, the performance of such atree would be reduced. However, since we are combining several trees in agradient-boosting, we can add more estimators to overcomethis issue.This algorithm is called `HistGradientBoostingClassifier` and`HistGradientBoostingRegressor`. Each feature in the dataset `X` is firstbinned by computing histograms which are later used to evaluate the potentialsplits. The numberof splits to evaluate is then much smaller. This algorithm becomes much moreefficient than gradient bossting when the dataset has 10,000+ samples.Below we will give an example of a large dataset and we can comparecomputation time with the earlier experiment in the previous section.<jupyter_code>from sklearn.experimental import enable_hist_gradient_boosting
from sklearn.ensemble import HistGradientBoostingRegressor

histogram_gradient_boosting = HistGradientBoostingRegressor(
    max_iter=200, random_state=0,
)
start_time = time()
histogram_gradient_boosting.fit(X_train, y_train)
fit_time_histogram_gradient_boosting = time() - start_time

print(
    f"The performance of histogram gradient-boosting are:"
    f"{histogram_gradient_boosting.score(X_test, y_test):.3f}"
)
print(f"Fitting time took: {fit_time_histogram_gradient_boosting:.2f} seconds")<jupyter_output>The performance of histogram gradient-boosting are:0.845
Fitting time took: 2.78 seconds
{'fraction_non_alphanumeric': 0.08272654281680002, 'fraction_numerical': 0.012929483097538691, 'mean_word_length': 5.315789473684211, 'pattern_counts': {'":': 8, '<': 122, '<?xml version=': 0, '>': 122, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '33', 'dir': 'jupyter-structured-clean-dedup', 'id': '563614', 'license': 'CC-BY-4.0', 'n_tokens_mistral': 11366, 'n_tokens_neox': 10721, 'n_words': 3938, 'path': 'rendered_notebooks/ensemble.ipynb', 'repo_name': 'jeremiedbb/scikit-learn-mooc', 'size': '36119'}