Scikit-Learn Regression

Real-World Scenario

A property valuation company needs to predict house prices from 80 features including location, square footage, age, and neighborhood demographics. Linear models explain which features drive price (interpretable for regulators). Gradient Boosting achieves the lowest RMSE for the automated valuation model. Both have their place in production.

Linear Regression and Regularization

from sklearn.datasets import fetch_california_housing
from sklearn.linear_model import LinearRegression, Ridge, Lasso, ElasticNet
from sklearn.model_selection import train_test_split, cross_val_score
from sklearn.preprocessing import StandardScaler
from sklearn.metrics import mean_absolute_error, root_mean_squared_error, r2_score
import numpy as np
import pandas as pd

X, y = fetch_california_housing(return_X_y=True, as_frame=True)
feature_names = X.columns.tolist()

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

scaler = StandardScaler()
X_train_s = scaler.fit_transform(X_train)
X_test_s  = scaler.transform(X_test)

models = {
    "Linear Regression": LinearRegression(),
    "Ridge (L2)":        Ridge(alpha=1.0),
    "Lasso (L1)":        Lasso(alpha=0.01, max_iter=10000),
    "Elastic Net":       ElasticNet(alpha=0.01, l1_ratio=0.5, max_iter=10000),
}

results = []
for name, model in models.items():
    model.fit(X_train_s, y_train)
    y_pred = model.predict(X_test_s)
    results.append({
        "Model": name,
        "MAE":   mean_absolute_error(y_test, y_pred),
        "RMSE":  root_mean_squared_error(y_test, y_pred),
        "R²":    r2_score(y_test, y_pred),
    })

print(pd.DataFrame(results).round(4).to_string(index=False))

# Lasso coefficients — feature selection via zero coefficients
lasso = Lasso(alpha=0.01, max_iter=10000).fit(X_train_s, y_train)
coef_df = pd.DataFrame({"feature": feature_names, "coefficient": lasso.coef_})
coef_df = coef_df[coef_df["coefficient"] != 0].sort_values("coefficient", key=abs, ascending=False)
print("\nNon-zero Lasso coefficients:")
print(coef_df.to_string(index=False))

Finding the Best Regularization Strength

from sklearn.linear_model import RidgeCV, LassoCV
from sklearn.datasets import fetch_california_housing
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
from sklearn.metrics import root_mean_squared_error
import numpy as np

X, y = fetch_california_housing(return_X_y=True)
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)

scaler = StandardScaler()
X_train_s = scaler.fit_transform(X_train)
X_test_s  = scaler.transform(X_test)

# RidgeCV finds optimal alpha via cross-validation automatically
alphas = np.logspace(-3, 4, 50)   # try 50 values from 0.001 to 10000
ridge_cv = RidgeCV(alphas=alphas, cv=5, scoring="neg_root_mean_squared_error")
ridge_cv.fit(X_train_s, y_train)
print(f"Ridge best alpha: {ridge_cv.alpha_:.4f}")
print(f"Ridge RMSE: {root_mean_squared_error(y_test, ridge_cv.predict(X_test_s)):.4f}")

# LassoCV with path regularization
lasso_cv = LassoCV(cv=5, max_iter=10000, random_state=42)
lasso_cv.fit(X_train_s, y_train)
print(f"Lasso best alpha: {lasso_cv.alpha_:.6f}")
print(f"Lasso RMSE: {root_mean_squared_error(y_test, lasso_cv.predict(X_test_s)):.4f}")

Gradient Boosting for Regression

from sklearn.datasets import fetch_california_housing
from sklearn.ensemble import HistGradientBoostingRegressor, GradientBoostingRegressor
from sklearn.model_selection import train_test_split, cross_val_score
from sklearn.metrics import mean_absolute_error, root_mean_squared_error, r2_score
import numpy as np

X, y = fetch_california_housing(return_X_y=True)
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)

# HistGradientBoosting — faster, handles missing values, state-of-the-art
model = HistGradientBoostingRegressor(
    max_iter=500,
    learning_rate=0.05,
    max_depth=4,
    min_samples_leaf=20,
    l2_regularization=0.1,
    early_stopping=True,
    validation_fraction=0.1,
    n_iter_no_change=20,
    random_state=42,
)
model.fit(X_train, y_train)

y_pred = model.predict(X_test)
print(f"MAE:  {mean_absolute_error(y_test, y_pred) * 100_000:.0f} USD")
print(f"RMSE: {root_mean_squared_error(y_test, y_pred) * 100_000:.0f} USD")
print(f"R²:   {r2_score(y_test, y_pred):.4f}")
print(f"Stopped at: {model.n_iter_} iterations")

Prediction Intervals with Quantile Regression

from sklearn.ensemble import GradientBoostingRegressor
from sklearn.datasets import fetch_california_housing
from sklearn.model_selection import train_test_split
import numpy as np

X, y = fetch_california_housing(return_X_y=True)
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)

# Train three models: lower bound, median, upper bound
quantile_models = {}
for q, name in [(0.1, "lower_10"), (0.5, "median"), (0.9, "upper_90")]:
    model = GradientBoostingRegressor(
        loss="quantile", alpha=q,
        n_estimators=200, learning_rate=0.05, random_state=42
    )
    model.fit(X_train, y_train)
    quantile_models[name] = model

# Predictions with 80% prediction interval
preds_lower  = quantile_models["lower_10"].predict(X_test[:10])
preds_median = quantile_models["median"].predict(X_test[:10])
preds_upper  = quantile_models["upper_90"].predict(X_test[:10])
actuals      = y_test[:10]

print("Sample predictions (housing price in $100k):")
print(f"{'Actual':>8} {'Lower':>8} {'Median':>8} {'Upper':>8} {'Covered':>8}")
for actual, lo, med, hi in zip(actuals, preds_lower, preds_median, preds_upper):
    covered = "✓" if lo <= actual <= hi else "✗"
    print(f"{actual:8.2f} {lo:8.2f} {med:8.2f} {hi:8.2f} {covered:>8}")

# Coverage rate should be ~80% (matching the 10th–90th percentile interval)
covered = np.mean(
    (quantile_models["lower_10"].predict(X_test) <= y_test) &
    (y_test <= quantile_models["upper_90"].predict(X_test))
)
print(f"\nActual coverage: {covered:.1%} (target: 80%)")

Residual Analysis

from sklearn.linear_model import Ridge
from sklearn.datasets import fetch_california_housing
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
import numpy as np

X, y = fetch_california_housing(return_X_y=True)
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)

scaler = StandardScaler()
X_train_s = scaler.fit_transform(X_train)
X_test_s  = scaler.transform(X_test)

model = Ridge(alpha=1.0)
model.fit(X_train_s, y_train)
y_pred = model.predict(X_test_s)

residuals = y_test - y_pred

# Residuals should be approximately normally distributed with mean ~0
print(f"Residual mean:     {residuals.mean():.4f}  (want: ~0)")
print(f"Residual std:      {residuals.std():.4f}")
print(f"Max over-predict:  {residuals.min():.2f}")  # model predicted too high
print(f"Max under-predict: {residuals.max():.2f}")  # model predicted too low

# Check for heteroscedasticity: residuals vs predicted values
# If residuals increase with predicted value, variance is non-constant
bins = np.percentile(y_pred, [0, 25, 50, 75, 100])
for i in range(len(bins) - 1):
    mask = (y_pred >= bins[i]) & (y_pred < bins[i+1])
    print(f"Predicted quartile {i+1}: residual std = {residuals[mask].std():.3f}")