<jupyter_start><jupyter_text>Housing Price Prediction Project<NAME> Working through Ch2 on my own.<jupyter_code>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()
    print("Fetched")<jupyter_output><empty_output><jupyter_text>Import Data<jupyter_code>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)
fetch_housing_data()
housing = load_housing_data()
housing.head()
housing.info()
housing.describe()
housing["ocean_proximity"].value_counts()
%matplotlib inline #ipython magic function
import matplotlib.pyplot as plt
housing.hist(bins=50, figsize=(20,15))
plt.show()<jupyter_output><empty_output><jupyter_text>Set aside the test setAnd don't look at it again until evaluation time. (Avoid the data snooping bias)This shuffle method works ok for now, but if the data set is evolving over time we can't adapt to the new data and keep our existing test/train sets. A better way might be to use last digit of the hash of an ID column as a definition of test vs train membership. Then sort rows based on this membership test.<jupyter_code>import numpy as np
def split_train_test(data, test_ratio):
    np.random.seed(42) #so we shuffle in a known repeatable way.
    shuffled_indices= np.random.permutation(len(data))
    test_set_size = int(len(data) * test_ratio)
    test_indices = shuffled_indices[:test_set_size]
    train_indices = shuffled_indices[test_set_size]
    return data.iloc[train_indices], data.iloc[test_indices]
random_train_set, random_test_set = split_train_test(housing, 0.2)
#Or we can use sci kit to just do this for us.
from sklearn.model_selection import train_test_split
train_set_sk, test_set_sk = train_test_split(housing, test_size=0.2, random_state=42)<jupyter_output><empty_output><jupyter_text>Avoiding Sampling BiasIf you randomly sample, you may not be getting a representive slice of the data in your test and training sets. Use "Stratified Sampling" to avoid this.Don't have too many strata, and make sure each is large enough.<jupyter_code>#Limit and discretize income categories
housing["income_cat"] = np.ceil(housing["median_income"] / 1.5) 
# Label those above 5 as 5
# Where the predicate test is NOT true, then we replace with the value. Confusing.
housing["income_cat"].where(housing["income_cat"] < 5, 5.0, inplace=True)

#Compare input vs output after we setup our strata
housing["income_cat"].value_counts()
#Compare input vs output after we setup our strata
housing["median_income"].hist()
housing["income_cat"].hist()
#time to do the stratified sampling
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]<jupyter_output><empty_output><jupyter_text>Compare sampling error of Strat vs Random<jupyter_code>#how did the strat split do, vs random.
sample_stats = pd.DataFrame()
sample_stats["overall"] = housing["income_cat"].value_counts() / len(housing)
sample_stats["strat"] = strat_train_set["income_cat"].value_counts() / len(strat_train_set)

#Rerun random sample, now that housing has income_cat
random_train_set, random_test_set = train_test_split(housing, test_size=0.2, random_state=42)
sample_stats["random"] = random_train_set["income_cat"].value_counts() / len(random_train_set)

sample_stats.sort_index()
#Compute error percentages
sample_stats["Random  % Error"] = (sample_stats["overall"] - sample_stats["random"]) / sample_stats["overall"] * 100
sample_stats["Strat  % Error"] = (sample_stats["overall"] - sample_stats["strat"]) / sample_stats["overall"] * 100

sample_stats<jupyter_output><empty_output><jupyter_text>Clearly, Strat has lower error percentages, so it is a better representation of the whole set.<jupyter_code>#Finally, cleanup income_cat off the test/train data.
for set_ in (strat_train_set, strat_test_set):
    set_.drop("income_cat", axis=1, inplace=True)<jupyter_output><empty_output><jupyter_text>Discover and Visualize the Data to Gain InsightsNow that we've prepped the test and train sets, we can dig into it to see what we can learn.At this point you are only to look at the training set (ignore the test set!). If the training set is large, you can sample it to make exploration more efficient. Ours is small though.<jupyter_code>#make a copy of train set so we don't mangle it.
housing = strat_train_set.copy()
#rough plot
housing.plot(kind="scatter", x="longitude", y="latitude")
#Better plot with alpha, to visualize density
housing.plot(kind="scatter", x="longitude", y="latitude", alpha=0.1)<jupyter_output><empty_output><jupyter_text>Keep improving the visualization, maybe patterns will emerge<jupyter_code>#s is radius
#c is color
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)<jupyter_output><empty_output><jupyter_text>Looking For Correlations<jupyter_code>corr_matrix = housing.corr() #Pearson's r, or standard correlation coeff
corr_matrix["median_house_value"].sort_values(ascending=False)<jupyter_output><empty_output><jupyter_text>Close to 1 means strongly correlated, (as opposed to -1)<jupyter_code>#using scatter matrix plot
from pandas.plotting import scatter_matrix

attributes = ["median_house_value", "median_income", "total_rooms", "housing_median_age"]
scatter_matrix(housing[attributes], figsize=(12,8))<jupyter_output><empty_output><jupyter_text>Median income seems like the most promising correlation, so lets zoom in on that<jupyter_code>housing.plot(kind="scatter", x="median_income", y="median_house_value", alpha=0.1)<jupyter_output><empty_output><jupyter_text>Attribute CombinationsRoom numbers may not be useful unless you compare it to number of rooms. You could combine the too attributes into a new one that is more useful.<jupyter_code>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"]
#recompute correlation matrix
corr_matrix = housing.corr()
corr_matrix["median_house_value"].sort_values(ascending=False)<jupyter_output><empty_output><jupyter_text>Prepare the Data for ML AlgorithmsYou should not do this manually, write functions instead.<jupyter_code>#Start with fresh data set again for this section
#Separate the labels from the data (the label is the median house value, aka what we are trying to predict)
housing =  strat_train_set.drop("median_house_value", axis=1)
housing_labels = strat_train_set["median_house_value"].copy()<jupyter_output><empty_output><jupyter_text>Data CleaningML hates missing features (aka incomplete data).A few options:* get rid of rows with missing features.* get rid of the attr that has missing data.* set missing values to some value (zero, median, mean, something not harmful).<jupyter_code>#housing.dropna(subset=["total_bedrooms"]) #drop rows with the attr missing.
#housing.drop("total_bedrooms", axis=1) #drop whole attr
median = housing["total_bedrooms"].median() #be sure to save the median for future use.
housing["total_bedrooms"].fillna(median, inplace=True) #option 3.
#using Imputer from scikit
from sklearn.preprocessing import Imputer
imputer = Imputer(strategy="median")
housing_num = housing.drop("ocean_proximity", axis=1) #imputer needs numeric only data to work.
imputer.fit(housing_num)
imputer.statistics_
housing_num.median().values
#Use the imputer to transform training set to replace all missing values
X = imputer.transform(housing_num)
housing_tr = pd.DataFrame(X, columns=housing_num.columns)<jupyter_output><empty_output><jupyter_text>Handling Text and categorical attributes<jupyter_code>housing_cat = housing["ocean_proximity"]
housing_cat.head(10)
#ML prefers numbers over text, so we can convert text into numbers.
housing_cat_encoded, housing_categories = housing_cat.factorize()
housing_cat_encoded[:10]
housing_categories<jupyter_output><empty_output><jupyter_text>Some ML algorithms assume that numbers that are close together have significance, so this could present problems for our numerical encoding here,  we don't want to imply similarity where there is none.A fix is to use "one hot" binary encoding<jupyter_code>from sklearn.preprocessing import OneHotEncoder
encoder = OneHotEncoder()
housing_cat_1hot = encoder.fit_transform(housing_cat_encoded.reshape(-1, 1))
housing_cat_1hot #efficient sparse matrix format.
housing_cat_encoded.reshape(-1, 1) 
#fit_transform expected a 2d array, so we transposed the 1d array using reshape.<jupyter_output><empty_output><jupyter_text>Newer versions of Scikit-learn have a `CategoricalEncoder` class, which could simplify this one hot encoding.If you find the number of categories you have is really large, this could be detrimental to training and performance. Look into "embeddings", chapter 14.  SKLearn custom transformersYou can write your own data transformers. Following example adds in the combined attributes that we computed manually earlier.<jupyter_code>from sklearn.base import BaseEstimator, TransformerMixin

rooms_ix, bedrooms_ix, population_ix, household_ix = 3,4,5,6

class CombinedAttributesAdder(BaseEstimator, TransformerMixin):
    def __init__(self, add_bedrooms_per_room = True): #don't use *args or **kargs
        #This variable is a hyperparameter of the ML algorithm.
        self.add_bedrooms_per_room = add_bedrooms_per_room
        
    def fit(self, X, y=None):
        return self
    
    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]
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()<jupyter_output><empty_output><jupyter_text>Feature Scaling and Transformation Pipelines.ML Algorithms usually don't perform well when the attributes have vastly different scales on the numerical attributes. Use min-max scaling (normalization), or standardization. SK learn has builtins for these.Use Transformation Pipelines to do this processing.<jupyter_code>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)<jupyter_output><empty_output><jupyter_text>It would be great to add a Pandas DataFrame directly into the pipeline. To do that we can create a custom transformer for the start of the pipe.<jupyter_code>from sklearn.base import BaseEstimator, TransformerMixin

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, y=None):
        return X[self.attribute_names].values<jupyter_output><empty_output><jupyter_text>Now we show using multiple pipelines and then (one for numeric and one for category data) and then how to join them together.<jupyter_code>num_attribs = list(housing_num)
cat_attribs = ["ocean_proximity"]
num_attribs #list of attributes (columns) from DataFrame
num_pipeline = Pipeline([
    ('selector', DataFrameSelector(num_attribs)),
    ('imputer', Imputer(strategy="median")),
    ('attribs_adder', CombinedAttributesAdder()),
    ('std_scaler', StandardScaler())
])

from future_encoders import OneHotEncoder
cat_pipeline = Pipeline([
    ('selector', DataFrameSelector(cat_attribs)),
    ('onehot_encoder', OneHotEncoder())
])<jupyter_output><empty_output><jupyter_text>Warning: earlier versions of the book used the LabelBinarizer or CategoricalEncoder classes to convert each categorical value to a one-hot vector. It is now preferable to use the OneHotEncoder class. Right now it can only handle integer categorical inputs, but in Scikit-Learn 0.20 it will also handle string categorical inputs (see PR 10521). So for now we import it from future_encoders.py, but when Scikit-Learn 0.20 is released, you can import it from sklearn.preprocessing instead:<jupyter_code>#Join pipes using FeatureUnion
from sklearn.pipeline import FeatureUnion

full_pipeline = FeatureUnion(transformer_list=[
    ('num_pipeline', num_pipeline),
    ('cat_pipeline', cat_pipeline)                                    
])

housing_prepared = full_pipeline.fit_transform(housing)
housing_prepared.shape<jupyter_output><empty_output><jupyter_text>Select and Train a ModelFinally! We've done the prep work to make the rest of this go as smoothly as possible.<jupyter_code>#Linear Regression model.
from sklearn.linear_model import LinearRegression

lin_reg = LinearRegression()
lin_reg.fit(housing_prepared, housing_labels)
#Try out linear model on a some data from the training set.
some_data = housing.iloc[:5]
some_labels= housing_labels.iloc[:5]

#yay pipeline reuse
some_data_prepared = full_pipeline.transform(some_data)

#predict
print("Predictions:", lin_reg.predict(some_data_prepared))
print("labels:", list(some_labels))
#Measure RMSE to see how the model performed.
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<jupyter_output><empty_output><jupyter_text>This means typical error is $68K. not great. Model might be underfit to training data. Perhaps the features are not good for predictions, or the model is not powerful enough.Lets try a more powerful model than Linear next, Decision Tree<jupyter_code>from sklearn.tree import DecisionTreeRegressor

tree_reg = DecisionTreeRegressor()
tree_reg.fit(housing_prepared, housing_labels)
housing_predictions = tree_reg.predict(housing_prepared)

tree_mse = mean_squared_error(housing_labels, housing_predictions)
tree_rmse = np.sqrt(tree_mse)
tree_rmse<jupyter_output><empty_output><jupyter_text>That can't be right. Perhaps we are badly overfit. But still *DON'T TOUCH TEST SET UNTIL LAUNCH*We can use part of the training data for Cross Validation Cross Validation for Better Evaluation<jupyter_code>#Randomly split the training set into 10 distinct subsets (called folds)
#Trains and evals the models 10 times and cross validates.
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)

#SK learn cross validation expects utility function (great is better) compared to cost function (lower is better)
#thats why we use -scores before calculating sqrt
tree_rmse_scores = np.sqrt(-scores)

def display_scores(scores):
    print('Scores:', scores)
    print('Mean:', scores.mean())
    print('Std Dev:', scores.std())
    
display_scores(tree_rmse_scores)<jupyter_output>Scores: [69316.02634772 65498.84994772 71404.25935862 69098.46240168
 70580.30735263 75540.88413124 69717.93143674 70428.42648461
 75888.17618283 68976.12268448]
Mean: 70644.94463282847
Std Dev: 2938.9378926251234<jupyter_text>These metrics seem more realistic than 0 error after doing cross validation. Seems worse than Linear regression, but lets use cross validation on Linear to have apples to apples comparison.<jupyter_code>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)<jupyter_output>Scores: [66782.7384406  66960.1180821  70347.9524213  74739.57054676
 68031.13387002 71193.84181749 64969.63056717 68281.61139934
 71552.91569222 67665.10078121]
Mean: 69052.46136182113
Std Dev: 2731.6740070650844<jupyter_text>This shows that our decision tree is overfitting and is even worse than a simple linear regression. Lets try one more model, Random Forest.<jupyter_code>from sklearn.ensemble import RandomForestRegressor
forest_reg = RandomForestRegressor(random_state=42)
forest_reg.fit(housing_prepared, housing_labels)
#Full Set
housing_predictions = forest_reg.predict(housing_prepared)
forest_mse = mean_squared_error(housing_labels, housing_predictions)
forest_rmse = np.sqrt(forest_mse)
forest_rmse
#Cross Val for forest
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_scores(forest_rmse_scores)<jupyter_output>Scores: [51650.94405471 48920.80645498 52979.16096752 54412.74042021
 50861.29381163 56488.55699727 51866.90120786 49752.24599537
 55399.50713191 53309.74548294]
Mean: 52564.19025244012
Std Dev: 2301.873803919754<jupyter_text>We see that the cross val scores for Forest show worse performance than the RMSE for the whole training set. This means we may still be overfitting the data using this model.  Fine Tuning the ModelAssume you have a short list of a few promising models. now you want to fine tune them.<jupyter_code>#Using Grid search to tune hyperparameters for Random Forest
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()
grid_search = GridSearchCV(forest_reg, param_grid, cv=5, scoring='neg_mean_squared_error')
grid_search.fit(housing_prepared, housing_labels)<jupyter_output><empty_output><jupyter_text>This grid search will check all combinations of parameters (in this case, 18 combos in all). Since CV is 5, it means there will be 90 rounds of training against the training set (for cross validation). Might take a while.<jupyter_code>#Once grid search is done, see what it found.
grid_search.best_params_
#get best model directly
grid_search.best_estimator_
#get scores
cvres = grid_search.cv_results_
for mean_score, params in sorted(zip(cvres["mean_test_score"], cvres["params"])):
    print(np.sqrt(-mean_score), params)<jupyter_output>64583.23831343775 {'max_features': 2, 'n_estimators': 3}
62470.342102357106 {'bootstrap': False, 'max_features': 2, 'n_estimators': 3}
60717.61577683773 {'max_features': 4, 'n_estimators': 3}
60353.800853111454 {'bootstrap': False, 'max_features': 3, 'n_estimators': 3}
59152.22141463088 {'bootstrap': False, 'max_features': 4, 'n_estimators': 3}
58650.6191468369 {'max_features': 8, 'n_estimators': 3}
58493.21467912166 {'max_features': 6, 'n_estimators': 3}
55576.873722117016 {'max_features': 2, 'n_estimators': 10}
54496.78629093729 {'bootstrap': False, 'max_features': 2, 'n_estimators': 10}
52956.30662941401 {'max_features': 2, 'n_estimators': 30}
52634.105267438536 {'max_features': 4, 'n_estimators': 10}
52397.67552164831 {'bootstrap': False, 'max_features': 3, 'n_estimators': 10}
52369.18784516499 {'max_features': 8, 'n_estimators': 10}
52208.348666826736 {'max_features': 6, 'n_estimators': 10}
51855.47219074427 {'bootstrap': False, 'max_features': 4, 'n_estimators': 10}
50528.7451100[...]<jupyter_text>Randomized SearchGrid search works ok if the state space of hyperparameters is not too large. Otherwise using RandomizedSearchCV might be more tractable  Analyze Best models and their Errors<jupyter_code>#RandomForest can tell you which attributes ended up being most important.
feature_importances = grid_search.best_estimator_.feature_importances_
feature_importances
#Display in context
extra_attribs = ["rooms_per_hhold", "pop_per_hhold", "bedrooms_per_room"]
cat_encoder = cat_pipeline.named_steps["onehot_encoder"]
cat_one_hot_attribs = list(cat_encoder.categories_[0])
attributes = num_attribs + extra_attribs + cat_one_hot_attribs
sorted(zip(feature_importances, attributes), reverse=True)<jupyter_output><empty_output><jupyter_text>Based on these results we could try dropping some of the less important attributes/categories, since they are not contributing much to the model.  Evaluate System on the Test Set<jupyter_code>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) #just transform, not fit!

final_predictions = final_model.predict(X_test_prepared)
final_mse = mean_squared_error(y_test, final_predictions)
final_rmse = np.sqrt(final_mse)
final_rmse<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.09138381201044386, 'fraction_numerical': 0.04154642480875819, 'mean_word_length': 6.638908327501749, 'pattern_counts': {'":': 0, '<': 119, '<?xml version=': 0, '>': 118, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '30', 'dir': 'jupyter-structured-clean-dedup', 'id': '493398', 'license': 'MIT', 'n_tokens_mistral': 7689, 'n_tokens_neox': 6892, 'n_words': 1925, 'path': 'handsOn/ch2/Housing.ipynb', 'repo_name': 'codyhanson/ml-toys', 'size': '21836'}