End-to-End ML Project

Before any code, three questions. They take 30 minutes to answer and save you weeks downstream.

Question 1: What is the business or research objective? "Predict median house price in California districts, accurate enough to inform an automated valuation pipeline that flags overpriced listings."

Question 2: How will success be measured? Median absolute error in dollars, on a held-out random sample of 2024-era California districts. The product team can tolerate ±$25k. Better than ±$25k is the goal.

Question 3: What is the current baseline? A team of human appraisers, manually entering prices, with a reported median absolute error of ±$30k. So beating $30k means we contributed something.

These three answers determine everything else. They tell you:

• The task is supervised regression (continuous target).
• The metric is median absolute error — the median of the per-district absolute errors, which a few wildly-mispredicted districts cannot drag upward, so it stays robust to outliers. (Note this is not MAE: MAE is the mean of those same absolute errors, which the capped, mispredicted districts can pull upward. The two are distinct metrics, and the code in §9 computes both.)
• The baseline you must beat is the human appraisal at ±$30k.

If you cannot answer these three questions about your own project, stop coding and answer them first. They are the difference between an ML project and a programming exercise.

from pathlib import Path
import pandas as pd
import urllib.request
import tarfile

def fetch_housing(data_root: Path = Path("datasets/housing")) -> pd.DataFrame:
    """Download and extract the California Housing CSV."""
    data_root.mkdir(parents=True, exist_ok=True)
    tgz = data_root / "housing.tgz"
    if not tgz.exists():
        url = "https://github.com/ageron/data/raw/main/housing.tgz"
        urllib.request.urlretrieve(url, tgz)
        with tarfile.open(tgz) as tarball:
            tarball.extractall(path=data_root)
    return pd.read_csv(data_root / "housing" / "housing.csv")

housing = fetch_housing()
print(housing.shape)            # (20640, 10)
print(housing.head())

The columns:

Five things to do, in order. They take 15 minutes total. They surface 90% of the dataset's quirks.

3a. Summary statistics and dtypes

housing.info()

What you learn: total_bedrooms has 207 missing values out of 20,640. ocean_proximity is object-typed (categorical, not numeric). Everything else is float.

housing.describe()

What you learn: housing_median_age and median_house_value are capped (the max for both looks suspiciously round at 52.0 and 500001.0). The cap means the model will systematically under-predict expensive houses.

3b. Histograms

import matplotlib.pyplot as plt
housing.hist(bins=50, figsize=(12, 8))
plt.tight_layout()
plt.show()

What you learn: most features are heavily right-skewed (total_rooms, total_bedrooms, population, households). The income distribution is bounded above by capping. Many distributions look like they want a log transform.

3c. Geographic scatter

housing.plot(kind="scatter", x="longitude", y="latitude",
             c="median_house_value", cmap="jet",
             s=housing["population"] / 100,
             alpha=0.4, figsize=(10, 7))
plt.show()

What you learn: California shape becomes visible. Coastal districts are more expensive. The Bay Area and LA cluster around the high prices. Geography matters.

3d. Correlations with the target

corr = housing.select_dtypes(include="number").corr()
print(corr["median_house_value"].sort_values(ascending=False))

What you learn: median_income correlates strongly (~0.69) with house price. The rest are weak (≤0.13 in magnitude). This is the most important diagnostic in the whole EDA: income dominates the signal.

3e. Engineered ratio features

housing["rooms_per_household"] = housing["total_rooms"] / housing["households"]
housing["bedrooms_per_room"] = housing["total_bedrooms"] / housing["total_rooms"]
housing["population_per_household"] = housing["population"] / housing["households"]
corr = housing.select_dtypes(include="number").corr()
print(corr["median_house_value"].sort_values(ascending=False))

What you learn: bedrooms_per_room has a stronger correlation (-0.26) than any of the underlying counts. Engineered features matter. This will become a theme.

Random splitting is fine when the dataset is large and well-mixed. For California Housing, income is so predictive that an income-imbalanced train/test split would bias the evaluation. The fix: stratify on a binned version of the income variable. Stratifying means drawing the test set so each income bracket keeps the same proportion it has in the full dataset, instead of leaving that to chance; pd.cut is what creates those brackets, chopping continuous income into five labeled bins, and StratifiedShuffleSplit samples the held-out 20% from within each bin.

import numpy as np
from sklearn.model_selection import StratifiedShuffleSplit

# Discretize median_income into 5 buckets
housing["income_cat"] = pd.cut(housing["median_income"],
                                bins=[0., 1.5, 3.0, 4.5, 6., np.inf],
                                labels=[1, 2, 3, 4, 5])

split = StratifiedShuffleSplit(n_splits=1, test_size=0.2, random_state=42)
for train_idx, test_idx in split.split(housing, housing["income_cat"]):
    strat_train = housing.iloc[train_idx]
    strat_test = housing.iloc[test_idx]

# Verify the bucket proportions match
print(strat_train["income_cat"].value_counts(normalize=True).sort_index())
print(strat_test["income_cat"].value_counts(normalize=True).sort_index())

# Drop the helper column
for s in (strat_train, strat_test):
    s.drop("income_cat", axis=1, inplace=True)

The stratification matters most for small test sets. The general rule: if the dataset has a feature so predictive that imbalance in that feature would bias the evaluation, stratify on a binned version of it. Time series get a different rule (chronological split, never random).

This is the section where every bug you saw in another chapter (leakage, pipeline contamination, fit-on-test) gets prevented architecturally. The tool is sklearn.pipeline.Pipeline composed with sklearn.compose.ColumnTransformer — the latter routes different columns to different transformers (numeric ones to one sub-pipeline, the categorical one to another) and stitches the outputs back into a single matrix.

Three preprocessing operations on the numeric columns:

1. Impute missing values with the column median (we saw total_bedrooms has NaNs). SimpleImputer(strategy="median") is not just a one-off fillna: during fit it records each column's median as a learned statistic, then transform fills NaNs with that stored value — so the test set is patched with the training median, never its own.
2. Add the three engineered ratio features.
3. Scale every numeric column to zero mean and unit variance.

One preprocessing operation on the categorical column:

1. One-hot encode ocean_proximity into five binary columns — one column per category, with a 1 in the column for that row's value and 0 elsewhere. The alternative, mapping the five categories to integers 0–4, would falsely tell the model that NEAR BAY (say, 3) is "more than" INLAND (1); one-hot avoids inventing an order that isn't there.

from sklearn.compose import ColumnTransformer
from sklearn.pipeline import Pipeline
from sklearn.impute import SimpleImputer
from sklearn.preprocessing import StandardScaler, OneHotEncoder
from sklearn.base import BaseEstimator, TransformerMixin

class RatioFeatures(BaseEstimator, TransformerMixin):
    """Append rooms_per_household, bedrooms_per_room, population_per_household.

    Resolve indices by name from the column order passed in, so the ratios
    can't silently point at the wrong columns if num_columns is reordered.
    """
    def __init__(self, columns):
        self.columns = columns
    def fit(self, X, y=None):
        return self
    def transform(self, X):
        rooms = self.columns.index("total_rooms")
        bedrooms = self.columns.index("total_bedrooms")
        pop = self.columns.index("population")
        household = self.columns.index("households")
        rooms_per_household = X[:, rooms] / X[:, household]
        bedrooms_per_room = X[:, bedrooms] / X[:, rooms]
        pop_per_household = X[:, pop] / X[:, household]
        return np.c_[X, rooms_per_household, bedrooms_per_room, pop_per_household]

num_columns = ["longitude", "latitude", "housing_median_age", "total_rooms",
               "total_bedrooms", "population", "households", "median_income"]
cat_columns = ["ocean_proximity"]

num_pipeline = Pipeline([
    ("impute", SimpleImputer(strategy="median")),
    ("ratios", RatioFeatures(num_columns)),
    ("scale", StandardScaler()),
])

full_pipeline = ColumnTransformer([
    ("num", num_pipeline, num_columns),
    ("cat", OneHotEncoder(handle_unknown="ignore"), cat_columns),
])

# Separate target before fitting
X_train = strat_train.drop("median_house_value", axis=1)
y_train = strat_train["median_house_value"].copy()
X_test = strat_test.drop("median_house_value", axis=1)
y_test = strat_test["median_house_value"].copy()

X_train_prepared = full_pipeline.fit_transform(X_train)
X_test_prepared = full_pipeline.transform(X_test)   # transform, NOT fit_transform

print(X_train_prepared.shape)

The critical pattern: fit_transform on the train set, transform on the test set. The pipeline has learned the median for imputation and the mean/std for scaling from the training data only. The test set never contributes to those statistics.

This is the structural fix to data leakage. Every preprocessing step that learns something (mean, std, vocabulary, PCA basis) becomes a pipeline stage. The pipeline composes them and exposes a single fit_transform/transform interface. You cannot accidentally leak. The companion notebook makes this visceral: it fits a feature selector on the whole dataset before cross-validating against a pure-noise target and watches a clearly positive $R^2$ appear — the coefficient of determination $R^2 = 1 - \mathrm{SS}_\text{res}/\mathrm{SS}_\text{tot}$, the fraction of the target's variance the model explains, where guessing the mean scores 0 and a perfect fit scores 1, so on noise the honest answer is about 0. Moving the selector inside the pipeline collapses that fabricated score back to zero.

RatioFeatures above is your first custom transformer: subclass BaseEstimator, TransformerMixin, implement fit (learn any statistics, here none, so it just returns self) and transform (return the modified array), and it slots into a pipeline exactly like a built-in. TransformerMixin is what gives you fit_transform for free.

The first version of every step is a baseline. Train three model families. Compare on held-out validation (or via cross-validation). Two of them are new here. A decision tree repeatedly splits the data on one feature at a time ("median_income < 3.1?"), funneling each district down to a leaf that predicts the average price of the training rows that landed there; left unconstrained it grows one leaf per training row and memorizes the set. A random forest trains many such trees, each on a bootstrap sample (a random draw of rows with replacement) and a random subset of features, then averages their predictions — the per-tree errors are largely independent, so averaging cancels most of them and the variance drops. (Both get their own chapter later; here you just need to read the scores.)

from sklearn.linear_model import LinearRegression
from sklearn.tree import DecisionTreeRegressor
from sklearn.ensemble import RandomForestRegressor
from sklearn.model_selection import cross_val_score
import numpy as np

def evaluate_model(model, X, y, name: str, cv: int = 5):
    scores = cross_val_score(model, X, y, scoring="neg_mean_squared_error", cv=cv)
    rmse_scores = np.sqrt(-scores)
    print(f"{name:25s}  RMSE = {rmse_scores.mean():.0f}  ± {rmse_scores.std():.0f}")
    return rmse_scores

results = {}
for name, model in [
    ("LinearRegression", LinearRegression()),
    ("DecisionTreeRegressor", DecisionTreeRegressor(random_state=42)),
    ("RandomForestRegressor", RandomForestRegressor(n_estimators=100, random_state=42, n_jobs=-1)),
]:
    results[name] = evaluate_model(model, X_train_prepared, y_train, name)

Typical output on California Housing:

LinearRegression           RMSE = 68627  ± 2280
DecisionTreeRegressor      RMSE = 71407  ± 2300
RandomForestRegressor      RMSE = 50182  ± 1700

RMSE is the square root of the mean squared error; taking the root puts the number back in the target's own units, so an RMSE of 50,182 reads directly as "typically off by about $50k." (Squaring before averaging makes it weight large errors more heavily than the median AE does, which is why the two metrics disagree later.)

The random forest wins by ≈30%. Linear regression is competitive but limited. The single decision tree is overfit (it gets near-zero error on the training set but a much worse cross-validation RMSE — it has effectively memorized the training districts). The random forest's averaging over independent trees is exactly the variance-reduction that fixes that.

The lesson: try multiple model families before tuning. The difference between the best and worst model class is usually 2-10x the difference between a tuned and an untuned best model.

Random forest is the front-runner. Now tune its hyperparameters.

Two approaches. Grid search exhaustively enumerates combinations. Randomized search samples from distributions and is more efficient when you don't know the right ranges.

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

param_dist = {
    "n_estimators": randint(50, 300),
    "max_features": randint(2, 12),
    "max_depth": [None, 10, 20, 30],
    "min_samples_leaf": randint(1, 5),
}

rf = RandomForestRegressor(random_state=42, n_jobs=-1)
search = RandomizedSearchCV(
    rf, param_dist, n_iter=30,
    cv=5, scoring="neg_mean_squared_error",
    random_state=42, n_jobs=-1,
)
search.fit(X_train_prepared, y_train)

print("Best params:", search.best_params_)
print(f"Best CV RMSE: {np.sqrt(-search.best_score_):.0f}")

A typical result drops RMSE from ≈50k to ≈47k. The gains from tuning are real but small. The gains from picking the right model family in step 6 were five times larger. Internalize this. Most ML practitioners spend too much time on hyperparameter search and not enough on feature engineering and model selection.

Feature importance. A random forest scores how much each feature contributed to its splits across all the trees, exposed as feature_importances_ (the scores sum to 1, so each is a fraction of the total). You get it for free:

feature_importance = search.best_estimator_.feature_importances_
# Combine with feature names (you have to assemble them from the pipeline)
num_features = num_columns + ["rooms_per_household", "bedrooms_per_room", "population_per_household"]
cat_features = list(full_pipeline.named_transformers_["cat"].get_feature_names_out(cat_columns))
all_features = num_features + cat_features

importance_pairs = sorted(zip(feature_importance, all_features), reverse=True)
for imp, name in importance_pairs[:10]:
    print(f"  {imp:.3f}  {name}")

You will see median_income at ≈0.35, latitude/longitude clustered at ≈0.15 each, and the engineered ratios in the next tier. Three or four features carry most of the signal.

Now the part most ML tutorials skip: look at the predictions. Where does the model fail?

final_model = search.best_estimator_

# Predictions on the train set
y_train_pred = final_model.predict(X_train_prepared)
residuals = y_train - y_train_pred

# Worst predictions
worst_idx = np.argsort(np.abs(residuals))[-10:]
print("Worst predictions:")
print(pd.DataFrame({
    "actual": y_train.iloc[worst_idx].values,
    "predicted": y_train_pred[worst_idx],
    "error": residuals.iloc[worst_idx].values,
}).round(0))

# Plot residuals vs actual
plt.figure(figsize=(8, 6))
plt.scatter(y_train, residuals, alpha=0.3, s=5)
plt.axhline(0, color="red")
plt.xlabel("actual")
plt.ylabel("residual")
plt.show()

# Plot residuals vs each feature
for feat in ["median_income", "housing_median_age", "longitude"]:
    plt.figure()
    plt.scatter(X_train[feat], residuals, alpha=0.3, s=5)
    plt.axhline(0, color="red")
    plt.xlabel(feat)
    plt.ylabel("residual")
    plt.show()

Three patterns to look for:

1. Heteroscedasticity (error spread that isn't constant across the range): residuals fan out at high predicted values. The model is less precise for expensive houses. Often related to the $500,001 cap in the target.
2. Systematic over/under prediction: residuals are systematically positive (or negative) in certain feature ranges. Suggests a missing nonlinearity or interaction.
3. Outlier districts: a handful of districts with huge errors. Worth looking at individually. They often reveal data quality issues (a district with total_bedrooms = 0, or a population of 10).

This is where you go back to feature engineering. The "iterate" arrow in the workflow points back from here.

This is the moment the test set gets touched. Once. After this, every additional change to the pipeline invalidates the test number.

X_test = strat_test.drop("median_house_value", axis=1)
y_test = strat_test["median_house_value"].copy()
X_test_prepared = full_pipeline.transform(X_test)

y_test_pred = final_model.predict(X_test_prepared)
test_rmse = np.sqrt(np.mean((y_test - y_test_pred) ** 2))
test_mae = np.mean(np.abs(y_test - y_test_pred))
test_median_ae = np.median(np.abs(y_test - y_test_pred))

print(f"Test RMSE        = {test_rmse:.0f}")
print(f"Test MAE         = {test_mae:.0f}")
print(f"Test median AE   = {test_median_ae:.0f}")

# 95% confidence interval on the RMSE, by bootstrap (resample the errors with replacement)
from scipy import stats
def rmse(squared_errors):
    return np.sqrt(np.mean(squared_errors))
squared_errors = (y_test - y_test_pred) ** 2
boot = stats.bootstrap([squared_errors], rmse, confidence_level=0.95, random_state=42)
ci_lo, ci_hi = boot.confidence_interval
print(f"95% CI for test RMSE: [{ci_lo:.0f}, {ci_hi:.0f}]")

Reporting a confidence interval is the second most-skipped step in ML practice (after error analysis). A test RMSE of 47,500 means nothing without "± what?". The 95% interval is the range you'd expect the true RMSE to fall in if you re-ran this on fresh test sets from the same distribution; the bootstrap above gets it by resampling the errors with replacement thousands of times and recomputing the RMSE each time, which assumes nothing about how the errors are distributed. For California Housing with 4,128 test samples, the 95% CI is typically [44,500, 50,500]. (The companion notebook loads the same data via fetch_california_housing, whose target is in units of $100k, so it prints this same RMSE as ≈ 0.47 — the dollar figure divided by 100,000.) Now you know the model beats the $30k baseline... or do you? Re-read the framing: the baseline was stated as a median absolute error of $30k, so the only number that legitimately compares against it is the model's median AE — printed above as test_median_ae, not test_rmse or test_mae. RMSE, MAE, and median AE are three different numbers — RMSE's squaring inflates it well above the median error, and MAE (the mean) sits between the two — so compare against the baseline on its own metric.

A trained sklearn pipeline contains the model and all the preprocessing. You save them together. You load them together. You never re-implement the preprocessing at inference time, because the bugs are guaranteed.

import joblib

# Save: bundle preprocessing pipeline + model
final_pipeline = Pipeline([
    ("preprocess", full_pipeline),
    ("model", final_model),
])
final_pipeline.fit(strat_train.drop("median_house_value", axis=1), y_train)
joblib.dump(final_pipeline, "housing_model.joblib")

# Load and use:
loaded = joblib.load("housing_model.joblib")
new_district = pd.DataFrame([{
    "longitude": -118.5,
    "latitude": 34.2,
    "housing_median_age": 25,
    "total_rooms": 1500,
    "total_bedrooms": 300,
    "population": 800,
    "households": 280,
    "median_income": 5.0,
    "ocean_proximity": "INLAND",
}])
predicted_price = loaded.predict(new_district)[0]
print(f"Predicted: ${predicted_price:,.0f}")

joblib.dump is the sklearn convention; pickle works too. For very large models, use joblib.dump(..., compress=3). For interop with other languages, export to ONNX (skl2onnx). For LLMs, use safetensors (another chapter).

Shipping is not done. After deployment, you monitor four things:

1. Input distribution: are the features the model sees in production drawn from the same distribution as training? Kolmogorov-Smirnov test per feature, weekly. Alert on shift.
2. Output distribution: are the predictions you produce drifting? Compare current week's prediction histogram to the test-set prediction histogram.
3. Performance metric: when ground truth becomes available (delayed labels), compute the production performance metric. Compare to the test-set CI from step 9. If you fall outside, retrain.
4. Latency and error rates: are predictions returning within SLA? Are there exceptions in the pipeline? Standard service-monitoring rules.

The Made-With-ML monitoring chapter (06-practice/madewithml-mlops-monitoring) walks through the tooling: WhyLabs, Evidently AI, custom dashboards in Grafana. The principles are the same regardless of tool: measure shift, alert on shift, retrain on shift.

A typical retraining cadence for a stable problem like housing prices is monthly or quarterly. For a fast-shifting problem (fraud, recommendations), it can be hourly. The cadence is a business decision informed by how fast the world changes.

For reference, here is the entire end-to-end pipeline in one place. Open this as notebook.ipynb in Colab, run cell by cell, modify, iterate. The whole thing is ≈80 lines.

# Cell 1: imports
import urllib.request, tarfile
from pathlib import Path
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.compose import ColumnTransformer
from sklearn.pipeline import Pipeline
from sklearn.impute import SimpleImputer
from sklearn.preprocessing import StandardScaler, OneHotEncoder
from sklearn.model_selection import StratifiedShuffleSplit, RandomizedSearchCV
from sklearn.ensemble import RandomForestRegressor
from scipy.stats import randint
import joblib

# Cell 2: data
data_root = Path("datasets/housing")
data_root.mkdir(parents=True, exist_ok=True)
tgz = data_root / "housing.tgz"
if not tgz.exists():
    urllib.request.urlretrieve("https://github.com/ageron/data/raw/main/housing.tgz", tgz)
    with tarfile.open(tgz) as t:
        t.extractall(path=data_root)
housing = pd.read_csv(data_root / "housing" / "housing.csv")

# Cell 3: stratified split
housing["income_cat"] = pd.cut(housing["median_income"],
                                bins=[0., 1.5, 3.0, 4.5, 6., np.inf],
                                labels=[1, 2, 3, 4, 5])
split = StratifiedShuffleSplit(n_splits=1, test_size=0.2, random_state=42)
for train_idx, test_idx in split.split(housing, housing["income_cat"]):
    strat_train, strat_test = housing.iloc[train_idx], housing.iloc[test_idx]
for s in (strat_train, strat_test):
    s.drop("income_cat", axis=1, inplace=True)

# Cell 4: pipeline
X_train = strat_train.drop("median_house_value", axis=1)
y_train = strat_train["median_house_value"].copy()
num_columns = ["longitude", "latitude", "housing_median_age", "total_rooms",
               "total_bedrooms", "population", "households", "median_income"]
cat_columns = ["ocean_proximity"]
num_pipeline = Pipeline([("impute", SimpleImputer(strategy="median")), ("scale", StandardScaler())])
full_pipeline = ColumnTransformer([
    ("num", num_pipeline, num_columns),
    ("cat", OneHotEncoder(handle_unknown="ignore"), cat_columns),
])

# Cell 5: train + tune
param_dist = {"n_estimators": randint(50, 300), "max_features": randint(2, 12)}
rf = RandomForestRegressor(random_state=42, n_jobs=-1)
final_pipeline = Pipeline([("prep", full_pipeline), ("model", rf)])
search = RandomizedSearchCV(
    final_pipeline,
    {"model__" + k: v for k, v in param_dist.items()},
    n_iter=20, cv=5, scoring="neg_mean_squared_error",
    random_state=42, n_jobs=-1,
)
search.fit(X_train, y_train)
print(f"Best CV RMSE: {np.sqrt(-search.best_score_):.0f}")

# Cell 6: test (once!)
X_test = strat_test.drop("median_house_value", axis=1)
y_test = strat_test["median_house_value"].copy()
y_pred = search.best_estimator_.predict(X_test)
print(f"Test RMSE: {np.sqrt(np.mean((y_test - y_pred)**2)):.0f}")

# Cell 7: save
joblib.dump(search.best_estimator_, "housing_model.joblib")

That is the entire chapter, runnable in three minutes on a free Colab. The rest is interpretation, intuition, and iteration.