Scikit-Learn Regression

Regression predicts a number. Every example here ran with scikit-learn 1.9.1 on the California housing data (20,640 districts, target = median house value in $100,000s), which fetch_california_housing downloads once and then caches.

Baseline First

Compare every model against a DummyRegressor that predicts the training mean.

from sklearn.datasets import fetch_california_housing
from sklearn.dummy import DummyRegressor
from sklearn.ensemble import HistGradientBoostingRegressor
from sklearn.linear_model import LinearRegression
from sklearn.metrics import r2_score, root_mean_squared_error
from sklearn.model_selection import train_test_split
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler

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)

models = {"mean baseline": DummyRegressor(),
          "linear": make_pipeline(StandardScaler(), LinearRegression()),
          "hist_gb": HistGradientBoostingRegressor(random_state=0)}
for name, model in models.items():
    pred = model.fit(X_train, y_train).predict(X_test)
    print(f"{name:13s} RMSE {root_mean_squared_error(y_test, pred):.3f}"
          f"  R2 {r2_score(y_test, pred):6.3f}")

$ python compare.py
mean baseline RMSE 1.145  R2 -0.000
linear        RMSE 0.746  R2  0.576
hist_gb       RMSE 0.462  R2  0.837

The mean baseline scores R² of 0 by definition, and its RMSE of 1.145 ($114,500) is the yardstick. The linear model cuts the error by a third. Gradient boosting cuts it by more than half, because house value depends on location in ways a straight line can’t express.

Correlated Features Break Plain Least Squares

When two features carry the same information, ordinary least squares can’t decide how to split the weight between them. Here x2 is x1 plus a tiny amount of noise.

import numpy as np
from sklearn.linear_model import LinearRegression, Ridge

rng = np.random.default_rng(0)
for trial in range(3):
    x1 = rng.normal(size=100)
    x2 = x1 + rng.normal(scale=0.01, size=100)        # near-duplicate
    y = 3 * x1 + rng.normal(scale=0.5, size=100)
    X = np.column_stack([x1, x2])
    ols = LinearRegression().fit(X, y).coef_.round(1)
    ridge = Ridge(alpha=1.0).fit(X, y).coef_.round(2)
    print(f"sample {trial}: OLS coefs {ols}   Ridge coefs {ridge}")

$ python collinear.py
sample 0: OLS coefs [1.2 1.8]   Ridge coefs [1.49 1.49]
sample 1: OLS coefs [2.1 0.9]   Ridge coefs [1.51 1.51]
sample 2: OLS coefs [1.9 1.1]   Ridge coefs [1.49 1.49]

The true model is 3 * x1, and the two features share that weight of 3. Plain least squares splits it differently for every sample (1.2/1.8, then 2.1/0.9), because almost any split fits the data equally well. The coefficients are unstable even though predictions barely change. Ridge’s penalty prefers the smallest coefficients that fit, which is an even 1.5/1.5 split every time.

Lasso Zeroes Out Features

Lasso’s L1 penalty drives weak coefficients to exactly zero, and a larger alpha removes more of them.

from sklearn.datasets import fetch_california_housing
from sklearn.linear_model import Lasso
from sklearn.preprocessing import StandardScaler

data = fetch_california_housing()
X = StandardScaler().fit_transform(data.data)

for alpha in [0.001, 0.05, 0.2, 0.5]:
    coef = Lasso(alpha=alpha).fit(X, data.target).coef_
    kept = [n for n, c in zip(data.feature_names, coef) if c != 0]
    print(f"alpha={alpha:<5}  {len(kept)} kept: {kept}")

$ python lasso_path.py
alpha=0.001  8 kept: ['MedInc', 'HouseAge', 'AveRooms', 'AveBedrms', 'Population', 'AveOccup', 'Latitude', 'Longitude']
alpha=0.05   4 kept: ['MedInc', 'HouseAge', 'Latitude', 'Longitude']
alpha=0.2    1 kept: ['MedInc']
alpha=0.5    1 kept: ['MedInc']

At a tiny alpha Lasso keeps all eight features, like ordinary regression. As alpha grows, weaker features drop out entirely, and by 0.2 only median income is left, the single strongest predictor. Lasso gives a ranked shortlist of features along with the model.

Choosing alpha by Cross-Validation

RidgeCV and LassoCV search the penalty strength internally.

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

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)

for est in [RidgeCV(alphas=np.logspace(-3, 3, 13)), LassoCV(cv=5, random_state=0)]:
    pipe = make_pipeline(StandardScaler(), est).fit(X_train, y_train)
    rmse = root_mean_squared_error(y_test, pipe.predict(X_test))
    print(f"{type(est).__name__:8s} alpha {pipe[-1].alpha_:.4g}   test RMSE {rmse:.4f}")

$ python alpha_cv.py
RidgeCV  alpha 3.162   test RMSE 0.7455
LassoCV  alpha 0.0007985   test RMSE 0.7448

Both pick a small penalty and land on the same test error as the unregularised linear model (0.746). With 16,000 rows and only 8 features there is little overfitting for regularisation to fix. Its payoff comes with many features, few rows, or correlated inputs.

Prediction Intervals with Quantile Loss

A point prediction hides uncertainty. Two quantile models give a range instead.

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

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)

lo, hi = (HistGradientBoostingRegressor(loss="quantile", quantile=q, random_state=0)
          .fit(X_train, y_train).predict(X_test) for q in (0.1, 0.9))

print("first 4 intervals ($100k):")
for t, a, b in zip(y_test[:4], lo[:4], hi[:4]):
    print(f"  true {t:.2f}  interval [{a:.2f}, {b:.2f}]")
print(f"coverage on {len(y_test)} test rows: {np.mean((lo <= y_test) & (y_test <= hi)):.1%}")

$ python quantiles.py
first 4 intervals ($100k):
  true 0.48  interval [0.45, 0.83]
  true 0.46  interval [0.64, 1.56]
  true 5.00  interval [2.53, 5.29]
  true 2.19  interval [2.06, 3.02]
coverage on 4128 test rows: 73.1%

The two models bound a nominal 80% interval, but only 73% of test values actually fall inside it, so these intervals are too narrow. Quantile models trained this way often under-cover, which is why coverage must be measured, not assumed. The third row is a district at the dataset’s $500,000 cap, which the next section explains.

Residuals and a Capped Target

Residuals are true minus predicted values. Grouping them by the size of the true value shows where the model is biased.

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

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)
resid = y_test - HistGradientBoostingRegressor(random_state=0).fit(X_train, y_train).predict(X_test)

print(f"max target {y.max():.5f}, rows at that max: {(y == y.max()).sum()}")
print(f"mean residual overall: {resid.mean():+.3f}")
for lo, hi in [(0, 1), (1, 2), (2, 3), (3, 4.9), (4.9, 6)]:
    m = (y_test >= lo) & (y_test < hi)
    print(f"true value {lo}-{hi}: n={m.sum():4d}  mean residual {resid[m].mean():+.3f}")

$ python residuals.py
max target 5.00001, rows at that max: 965
mean residual overall: -0.005
true value 0-1: n= 730  mean residual -0.186
true value 1-2: n=1684  mean residual -0.092
true value 2-3: n= 956  mean residual +0.008
true value 3-4.9: n= 568  mean residual +0.258
true value 4.9-6: n= 190  mean residual +0.603

The average residual is near zero, yet the model over-predicts cheap districts (negative residuals: the prediction is above the true value) and under-predicts expensive ones. Some of that pattern appears whenever you group by the true value, because predictions pull toward the middle. The top bucket is extreme for another reason: 965 districts are recorded at exactly 5.00001, a cap in the source data, so their real values are higher than any label says. No model can learn values the data doesn’t contain; treat capped rows separately or drop them.

Practice

from sklearn.datasets import fetch_california_housing
from sklearn.ensemble import RandomForestRegressor
from sklearn.metrics import r2_score, root_mean_squared_error
from sklearn.model_selection import train_test_split

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)
pred = RandomForestRegressor(n_estimators=100, n_jobs=-1, random_state=0).fit(
    X_train, y_train).predict(X_test)
print(f"forest RMSE {root_mean_squared_error(y_test, pred):.3f}  R2 {r2_score(y_test, pred):.3f}")

$ python p1_forest.py
forest RMSE 0.504  R2 0.806

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

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)
lo, hi = (HistGradientBoostingRegressor(loss="quantile", quantile=q, random_state=0)
          .fit(X_train, y_train).predict(X_test) for q in (0.05, 0.95))
print(f"90% interval coverage: {np.mean((lo <= y_test) & (y_test <= hi)):.1%}")
print(f"mean width: {np.mean(hi - lo):.2f} ($100k)")

$ python p2_wider.py
90% interval coverage: 84.1%
mean width: 1.36 ($100k)

from sklearn.datasets import fetch_california_housing
from sklearn.linear_model import Ridge
from sklearn.preprocessing import StandardScaler

data = fetch_california_housing()
X = StandardScaler().fit_transform(data.data)
for alpha in [1, 10_000]:
    coef = Ridge(alpha=alpha).fit(X, data.target).coef_
    print(f"alpha={alpha:<6}", " ".join(f"{c:+.2f}" for c in coef))
print("features:   ", " ".join(n[:5] for n in data.feature_names))

$ python p3_ridge_shrink.py
alpha=1      +0.83 +0.12 -0.27 +0.31 -0.00 -0.04 -0.90 -0.87
alpha=10000  +0.53 +0.12 +0.03 -0.02 +0.00 -0.03 -0.16 -0.12
features:    MedIn House AveRo AveBe Popul AveOc Latit Longi