The Project | About Us | Contribute | Donations | License

HOME

Eight modules with MercaFresh have shown you the way; now it is time to walk it on your own. In this first project we temporarily leave the supermarket behind to solve a classic end-to-end regression problem: predicting the median house price in California districts from census data. You will practice the full 7-phase workflow you saw in 01-05 "The Machine Learning project workflow": definition, data, exploration, preprocessing, modeling, evaluation and conclusions. We will use the California Housing dataset, bundled with scikit-learn (fetch_california_housing), so you need no external accounts or Kaggle downloads: everything runs in your usual Python environment.

Problem definition and target metric
Getting and loading the data
Exploratory data analysis (EDA)
Preprocessing with a Pipeline
Modeling: from the baseline to ensembles
Honest comparison with cross-validation
Error analysis
Final evaluation on the test set
Project conclusions

Problem definition and target metric

Before writing a single line of code, we define the problem as you learned in 01-05:

Business question: given a California district described by 8 census variables, estimate the median value of its homes.
Problem type: supervised regression (the target variable is continuous), as in 04-01 "Linear regression".
Primary metric: RMSE (root mean squared error), because it penalizes large errors more heavily and is expressed in the same units as the price. We will pair it with MAE (more robust to outliers) and R² (proportion of variance explained), all three introduced in 06-02 "Evaluation metrics".
Success criterion: clearly beat a trivial baseline and achieve a competitive R² (for this dataset, anything above 0.80 is a good result).

The target price is expressed in hundreds of thousands of dollars (a value of 2.5 means $250,000). One important detail we will uncover during EDA: prices are censored at 5.0 (every district above $500,000 shows up as 5.00001).

Getting and loading the data

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.datasets import fetch_california_housing

data = fetch_california_housing(as_frame=True)
df = data.frame  # DataFrame with features + target
print(df.shape)  # (20640, 9)
df.head()

With as_frame=True we get a pandas DataFrame directly. The data dictionary:

Variable	Meaning	Unit
`MedInc`	Median income of the district	Tens of thousands of $
`HouseAge`	Median age of the houses	Years
`AveRooms`	Average rooms per household	Rooms
`AveBedrms`	Average bedrooms per household	Bedrooms
`Population`	District population	People
`AveOccup`	Average occupants per household	People
`Latitude`	District latitude	Degrees
`Longitude`	District longitude	Degrees
`MedHouseVal`	Target: median price	Hundreds of thousands of $

First things first — as always since 06-01 "Data splitting: training, validation and test" — we set the test set aside before exploring anything, so that no decision gets contaminated by information we should not see:

from sklearn.model_selection import train_test_split

X = df.drop(columns="MedHouseVal")
y = df["MedHouseVal"]

X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.2, random_state=42
)
print(len(X_train), len(X_test))  # 16512 / 4128

From here on, all EDA is done on the training set only. The test set stays under lock and key until the very end.

Exploratory data analysis (EDA)

Distributions

train = X_train.copy()
train["MedHouseVal"] = y_train

train.hist(bins=50, figsize=(14, 9))
plt.tight_layout()
plt.show()

print(train.describe().round(2))
print("Nulls per column:\n", train.isna().sum())

Key observations (compare them with what you see on screen):

There are no missing values: we are spared the strategies from 03-02 "Handling missing data" — a rare luxury in real-world data.
MedHouseVal has an artificial spike at 5.0: this is the censoring mentioned earlier. Around 800 districts in the training set are worth "$500,000 or more"; the model will not be able to distinguish anything within that group.
MedInc and the target are positively skewed (right-tailed), the pattern you learned to handle with logarithms in 03-03 "Data transformation". Here the skew of the price is moderate, so we will try with and without the transformation in the challenges; in the body of the project we work with the price on its original scale so we can interpret errors in dollars.
AveRooms, AveBedrms, Population and AveOccup have extreme outliers (districts with group housing, hotels, and so on): an AveOccup of 1,200 people per household is not a household. This will affect linear models more than tree-based ones, as you reasoned in 04-03 "Decision trees".

Correlations

As in 02-03 "Correlation and covariance", we measure the linear relationship between each feature and the target:

corr = train.corr(numeric_only=True)["MedHouseVal"].sort_values(ascending=False)
print(corr.round(3))

Typical result: MedInc dominates (≈ 0.69), followed at a great distance by AveRooms (≈ 0.15); Latitude and Longitude have small, negative linear correlations, but that does not mean they don't matter: their relationship with price is not linear, it is geographic.

The map: latitude and longitude colored by price

plt.figure(figsize=(9, 7))
sc = plt.scatter(
    train["Longitude"], train["Latitude"],
    c=train["MedHouseVal"], cmap="viridis",
    s=8, alpha=0.4
)
plt.colorbar(sc, label="Median price (x$100,000)")
plt.xlabel("Longitude")
plt.ylabel("Latitude")
plt.title("California: price is geography")
plt.show()

The picture that emerges is a map of California: high prices cluster along the coast (San Francisco Bay, Los Angeles, San Diego) and fall away inland. This spatial structure is exactly the kind of non-linear, interactive pattern (what matters is the latitude-longitude combination, not each one separately) where the trees and ensembles from module 7 should outperform linear regression. We have just formulated a testable hypothesis — that is EDA with a purpose.

Preprocessing with a Pipeline

All variables are numeric, so preprocessing is straightforward: standardize for the scale-sensitive models (linear, Ridge, Lasso — remember 03-05 "Normalization and standardization") and leave the data untouched for the tree-based ones, which are invariant to monotonic transformations. We wrap it all in a Pipeline as you learned at the end of module 3, so that scaling is fitted only on the training portion of each split, avoiding the data leakage from 06-01:

from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LinearRegression, Ridge, Lasso
from sklearn.ensemble import RandomForestRegressor, HistGradientBoostingRegressor
from sklearn.dummy import DummyRegressor

models = {
    "Dummy (mean)": DummyRegressor(strategy="mean"),
    "Linear": Pipeline([("scale", StandardScaler()), ("m", LinearRegression())]),
    "Ridge": Pipeline([("scale", StandardScaler()), ("m", Ridge(alpha=1.0))]),
    "Lasso": Pipeline([("scale", StandardScaler()), ("m", Lasso(alpha=0.001))]),
    "Random Forest": RandomForestRegressor(
        n_estimators=200, n_jobs=-1, random_state=42
    ),
    "HistGradientBoosting": HistGradientBoostingRegressor(random_state=42),
}

Notice how the dictionary already tells the story of the course: the trivial baseline, the linear regression from 04-01, the regularized versions from 07-01 "Regularization: Ridge, Lasso and Elastic Net", the Random Forest from 07-02 "Ensemble Learning" and the histogram-based gradient boosting from 07-03 "Gradient Boosting".

Modeling: from the baseline to ensembles

Why start with DummyRegressor

The baseline answers an essential question: how good is "knowing nothing"? DummyRegressor(strategy="mean") always predicts the mean price. Its RMSE is, essentially, the standard deviation of the target (≈ 1.15, that is, ±$115,000). Any model that does not beat it comfortably has learned nothing; and its figure gives scale to all the others.

Honest comparison with cross-validation

We compare all the models with 5-fold cross-validation on the training set, as you learned in 06-03 "Cross-validation": no model ever sees the test set, and every figure comes with its own variability.

from sklearn.model_selection import cross_validate, KFold

cv = KFold(n_splits=5, shuffle=True, random_state=42)
metrics = {
    "RMSE": "neg_root_mean_squared_error",
    "MAE": "neg_mean_absolute_error",
    "R2": "r2",
}

rows = []
for name, model in models.items():
    res = cross_validate(model, X_train, y_train, cv=cv, scoring=metrics)
    rows.append({
        "Model": name,
        "RMSE": -res["test_RMSE"].mean(),
        "RMSE_std": res["test_RMSE"].std(),
        "MAE": -res["test_MAE"].mean(),
        "R2": res["test_R2"].mean(),
    })

table = pd.DataFrame(rows).sort_values("RMSE").round(3)
print(table.to_string(index=False))

Ballpark results (yours will vary slightly in the decimals):

Model	RMSE (CV)	MAE (CV)	R² (CV)	Reading
HistGradientBoosting	≈ 0.47	≈ 0.31	≈ 0.83	The best: it captures the geographic non-linearity
Random Forest	≈ 0.51	≈ 0.33	≈ 0.80	Very close, slower and heavier
Linear	≈ 0.72	≈ 0.53	≈ 0.61	Respectable, but blind to interactions
Ridge	≈ 0.72	≈ 0.53	≈ 0.61	Regularizing doesn't help: there was no linear overfitting
Lasso	≈ 0.72	≈ 0.53	≈ 0.61	Ditto; with a high alpha it starts zeroing out coefficients
Dummy (mean)	≈ 1.15	≈ 0.91	≈ 0.00	The minimum bar

Three honest readings, in the spirit of 06-05 "Overfitting and underfitting":

The linear models suffer from underfitting: with 8 features and 16,500 rows, the problem is not variance but bias — the price-geography relationship is not linear. That is why Ridge and Lasso, which fight overfitting, contribute nothing here.
The EDA hypothesis is confirmed: tree-based models, able to partition the latitude-longitude space into regions, cut the RMSE by 35%.
Standard deviations matter: if two models differ by less than the variability across folds, declaring a "winner" is noise.

Error analysis

Before signing off on the best model, we look at where it fails. We train the HistGradientBoosting on the whole training set and analyze its residuals through internal validation:

from sklearn.model_selection import cross_val_predict

best = HistGradientBoostingRegressor(random_state=42)
y_pred_cv = cross_val_predict(best, X_train, y_train, cv=cv)
residuals = y_train - y_pred_cv

fig, axes = plt.subplots(1, 2, figsize=(13, 5))
axes[0].scatter(y_pred_cv, residuals, s=5, alpha=0.3)
axes[0].axhline(0, color="red", lw=1)
axes[0].set_xlabel("Prediction"); axes[0].set_ylabel("Residual")
axes[0].set_title("Residuals vs. prediction")

axes[1].scatter(X_train["Longitude"], X_train["Latitude"],
                c=residuals.abs(), cmap="Reds", s=8, alpha=0.4)
axes[1].set_title("Absolute error on the map")
plt.show()

What to look for and what usually shows up:

A diagonal band of positive residuals on the right of the first plot: these are the districts censored at 5.0. The model predicts 3.5–4.5 for districts that are "worth 5+", and it cannot do any better because the label has been clipped. That is a limit of the data, not of the model.
On the error map, the big misses concentrate in expensive coastal areas and in atypical districts (few homes, non-residential uses). Knowing where the model fails is worth more than a tenth of RMSE: it tells you what additional data to ask for.

Final evaluation on the test set

Only now, with the model chosen and no decisions left to make, do we open the test set. One single time, as you swore in 06-01:

from sklearn.metrics import root_mean_squared_error, mean_absolute_error, r2_score

best.fit(X_train, y_train)
y_pred = best.predict(X_test)

print("RMSE test:", round(root_mean_squared_error(y_test, y_pred), 3))
print("MAE  test:", round(mean_absolute_error(y_test, y_pred), 3))
print("R²   test:", round(r2_score(y_test, y_pred), 3))

You should get figures very close to the cross-validation ones (RMSE ≈ 0.46–0.48, R² ≈ 0.83). That agreement is the headline: it means the protocol was clean and the CV estimate was trustworthy. If the test result had come out clearly worse, the suspects would be data leakage or decisions overfitted to the validation.

Translated into business terms: a MAE of ≈ 0.31 means the model is off, per median district, by about $31,000 — versus the $91,000 of the baseline.

Common Mistakes and Tips

Exploring before splitting. If you compute correlations or draw the map with the full dataset, your decisions have already "seen" the test set. Split first, always.
Scaling outside the Pipeline. Fitting StandardScaler on all the data and then running CV is the classic leak from 06-01. The Pipeline does it correctly, automatically.
Ignoring the censoring of the target. The spike at 5.0 is not a market cluster: it is a dataset clipping. In a real report you would have to disclose it (and consider excluding those districts or modeling them separately).
Choosing a model on a single figure. Compare RMSE, MAE and R² together with their deviation across folds; two models that "differ" in the third decimal are the same model.
Reading R² as a hit rate. R² = 0.83 does not mean "it gets 83% right": it means it explains 83% of the variance in the price. The dollar figure (MAE/RMSE) is the one the business understands.

Challenges to go further

Interaction features. Create new variables in the style of 03-06 "Feature engineering": rooms_per_occupant = AveRooms / AveOccup, bedrooms_ratio = AveBedrms / AveRooms, or the distance from each district to the coordinates of San Francisco and Los Angeles. Hint: measure their effect on the linear regression (where they should help most) and on the boosting (which may have been discovering them on its own already).
XGBoost versus HistGradientBoosting. Install xgboost and pit it against the winner with default hyperparameters; then tune both with RandomizedSearchCV as in 07-05 "Hyperparameter optimization". Hint: start with n_estimators, learning_rate and max_depth.
Prediction intervals. A point price with no uncertainty is half an answer. Look into HistGradientBoostingRegressor(loss="quantile", quantile=0.05) and its twin with 0.95 to build a 90% interval. Hint: check what fraction of the real test prices falls inside the interval — it should be around 90%.

Conclusion

First project completed: you have walked the entire workflow that MercaFresh taught you step by step, this time on your own — defining the problem and its metric, setting the test set aside before looking at anything, exploring with hypotheses (geography rules), preprocessing inside a Pipeline, climbing from a trivial baseline to gradient boosting, comparing with honest cross-validation, analyzing residuals to learn where and why the model fails, and confirming on the test set that nothing had been contaminated. The underlying lesson: the best model did not win by magic, but because the EDA had already announced a non-linear structure that trees capture and a straight line cannot. In the next project we change terrain completely: from tables of numbers to images, where you will experience first-hand why classical methods fall short and the convolutional networks from 07-04 take over.

Project 1: Housing price prediction

Contents

Problem definition and target metric

Getting and loading the data

Exploratory data analysis (EDA)

Distributions

Correlations

The map: latitude and longitude colored by price

Preprocessing with a Pipeline

Modeling: from the baseline to ensembles

Why start with DummyRegressor

Honest comparison with cross-validation

Error analysis

Final evaluation on the test set

Common Mistakes and Tips

Challenges to go further

Conclusion

Machine Learning Course

Module 1: Introduction to Machine Learning

Module 2: Foundations of Statistics and Probability

Module 3: Data Preprocessing

Module 4: Supervised Machine Learning Algorithms

Module 5: Unsupervised Machine Learning Algorithms

Module 6: Model Evaluation and Validation

Module 7: Advanced Techniques and Optimization

Module 8: Model Implementation and Deployment

Module 9: Hands-On Projects

Module 10: Additional Resources