XGBoost: Implementation and Final Comparison

import xgboost as xgb
from xgboost import XGBRegressor, XGBClassifier
from sklearn.datasets import fetch_california_housing, load_breast_cancer
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_squared_error, r2_score, roc_auc_score
import numpy as np

# Regression
ch = fetch_california_housing()
X_r, y_r = ch.data, ch.target
Xr_tr, Xr_te, yr_tr, yr_te = train_test_split(X_r, y_r, test_size=0.2, random_state=42)

# Classification
bc = load_breast_cancer()
X_c, y_c = bc.data, bc.target
Xc_tr, Xc_te, yc_tr, yc_te = train_test_split(X_c, y_c, test_size=0.2, random_state=42, stratify=y_c)

xgb_reg = XGBRegressor(
    n_estimators=500,
    learning_rate=0.05,
    max_depth=4,
    subsample=0.8,
    colsample_bytree=0.8,
    reg_alpha=0.1,       # L1 regularization
    reg_lambda=1.0,      # L2 regularization
    min_child_weight=5,
    random_state=42,
    n_jobs=-1,
    tree_method='hist'   # histogram-based split finding
)
xgb_reg.fit(Xr_tr, yr_tr, eval_set=[(Xr_te, yr_te)], verbose=False)
y_pred_r = xgb_reg.predict(Xr_te)
print(f"XGB Regressor: RMSE={np.sqrt(mean_squared_error(yr_te, y_pred_r)):.4f}, R²={r2_score(yr_te, y_pred_r):.4f}")

xgb_es = XGBRegressor(
    n_estimators=2000,           # upper bound — won't reach this
    learning_rate=0.05,
    max_depth=4,
    subsample=0.8,
    colsample_bytree=0.8,
    reg_alpha=0.1,
    reg_lambda=1.0,
    random_state=42,
    n_jobs=-1,
    tree_method='hist',
    early_stopping_rounds=50     # stop after 50 rounds without improvement
)
xgb_es.fit(
    Xr_tr, yr_tr,
    eval_set=[(Xr_te, yr_te)],
    verbose=100
)
print(f"Best iteration: {xgb_es.best_iteration}")
print(f"Best RMSE:      {xgb_es.best_score:.4f}")

dtrain = xgb.DMatrix(Xr_tr, label=yr_tr, feature_names=ch.feature_names.tolist())
dtest  = xgb.DMatrix(Xr_te, label=yr_te, feature_names=ch.feature_names.tolist())

params = {
    'objective':        'reg:squarederror',
    'learning_rate':    0.05,
    'max_depth':        4,
    'subsample':        0.8,
    'colsample_bytree': 0.8,
    'reg_alpha':        0.1,
    'reg_lambda':       1.0,
    'min_child_weight': 5,
    'seed':             42,
    'tree_method':      'hist'
}

evals_result = {}
model_native = xgb.train(
    params,
    dtrain,
    num_boost_round=500,
    evals=[(dtrain, 'train'), (dtest, 'validation')],
    early_stopping_rounds=50,
    evals_result=evals_result,
    verbose_eval=False
)
print(f"Best iteration: {model_native.best_iteration}")

print(f"{'colsample_bytree':>18} | {'RMSE':>8}")
for col_tree in [0.5, 0.7, 0.8, 1.0]:
    xgb_c = XGBRegressor(n_estimators=300, learning_rate=0.1, max_depth=4,
                          colsample_bytree=col_tree, random_state=42, n_jobs=-1)
    xgb_c.fit(Xr_tr, yr_tr)
    rmse = np.sqrt(mean_squared_error(yr_te, xgb_c.predict(Xr_te)))
    print(f"{col_tree:>18} | {rmse:>8.4f}")

from sklearn.model_selection import RandomizedSearchCV
from scipy.stats import uniform

param_dist = {
    'n_estimators':     [200, 300, 500],
    'max_depth':        [3, 4, 5, 6],
    'learning_rate':    uniform(0.01, 0.2),
    'subsample':        uniform(0.6, 0.4),
    'colsample_bytree': uniform(0.6, 0.4),
    'min_child_weight': [1, 3, 5],
    'reg_alpha':        [0, 0.1, 0.5],
}
rs = RandomizedSearchCV(
    XGBRegressor(tree_method='hist', random_state=42, n_jobs=-1),
    param_dist, n_iter=30, cv=5,
    scoring='neg_root_mean_squared_error', random_state=42
)
rs.fit(Xr_tr, yr_tr)
print(f"Best params: {rs.best_params_}")
print(f"Best CV RMSE: {-rs.best_score_:.4f}")
print(f"Test RMSE:    {np.sqrt(mean_squared_error(yr_te, rs.best_estimator_.predict(Xr_te))):.4f}")

xgb_clf = XGBClassifier(
    n_estimators=300,
    learning_rate=0.1,
    max_depth=4,
    subsample=0.8,
    colsample_bytree=0.8,
    eval_metric='logloss',
    random_state=42,
    n_jobs=-1
)
xgb_clf.fit(Xc_tr, yc_tr, eval_set=[(Xc_te, yc_te)], verbose=False)
y_prob_c = xgb_clf.predict_proba(Xc_te)[:, 1]

from sklearn.metrics import classification_report, confusion_matrix
print(f"Test Accuracy: {xgb_clf.score(Xc_te, yc_te):.4f}")
print(f"AUC-ROC:       {roc_auc_score(yc_te, y_prob_c):.4f}")
print(confusion_matrix(yc_te, xgb_clf.predict(Xc_te)))

# Gain = average reduction in loss per split using this feature (most informative)
imp_gain   = xgb_reg.get_booster().get_score(importance_type='gain')
# Weight = number of times a feature appears in any split
imp_weight = xgb_reg.get_booster().get_score(importance_type='weight')
# Cover = average number of samples affected per split
imp_cover  = xgb_reg.get_booster().get_score(importance_type='cover')

import pandas as pd
feat_imp = pd.DataFrame({
    'Feature': list(imp_gain.keys()),
    'Gain':    list(imp_gain.values())
}).sort_values('Gain', ascending=False)
print(feat_imp.head(5))

from sklearn.ensemble import RandomForestRegressor, GradientBoostingRegressor, AdaBoostRegressor
from sklearn.tree import DecisionTreeRegressor

methods = [
    ('Decision Tree',    DecisionTreeRegressor(max_depth=5, random_state=42)),
    ('Random Forest',    RandomForestRegressor(n_estimators=200, random_state=42, n_jobs=-1)),
    ('AdaBoost',         AdaBoostRegressor(n_estimators=200, random_state=42)),
    ('Gradient Boosting', GradientBoostingRegressor(n_estimators=200, learning_rate=0.1, max_depth=3, random_state=42)),
    ('XGBoost',          XGBRegressor(n_estimators=500, learning_rate=0.05, max_depth=4,
                                       random_state=42, n_jobs=-1, tree_method='hist')),
]

print(f"{'Method':25s} | {'Train RMSE':>10} | {'Test RMSE':>10} | {'R²':>8}")
print("-" * 65)
for name, model in methods:
    model.fit(Xr_tr, yr_tr)
    tr_rmse = np.sqrt(mean_squared_error(yr_tr, model.predict(Xr_tr)))
    te_rmse = np.sqrt(mean_squared_error(yr_te, model.predict(Xr_te)))
    r2      = r2_score(yr_te, model.predict(Xr_te))
    print(f"{name:25s} | {tr_rmse:>10.4f} | {te_rmse:>10.4f} | {r2:>8.4f}")