Real Estate Price Prediction with ML: Python Tutorial (sklearn + XGBoost)

Tested with: Python 3.12.3, scikit-learn 1.9.1, XGBoost 3.4.1, pandas 3.0.6, NumPy 2.5.3, matplotlib 3.11.2, SHAP 0.52.0 (CPU only, 4 cores). Last run 2026-09-27.

This tutorial walks through a complete ML pipeline for predicting house prices: data loading, exploratory analysis, feature engineering, training three different models, and comparing them on held-out data. We use the California Housing dataset that ships with scikit-learn (fetch_california_housing downloads it once and caches it under ~/scikit_learn_data), so there is no CSV to find. Every output below is from an actual run with the versions listed above and fixed random_state values.

pip install scikit-learn xgboost pandas numpy matplotlib shap

---

The Problem

Given features about a housing district in California (median income, house age, location, population), predict the median house value. It is a regression problem, and the data is messy in ways real data usually is: the target is right-skewed and capped at 5.0, several features overlap, and location drives price in a way no single column captures.

---

Step 1: Load the Data

import pandas as pd
import numpy as np
from sklearn.datasets import fetch_california_housing

housing = fetch_california_housing(as_frame=True)
df = housing.frame

print(df.shape)          # (20640, 9)
print(df.dtypes)       # all 9 columns are float64
print(df.describe())

Features:

• MedInc: median income in block group (in $10k)
• HouseAge: median house age
• AveRooms, AveBedrms: average rooms/bedrooms per household
• Population, AveOccup: block group population and average occupancy
• Latitude, Longitude: geographic coordinates
• MedHouseVal: target (median house value in $100k, capped at 5.0)

---

Step 2: Exploratory Data Analysis

import matplotlib.pyplot as plt

# Target distribution
plt.figure(figsize=(8, 4))
plt.hist(df["MedHouseVal"], bins=50, edgecolor="k")
plt.title("Distribution of Median House Values")
plt.xlabel("Median House Value ($100k)")
plt.ylabel("Count")
plt.tight_layout()
plt.savefig("target_distribution.png", dpi=100)

# Correlation with target
print(df.corr()["MedHouseVal"].sort_values(ascending=False))

Output:

MedHouseVal    1.000000
MedInc         0.688075
AveRooms       0.151948
HouseAge       0.105623
AveOccup      -0.023737
Population    -0.024650
Longitude     -0.045967
AveBedrms     -0.046701
Latitude      -0.144160
Name: MedHouseVal, dtype: float64

Key finding: MedInc has by far the strongest correlation with price (0.69). Latitude (-0.14) and longitude (-0.05) look weak here, but Pearson correlation only measures a straight-line relationship. Prices peak around the Bay Area and Los Angeles, not at one end of the state, so a linear coefficient misses most of the geographic signal. That is why we engineer a distance feature below.

---

Step 3: Feature Engineering

The raw columns describe households and coordinates, but price depends more on ratios and on where a district sits relative to expensive areas. We'll add four derived features:

def engineer_features(df):
    df = df.copy()

    # Rooms and bedrooms per person are more informative than totals
    df["RoomsPerPerson"] = df["AveRooms"] / df["AveOccup"]
    df["BedrmsPerRoom"] = df["AveBedrms"] / df["AveRooms"]

    # High income flag (top quartile)
    df["HighIncome"] = (df["MedInc"] > df["MedInc"].quantile(0.75)).astype(int)

    # Distance from San Francisco (a major price driver)
    sf_lat, sf_lon = 37.7749, -122.4194
    df["DistFromSF"] = np.sqrt(
        (df["Latitude"] - sf_lat) ** 2 + (df["Longitude"] - sf_lon) ** 2
    )

    return df

# Note: the HighIncome threshold is computed on the full dataset, so it sees
# the test rows. The effect here is tiny, but in production compute it on
# the training split only.
df_feat = engineer_features(df)
print(df_feat.shape)  # (20640, 13)

---

Step 4: Train/Test Split and Scaling

from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler

TARGET = "MedHouseVal"
features = [c for c in df_feat.columns if c != TARGET]

X = df_feat[features].values
y = df_feat[TARGET].values

X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.2, random_state=42
)

scaler = StandardScaler()
X_train_scaled = scaler.fit_transform(X_train)
X_test_scaled = scaler.transform(X_test)

print(f"Train: {X_train.shape}, Test: {X_test.shape}")
# Train: (16512, 12), Test: (4128, 12)

---

Step 5: Train Three Models

from sklearn.linear_model import Ridge
from sklearn.ensemble import RandomForestRegressor
import xgboost as xgb

models = {
    "Ridge Regression": Ridge(alpha=1.0),
    "Random Forest": RandomForestRegressor(
        n_estimators=200, max_features=0.5, n_jobs=-1, random_state=42
    ),
    "XGBoost": xgb.XGBRegressor(
        n_estimators=500,
        learning_rate=0.05,
        max_depth=6,
        subsample=0.8,
        colsample_bytree=0.8,
        random_state=42,
        n_jobs=-1,
    ),
}

# Ridge uses scaled features; tree models don't need scaling
fit_data = {
    "Ridge Regression": (X_train_scaled, X_test_scaled),
    "Random Forest": (X_train, X_test),
    "XGBoost": (X_train, X_test),
}

for name, model in models.items():
    X_tr, _ = fit_data[name]
    model.fit(X_tr, y_train)
    print(f"Trained {name}")

---

Step 6: Evaluate and Compare

from sklearn.metrics import mean_squared_error, r2_score

print(f"\n{'Model':<22} {'RMSE':>8} {'R²':>8}")
print("-" * 42)

results = {}
for name, model in models.items():
    _, X_te = fit_data[name]
    preds = model.predict(X_te)
    rmse = np.sqrt(mean_squared_error(y_test, preds))
    r2 = r2_score(y_test, preds)
    results[name] = {"rmse": rmse, "r2": r2}
    print(f"{name:<22} {rmse:>8.3f} {r2:>8.3f}")

Output:

Model                      RMSE       R²
------------------------------------------
Ridge Regression          0.674    0.653
Random Forest             0.486    0.820
XGBoost                   0.440    0.852

RMSE is in units of $100k, so XGBoost's error of about $44k is meaningfully better than Ridge's $67k.

---

Step 7: Feature Importances (XGBoost)

xgb_model = models["XGBoost"]
importances = sorted(
    zip(features, xgb_model.feature_importances_),
    key=lambda x: -x[1]
)

print("\nTop features (XGBoost):")
for name, imp in importances[:6]:
    bar = "█" * int(imp * 100)
    print(f"  {name:<20} {imp:.3f}  {bar}")

Output:

Top features (XGBoost):
  HighIncome           0.413  █████████████████████████████████████████
  MedInc               0.248  ████████████████████████
  AveOccup             0.073  ███████
  DistFromSF           0.061  ██████
  Longitude            0.047  ████
  Latitude             0.046  ████

Income dominates, but notice how it is split. feature_importances_ on an XGBRegressor reports gain by default, and HighIncome is just a thresholded copy of MedInc. The early splits that separate high-income districts happen to land on the binary flag, so it collects the largest share of gain even though it adds almost no new information. Read the two rows together as "income". The engineered DistFromSF made the top 4, ahead of raw latitude and longitude, which supports the geographic feature. RoomsPerPerson did not make the top 6. Gain importances are easily skewed by correlated features like this, which is one more reason to check SHAP values below.

---

Which Model to Use and Why

Model	RMSE	R²	Training time (4-core CPU)	Explainability
Ridge Regression	0.674	0.65	<0.01s	High (coefficients)
Random Forest	0.486	0.82	~3.5s	Medium (importances)
XGBoost	0.440	0.85	~0.8s	Medium (SHAP values)

Use Ridge when you need a fast baseline or need to explain the model coefficient-by-coefficient to a stakeholder.

Use Random Forest when you want solid accuracy without hyperparameter tuning.

Use XGBoost when accuracy matters most and you're willing to tune learning rate, depth, and subsampling. Add SHAP for explainability:

import shap

explainer = shap.TreeExplainer(xgb_model)
shap_values = explainer.shap_values(X_test[:100])
shap.summary_plot(shap_values, X_test[:100], feature_names=features)

---

Related articles

• Top 10 Machine Learning Algorithms Every Practitioner Should Know
• Fix sklearn's "Input Contains NaN" in Pipelines
• Fixing pandas Datetime Parsing Errors