Dataset
   ↓
Data Exploration
   ↓
Feature Engineering
   ↓
Train/Test Split
   ↓
Preprocessing Pipeline
   ↓
Multiple Regression Models
   ↓
Model Comparison
   ↓
Hyperparameter Tuning
   ↓
Final Evaluation
   ↓
Evaluation Plots
   ↓
Save Model

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

from sklearn.datasets import fetch_california_housing

from sklearn.model_selection import train_test_split, GridSearchCV

from sklearn.pipeline import Pipeline
from sklearn.compose import ColumnTransformer

from sklearn.preprocessing import StandardScaler
from sklearn.impute import SimpleImputer

from sklearn.base import BaseEstimator, TransformerMixin

from sklearn.linear_model import LinearRegression
from sklearn.ensemble import (
    RandomForestRegressor,
    GradientBoostingRegressor
)

from sklearn.metrics import (
    mean_absolute_error,
    mean_squared_error,
    r2_score
)

import joblib

housing = fetch_california_housing(as_frame=True)

df = housing.frame

print("Dataset shape:", df.shape)

Dataset shape: (20640, 9)

df.head()

df.info()

Rows
Columns
Data types
Missing values

df.describe()

print(df.isnull().sum())

SimpleImputer(strategy="median")

X = df.drop("MedHouseVal", axis=1)

y = df["MedHouseVal"]

print("Features shape:", X.shape)
print("Target shape:", y.shape)

X → Features / Inputs
y → Target / Output

X = 8 features
y = MedHouseVal

plt.figure(figsize=(8, 5))

plt.hist(y, bins=50)

plt.xlabel("Median House Value")
plt.ylabel("Frequency")
plt.title("Distribution of House Values")

plt.show()

RoomsPerOccupant
BedroomsPerRoom
PopulationPerRoom

class FeatureEngineering(BaseEstimator, TransformerMixin):

    def fit(self, X, y=None):
        return self

    def transform(self, X):
        X = X.copy()

        X["RoomsPerOccupant"] = (
            X["AveRooms"] / X["AveOccup"]
        )

        X["BedroomsPerRoom"] = (
            X["AveBedrms"] / X["AveRooms"]
        )

        X["PopulationPerRoom"] = (
            X["Population"] / X["AveRooms"]
        )

        return X

Raw Data
   ↓
Feature Engineering
   ↓
Preprocessing
   ↓
Model

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

print("Training data:", X_train.shape)
print("Testing data:", X_test.shape)

80% → Training
20% → Testing

feature_engineering = FeatureEngineering()

X_train_engineered = feature_engineering.fit_transform(
    X_train
)

print(X_train_engineered.columns.tolist())

numeric_features = X_train_engineered.columns.tolist()

print(numeric_features)

MedInc
HouseAge
AveRooms
AveBedrms
Population
AveOccup
Latitude
Longitude
RoomsPerOccupant
BedroomsPerRoom
PopulationPerRoom

numeric_pipeline = Pipeline([
    ("imputer", SimpleImputer(strategy="median")),
    ("scaler", StandardScaler())
])

Missing Values
      ↓
Median Imputation
      ↓
StandardScaler

preprocessor = ColumnTransformer([
    (
        "numeric",
        numeric_pipeline,
        numeric_features
    )
])

linear_pipeline = Pipeline([
    ("feature_engineering", FeatureEngineering()),
    ("preprocessing", preprocessor),
    ("model", LinearRegression())
])

Raw Data
   ↓
Feature Engineering
   ↓
Preprocessing
   ↓
Linear Regression

linear_pipeline.fit(
    X_train,
    y_train
)

Housing Features → House Value

linear_predictions = linear_pipeline.predict(
    X_test
)

y_test

linear_mae = mean_absolute_error(
    y_test,
    linear_predictions
)

linear_rmse = np.sqrt(
    mean_squared_error(
        y_test,
        linear_predictions
    )
)

linear_r2 = r2_score(
    y_test,
    linear_predictions
)

print("Linear Regression")
print("-----------------")
print("MAE :", linear_mae)
print("RMSE:", linear_rmse)
print("R²  :", linear_r2)

rf_pipeline = Pipeline([
    ("feature_engineering", FeatureEngineering()),
    ("preprocessing", preprocessor),
    (
        "model",
        RandomForestRegressor(
            n_estimators=100,
            random_state=42,
            n_jobs=-1
        )
    )
])

rf_pipeline.fit(
    X_train,
    y_train
)

rf_predictions = rf_pipeline.predict(
    X_test
)

rf_mae = mean_absolute_error(
    y_test,
    rf_predictions
)

rf_rmse = np.sqrt(
    mean_squared_error(
        y_test,
        rf_predictions
    )
)

rf_r2 = r2_score(
    y_test,
    rf_predictions
)

print("Random Forest")
print("-------------")
print("MAE :", rf_mae)
print("RMSE:", rf_rmse)
print("R²  :", rf_r2)

gradient_pipeline = Pipeline([
    ("feature_engineering", FeatureEngineering()),
    ("preprocessing", preprocessor),
    (
        "model",
        GradientBoostingRegressor(
            random_state=42
        )
    )
])

gradient_pipeline.fit(
    X_train,
    y_train
)

gradient_predictions = gradient_pipeline.predict(
    X_test
)

gradient_mae = mean_absolute_error(
    y_test,
    gradient_predictions
)

gradient_rmse = np.sqrt(
    mean_squared_error(
        y_test,
        gradient_predictions
    )
)

gradient_r2 = r2_score(
    y_test,
    gradient_predictions
)

print("Gradient Boosting")
print("-----------------")
print("MAE :", gradient_mae)
print("RMSE:", gradient_rmse)
print("R²  :", gradient_r2)

results = pd.DataFrame({
    "Model": [
        "Linear Regression",
        "Random Forest",
        "Gradient Boosting"
    ],
    "MAE": [
        linear_mae,
        rf_mae,
        gradient_mae
    ],
    "RMSE": [
        linear_rmse,
        rf_rmse,
        gradient_rmse
    ],
    "R2": [
        linear_r2,
        rf_r2,
        gradient_r2
    ]
})

results

MAE  → Lower is better
RMSE → Lower is better
R²   → Higher is better

Random Forest

param_grid = {
    "model__n_estimators": [
        100,
        200
    ],
    "model__max_depth": [
        10,
        20,
        None
    ],
    "model__min_samples_split": [
        2,
        5
    ]
}

model__parameter

model__n_estimators

Pipeline
   ↓
model
   ↓
n_estimators

grid_search = GridSearchCV(
    estimator=rf_pipeline,
    param_grid=param_grid,
    cv=5,
    scoring="neg_root_mean_squared_error",
    n_jobs=-1,
    verbose=1
)

grid_search.fit(
    X_train,
    y_train
)

Parameter Combination
        ↓
5-Fold Cross Validation
        ↓
Calculate Score
        ↓
Compare Scores
        ↓
Select Best Combination

print("Best Parameters:")
print(grid_search.best_params_)

{
    'model__max_depth': 20,
    'model__min_samples_split': 2,
    'model__n_estimators': 200
}

print(
    "Best CV RMSE:",
    -grid_search.best_score_
)

best_model = grid_search.best_estimator_

print(best_model)

Feature Engineering
        ↓
Preprocessing
        ↓
Tuned Random Forest

final_predictions = best_model.predict(
    X_test
)

final_mae = mean_absolute_error(
    y_test,
    final_predictions
)

print("Final MAE:", final_mae)

final_rmse = np.sqrt(
    mean_squared_error(
        y_test,
        final_predictions
    )
)

final_r2 = r2_score(
    y_test,
    final_predictions
)

print("Final RMSE:", final_rmse)
print("Final R²:", final_r2)

print("Final Model Performance")
print("=======================")
print(f"MAE : {final_mae:.4f}")
print(f"RMSE: {final_rmse:.4f}")
print(f"R²  : {final_r2:.4f}")

MAE  → Average absolute error
RMSE → Penalizes larger errors
R²   → Explained variance

plt.figure(figsize=(8, 6))

plt.scatter(
    y_test,
    final_predictions,
    alpha=0.4
)

plt.plot(
    [y_test.min(), y_test.max()],
    [y_test.min(), y_test.max()],
    linestyle="--"
)

plt.xlabel("Actual House Value")
plt.ylabel("Predicted House Value")
plt.title("Actual vs Predicted House Values")

plt.tight_layout()

plt.savefig(
    "actual_vs_predicted.png",
    dpi=300
)

plt.show()

residuals = y_test - final_predictions

print(residuals.head())

Residual = Actual - Predicted

plt.figure(figsize=(8, 6))

plt.scatter(
    final_predictions,
    residuals,
    alpha=0.4
)

plt.axhline(
    0,
    linestyle="--"
)

plt.xlabel("Predicted House Value")
plt.ylabel("Residual")
plt.title("Residual Plot")

plt.tight_layout()

plt.savefig(
    "residual_plot.png",
    dpi=300
)

plt.show()

Residual
   |
 + | .   .    .
 0 |----------------
 - |   .   .    .
   |
   +----------------
       Predicted

plt.figure(figsize=(9, 6))

plt.bar(
    results["Model"],
    results["RMSE"]
)

plt.xlabel("Model")
plt.ylabel("RMSE")
plt.title("Model Comparison - RMSE")

plt.xticks(rotation=15)

plt.tight_layout()

plt.savefig(
    "model_comparison_rmse.png",
    dpi=300
)

plt.show()

plt.figure(figsize=(9, 6))

plt.bar(
    results["Model"],
    results["R2"]
)

plt.xlabel("Model")
plt.ylabel("R²")
plt.title("Model Comparison - R²")

plt.xticks(rotation=15)

plt.tight_layout()

plt.savefig(
    "model_comparison_r2.png",
    dpi=300
)

plt.show()

joblib.dump(
    best_model,
    "housing_model.pkl"
)

print("Model saved successfully.")

housing_model.pkl

loaded_model = joblib.load(
    "housing_model.pkl"
)

print("Model loaded successfully.")

loaded_predictions = loaded_model.predict(
    X_test
)

loaded_r2 = r2_score(
    y_test,
    loaded_predictions
)

print("Loaded Model R²:", loaded_r2)

print(
    np.allclose(
        final_predictions,
        loaded_predictions
    )
)

True

final_results = pd.DataFrame({
    "Metric": [
        "MAE",
        "RMSE",
        "R²"
    ],
    "Score": [
        final_mae,
        final_rmse,
        final_r2
    ]
})

final_results

## Conclusion

This project implemented an end-to-end machine learning workflow for California house value prediction.

### Workflow

1. Loaded the California Housing dataset.
2. Explored the dataset and checked for missing values.
3. Created additional features using feature engineering.
4. Split the data into training and testing sets.
5. Built preprocessing pipelines using imputation and standardization.
6. Trained Linear Regression, Random Forest, and Gradient Boosting models.
7. Compared the models using MAE, RMSE, and R².
8. Selected Random Forest as the candidate model for tuning.
9. Used GridSearchCV with 5-fold cross-validation to find better hyperparameters.
10. Evaluated the final model on the unseen test set.
11. Created actual-vs-predicted and residual plots.
12. Saved the complete trained pipeline using Joblib.

### Final Model

The final model is saved as:

`housing_model.pkl`

The saved pipeline contains feature engineering, preprocessing, and the trained regression model, allowing the complete workflow to be reused for future predictions.

housing_pipeline.ipynb
housing_model.pkl
actual_vs_predicted.png
residual_plot.png
model_comparison_rmse.png
model_comparison_r2.png

pip install numpy pandas matplotlib scikit-learn joblib

jupyter notebook housing_pipeline.ipynb