Case studies

import numpy as np
import pandas as pd

df = pd.read_csv('housing.csv')
df['value'] = df.median_house_value / 100_000            # in $100,000s
rooms = df.total_rooms / df.households

n_rows = len(df)
at_value_cap = (df.median_house_value >= 500_001).sum()   # top-coded at $500,001
at_age_cap = (df.housing_median_age == 52).sum()          # top-coded at 52 years
at_income_cap = (df.median_income >= 15).sum()            # top-coded at 15 ($150,000)
print(n_rows, 'districts')
print('at the value cap:', at_value_cap, f'({at_value_cap / n_rows:.1%})')
print('at the age cap:', at_age_cap, '| at the income cap:', at_income_cap)
print('rooms per household: median', round(rooms.median(), 1), '| max', round(rooms.max(), 1))

from sklearn.model_selection import KFold, PredefinedSplit, cross_val_score
from sklearn.linear_model import LinearRegression
from sklearn.neighbors import KNeighborsRegressor

X = pd.DataFrame({
    'income': df.median_income, 'age': df.housing_median_age, 'rooms': rooms,
    'bedroom_share': df.total_bedrooms / df.total_rooms,
    'people_per_household': df.population / df.households,
    'lat': df.latitude, 'lon': df.longitude})
jitter = np.random.default_rng(0).normal(0, 1e-6, (len(X), 2))
X[['lat', 'lon']] += jitter      # districts that share coordinates: makes "nearest" unique
y = df.value
baseline_mae = (y - y.mean()).abs().mean()

block = np.floor(df.latitude).astype(int) * 1000 + np.floor(df.longitude).astype(int)
sizes = block.value_counts().sort_index().sort_values(ascending=False, kind='stable')
n_blocks = len(sizes)
load, fold = np.zeros(5), {}
for b, n in sizes.items():       # largest region first, each to the emptiest fold
    k = int(load.argmin()); fold[b] = k; load[k] += n
random_cv = KFold(5, shuffle=True, random_state=0)
region_cv = PredefinedSplit(block.map(fold))

def r2(model, cols):
    a = cross_val_score(model, X[cols], y, cv=random_cv, scoring='r2').mean()
    b = cross_val_score(model, X[cols], y, cv=region_cv, scoring='r2').mean()
    return round(a, 3), round(b, 3)

features = ['income', 'age', 'rooms', 'bedroom_share', 'people_per_household']
linear = r2(LinearRegression(), features)
knn = r2(KNeighborsRegressor(10), ['lat', 'lon'])
print('baseline MAE:', round(baseline_mae, 3))
print('R² (random folds, new regions)')
print('  linear, no location:', linear)
print('  10 nearest districts:', knn)

from sklearn.ensemble import HistGradientBoostingRegressor

boost = r2(HistGradientBoostingRegressor(random_state=0), list(X.columns))
print('gradient boosting, all features:', boost)

from sklearn.model_selection import train_test_split

X_tr, X_te, y_tr, y_te = train_test_split(X, y, test_size=0.2, random_state=0)
model = HistGradientBoostingRegressor(random_state=0).fit(X_tr, y_tr)
pred = pd.Series(model.predict(X_te), index=y_te.index)
err = (pred - y_te).abs()

capped = y_te >= 5.00001
worst = err >= err.quantile(0.95)                         # the 5% largest errors
print('test MAE:', round(err.mean(), 3), '| capped districts in test:', capped.sum())
print('mean prediction for capped districts:', round(pred[capped].mean(), 2))
print('share of capped districts among the worst 5%:',
      round((capped & worst).sum() / worst.sum(), 3))
print('highest prediction:', round(pred.max(), 2))