<jupyter_start><jupyter_text>Gradient boosting with mixed estimators Why only use decision trees?  imports<jupyter_code>from itertools import cycle, islice
from sklearn.base import BaseEstimator, RegressorMixin, clone
from sklearn.metrics import mean_squared_error
from sklearn.tree import DecisionTreeRegressor
from sklearn.linear_model import LinearRegression
from sklearn.neighbors import KNeighborsRegressor
from sklearn.datasets import make_regression
from sklearn.datasets import fetch_california_housing
from sklearn.model_selection import train_test_split, cross_val_score
from sklearn.ensemble._gb_losses import LeastSquaresError
from sklearn.ensemble import GradientBoostingRegressor, RandomForestRegressor
from sklearn.dummy import DummyRegressor
from sklearn.ensemble import BaggingRegressor
import numpy as np<jupyter_output><empty_output><jupyter_text>load data<jupyter_code># fake data
X, y = make_regression(1000, 50, n_informative=40, random_state=0, noise=1)
# or california housing data
X, y = fetch_california_housing(return_X_y=True)
X_train, X_valid, y_train, y_valid = train_test_split(X, y, random_state=0)
X.std(0)
X.min(0)
X.max(0)
y.min(), y.max(), y.mean(), y.std()<jupyter_output><empty_output><jupyter_text>Purposefully no feature engineering  custom code no `set_params` on estimators implemented yet  light weight implementation<jupyter_code>class ArbitraryBoostingRegressor(BaseEstimator, RegressorMixin):
    def __init__(self, estimators, n_estimators=10, learning_rate=0.1):
        self.estimators = estimators
        self.n_estimators = n_estimators
        self.learning_rate = learning_rate
        
        assert self.estimators
        
    def _estimator_cycle(self):
        for estimator in islice(cycle(self.estimators), self.n_estimators):
            yield clone(estimator)
            
    def fit(self, X, y, **fit_params):
        self.loss_ = LeastSquaresError(1)
        self.estimators_ = []
        self.errors_ = []
        y_resid = y

        self.init_ = self.loss_.init_estimator()
        y_resid = y - self.init_.fit(X, y, **fit_params).predict(X)

        for estimator in self._estimator_cycle():
            y_pred = estimator.fit(X, y_resid, **fit_params).predict(X)
            self.estimators_.append(estimator)
            y_resid = y_resid - y_pred * self.learning_rate
            self.errors_.append(mean_squared_error(y_resid, y_pred))
        return self
    
    def predict(self, X):
        y_pred = self.init_.predict(X)
        for estimator in self.estimators_:
            y_pred += self.learning_rate * estimator.predict(X)
        return y_pred<jupyter_output><empty_output><jupyter_text>inheriting from GradientBoostingRegressor<jupyter_code>class TreeWrapper:
    def __init__(self, obj):
        self.obj = obj
        
    @property
    def tree_(self):
        return self.obj
    
    def __getattr__(self, attr):
        return getattr(self.obj, attr)


class ArbitraryBoostingRegressor2(GradientBoostingRegressor):
    def __init__(
        self,
        estimators,
        *,
        loss='ls',
        learning_rate=0.1,
        n_estimators=10,
        subsample=1.0,
        criterion='friedman_mse',
        min_samples_split=2,
        min_samples_leaf=1,
        min_weight_fraction_leaf=0.0,
        max_depth=3,
        min_impurity_decrease=0.0,
        min_impurity_split=None,
        init=None,
        random_state=None,
        max_features=None,
        alpha=0.9,
        verbose=0,
        max_leaf_nodes=None,
        warm_start=False,
        presort='deprecated',
        validation_fraction=0.1,
        n_iter_no_change=None,
        tol=0.0001,
        ccp_alpha=0.0,
    ):
        super().__init__(
            loss=loss,
            learning_rate=learning_rate,
            n_estimators=n_estimators,
            subsample=subsample,
            criterion=criterion,
            min_samples_split=min_samples_split,
            min_samples_leaf=min_samples_leaf,
            min_weight_fraction_leaf=min_weight_fraction_leaf,
            max_depth=max_depth,
            min_impurity_decrease=min_impurity_decrease,
            min_impurity_split=min_impurity_split,
            init=init,
            random_state=random_state,
            max_features=max_features,
            alpha=alpha,
            verbose=verbose,
            max_leaf_nodes=max_leaf_nodes,
            warm_start=warm_start,
            presort=presort,
            validation_fraction=validation_fraction,
            n_iter_no_change=n_iter_no_change,
            tol=tol,
            ccp_alpha=ccp_alpha,
        )
        self.estimators = estimators

    def _fit_stage(
        self, i, X, y, raw_predictions, sample_weight, sample_mask,
        random_state, X_idx_sorted, X_csc=None, X_csr=None):
        """Fit another stage of ``n_classes_`` trees to the boosting model. """

        assert sample_mask.dtype == np.bool
        loss = self.loss_
        original_y = y

        # Need to pass a copy of raw_predictions to negative_gradient()
        # because raw_predictions is partially updated at the end of the loop
        # in update_terminal_regions(), and gradients need to be evaluated at
        # iteration i - 1.
        raw_predictions_copy = raw_predictions.copy()

        for k in range(loss.K):
            if loss.is_multi_class:
                y = np.array(original_y == k, dtype=np.float64)

            residual = loss.negative_gradient(y, raw_predictions_copy, k=k,
                                              sample_weight=sample_weight)

            # induce regression tree on residuals
            tree = TreeWrapper(clone(next(islice(cycle(self.estimators), i, i+1))))

            if self.subsample < 1.0:
                # no inplace multiplication!
                sample_weight = sample_weight * sample_mask.astype(np.float64)

            X = X_csr if X_csr is not None else X
            tree.fit(X, residual)  # remaining arguments not supported

            # update tree leaves
            loss.update_terminal_regions(
                tree, X, y, residual, raw_predictions, sample_weight,
                sample_mask, learning_rate=self.learning_rate, k=k)

            # add tree to ensemble
            self.estimators_[i, k] = tree

        return raw_predictions

    def _raw_predict_init(self, X):
        """Check input and compute raw predictions of the init estimator."""
        self._check_initialized()
        if X.shape[1] != self.n_features_:
            raise ValueError("X.shape[1] should be {0:d}, not {1:d}.".format(
                self.n_features_, X.shape[1]))
        if self.init_ == 'zero':
            raw_predictions = np.zeros(shape=(X.shape[0], self.loss_.K),
                                       dtype=np.float64)
        else:
            raw_predictions = self.loss_.get_init_raw_predictions(
                X, self.init_).astype(np.float64)
        return raw_predictions

    def _raw_predict(self, X):
        """Return the sum of the trees raw predictions (+ init estimator)."""
        raw_predictions = self._raw_predict_init(X)
        for i in range(self.n_estimators):
            regr = self.estimators_[i][0]
            raw_predictions += self.learning_rate * regr.predict(X).reshape(-1, 1)
        return raw_predictions<jupyter_output><empty_output><jupyter_text>experiments  fit a single model  linear regression<jupyter_code>abr = ArbitraryBoostingRegressor([LinearRegression()])
abr.fit(X_train, y_train)
mean_squared_error(y_valid, abr.predict(X_valid))
abr.errors_<jupyter_output><empty_output><jupyter_text>linear + dt + knn<jupyter_code>regr = ArbitraryBoostingRegressor([LinearRegression(), DecisionTreeRegressor(), KNeighborsRegressor()])
abr.fit(X_train, y_train)
mean_squared_error(y_valid, abr.predict(X_valid))
abr.errors_<jupyter_output><empty_output><jupyter_text>cross validate approach compare the 2 implementations  dt<jupyter_code>%%time
regr = ArbitraryBoostingRegressor([DecisionTreeRegressor()])
scores = cross_val_score(regr, X, y, cv=5, scoring='neg_mean_squared_error')
print(-np.mean(scores))
print(scores)
%%time
regr = ArbitraryBoostingRegressor2([DecisionTreeRegressor()])
scores = cross_val_score(regr, X, y, cv=5, scoring='neg_mean_squared_error')
print(-np.mean(scores))
print(scores)<jupyter_output>0.6722669141379157
[-0.58848153 -0.59639745 -0.70226855 -0.71117485 -0.76301219]
CPU times: user 8.53 s, sys: 43.8 ms, total: 8.57 s
Wall time: 8.57 s<jupyter_text>use default GradientBoostingRegressor hyper params<jupyter_code>%%time
dt = DecisionTreeRegressor(
    criterion='friedman_mse',
    min_samples_split=2,
    min_samples_leaf=1,
    min_weight_fraction_leaf=0.0,
    max_depth=3,
    min_impurity_decrease=0.0,
    min_impurity_split=None,
)
regr = ArbitraryBoostingRegressor2([dt])
scores = cross_val_score(regr, X, y, cv=5, scoring='neg_mean_squared_error')
print(-np.mean(scores))
print(scores)<jupyter_output>0.7011773270758519
[-0.62248244 -0.63549276 -0.74758768 -0.78483192 -0.71549184]
CPU times: user 1.74 s, sys: 40 µs, total: 1.74 s
Wall time: 1.74 s<jupyter_text>lr<jupyter_code>%%time
regr = ArbitraryBoostingRegressor([LinearRegression()])
scores = cross_val_score(regr, X, y, cv=5, scoring='neg_mean_squared_error')
print(-np.mean(scores))
print(scores)
%%time
regr = ArbitraryBoostingRegressor2([LinearRegression()])
scores = cross_val_score(regr, X, y, cv=5, scoring='neg_mean_squared_error')
print(-np.mean(scores))
print(scores)<jupyter_output>0.6535163904654909
[-0.5680934  -0.69759844 -0.75928304 -0.65383644 -0.58877063]
CPU times: user 536 ms, sys: 7.73 ms, total: 544 ms
Wall time: 137 ms<jupyter_text>knn<jupyter_code>%%time
regr = ArbitraryBoostingRegressor([KNeighborsRegressor()])
scores = cross_val_score(regr, X, y, cv=5, scoring='neg_mean_squared_error')
print(-np.mean(scores))
print(scores)
%%time
regr = ArbitraryBoostingRegressor2([KNeighborsRegressor()])
scores = cross_val_score(regr, X, y, cv=5, scoring='neg_mean_squared_error')
print(-np.mean(scores))
print(scores)<jupyter_output>1.2037599868343167
[-1.09356936 -1.03899702 -1.35285596 -1.18142013 -1.35195747]
CPU times: user 5.4 s, sys: 32 ms, total: 5.43 s
Wall time: 5.43 s<jupyter_text>lr + dt<jupyter_code>%%time
regr = ArbitraryBoostingRegressor([LinearRegression(), dt])
scores = cross_val_score(regr, X, y, cv=5, scoring='neg_mean_squared_error')
print(-np.mean(scores))
print(scores)
%%time
regr = ArbitraryBoostingRegressor2([LinearRegression(), dt])
scores = cross_val_score(regr, X, y, cv=5, scoring='neg_mean_squared_error')
print(-np.mean(scores))
print(scores)<jupyter_output>0.6334443270638227
[-0.55734265 -0.63247685 -0.70368564 -0.67486366 -0.59885283]
CPU times: user 3.79 s, sys: 48 ms, total: 3.84 s
Wall time: 960 ms<jupyter_text>lr + dt + knn<jupyter_code>%%time
regr = ArbitraryBoostingRegressor([LinearRegression(), dt, KNeighborsRegressor()])
scores = cross_val_score(regr, X, y, cv=5, scoring='neg_mean_squared_error')
print(-np.mean(scores))
print(scores)
%%time
regr = ArbitraryBoostingRegressor2([LinearRegression(), dt, KNeighborsRegressor()])
scores = cross_val_score(regr, X, y, cv=5, scoring='neg_mean_squared_error')
print(-np.mean(scores))
print(scores)<jupyter_output>0.7240273214247799
[-0.63529467 -0.70157353 -0.81329003 -0.7485971  -0.72138127]
CPU times: user 9.85 s, sys: 132 ms, total: 9.98 s
Wall time: 2.49 s<jupyter_text>benchmarks  dummy<jupyter_code>%%time
regr = DummyRegressor()
scores = cross_val_score(regr, X, y, cv=5, scoring='neg_mean_squared_error')
print(-np.mean(scores))
print(scores)<jupyter_output>1.3702550894247665
[-1.3064619  -1.19893373 -1.57273505 -1.25865506 -1.5144897 ]
CPU times: user 42.7 ms, sys: 34 µs, total: 42.8 ms
Wall time: 10.9 ms<jupyter_text>vanilla gradient boosting<jupyter_code>%%time
regr = GradientBoostingRegressor(n_estimators=10)
scores = cross_val_score(regr, X, y, cv=5, scoring='neg_mean_squared_error')
print(-np.mean(scores))
print(scores)<jupyter_output>0.7011773270758519
[-0.62248244 -0.63549276 -0.74758768 -0.78483192 -0.71549184]
CPU times: user 2.26 s, sys: 7.96 ms, total: 2.26 s
Wall time: 1.75 s<jupyter_text>random forest<jupyter_code>%%time
regr = RandomForestRegressor(n_estimators=10)
scores = cross_val_score(regr, X, y, cv=5, scoring='neg_mean_squared_error')
print(-np.mean(scores))
print(scores)<jupyter_output>0.4847291292478402
[-0.62378846 -0.38590051 -0.41904684 -0.48970468 -0.50520516]
CPU times: user 5.16 s, sys: 12 ms, total: 5.17 s
Wall time: 5.17 s<jupyter_text>bagging + lr<jupyter_code>%%time
regr = BaggingRegressor(LinearRegression(), n_estimators=10)
scores = cross_val_score(regr, X, y, cv=5, scoring='neg_mean_squared_error')
print(-np.mean(scores))
print(scores)<jupyter_output>0.5587854379934888
[-0.49734083 -0.62139618 -0.63548575 -0.54329421 -0.49641022]
CPU times: user 962 ms, sys: 8.06 ms, total: 970 ms
Wall time: 249 ms<jupyter_text>bagging + dt should roughly equal rf<jupyter_code>%%time
regr = BaggingRegressor(DecisionTreeRegressor(), n_estimators=10)
scores = cross_val_score(regr, X, y, cv=5, scoring='neg_mean_squared_error')
print(-np.mean(scores))
print(scores)<jupyter_output>0.46562762669283
[-0.51150194 -0.37253359 -0.42157403 -0.51776368 -0.5047649 ]
CPU times: user 5.65 s, sys: 56 ms, total: 5.71 s
Wall time: 5.13 s
{'fraction_non_alphanumeric': 0.109107303877367, 'fraction_numerical': 0.06792906522392546, 'mean_word_length': 3.778815080789946, 'pattern_counts': {'":': 0, '<': 59, '<?xml version=': 0, '>': 58, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '17', 'dir': 'jupyter-structured-clean-dedup', 'id': '645902', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 5382, 'n_tokens_neox': 4665, 'n_words': 912, 'path': 'gradient_boosting_mixed_estimators.ipynb', 'repo_name': 'BenjaminBossan/playground', 'size': '13308'}