Predict Housing Prices¶
The goal here is to build a machine learning model to predict housing prices in California using the California Census Data. The data has metrics such as population, median income, median housing prices, and so on. The data was mostly raw and lacked any form of preprocessing. Hence a lot of time was invested for processing and cleaning the data thus getting it ready for a mathematical model to gain insights. The data was tested with linear-regression, decision-tree-regressor and random-forest-regressor.¶
Importing few common libraries, ensure matplotlib plots figures inline and prepare a function to save the figures:¶
In [1]:
import numpy as np
import os
np.random.seed(42)
%matplotlib inline
import matplotlib
import matplotlib.pyplot as plt
plt.rcParams['axes.labelsize'] = 14
plt.rcParams['xtick.labelsize'] = 12
plt.rcParams['ytick.labelsize'] = 12
# Where to save the figures
PROJECT_ROOT_DIR = "."
CHAPTER_ID = "end_to_end_project"
IMAGES_PATH = os.path.join(PROJECT_ROOT_DIR, "images", CHAPTER_ID)
def save_fig(fig_id, tight_layout=True, fig_extension="png", resolution=300):
path = os.path.join(IMAGES_PATH, fig_id + "." + fig_extension)
print("Saving figure", fig_id)
if tight_layout:
plt.tight_layout()
plt.savefig(path, format=fig_extension, dpi=resolution)
import warnings
warnings.filterwarnings(action="ignore", module="scipy", message="^internal gelsd")
Function to get the data from the server.¶
This function should only be called when new data is required to be fetched from the server
In [2]:
import os
import tarfile
from six.moves import urllib
DOWNLOAD_ROOT = "https://raw.githubusercontent.com/ageron/handson-ml/master/"
HOUSING_PATH = os.path.join("datasets", "housing")
HOUSING_URL = DOWNLOAD_ROOT + "datasets/housing/housing.tgz"
def fetch_housing_data(housing_url= HOUSING_URL, housing_path= HOUSING_PATH):
if not os.path.isdir(housing_path):
os.makedirs(housing_path)
tgz_path = os.path.join(housing_path, "housing.tgz")
urllib.request.urlretrieve(housing_url, tgz_path)
housing_tgz = tarfile.open(tgz_path)
housing_tgz.extractall(path= housing_path)
housing_tgz.close()
In [3]:
#fetch_housing_data()
Reading the raw data¶
Read and study the unprocessed data. By plotting various graphs of the data, we can gain some insights. Plotting various figures will help us decide which features have will an impact on the housing prices. This will further allow us to determine which features should be removed and which new features to create, using the existing features.
In [4]:
import pandas as pd
def load_housing_data(housing_path= HOUSING_PATH):
csv_path = os.path.join(housing_path, "housing.csv")
return pd.read_csv(csv_path)
In [5]:
housing = load_housing_data()
housing.head()
Out[5]:
In [6]:
housing.info()
In [7]:
housing["ocean_proximity"].value_counts()
Out[7]:
In [8]:
housing.describe()
Out[8]:
In [9]:
%matplotlib inline
import matplotlib.pyplot as plt
housing.hist(bins= 50, figsize=(20,15))
save_fig("attribute_histogram_plots")
plt.show()
In [10]:
from sklearn.model_selection import train_test_split
train_set, test_set = train_test_split(housing, test_size=0.2, random_state=42)
In [11]:
test_set.head()
Out[11]:
In [12]:
housing["median_income"].hist()
Out[12]:
In [13]:
#Divide by 1.5 to limit the no. of income categories
housing["income_cat"] = np.ceil(housing["median_income"] / 1.5)
#Label those above 5 as 5
housing["income_cat"].where(housing["income_cat"] < 5, 5.0, inplace=True) #Where cond is True, keep the original value. Where False, replace with corresponding value from other.
In [14]:
housing["income_cat"].value_counts()
Out[14]:
In [15]:
housing["income_cat"].hist()
Out[15]:
In [16]:
from sklearn.model_selection import StratifiedShuffleSplit
split = StratifiedShuffleSplit(n_splits= 1, test_size= 0.2, random_state= 42)
for train_index, test_index in split.split(housing, housing["income_cat"]):
strat_train_set = housing.loc[train_index]
strat_test_set = housing.loc[test_index]
In [17]:
strat_test_set["income_cat"].value_counts() / len(strat_test_set)
Out[17]:
In [18]:
housing["income_cat"].value_counts() / len(housing)
Out[18]:
In [19]:
for set_ in (strat_train_set, strat_test_set):
set_.drop("income_cat",axis=1, inplace=True)
Discover and Visualize the data to gain insights¶
In [20]:
housing = strat_train_set.copy() # make a copy to be on safe side
In [21]:
housing.plot(kind= "scatter", x="longitude", y="latitude", alpha=0.1)
save_fig("better_visualization_plot")
In [22]:
housing.plot(kind="scatter", x="longitude", y="latitude", alpha=0.4,
s=housing["population"]/100, label="population", figsize=(10,7),
c="median_house_value", cmap=plt.get_cmap("jet"), colorbar=True,
sharex=False)
plt.legend()
save_fig("housing_prices_scatterplot")
In [23]:
import matplotlib.image as mpimg
california_img=mpimg.imread(PROJECT_ROOT_DIR + '/images/end_to_end_project/california.png')
ax = housing.plot(kind="scatter", x="longitude", y="latitude", figsize=(10,7),
s=housing['population']/100, label="Population",
c="median_house_value", cmap=plt.get_cmap("jet"),
colorbar=False, alpha=0.4,
)
plt.imshow(california_img, extent=[-124.55, -113.80, 32.45, 42.05], alpha=0.5,
cmap=plt.get_cmap("jet"))
plt.ylabel("Latitude", fontsize=14)
plt.xlabel("Longitude", fontsize=14)
prices = housing["median_house_value"]
tick_values = np.linspace(prices.min(), prices.max(), 11)
cbar = plt.colorbar()
cbar.ax.set_yticklabels(["$%dk"%(round(v/1000)) for v in tick_values], fontsize=14)
cbar.set_label('Median House Value', fontsize=16)
plt.legend(fontsize=16)
save_fig("california_housing_prices_plot")
plt.show()
Looking for Correlations¶
In [24]:
corr_matrix = housing.corr()
In [25]:
corr_matrix["median_house_value"].sort_values(ascending=False)
Out[25]:
In [26]:
from pandas.plotting import scatter_matrix
attributes = ["median_house_value", "median_income", "total_rooms",
"housing_median_age"]
scatter_matrix(housing[attributes], figsize=(12, 8))
save_fig("scatter_matrix_plot")
In [27]:
housing.plot(kind="scatter", x="median_income", y="median_house_value",
alpha=0.1)
plt.axis([0, 16, 0, 550000])
save_fig("income_vs_house_value_scatterplot")
Experimenting with Attribute Combinations¶
Drawing out new features from the existing features may help us to create new meaningful data
In [28]:
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"]
In [29]:
corr_matrix = housing.corr()
corr_matrix["median_house_value"].sort_values(ascending=False)
Out[29]:
rooms_per_household:
In [30]:
housing.plot(kind="scatter", x="rooms_per_household", y="median_house_value",
alpha=0.2)
plt.axis([0, 5, 0, 520000])
plt.axis()
Out[30]:
bedrooms_per_room:¶
In [31]:
housing.plot(kind="scatter", x="bedrooms_per_room", y="median_house_value",
alpha=0.2)
plt.axis([0, 5, 0, 520000])
plt.axis()
Out[31]:
population_per_household¶
In [32]:
housing.plot(kind="scatter", x="population_per_household", y="median_house_value",
alpha=0.2)
plt.axis([0, 5, 0, 520000])
plt.axis()
Out[32]:
In [33]:
housing.describe()
Out[33]:
Prepare the data for ML Algorithms¶
In [34]:
housing = strat_train_set.drop("median_house_value", axis=1) # drop labels for training data
housing_labels = strat_train_set["median_house_value"].copy()
In [35]:
housing_labels.head()
Out[35]:
In [36]:
housing.head()
Out[36]:
In [37]:
from sklearn.preprocessing import Imputer
imputer = Imputer(strategy="median")
In [38]:
housing_num = housing.drop("ocean_proximity", axis=1)
housing_num.head()
Out[38]:
In [39]:
imputer.fit(housing_num)
Out[39]:
In [40]:
imputer.statistics_
Out[40]:
In [41]:
housing_num.median().values
Out[41]:
Transforming the training set
In [42]:
X = imputer.transform(housing_num)
In [43]:
housing_tr = pd.DataFrame(X, columns=housing_num.columns)
housing_tr.head()
Out[43]:
Handling Text and Categorical Attributes¶
In [44]:
housing_cat = housing["ocean_proximity"]
housing_cat.head(10)
Out[44]:
In [45]:
from sklearn.preprocessing import LabelEncoder
encoder = LabelEncoder()
housing_cat_encoded = encoder.fit_transform(housing_cat)
housing_cat_encoded[:10]
Out[45]:
In [46]:
from sklearn.preprocessing import LabelBinarizer
encoder = LabelBinarizer()
housing_cat_1hot = encoder.fit_transform(housing_cat)
housing_cat_1hot
Out[46]:
Custom Transformers¶
In [47]:
from sklearn.base import BaseEstimator, TransformerMixin
#column index
rooms_ix, bedrooms_ix, population_ix, household_ix = 3, 4, 5, 6
class CombinedAttributesAdder(BaseEstimator, TransformerMixin):
def __init__(self, add_bedrooms_per_room = True):
self.add_bedrooms_per_room = add_bedrooms_per_room
def fit(self, X, y=None):
return self #nothing to do
def transform(self, X, y=None):
rooms_per_household = X[:, rooms_ix] / X[:, household_ix]
population_per_household = X[:, population_ix] / X[:, household_ix]
if self.add_bedrooms_per_room:
bedrooms_per_room = X[:, bedrooms_ix] / X[:, rooms_ix]
return np.c_[X, rooms_per_household, population_per_household,
bedrooms_per_room]
else:
return np.c_[X, rooms_per_household, population_per_household]
In [48]:
attr_adder = CombinedAttributesAdder(add_bedrooms_per_room=False)
housing_extra_attribs = attr_adder.transform(housing.values)
housing_extra_attribs = pd.DataFrame(
housing_extra_attribs,
columns=list(housing.columns)+["rooms_per_household","population_per_household"])
housing_extra_attribs.head()
Out[48]:
Building a Pipeline for preprocessing the numerical attributes:
In [49]:
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler
num_pipeline = Pipeline([
('imputer', Imputer(strategy="median")),
('attribs_adder', CombinedAttributesAdder()),
('std_scaler', StandardScaler()),
])
housing_num_tr = num_pipeline.fit_transform(housing_num)
In [50]:
housing_num_tr
Out[50]:
In [51]:
#Create a class to select numerical or categorical columns
#since scikit-learn doesn't handle Dataframes yet
class DataFrameSelector(BaseEstimator, TransformerMixin):
def __init__(self, attribute_names):
self.attribute_names = attribute_names
def fit(self, X, y=None):
return self
def transform(self, X):
return X[self.attribute_names].values
Note: The following code was copied from https://github.com/scikit-learn/scikit-learn/pull/7375/files#diff-1e175ddb0d84aad0a578d34553f6f9c6 since LabelBinarizer only take 2 arguments (self and labels), but we need to send 3 arguments (self, features and labels)¶
There are 3 ways to solve this issue:¶
1) use scikit-learn 0.20 OneHotEncoder which supports string features (still in development at the time of writing this code)¶
2) use scikit-learn 0.18¶
3) use custom LabelBinarizer which can handle 3 arguments¶
I have gone with the 3rd option, as the code was readily available on the scikit-learn github repo.¶
In [52]:
class LabelBinarizerPipelineFriendly(LabelBinarizer):
def fit(self, X, y=None):
"""this would allow us to fit the model based on the X input."""
super(LabelBinarizerPipelineFriendly, self).fit(X)
def transform(self, X, y=None):
return super(LabelBinarizerPipelineFriendly, self).transform(X)
def fit_transform(self, X, y=None):
return super(LabelBinarizerPipelineFriendly, self).fit(X).transform(X)
In [53]:
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler
num_attribs = list(housing_num)
cat_attribs = ["ocean_proximity"]
num_pipeline = Pipeline([
('selector', DataFrameSelector(num_attribs)),
('imputer', Imputer(strategy="median")),
('attribs_adder', CombinedAttributesAdder()),
('std_scaler', StandardScaler()),
])
cat_pipeline = Pipeline([
('selector', DataFrameSelector(cat_attribs)),
('label_binarizer_pipeline_friendly', LabelBinarizerPipelineFriendly()),
])
In [54]:
from sklearn.pipeline import FeatureUnion
full_pipeline = FeatureUnion(transformer_list = [
("num_pipeline", num_pipeline),
("cat_pipeline", cat_pipeline),
])
In [55]:
housing_prepared = full_pipeline.fit_transform(housing)
housing_prepared
Out[55]:
In [56]:
housing_prepared.shape
Out[56]:
Select and train a model¶
Trying Linear Regression without cross-validation¶
In [57]:
#training
from sklearn.linear_model import LinearRegression
lin_reg = LinearRegression()
lin_reg.fit(housing_prepared, housing_labels)
Out[57]:
In [58]:
# trying full pipeline on few training instances
some_data = housing.iloc[:5]
some_labels = housing_labels.iloc[:5]
some_data_prepared = full_pipeline.transform(some_data)
print("Prediction:", lin_reg.predict(some_data_prepared))
Compare aganist the actual value
In [59]:
print("Labels:", list(some_labels))
In [60]:
some_data_prepared
Out[60]:
In [61]:
# Calculating RMS for the entire dataset for Linear Regression
from sklearn.metrics import mean_squared_error
housing_predictions = lin_reg.predict(housing_prepared)
lin_mse = mean_squared_error(housing_labels, housing_predictions)
lin_rmse = np.sqrt(lin_mse)
lin_rmse
Out[61]:
Trying DecisionTreeRegressor without cross-validation¶
In [62]:
from sklearn.tree import DecisionTreeRegressor
tree_reg = DecisionTreeRegressor()
tree_reg.fit(housing_prepared, housing_labels)
Out[62]:
In [63]:
housing_predictions = tree_reg.predict(housing_prepared)
tree_mse = mean_squared_error(housing_labels, housing_predictions)
tree_rmse = np.sqrt(tree_mse)
tree_rmse
Out[63]:
Clearly the DecisionTreeRegressor has badly overfitted the data. Hence, now I will be using cross-validation and again use DecisionTreeRegressor
Trying DecisionTreeRegressor with cross-validation¶
In [64]:
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)
In [65]:
def display_scores(scores):
print("Scores:", scores)
print("Mean:", scores.mean())
print("Standard Deviation:", scores.std())
display_scores(tree_rmse_scores)
Trying LinearRegression with cross-validation¶
In [66]:
lin_scores = cross_val_score(lin_reg, housing_prepared, housing_labels,
scoring="neg_mean_squared_error", cv=10)
lin_rmse_scores = np.sqrt(-lin_scores)
display_scores(lin_rmse_scores)
Trying RandomForestRegressor with cross-validation¶
In [67]:
from sklearn.ensemble import RandomForestRegressor
forest_reg = RandomForestRegressor()
forest_reg.fit(housing_prepared, housing_labels)
Out[67]:
In [68]:
forest_scores = cross_val_score(forest_reg, housing_prepared, housing_labels,
scoring = "neg_mean_squared_error", cv=10)
forest_rmse_scores = np.sqrt(-forest_scores)
display(forest_rmse_scores)
Trying Support Vector Regressor without cross-validation¶
In [69]:
from sklearn.svm import SVR
svm_reg = SVR(kernel="linear")
svm_reg.fit(housing_prepared, housing_labels)
housing_predictions = svm_reg.predict(housing_prepared)
svm_mse = mean_squared_error(housing_labels, housing_predictions)
svm_rmse = np.sqrt(svm_mse)
svm_rmse
Out[69]:
Fine tuning the models¶
After that we have shortlisted some of the models, the next goal is to fine-tune them. One way to do that is to fiddle with the hyperparameters manually until we find a great combination of hyperparameter values. But that is a tedious and time consuming process. Instead we can take help of Scikit-Learn's:
1) GridSearchCV
2) RandomizedSearchCV
Grid Search¶
In [70]:
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)
Out[70]:
The best hyperparameter combination found:¶
In [71]:
grid_search.best_params_
Out[71]:
In [72]:
grid_search.best_estimator_
Out[72]:
Looking at the score of each hyperparameter combination tested during the grid search:¶
In [73]:
cvres = grid_search.cv_results_
for mean_score, params in zip(cvres["mean_test_score"], cvres["params"]):
print(np.sqrt(-mean_score), params)
In [74]:
pd.DataFrame(grid_search.cv_results_)
Out[74]:
Randomized Search¶
In [75]:
from sklearn.model_selection import RandomizedSearchCV
from scipy.stats import randint
param_distribs = {
'n_estimators': randint(low=1, high=200),
'max_features': randint(low=1, high=8),
}
forest_reg = RandomForestRegressor(random_state=42)
rnd_search = RandomizedSearchCV(forest_reg, param_distributions=param_distribs,
n_iter=10, cv=5, scoring='neg_mean_squared_error', random_state=42)
rnd_search.fit(housing_prepared, housing_labels)
Out[75]:
In [76]:
cvres = rnd_search.cv_results_
for mean_score, params in zip(cvres["mean_test_score"], cvres["params"]):
print(np.sqrt(-mean_score), params)
Analysis the Best Models and Their Errors¶
In [77]:
feature_importances = grid_search.best_estimator_.feature_importances_
feature_importances
Out[77]:
Evaluate Our System on the Test Set¶
In [78]:
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_mse = mean_squared_error(y_test, final_predictions)
final_rmse = np.sqrt(final_mse)
In [79]:
final_rmse # Final RMS Error
Out[79]: