What you’ll learn

• the full roadmap of a real ML project, start to finish, in one picture
• how the previous pages in this phase fit together into housing_preparedhousing_prepared
• training and comparing several models: Linear Regression, Decision Tree, Random Forest
• why a model with zero training error is a red flag, not a win
• evaluating models honestly with k-fold cross-validation
• fine-tuning hyperparameters automatically with GridSearchCV
• reading feature_importances_feature_importances_ to see what the best model actually relies on
• the one and only time you’re allowed to touch the test set, plus a confidence interval for the final RMSE

The whole project, in one picture

Every page in this phase has been one stage of a single project. This page is the capstone: it strings every stage together and adds the two stages that come after data preparation — training models and fine-tuning them — finishing with the one evaluation that actually counts.

  flowchart LR
  A["Frame the problem"] --> B["Get the data"]
  B --> C["Explore & visualize"]
  C --> D["Prepare the data
(pipelines)"]
  D --> E["Select & train models"]
  E --> F["Fine-tune
(cross-val + grid search)"]
  F --> G["Evaluate on test set"]
  G --> H["Launch, monitor,
maintain"]

Recap: how we got to housing_preparedhousing_prepared

If you’ve worked through this phase in order, you already have everything you need sitting in two variables. Quick recap of where they came from:

1. Frame the problem & get the data — a supervised, batch, multiple regression problem: predict median_house_valuemedian_house_value from census features, measured with RMSE.
2. Create a test set — StratifiedShuffleSplitStratifiedShuffleSplit on income brackets carves off strat_test_setstrat_test_set immediately, and it’s locked away until the very end.
3. Explore & find correlations — median_incomemedian_income is by far the strongest predictor; engineered ratios like bedrooms_per_roombedrooms_per_room beat the raw counts.
4. Clean, encode, scale — missing total_bedroomstotal_bedrooms values imputed, the ocean_proximityocean_proximity text column one-hot encoded, every numeric column scaled.
5. Pipelines — all of the above chained into one full_pipelinefull_pipeline, a ColumnTransformerColumnTransformer wrapping a numeric PipelinePipeline and a categorical encoder.

housing = strat_train_set.drop("median_house_value", axis=1)
housing_labels = strat_train_set["median_house_value"].copy()
 
housing_prepared = full_pipeline.fit_transform(housing)
print(housing_prepared.shape)

housing = strat_train_set.drop("median_house_value", axis=1)
housing_labels = strat_train_set["median_house_value"].copy()
 
housing_prepared = full_pipeline.fit_transform(housing)
print(housing_prepared.shape)

housing_preparedhousing_prepared is a plain numeric matrix — every row model-ready. From here on, it’s finally time to train something.

Select and train a model

The good news: after all that prep, training itself is almost anticlimactic. Start with the simplest reasonable model, Linear Regression:

from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error
import numpy as np
 
lin_reg = LinearRegression()
lin_reg.fit(housing_prepared, housing_labels)
 
housing_predictions = lin_reg.predict(housing_prepared)
lin_rmse = np.sqrt(mean_squared_error(housing_labels, housing_predictions))
print(round(lin_rmse, 0))
# 68628.0

from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error
import numpy as np
 
lin_reg = LinearRegression()
lin_reg.fit(housing_prepared, housing_labels)
 
housing_predictions = lin_reg.predict(housing_prepared)
lin_rmse = np.sqrt(mean_squared_error(housing_labels, housing_predictions))
print(round(lin_rmse, 0))
# 68628.0

Most districts’ median_house_valuemedian_house_value sits between 120,000 and 120,000 and 265,000, so a typical error of ~$68,628 isn’t great — a sign of underfitting: the model isn’t powerful enough (or the features aren’t informative enough) to capture the pattern. Let’s try something more expressive, a Decision Tree:

from sklearn.tree import DecisionTreeRegressor
 
tree_reg = DecisionTreeRegressor(random_state=42)
tree_reg.fit(housing_prepared, housing_labels)
 
housing_predictions = tree_reg.predict(housing_prepared)
tree_rmse = np.sqrt(mean_squared_error(housing_labels, housing_predictions))
print(tree_rmse)
# 0.0

from sklearn.tree import DecisionTreeRegressor
 
tree_reg = DecisionTreeRegressor(random_state=42)
tree_reg.fit(housing_prepared, housing_labels)
 
housing_predictions = tree_reg.predict(housing_prepared)
tree_rmse = np.sqrt(mean_squared_error(housing_labels, housing_predictions))
print(tree_rmse)
# 0.0

Zero error sounds like a win — it isn’t. A Decision Tree can carve the training data into arbitrarily small, pure regions until every training row is predicted perfectly. That’s overfitting, not skill, and you can’t tell the difference by looking at training error alone. You need to evaluate on data the model hasn’t memorized.

Better evaluation using cross-validation

You could carve out a separate validation set with train_test_splittrain_test_split, but k-fold cross-validation gets more out of the same data: it splits the training set into kk folds, trains on k - 1k - 1 of them, evaluates on the fold left out, and repeats until every fold has been the evaluation fold exactly once.

from sklearn.model_selection import cross_val_score
 
scores = cross_val_score(
    tree_reg, housing_prepared, housing_labels,
    scoring="neg_mean_squared_error", cv=10,
)
tree_rmse_scores = np.sqrt(-scores)
 
def display_scores(scores):
    print("Scores:", scores)
    print("Mean:", scores.mean())
    print("Standard deviation:", scores.std())
 
display_scores(tree_rmse_scores)
# Mean: ~71407   Standard deviation: ~2439

from sklearn.model_selection import cross_val_score
 
scores = cross_val_score(
    tree_reg, housing_prepared, housing_labels,
    scoring="neg_mean_squared_error", cv=10,
)
tree_rmse_scores = np.sqrt(-scores)
 
def display_scores(scores):
    print("Scores:", scores)
    print("Mean:", scores.mean())
    print("Standard deviation:", scores.std())
 
display_scores(tree_rmse_scores)
# Mean: ~71407   Standard deviation: ~2439

Scikit-learn’s cross-validation tools expect a utility function (bigger is better), not a cost function — that’s why scoring="neg_mean_squared_error"scoring="neg_mean_squared_error" returns negative values, and the code flips the sign with -scores-scores before taking the square root.

Once you cross-validate honestly, the Decision Tree’s real performance shows up: it’s actually worse than plain Linear Regression (mean RMSE ≈ 71,407 vs. ≈ 69,052) — it was simply overfitting the training set, not learning it. Cross- validation gives you both a performance estimate and a sense of how precise that estimate is (the standard deviation), something a single train/test split can’t tell you.

Try one more model, a Random Forest — many Decision Trees trained on random subsets of the data and features, with their predictions averaged (an idea called Ensemble Learning, covered in depth in Phase 5):

from sklearn.ensemble import RandomForestRegressor
 
forest_reg = RandomForestRegressor(n_estimators=100, random_state=42)
forest_scores = cross_val_score(
    forest_reg, housing_prepared, housing_labels,
    scoring="neg_mean_squared_error", cv=10,
)
display_scores(np.sqrt(-forest_scores))
# Mean: ~50182   Standard deviation: ~2097

from sklearn.ensemble import RandomForestRegressor
 
forest_reg = RandomForestRegressor(n_estimators=100, random_state=42)
forest_scores = cross_val_score(
    forest_reg, housing_prepared, housing_labels,
    scoring="neg_mean_squared_error", cv=10,
)
display_scores(np.sqrt(-forest_scores))
# Mean: ~50182   Standard deviation: ~2097

Random Forests are the clear winner of the three so far — but notice that its training-set error is still much lower than its cross-validated error, meaning it’s still overfitting somewhat too. The fix isn’t to panic; it’s to shortlist a couple of promising models (don’t over-tune any single one yet) and move on to fine-tuning.

Fine-tune with GridSearchCV

Manually trying hyperparameter combinations by hand is tedious and unsystematic. GridSearchCVGridSearchCV does it for you: give it a grid of values per hyperparameter and it cross-validates every combination, then hands back the best one.

from sklearn.model_selection import GridSearchCV
 
param_grid = [
    {"n_estimators": [3, 10, 30], "max_features": [2, 4, 6, 8]},
    {"bootstrap": [False], "n_estimators": [3, 10], "max_features": [2, 3, 4]},
]
 
forest_reg = RandomForestRegressor(random_state=42)
grid_search = GridSearchCV(
    forest_reg, param_grid, cv=5,
    scoring="neg_mean_squared_error", return_train_score=True,
)
grid_search.fit(housing_prepared, housing_labels)
 
print(grid_search.best_params_)
# {'max_features': 8, 'n_estimators': 30}

from sklearn.model_selection import GridSearchCV
 
param_grid = [
    {"n_estimators": [3, 10, 30], "max_features": [2, 4, 6, 8]},
    {"bootstrap": [False], "n_estimators": [3, 10], "max_features": [2, 3, 4]},
]
 
forest_reg = RandomForestRegressor(random_state=42)
grid_search = GridSearchCV(
    forest_reg, param_grid, cv=5,
    scoring="neg_mean_squared_error", return_train_score=True,
)
grid_search.fit(housing_prepared, housing_labels)
 
print(grid_search.best_params_)
# {'max_features': 8, 'n_estimators': 30}

The first dict tries 3 x 4 = 123 x 4 = 12 combinations; the second tries 2 x 3 = 62 x 3 = 6 more with bootstrapbootstrap forced to FalseFalse — 18 combinations total, each cross-validated 5 times, for 90 training runs in one .fit().fit() call. GridSearchCVGridSearchCV defaults to refit=Truerefit=True, so once it finds the winner it automatically retrains it on the whole training set — grid_search.best_estimator_grid_search.best_estimator_ is ready to use immediately.

You can also treat data-prep choices as hyperparameters — for example, searching over whether add_bedrooms_per_roomadd_bedrooms_per_room should be TrueTrue or FalseFalse inside a custom transformer, right alongside n_estimatorsn_estimators. If the search space gets too big for an exhaustive grid, RandomizedSearchCVRandomizedSearchCV samples a fixed number of random combinations instead of trying every single one.

Analyze the best model and its errors

Before you touch the test set, it’s worth digging into why the best model works. A fitted RandomForestRegressorRandomForestRegressor exposes feature_importances_feature_importances_ — a relative score of how much each input contributed to reducing error across the forest:

feature_importances = grid_search.best_estimator_.feature_importances_
 
extra_attribs = ["rooms_per_hhold", "pop_per_hhold", "bedrooms_per_room"]
cat_encoder = full_pipeline.named_transformers_["cat"]
cat_one_hot_attribs = list(cat_encoder.categories_[0])
attributes = num_attribs + extra_attribs + cat_one_hot_attribs
 
print(sorted(zip(feature_importances, attributes), reverse=True))
# [(0.366, 'median_income'), (0.165, 'INLAND'), (0.109, 'pop_per_hhold'), ...]

feature_importances = grid_search.best_estimator_.feature_importances_
 
extra_attribs = ["rooms_per_hhold", "pop_per_hhold", "bedrooms_per_room"]
cat_encoder = full_pipeline.named_transformers_["cat"]
cat_one_hot_attribs = list(cat_encoder.categories_[0])
attributes = num_attribs + extra_attribs + cat_one_hot_attribs
 
print(sorted(zip(feature_importances, attributes), reverse=True))
# [(0.366, 'median_income'), (0.165, 'INLAND'), (0.109, 'pop_per_hhold'), ...]

median_incomemedian_income dominates, just as the correlation analysis suggested back in Exploratory Data Analysis & Correlations — and the engineered pop_per_hholdpop_per_hhold ratio beats several raw columns. Meanwhile, ISLANDISLAND and NEAR BAYNEAR BAY contribute almost nothing. This kind of ranking tells you where to focus next: maybe drop the near-useless ocean_proximityocean_proximity categories, maybe engineer another ratio near the top of the list. You should also look at the specific districts your model gets most wrong and ask why — outliers, a missing feature, or noisy labels are all common culprits.

Evaluate your system on the test set

This is the moment the whole phase has been building toward: run the locked-away test set through the fitted full_pipelinefull_pipeline (.transform().transform(), never .fit_transform().fit_transform() — the test set must never influence fitting), and score the final model exactly once.

final_model = grid_search.best_estimator_
 
X_test = strat_test_set.drop("median_house_value", axis=1)
y_test = strat_test_set["median_house_value"].copy()
 
X_test_prepared = full_pipeline.transform(X_test)
final_predictions = final_model.predict(X_test_prepared)
 
final_rmse = np.sqrt(mean_squared_error(y_test, final_predictions))
print(round(final_rmse, 0))
# 47730.0

final_model = grid_search.best_estimator_
 
X_test = strat_test_set.drop("median_house_value", axis=1)
y_test = strat_test_set["median_house_value"].copy()
 
X_test_prepared = full_pipeline.transform(X_test)
final_predictions = final_model.predict(X_test_prepared)
 
final_rmse = np.sqrt(mean_squared_error(y_test, final_predictions))
print(round(final_rmse, 0))
# 47730.0

The test-set RMSE (~$47,730) came out slightly better than the cross-validated estimate here — but if yours comes out a bit worse, resist the urge to keep tweaking hyperparameters until the test-set number looks good. Any change made because of the test set defeats its entire purpose: it stops being an honest estimate of how the model performs on data it has never influenced.

A single RMSE number is a point estimate — it doesn’t say how precise that estimate is. scipy.statsscipy.stats can compute a 95% confidence interval for the generalization error from the per-instance squared errors:

from scipy import stats
 
confidence = 0.95
squared_errors = (final_predictions - y_test) ** 2
 
interval = np.sqrt(
    stats.t.interval(
        confidence, len(squared_errors) - 1,
        loc=squared_errors.mean(), scale=stats.sem(squared_errors),
    )
)
print(interval.round(0))
# [45685. 49691.]

from scipy import stats
 
confidence = 0.95
squared_errors = (final_predictions - y_test) ** 2
 
interval = np.sqrt(
    stats.t.interval(
        confidence, len(squared_errors) - 1,
        loc=squared_errors.mean(), scale=stats.sem(squared_errors),
    )
)
print(interval.round(0))
# [45685. 49691.]

That narrow band (roughly 45,685–45,685–49,691) is reassuring — it means the point estimate of ~$47,730 isn’t a fluke of which rows happened to land in the test set. A model whose confidence interval is wide relative to its RMSE deserves more scrutiny before you trust the headline number.

🧪 Try It Yourself

Exercise 1 – Fit a preprocessing + model pipeline

Exercise 2 – Cross-validate with cross_val_score

Exercise 3 – Find the best hyperparameters with GridSearchCV

Next

That’s the full Géron Chapter 2 roadmap, start to finish. Continue to Phase 3: Supervised Learning - Regression, where you’ll go much deeper into the models this page only introduced — starting with the math behind Linear Regression itself.