Advanced Regression Models

depths = [2, 5, 10]
for depth in depths:
    dt = DecisionTreeRegressor(max_depth=depth, random_state=42)
    dt.fit(X_train, y_train)
    
    train_r2 = r2_score(y_train, dt.predict(X_train))
    test_r2 = r2_score(y_test, dt.predict(X_test))
    
    print(f"Depth {depth:2d}: Train R2={train_r2:.4f}, Test R2={test_r2:.4f}")

Explanation: As depth increases, training R2 improves (approaching 1.0), but test R2 may decrease after a point. This is the classic overfitting pattern where the model memorizes training data but fails to generalize.

import matplotlib.pyplot as plt

alphas = np.linspace(0.0, 0.1, 20)
train_scores = []
test_scores = []

for alpha in alphas:
    dt = DecisionTreeRegressor(ccp_alpha=alpha, random_state=42)
    dt.fit(X_train, y_train)
    train_scores.append(r2_score(y_train, dt.predict(X_train)))
    test_scores.append(r2_score(y_test, dt.predict(X_test)))

plt.figure(figsize=(10, 6))
plt.plot(alphas, train_scores, 'b-', label='Training R2')
plt.plot(alphas, test_scores, 'r-', label='Test R2')
plt.xlabel('ccp_alpha')
plt.ylabel('R2 Score')
plt.legend()
plt.title('Cost Complexity Pruning')
plt.show()

Explanation: Cost complexity pruning adds a penalty term for tree complexity. As alpha increases, the tree becomes simpler. The optimal alpha is where test score is maximized before it starts declining significantly.

dt = DecisionTreeRegressor(max_depth=5, random_state=42)
dt.fit(X_train, y_train)

# Get feature importances
importances = dt.feature_importances_
feature_names = [f'Feature_{i}' for i in range(X.shape[1])]

# Sort by importance
sorted_indices = np.argsort(importances)[::-1]
for idx in sorted_indices:
    print(f"{feature_names[idx]}: {importances[idx]:.4f}")

Explanation: Feature importance in decision trees is based on how much each feature reduces impurity (variance for regression). Features used higher in the tree and in more splits have higher importance.

min_samples_values = [1, 5, 10, 20, 50]

print("min_samples_leaf | Train R2  | Test R2   | Gap")
print("-" * 55)
for min_samples in min_samples_values:
    dt = DecisionTreeRegressor(min_samples_leaf=min_samples, random_state=42)
    dt.fit(X_train, y_train)
    train_r2 = r2_score(y_train, dt.predict(X_train))
    test_r2 = r2_score(y_test, dt.predict(X_test))
    gap = train_r2 - test_r2
    print(f"{min_samples:16d} | {train_r2:.4f}   | {test_r2:.4f}   | {gap:.4f}")

Explanation: Higher min_samples_leaf forces larger leaf nodes, creating a smoother model. The gap between train and test R2 decreases as overfitting is reduced. Find the sweet spot where test R2 is maximized.

from sklearn.model_selection import cross_val_score
import matplotlib.pyplot as plt

depths = range(1, 21)
cv_means = []
cv_stds = []

for depth in depths:
    dt = DecisionTreeRegressor(max_depth=depth, random_state=42)
    scores = cross_val_score(dt, X_train, y_train, cv=5, scoring='r2')
    cv_means.append(scores.mean())
    cv_stds.append(scores.std())

# Plot results
plt.figure(figsize=(10, 6))
plt.errorbar(depths, cv_means, yerr=cv_stds, fmt='-o', capsize=3)
plt.xlabel('max_depth')
plt.ylabel('CV R2 Score')
plt.title('Finding Optimal Tree Depth with Cross-Validation')
plt.axhline(y=max(cv_means), color='r', linestyle='--', alpha=0.5)
best_depth = depths[np.argmax(cv_means)]
plt.axvline(x=best_depth, color='g', linestyle='--', alpha=0.5)
print(f"Optimal depth: {best_depth}")

Explanation: Cross-validation gives an unbiased estimate of test performance. The optimal depth is where CV score peaks. Beyond this, the tree overfits, and CV score decreases.

dt = DecisionTreeRegressor(max_depth=8, random_state=42)
dt.fit(X_train, y_train)
y_pred = dt.predict(X_test)

plt.figure(figsize=(8, 8))
plt.scatter(y_test, y_pred, alpha=0.5, edgecolors='none')

# Perfect prediction line
min_val = min(y_test.min(), y_pred.min())
max_val = max(y_test.max(), y_pred.max())
plt.plot([min_val, max_val], [min_val, max_val], 'r--', lw=2, label='Perfect prediction')

plt.xlabel('Actual Values')
plt.ylabel('Predicted Values')
plt.title(f'Decision Tree: Predicted vs Actual (R2={r2_score(y_test, y_pred):.4f})')
plt.legend()
plt.axis('equal')
plt.show()

# Calculate correlation
correlation = np.corrcoef(y_test, y_pred)[0, 1]
print(f"Pearson correlation: {correlation:.4f}")

Explanation: Points close to the diagonal indicate accurate predictions. Systematic deviations suggest model bias. Scatter around the line indicates variance. This visualization helps diagnose model performance.

dt = DecisionTreeRegressor(max_depth=8, random_state=42)
dt.fit(X_train, y_train)
y_pred = dt.predict(X_test)
residuals = y_test - y_pred

fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# Residual plot
axes[0].scatter(y_pred, residuals, alpha=0.5)
axes[0].axhline(y=0, color='r', linestyle='--')
axes[0].set_xlabel('Predicted Values')
axes[0].set_ylabel('Residuals')
axes[0].set_title('Residual Plot')

# Histogram of residuals
axes[1].hist(residuals, bins=30, edgecolor='black')
axes[1].axvline(x=0, color='r', linestyle='--')
axes[1].set_xlabel('Residual Value')
axes[1].set_ylabel('Frequency')
axes[1].set_title(f'Residual Distribution (Mean: {residuals.mean():.2f})')

plt.tight_layout()
plt.show()

Explanation: Good models have residuals randomly scattered around zero with no pattern. Patterns in the residual plot suggest the model is missing systematic relationships in the data.

max_features_options = ['sqrt', 'log2', 1.0]

for mf in max_features_options:
    rf_temp = RandomForestRegressor(n_estimators=100, max_features=mf, random_state=42)
    rf_temp.fit(X_train, y_train)
    r2 = r2_score(y_test, rf_temp.predict(X_test))
    print(f"max_features={mf}: R2={r2:.4f}")

Explanation: For regression, using all features (1.0) often works well. 'sqrt' and 'log2' increase diversity between trees but may reduce individual tree accuracy. The best choice depends on your specific data and feature set.

from sklearn.inspection import PartialDependenceDisplay

# Find most important feature
most_important_idx = np.argmax(rf.feature_importances_)

# Create partial dependence plot
fig, ax = plt.subplots(figsize=(8, 6))
PartialDependenceDisplay.from_estimator(
    rf, X_train, [most_important_idx],
    ax=ax, kind='average'
)
plt.title(f'Partial Dependence for Feature_{most_important_idx}')
plt.tight_layout()
plt.show()

Explanation: Partial dependence plots show how the predicted value changes as we vary one feature while keeping others at their average values. This reveals the marginal effect of that feature on predictions.

from sklearn.tree import DecisionTreeRegressor
from sklearn.ensemble import RandomForestRegressor

# Single tree
dt = DecisionTreeRegressor(max_depth=10, random_state=42)
dt.fit(X_train, y_train)
dt_pred = dt.predict(X_test)

# Random Forest
rf = RandomForestRegressor(n_estimators=100, max_depth=10, random_state=42)
rf.fit(X_train, y_train)
rf_pred = rf.predict(X_test)

print(f"Single Tree R2: {r2_score(y_test, dt_pred):.4f}")
print(f"Random Forest R2: {r2_score(y_test, rf_pred):.4f}")
print(f"\nPrediction variance (lower is more stable):")
print(f"Single Tree: {np.var(dt_pred - y_test):.4f}")
print(f"Random Forest: {np.var(rf_pred - y_test):.4f}")

Explanation: Random Forest typically achieves higher R2 and lower prediction variance. The ensemble averaging reduces the noise that single trees are prone to.

from sklearn.model_selection import cross_val_score

# OOB score
rf_oob = RandomForestRegressor(n_estimators=100, oob_score=True, random_state=42)
rf_oob.fit(X_train, y_train)
oob_score = rf_oob.oob_score_

# Cross-validation score
rf_cv = RandomForestRegressor(n_estimators=100, random_state=42)
cv_scores = cross_val_score(rf_cv, X_train, y_train, cv=5, scoring='r2')

print(f"OOB Score: {oob_score:.4f}")
print(f"CV Score: {cv_scores.mean():.4f} (+/- {cv_scores.std():.4f}")
print(f"\nDifference: {abs(oob_score - cv_scores.mean()):.4f}")

Explanation: OOB and CV scores are usually similar. OOB is "free" (computed during training) while CV requires extra training runs. Use OOB for quick estimates; use CV for more rigorous validation or when comparing with non-RF models.

rf = RandomForestRegressor(n_estimators=100, random_state=42)
rf.fit(X_train, y_train)

# Get predictions from all trees
all_tree_preds = np.array([tree.predict(X_test) for tree in rf.estimators_])

# Calculate percentiles for prediction intervals
lower = np.percentile(all_tree_preds, 5, axis=0)
upper = np.percentile(all_tree_preds, 95, axis=0)
mean_pred = rf.predict(X_test)

# Plot for first 20 samples
plt.figure(figsize=(12, 6))
x_range = range(20)
plt.errorbar(x_range, mean_pred[:20], 
             yerr=[mean_pred[:20]-lower[:20], upper[:20]-mean_pred[:20]],
             fmt='o', capsize=3, label='Prediction ± 90% interval')
plt.scatter(x_range, y_test[:20], color='red', marker='x', s=100, label='Actual')
plt.xlabel('Sample Index')
plt.ylabel('Value')
plt.legend()
plt.title('Random Forest Predictions with Uncertainty Intervals')
plt.show()

Explanation: The spread of individual tree predictions gives us an estimate of model uncertainty. Wider intervals indicate less confident predictions. This is useful for risk-aware decision making.

from sklearn.inspection import permutation_importance

rf = RandomForestRegressor(n_estimators=100, random_state=42)
rf.fit(X_train, y_train)

# Built-in importance
builtin_importance = rf.feature_importances_

# Permutation importance (on test set)
perm_importance = permutation_importance(rf, X_test, y_test, n_repeats=10, random_state=42)

# Compare
feature_names = [f'Feature_{i}' for i in range(X.shape[1])]
comparison = pd.DataFrame({
    'Feature': feature_names,
    'Built-in': builtin_importance,
    'Permutation': perm_importance.importances_mean
}).sort_values('Permutation', ascending=False)

print(comparison.to_string(index=False))

Explanation: Permutation importance measures how much performance drops when a feature's values are shuffled. It's computed on test data and doesn't suffer from the bias toward high-cardinality features that affects built-in importance.

# Using sklearn's HalvingRandomSearchCV as an efficient alternative
from sklearn.experimental import enable_halving_search_cv
from sklearn.model_selection import HalvingRandomSearchCV
from scipy.stats import randint, uniform

param_dist = {
    'max_depth': randint(5, 30),
    'min_samples_leaf': randint(1, 20),
    'max_features': uniform(0.3, 0.7)
}

rf = RandomForestRegressor(n_estimators=100, random_state=42, n_jobs=-1)

search = HalvingRandomSearchCV(
    rf, param_dist, n_candidates=50, cv=3, 
    scoring='r2', random_state=42, factor=2
)
search.fit(X_train, y_train)

print(f"Best parameters: {search.best_params_}")
print(f"Best CV R2: {search.best_score_:.4f}")
print(f"Test R2: {r2_score(y_test, search.predict(X_test)):.4f}")

Explanation: HalvingRandomSearchCV progressively eliminates poor candidates using increasing amounts of data, making it more efficient than standard random search. For true Bayesian optimization, install scikit-optimize or optuna.

learning_rates = [0.01, 0.1, 0.3]

for lr in learning_rates:
    model = xgb.XGBRegressor(n_estimators=100, learning_rate=lr, random_state=42)
    model.fit(X_train, y_train)
    r2 = r2_score(y_test, model.predict(X_test))
    print(f"Learning rate {lr}: R2={r2:.4f}")

Explanation: Lower learning rates require more trees but often achieve better results. With only 100 trees, a learning rate of 0.01 may underfit while 0.3 may overfit. The optimal depends on how many trees you are willing to train.

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

param_dist = {
    'learning_rate': uniform(0.01, 0.29),
    'max_depth': randint(2, 10),
    'reg_alpha': uniform(0, 1)
}

xgb_search = xgb.XGBRegressor(n_estimators=100, random_state=42)
random_search = RandomizedSearchCV(
    xgb_search, param_dist, n_iter=20, cv=5, scoring='r2', random_state=42
)
random_search.fit(X_train, y_train)

print(f"Best params: {random_search.best_params_}")
print(f"Best CV R2: {random_search.best_score_:.4f}")

Explanation: RandomizedSearchCV samples random combinations from the parameter distributions, which is more efficient than GridSearch for many hyperparameters. 20 iterations with 5-fold CV gives 100 model fits.

from sklearn.ensemble import GradientBoostingRegressor
import xgboost as xgb
import time

# sklearn Gradient Boosting
start = time.time()
gb_sklearn = GradientBoostingRegressor(n_estimators=100, max_depth=3, random_state=42)
gb_sklearn.fit(X_train, y_train)
sklearn_time = time.time() - start
sklearn_r2 = r2_score(y_test, gb_sklearn.predict(X_test))

# XGBoost
start = time.time()
gb_xgb = xgb.XGBRegressor(n_estimators=100, max_depth=3, random_state=42)
gb_xgb.fit(X_train, y_train)
xgb_time = time.time() - start
xgb_r2 = r2_score(y_test, gb_xgb.predict(X_test))

print(f"sklearn GB: R2={sklearn_r2:.4f}, Time={sklearn_time:.2f}s")
print(f"XGBoost:    R2={xgb_r2:.4f}, Time={xgb_time:.2f}s")

Explanation: XGBoost is typically faster due to optimized C++ implementation and parallel processing. Performance is usually similar, but XGBoost provides more regularization options and features like early stopping.

import matplotlib.pyplot as plt

# Train with evaluation
model = xgb.XGBRegressor(n_estimators=200, learning_rate=0.1, random_state=42)
eval_set = [(X_train, y_train), (X_test, y_test)]
model.fit(X_train, y_train, eval_set=eval_set, verbose=False)

# Get evaluation results
results = model.evals_result()

# Plot
plt.figure(figsize=(10, 6))
plt.plot(results['validation_0']['rmse'], label='Training RMSE')
plt.plot(results['validation_1']['rmse'], label='Validation RMSE')
plt.xlabel('Boosting Iteration')
plt.ylabel('RMSE')
plt.legend()
plt.title('XGBoost: Training Progress')
plt.show()

# Find best iteration
best_iter = np.argmin(results['validation_1']['rmse'])
print(f"Best iteration: {best_iter}")

Explanation: Training RMSE continuously decreases. Validation RMSE decreases initially then may increase (overfitting). The gap between them indicates overfitting severity. Early stopping would stop at the minimum validation RMSE.

regularization_settings = [
    {'reg_alpha': 0, 'reg_lambda': 0},      # No regularization
    {'reg_alpha': 1, 'reg_lambda': 0},      # L1 only
    {'reg_alpha': 0, 'reg_lambda': 1},      # L2 only
    {'reg_alpha': 0.5, 'reg_lambda': 0.5},  # Both
    {'reg_alpha': 2, 'reg_lambda': 2},      # Strong both
]

print("Regularization | Train R2 | Test R2 | Gap")
print("-" * 50)
for reg in regularization_settings:
    model = xgb.XGBRegressor(n_estimators=100, **reg, random_state=42)
    model.fit(X_train, y_train)
    train_r2 = r2_score(y_train, model.predict(X_train))
    test_r2 = r2_score(y_test, model.predict(X_test))
    print(f"α={reg['reg_alpha']:.1f}, λ={reg['reg_lambda']:.1f} | {train_r2:.4f}  | {test_r2:.4f} | {train_r2-test_r2:.4f}")

Explanation: Regularization reduces overfitting (smaller gap between train and test). L1 (alpha) encourages sparse features, L2 (lambda) keeps all weights small. The best setting depends on your data and feature structure.

import lightgbm as lgb
import xgboost as xgb
import time

# XGBoost
start = time.time()
xgb_model = xgb.XGBRegressor(n_estimators=100, max_depth=6, learning_rate=0.1, random_state=42)
xgb_model.fit(X_train, y_train)
xgb_time = time.time() - start
xgb_r2 = r2_score(y_test, xgb_model.predict(X_test))

# LightGBM
start = time.time()
lgb_model = lgb.LGBMRegressor(n_estimators=100, max_depth=6, learning_rate=0.1, random_state=42, verbose=-1)
lgb_model.fit(X_train, y_train)
lgb_time = time.time() - start
lgb_r2 = r2_score(y_test, lgb_model.predict(X_test))

print(f"XGBoost:  R2={xgb_r2:.4f}, Time={xgb_time:.3f}s")
print(f"LightGBM: R2={lgb_r2:.4f}, Time={lgb_time:.3f}s")
print(f"LightGBM is {xgb_time/lgb_time:.1f}x faster")

Explanation: LightGBM uses histogram-based algorithms and leaf-wise tree growth, often making it faster than XGBoost, especially on large datasets. Accuracy is typically comparable, but LightGBM may need different hyperparameter tuning.

# pip install shap
import shap

# Train model
model = xgb.XGBRegressor(n_estimators=100, random_state=42)
model.fit(X_train, y_train)

# Create SHAP explainer
explainer = shap.TreeExplainer(model)
shap_values = explainer.shap_values(X_test)

# Summary plot
feature_names = [f'Feature_{i}' for i in range(X.shape[1])]
shap.summary_plot(shap_values, X_test, feature_names=feature_names)

# For a single prediction
shap.waterfall_plot(shap.Explanation(
    values=shap_values[0],
    base_values=explainer.expected_value,
    data=X_test[0],
    feature_names=feature_names
))

Explanation: SHAP values show how each feature contributes to individual predictions. The summary plot shows global feature importance and the direction of effects. Waterfall plots explain single predictions - crucial for model interpretability.

import time
from sklearn.model_selection import train_test_split

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2)

print("Training Time Comparison:")
print("-" * 40)
for name, model in models.items():
    start = time.time()
    model.fit(X_train, y_train)
    duration = time.time() - start
    print(f"{name:20s}: {duration:.4f} seconds")

Explanation: Training time varies significantly between algorithms. SVR is typically slowest for larger datasets, while Random Forest benefits from parallelization. XGBoost is highly optimized and often fast despite being complex.

# 1. Split data (keep test set completely separate)
X_train_full, X_test_final, y_train_full, y_test_final = train_test_split(
    X, y, test_size=0.2, random_state=42
)

# 2. Compare models with CV on training data
best_score = -np.inf
best_model_name = None
for name, model in models.items():
    scores = cross_val_score(model, X_train_full, y_train_full, cv=5, scoring='r2')
    if scores.mean() > best_score:
        best_score = scores.mean()
        best_model_name = name

print(f"Best model: {best_model_name} (CV R2: {best_score:.4f})")

# 3. Fine-tune the best model
from sklearn.model_selection import RandomizedSearchCV

if 'XGBoost' in best_model_name:
    param_dist = {'learning_rate': uniform(0.01, 0.19), 'max_depth': randint(3, 8)}
    tuner = RandomizedSearchCV(
        xgb.XGBRegressor(n_estimators=100), param_dist, n_iter=10, cv=5
    )
    tuner.fit(X_train_full, y_train_full)
    final_model = tuner.best_estimator_
else:
    final_model = models[best_model_name]
    final_model.fit(X_train_full, y_train_full)

# 4. Final evaluation on held-out test set
final_pred = final_model.predict(X_test_final)
print(f"Final Test R2: {r2_score(y_test_final, final_pred):.4f}")

Explanation: This pipeline follows best practices: compare models with CV, tune the best one, and only touch the test set once at the very end. This gives an honest estimate of real-world performance.

import time

# Train all models first
trained_models = {}
for name, model in models.items():
    model.fit(X_train, y_train)
    trained_models[name] = model

# Measure prediction time
print("Prediction Time Comparison (1000 iterations):")
print("-" * 50)
for name, model in trained_models.items():
    start = time.time()
    for _ in range(1000):
        _ = model.predict(X_test)
    duration = time.time() - start
    print(f"{name:20s}: {duration:.4f} seconds")

Explanation: Decision Tree is typically fastest (just if-else statements). SVR can be slow with many support vectors. Random Forest and XGBoost are in between. Prediction speed matters for real-time applications.

from sklearn.model_selection import learning_curve

fig, axes = plt.subplots(1, 2, figsize=(14, 5))

models_to_plot = {
    'Random Forest': RandomForestRegressor(n_estimators=50, random_state=42),
    'XGBoost': xgb.XGBRegressor(n_estimators=50, random_state=42)
}

for ax, (name, model) in zip(axes, models_to_plot.items()):
    train_sizes, train_scores, val_scores = learning_curve(
        model, X, y, cv=5, train_sizes=np.linspace(0.1, 1.0, 10),
        scoring='r2', n_jobs=-1
    )
    
    ax.plot(train_sizes, train_scores.mean(axis=1), 'b-o', label='Training')
    ax.plot(train_sizes, val_scores.mean(axis=1), 'r-s', label='Validation')
    ax.fill_between(train_sizes, train_scores.mean(axis=1) - train_scores.std(axis=1),
                    train_scores.mean(axis=1) + train_scores.std(axis=1), alpha=0.1)
    ax.fill_between(train_sizes, val_scores.mean(axis=1) - val_scores.std(axis=1),
                    val_scores.mean(axis=1) + val_scores.std(axis=1), alpha=0.1)
    ax.set_xlabel('Training Set Size')
    ax.set_ylabel('R2 Score')
    ax.set_title(f'{name} Learning Curve')
    ax.legend()

plt.tight_layout()
plt.show()

Explanation: Learning curves show if more data would help. If validation score is still rising, more data will improve performance. If train and validation converge, the model may need more complexity.

from scipy.stats import ttest_rel
from sklearn.model_selection import cross_val_score

# Get CV scores for both models
rf = RandomForestRegressor(n_estimators=100, random_state=42)
xgb_model = xgb.XGBRegressor(n_estimators=100, random_state=42)

rf_scores = cross_val_score(rf, X, y, cv=10, scoring='r2')
xgb_scores = cross_val_score(xgb_model, X, y, cv=10, scoring='r2')

# Paired t-test
t_stat, p_value = ttest_rel(rf_scores, xgb_scores)

print(f"Random Forest mean R2: {rf_scores.mean():.4f}")
print(f"XGBoost mean R2: {xgb_scores.mean():.4f}")
print(f"\nPaired t-test:")
print(f"t-statistic: {t_stat:.4f}")
print(f"p-value: {p_value:.4f}")
print(f"\nDifference is {'significant' if p_value < 0.05 else 'NOT significant'} (α=0.05)")

Explanation: A paired t-test accounts for the fact that both models are evaluated on the same CV folds. A low p-value (< 0.05) indicates the difference is unlikely due to random chance.

from sklearn.ensemble import StackingRegressor
from sklearn.linear_model import RidgeCV

# Define base learners
base_learners = [
    ('rf', RandomForestRegressor(n_estimators=100, random_state=42)),
    ('xgb', xgb.XGBRegressor(n_estimators=100, random_state=42))
]

# Create stacking ensemble
stacking = StackingRegressor(
    estimators=base_learners,
    final_estimator=RidgeCV(),
    cv=5
)

# Compare with individual models
from sklearn.model_selection import cross_val_score

rf_scores = cross_val_score(RandomForestRegressor(n_estimators=100, random_state=42), 
                            X, y, cv=5, scoring='r2')
xgb_scores = cross_val_score(xgb.XGBRegressor(n_estimators=100, random_state=42), 
                             X, y, cv=5, scoring='r2')
stack_scores = cross_val_score(stacking, X, y, cv=5, scoring='r2')

print(f"Random Forest: {rf_scores.mean():.4f} (+/- {rf_scores.std():.4f})")
print(f"XGBoost:       {xgb_scores.mean():.4f} (+/- {xgb_scores.std():.4f})")
print(f"Stacking:      {stack_scores.mean():.4f} (+/- {stack_scores.std():.4f})")

Explanation: Stacking combines predictions from multiple models using a meta-learner. It often achieves slightly better results by leveraging the strengths of different algorithms. However, it's more complex and slower to train.

from sklearn.metrics import mean_absolute_error

# Train models
trained_models = {}
for name, model in models.items():
    model.fit(X_train, y_train)
    trained_models[name] = model

# Create test set with outliers
X_test_outliers = X_test.copy()
y_test_outliers = y_test.copy()

# Add 5% extreme outliers
n_outliers = int(0.05 * len(y_test))
outlier_indices = np.random.choice(len(y_test), n_outliers, replace=False)
y_test_outliers[outlier_indices] *= 10  # Make them 10x larger

print("Performance on Clean vs Outlier Data:")
print("-" * 60)
print(f"{'Model':<20} {'Clean MAE':<15} {'Outlier MAE':<15} {'Degradation':<15}")
print("-" * 60)
for name, model in trained_models.items():
    clean_mae = mean_absolute_error(y_test, model.predict(X_test))
    outlier_mae = mean_absolute_error(y_test_outliers, model.predict(X_test))
    degradation = (outlier_mae - clean_mae) / clean_mae * 100
    print(f"{name:<20} {clean_mae:<15.4f} {outlier_mae:<15.4f} {degradation:<14.1f}%")

Explanation: Tree-based models are generally robust to outliers because they use splitting thresholds, not the actual values. SVR with epsilon-insensitive loss is also robust. This test reveals which models degrade gracefully with noisy data.