<jupyter_start><jupyter_text>Most of the code in this notebook is from https://github.com/ageron/handson-ml. This is a great book - please buy the book to support the author.This notebook can serve as a good training tutorial for beginner data scientists.The 10-Step Machine Learning Project Workflow (My Version):1. Define business object2. Make sense of the data from a high level    - data types (number, text, object, etc.)    - continuous/discrete    - basic stats (min, max, std, median, etc.) using boxplot    - frequency via histogram    - scales and distributions of different features3. Create the traning and test sets using proper sampling methods, e.g., random vs. stratified4. Correlation analysis (pair-wise and attribute combinations)5. Data cleaning (missing data, outliers, data errors)6. Data transformation via pipelines (categorical text to number using one hot encoding, feature scaling via normalization/standardization, feature combinations)7. Train and cross validate different models and select the most promising one (Linear Regression, Decision Tree, and Random Forest were tried in this tutorial)8. Fine tune the model using trying different combinations of hyperparameters9. Evaluate the model with best estimators in the test set10. Launch, monitor, and refresh the model and system<jupyter_code># This Python 3 environment comes with many helpful analytics libraries installed
# It is defined by the kaggle/python docker image: https://github.com/kaggle/docker-python

import numpy as np
import pandas as pd

%matplotlib inline
import matplotlib.pyplot as plt
import seaborn as sns

# Input data files are available in the "../input/" directory.
# For example, running this (by clicking run or pressing Shift+Enter) will list the files in the input directory

import os
print(os.listdir("input"))

# Any results you write to the current directory are saved as output.
# loading data
data_path = "input/housing.csv"
housing = pd.read_csv(data_path)

# see the basic info
housing.info()<jupyter_output><class 'pandas.core.frame.DataFrame'>
RangeIndex: 20640 entries, 0 to 20639
Data columns (total 10 columns):
longitude             20640 non-null float64
latitude              20640 non-null float64
housing_median_age    20640 non-null float64
total_rooms           20640 non-null float64
total_bedrooms        20433 non-null float64
population            20640 non-null float64
households            20640 non-null float64
median_income         20640 non-null float64
median_house_value    20640 non-null float64
ocean_proximity       20640 non-null object
dtypes: float64(9), object(1)
memory usage: 1.6+ MB<jupyter_text>What information should we get from `info()`- total observations: 20640: make sure you understand what is each observation. In this case, each observation is the data about a district- total columns (features): 10- data type of each feature: 9 numbers and 1 object (will handle later)- meaning of each feature: it is very important to work with domain expert to fully understand each feature- any null values (e.g., total_bedrooms is 20433, which indicates null values - will handle later)<jupyter_code>housing.head(10)<jupyter_output><empty_output><jupyter_text>Data Types`head()` shows the valuse of top rows, which gives more idea on data types. Pandas guessed the data types of features when reading in the data, which may not always work. In this dataset, you can see ocean_proximity feature is text. Sometimes, the price feature may also be object type becaue the raw data has $ sign, in which case you need to convert the data type from object to float64 if you want to use this feature in the model.  Basic Stats`describe()` shows a summary of numerial features, which can be visualized using boxplots and histograms. `value_counts()` can be used to generate a summary of categorical features.<jupyter_code>housing.describe()
# boxplot 
housing.boxplot(['median_house_value'], figsize=(10, 10))
# histogram
housing.hist(bins=50, figsize=(15, 15))<jupyter_output><empty_output><jupyter_text>Given that `.botplot()` and `.hist()` only handle numerical features. We cannot forget ocean_proximity, which is object type (no need to change to string).<jupyter_code>housing['ocean_proximity'].value_counts()
op_count = housing['ocean_proximity'].value_counts()
plt.figure(figsize=(10,5))
sns.barplot(op_count.index, op_count.values, alpha=0.7)
plt.title('Ocean Proximity Summary')
plt.ylabel('Number of Occurrences', fontsize=12)
plt.xlabel('Ocean Proximity', fontsize=12)
plt.show()
# housing['ocean_proximity'].value_counts().hist()<jupyter_output><empty_output><jupyter_text>Make Sence of the DataWhat are the typical things we can learn from the basic statistics with visualizaitons? 1. Do the data make sence? scan each column and see whether the data make sense at a high level.     - longitude, latitude and house median age seem to be OK.     - total rooms and total bedrooms are in hundreds and thousands - this does make sense given each row is a district    - population seems to be OK but you want to know what's the unit for the number, e.g., thousands? millions? households are numbers of people living together, which is OK. Households mean is 499 and population mean is 1425, so you can tell that population is just the total number of people in each district not in thousands or millions.     - median income is apparently problematic - how can mean income be 3.87? Actually, this is because median income has been scaled and capped between 15.0001 (max) and  0.4999 (min). Preprocessed attributes like this are common. Knowning how the data is calculated is critical.    - median house value data is OK and this is our target variable, i.e., we want to build a model to predict this value. 2. Feature scaling: you have noticed that the features have very different scales, which we need to handle later3. Distribution: from the histograms, we can tell many of them are skewed, i.e., having a long tail on one side or the other. In many cases, we need to transform the data so that they have more bell-shaped distributions.  Create Training, Validation, and Test Data SetsTo avoid [data snooping bias](https://en.wikipedia.org/wiki/Data_dredging), you should set aside 20% of the data as the test set without further exploriation of the data.<jupyter_code>from sklearn.model_selection import train_test_split
train_set, test_set = train_test_split(housing, test_size=0.2, random_state=42)
print(len(train_set), "train + ", len(test_set), "test")
# check whether the test set is the same for every run
test_set.head(10)<jupyter_output><empty_output><jupyter_text>The above function will generate a different test set every time you run the program, which is undesirable. A few potential solutions:1. save the test set and training set for the first run into csv files2. set a fixed random number generator seed so that the random numbers are always the same (we use this approach in this notebook using `random_state=42`)1 and 2 won't work if new data is loaded. The book proposed a method of hashing the identifier of the data to get the same test set (see page 50 for details).   About SamplingIn order to make sure that the test set is a good representation of the overall population. We may want to consider different sampling techniques:- random sampling (what we used above): OK if the dataset is large enough (how large is enough?)- stratified sampling: the population is devided into homogeneous subgroups called *strata* and the correct instances are sampled from each *stratum*  (such as for a survey of 1000 people, given US population has 51.3% female and 48.7% male, 513 female and 487 male should be surveyed instead of pure random sampling) Suppose you learned that median income is very important for predicting median housing prices. You may want to use stratified sampling for the test set. To do that, you first need to change median income from a continuse attribute to a categorical attribute. As shown in the histogram below, we can see most of the income are around 2 and 5, and some are beyond.<jupyter_code>housing['median_income'].hist()<jupyter_output><empty_output><jupyter_text>We limit the number of categories by dividing the median income by 1.5 and merge all the income greater than 5 into 5. Then, we can use stratified sampling.<jupyter_code>housing['income_cat'] = np.ceil(housing['median_income']/1.5)
# DataFrame.where(cond, other=nan, inplace=False, axis=None, level=None, errors='raise', try_cast=False, raise_on_error=None)
# Where cond is True, keep the original value. Where False, replace with corresponding value from other
housing['income_cat'].where(housing['income_cat']<5, 5.0, inplace=True)
housing['income_cat'].hist()
# stratified sampling based on income categories
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]

strat_test_set.head(10)<jupyter_output><empty_output><jupyter_text>We can compare the different sampling results for the test set by comparing them to the overall population distribution as follows. As you can see, stratified sampling's distributions are much more similar to the overall distributions compared with random sampling<jupyter_code>housing['income_cat'].value_counts() / len(housing)
strat_test_set['income_cat'].value_counts() / len(strat_test_set)
# we need to do the random sampling again to include income_cat column
from sklearn.model_selection import train_test_split
train_set, test_set = train_test_split(housing, test_size=0.2, random_state=42)

test_set['income_cat'].value_counts() / len(test_set)
# drop the income_cat attributes
for set_ in (strat_train_set, strat_test_set):
    set_.drop("income_cat", axis=1, inplace=True)
# check the dropping result
strat_test_set.info()<jupyter_output><class 'pandas.core.frame.DataFrame'>
Int64Index: 4128 entries, 5241 to 2398
Data columns (total 10 columns):
longitude             4128 non-null float64
latitude              4128 non-null float64
housing_median_age    4128 non-null float64
total_rooms           4128 non-null float64
total_bedrooms        4079 non-null float64
population            4128 non-null float64
households            4128 non-null float64
median_income         4128 non-null float64
median_house_value    4128 non-null float64
ocean_proximity       4128 non-null object
dtypes: float64(9), object(1)
memory usage: 354.8+ KB<jupyter_text>Now you should set aside the test set and only use the training set. When the training set is too large, you can also create an exploration set out of the training set to make the initial analysis fast. we create a copy of the training set as follows:<jupyter_code>housing = strat_train_set.copy()
housing.info()<jupyter_output><class 'pandas.core.frame.DataFrame'>
Int64Index: 16512 entries, 17606 to 15775
Data columns (total 10 columns):
longitude             16512 non-null float64
latitude              16512 non-null float64
housing_median_age    16512 non-null float64
total_rooms           16512 non-null float64
total_bedrooms        16354 non-null float64
population            16512 non-null float64
households            16512 non-null float64
median_income         16512 non-null float64
median_house_value    16512 non-null float64
ocean_proximity       16512 non-null object
dtypes: float64(9), object(1)
memory usage: 1.4+ MB<jupyter_text>Additional Visualizations for Data ExplorationThe following geographical data visualizations show that the price is related to the location and population density.<jupyter_code>housing.plot(kind='scatter', x='longitude', y='latitude', alpha=0.1)
# option s: radius of each circle represent the population/100
# option c: color represents the median price
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>Correlation AnalysisWe want to further explore the data to look for correlations between different attributes. correlation coefficient is between -1 and 1, representing negative and positive correlations. 0 means there is no liner correlation. Correlation is said to be linear if the **ratio of change** is constant, otherwise is non-linear. Visual inspection of data is very important, which cannot be replaced by simply looking at the correlation coefficients, i.e., check out Anscombe's quartet: https://en.wikipedia.org/wiki/Anscombe%27s_quartet<jupyter_code># Anscombe's quartet: https://seaborn.pydata.org/examples/anscombes_quartet.html
sns.set(style="ticks")
anscombe = pd.read_csv("input/anscombe.csv")

# Show the results of a linear regression within each dataset
sns.lmplot(x="x", y="y", col="dataset", hue="dataset", data=anscombe,
           col_wrap=2, ci=None, palette="muted", size=4,
           scatter_kws={"s": 50, "alpha": 1})
# Pearson's r, aka, standard correlation coefficient for every pair
corr_matrix = housing.corr()
# Check the how much each attribute correlates with the median house value
corr_matrix['median_house_value'].sort_values(ascending=False)
from pandas.tools.plotting import scatter_matrix
attributes = ['median_house_value', 'median_income', 'total_rooms', 'housing_median_age']
scatter_matrix(housing[attributes], figsize=(12,12))<jupyter_output>/Users/harrywang/Documents/sandbox/data-analysis-tutorials/house-price-prediction/venv/lib/python3.6/site-packages/ipykernel_launcher.py:3: FutureWarning: 'pandas.tools.plotting.scatter_matrix' is deprecated, import 'pandas.plotting.scatter_matrix' instead.
  This is separate from the ipykernel package so we can avoid doing imports until<jupyter_text>We can see that median_income is a promising attribute to predict median_house_value. A close-up of the scatterplot is as follows:<jupyter_code>housing.plot(kind='scatter', x='median_income', y='median_house_value', alpha=0.2, figsize=(10,10))<jupyter_output><empty_output><jupyter_text>We can see there are a number of "horizontal lines" in the plot: one clear one at $500,000, one at $450,000, another one at $350,000, and a few other ones. Try to find out why that is happening. If you cannot figure out the reason, removing those data points (if not too many) might be a good idea before feeding the data to the algorithms.   Attribute CombinationsSometime, the combinations of attributes are more meaningful and interesting in terms of solving the business problems, e.g.,- rooms per household: total  of rooms per district is not useful but rooms per household may be interesting- bedroom/total room ratio- population per household<jupyter_code># calculated attributes
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']

# checkout the correlations again
corr_matrix = housing.corr()
corr_matrix['median_house_value'].sort_values(ascending=False)<jupyter_output><empty_output><jupyter_text>Two findings after combining attributes:1. rooms_per_household is slightly more correlated (0.146285) with house value than total_rooms (0.135097)2. bedrooms_per_room is much more correlated (-0.259984) than total_rooms (0.135097) and total_bedrooms (0.047689): houses with lower bedroom/room ratio is more expensive: this sort of make sense, more expensive houses may have more offices, dens, playrooms, etc.   Data Cleaning and TransformationTypically, data need to be cleaned and transformed before trying different ML algorithms. missing data in one attribute. Three ways to handle this:1.   remove the observations with missing values using `dropna()`; ```housing.dropna(subset=['total_bedrooms']```2.  remove the entire attribute using `drop()`;```housing.drop('total_bedrooms', axis=1)```3. set/impute the missing values using `fillna()````median = housing['total_bedroom'].median()housing['total_bedrooms'].fillna(median, inplace=True)```  Seperate the predictors (independent variables) and labels (target/dependent variables)We want to create a clean training set first.<jupyter_code>housing.info()
housing = strat_train_set.drop("median_house_value", axis=1) # drop target labels for training set
housing_labels = strat_train_set["median_house_value"].copy() # this is the target label vector
housing.info()
# using Scikit-Learn Imputer
from sklearn.preprocessing import Imputer
imputer = Imputer(strategy='median')

# remove non-numerical attributes for Imputer by making a copy of the dataframe
housing_num = housing.drop('ocean_proximity', axis=1)

imputer.fit(housing_num)  # this computes median for each attributes and store the result in statistics_ variable
imputer.statistics_  # same result as housing_num.median().values
# see attributes with missing values
housing_num.info()
x = imputer.transform(housing_num)  # this is a Numpy array
housing_tr = pd.DataFrame(x, columns=housing_num.columns)  # change a Numpy array to a DataFrame
housing_tr.info()  # no missing values<jupyter_output><class 'pandas.core.frame.DataFrame'>
RangeIndex: 16512 entries, 0 to 16511
Data columns (total 8 columns):
longitude             16512 non-null float64
latitude              16512 non-null float64
housing_median_age    16512 non-null float64
total_rooms           16512 non-null float64
total_bedrooms        16512 non-null float64
population            16512 non-null float64
households            16512 non-null float64
median_income         16512 non-null float64
dtypes: float64(8)
memory usage: 1.0 MB<jupyter_text>Text and Categorial AttributesMost ML algorithms work with numbers better. Therefore, we often need to convert text attributes into numerical attributes. For ocean_proximity, we have two ways to handle this problem:1. map each category to a number, such as "<1H OCEAN" is 0, "INLAND" is 1, 'NEAR OCEAN' is 4, etc. The problem with this solution is that ML algorithm may think 4 is greater than 0, which could cause problem.2. To address the problem in 1, we can also create a binary variable for each attribute, which is called one-hot encoding (only one is 1 hot, all others are 0 cold)<jupyter_code># Approach 1
from sklearn.preprocessing import LabelEncoder
encoder = LabelEncoder()
housing_cat = housing['ocean_proximity']
housing_cat.head()
housing_cat_encoded = encoder.fit_transform(housing_cat)
housing_cat_encoded
print(encoder.classes_)  # '<1H OCEAN' is 0, 'INLAND' is 1, etc.
# Approach 2
# reshape
from sklearn.preprocessing import OneHotEncoder
encoder = OneHotEncoder()  # don't forget the ()!!!
housing_cat_1hot = encoder.fit_transform(housing_cat_encoded.reshape(-1,1))  # this returns a sparse SciPy matrix
housing_cat_1hot.toarray()  # convert the sparse matrix to numpy array
# Combine Approch 1 and 2 in one shot
from sklearn.preprocessing import LabelBinarizer
encoder = LabelBinarizer()
housing_cat_1hot = encoder.fit_transform(housing_cat)
housing_cat_1hot<jupyter_output><empty_output><jupyter_text>Custom TransformersYou may need to develop custom transformers - you can just write a simple function for that or if you want your translormer work with Scikit-Learn, you need to develope the transformers as a class.<jupyter_code># A custom transformer
from sklearn.base import BaseEstimator, TransformerMixin
rooms_ix, bedrooms_ix, population_ix, household_ix = 3, 4, 5, 6  # hardcoded just for this dataset

class CombineAttributesAdder(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
    
    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]

            
housing.head()  # note that rooms_per_household, and population_per_household already calculated before
attr_adder = CombineAttributesAdder(add_bedrooms_per_room=False)  # add_bedrooms_per_room is called a hyperparameter
housing_extra_attribs = attr_adder.transform(housing.values)  # housing.values is the numpy N-array representation of the DataFrame
housing_extra_attribs
# check the stats of the training set for feature scaling
housing_tr.describe()<jupyter_output><empty_output><jupyter_text>Feature ScalingTypically, ML algorithms don't perform well when the input numerial attributes have very different scales. For example, in this housing dataset (shown above), you can see median_income ranges from 0.4999 to 15 while total rooms is between 6 and 39320. Note that scaling the target values is typically not required. Two major scaling methods (two different scalers):- normalization (aka, min-max scaling): values are rescaled to between 0 and 1 using (value-min)/(max-min)- standardization: values are rescaled to have unit variance: (value - mean)/varianceNormalization can be dramatically affected by outliner data while standardization handles outliers very well.<jupyter_code># Transformation Pipeline
# name/estimator pairs for pipeline steps
# each estimator except the last one must be transformers with fit_transform() method
# pipeline fit() calls each fit_transform() of each estimator and fit() for the last estimator
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler

num_pipeline = Pipeline([
    ('imputer',Imputer(strategy='median')),
    ('attribs_adder', CombineAttributesAdder()),
    ('std_scaler', StandardScaler()),
])

housing_num_tr = num_pipeline.fit_transform(housing_num)<jupyter_output><empty_output><jupyter_text>Scikit-Learn only handles Numpy arrays not Pandas Dataframes, we can create another transformer so that we can feed the pipeline a DataFrame<jupyter_code># this is the fix to the problem at https://stackoverflow.com/questions/46162855/fit-transform-takes-2-positional-arguments-but-3-were-given-with-labelbinarize
# CategoricalEncoder should be used instead of LabelEncoder in the latest version of Scikit-Learn
# Definition of the CategoricalEncoder class, copied from PR #9151.
# Just run this cell, or copy it to your code, do not try to understand it (yet).

from sklearn.base import BaseEstimator, TransformerMixin
from sklearn.utils import check_array
from sklearn.preprocessing import LabelEncoder
from scipy import sparse

class CategoricalEncoder(BaseEstimator, TransformerMixin):
    def __init__(self, encoding='onehot', categories='auto', dtype=np.float64,
                 handle_unknown='error'):
        self.encoding = encoding
        self.categories = categories
        self.dtype = dtype
        self.handle_unknown = handle_unknown

    def fit(self, X, y=None):
        """Fit the CategoricalEncoder to X.
        Parameters
        ----------
        X : array-like, shape [n_samples, n_feature]
            The data to determine the categories of each feature.
        Returns
        -------
        self
        """

        if self.encoding not in ['onehot', 'onehot-dense', 'ordinal']:
            template = ("encoding should be either 'onehot', 'onehot-dense' "
                        "or 'ordinal', got %s")
            raise ValueError(template % self.handle_unknown)

        if self.handle_unknown not in ['error', 'ignore']:
            template = ("handle_unknown should be either 'error' or "
                        "'ignore', got %s")
            raise ValueError(template % self.handle_unknown)

        if self.encoding == 'ordinal' and self.handle_unknown == 'ignore':
            raise ValueError("handle_unknown='ignore' is not supported for"
                             " encoding='ordinal'")

        X = check_array(X, dtype=np.object, accept_sparse='csc', copy=True)
        n_samples, n_features = X.shape

        self._label_encoders_ = [LabelEncoder() for _ in range(n_features)]

        for i in range(n_features):
            le = self._label_encoders_[i]
            Xi = X[:, i]
            if self.categories == 'auto':
                le.fit(Xi)
            else:
                valid_mask = np.in1d(Xi, self.categories[i])
                if not np.all(valid_mask):
                    if self.handle_unknown == 'error':
                        diff = np.unique(Xi[~valid_mask])
                        msg = ("Found unknown categories {0} in column {1}"
                               " during fit".format(diff, i))
                        raise ValueError(msg)
                le.classes_ = np.array(np.sort(self.categories[i]))

        self.categories_ = [le.classes_ for le in self._label_encoders_]

        return self

    def transform(self, X):
        """Transform X using one-hot encoding.
        Parameters
        ----------
        X : array-like, shape [n_samples, n_features]
            The data to encode.
        Returns
        -------
        X_out : sparse matrix or a 2-d array
            Transformed input.
        """
        X = check_array(X, accept_sparse='csc', dtype=np.object, copy=True)
        n_samples, n_features = X.shape
        X_int = np.zeros_like(X, dtype=np.int)
        X_mask = np.ones_like(X, dtype=np.bool)

        for i in range(n_features):
            valid_mask = np.in1d(X[:, i], self.categories_[i])

            if not np.all(valid_mask):
                if self.handle_unknown == 'error':
                    diff = np.unique(X[~valid_mask, i])
                    msg = ("Found unknown categories {0} in column {1}"
                           " during transform".format(diff, i))
                    raise ValueError(msg)
                else:
                    # Set the problematic rows to an acceptable value and
                    # continue `The rows are marked `X_mask` and will be
                    # removed later.
                    X_mask[:, i] = valid_mask
                    X[:, i][~valid_mask] = self.categories_[i][0]
            X_int[:, i] = self._label_encoders_[i].transform(X[:, i])

        if self.encoding == 'ordinal':
            return X_int.astype(self.dtype, copy=False)

        mask = X_mask.ravel()
        n_values = [cats.shape[0] for cats in self.categories_]
        n_values = np.array([0] + n_values)
        indices = np.cumsum(n_values)

        column_indices = (X_int + indices[:-1]).ravel()[mask]
        row_indices = np.repeat(np.arange(n_samples, dtype=np.int32),
                                n_features)[mask]
        data = np.ones(n_samples * n_features)[mask]

        out = sparse.csc_matrix((data, (row_indices, column_indices)),
                                shape=(n_samples, indices[-1]),
                                dtype=self.dtype).tocsr()
        if self.encoding == 'onehot-dense':
            return out.toarray()
        else:
            return out
# given a list of attributes names, this transformer converts the dataframe to a numpy array
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):
        return X[self.attribute_names].values
# create two pipelines and use feture union to join them
num_attribs = list(housing_num)
cat_attribs = ['ocean_proximity']

num_pipeline = Pipeline([
    ('selector', DataFrameSelector(num_attribs)),
    ('imputer', Imputer(strategy='median')),
    ('attribs_adder', CombineAttributesAdder()),
    ('std_scaler', StandardScaler()),
])

cat_pipeline = Pipeline([
    ('selector', DataFrameSelector(cat_attribs)),
    ('label_binarizer', CategoricalEncoder()),
    # ('label_binarizer', LabelBinarizer()),  # LabelBinarizer does not work this way with last Scikit-Learn
])

housing_num_tr = num_pipeline.fit_transform(housing)
housing_num_tr.shape
num_attribs
housing_cat_tr = cat_pipeline.fit_transform(housing)
housing_cat_tr
from sklearn.pipeline import FeatureUnion

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

# run the whole pipeline
housing_prepared = full_pipeline.fit_transform(housing)
housing_prepared.shape<jupyter_output><empty_output><jupyter_text>Select and Train a ModelWe are going to try Linear Regression, Decision Tree, and Random Forest models.<jupyter_code># Linear Regression
from sklearn.linear_model import LinearRegression

lin_reg = LinearRegression()
lin_reg.fit(housing_prepared, housing_labels)  # housing_prepared are independent variables and housing_labels are dependent variables
# test out the linear regression model
some_data = housing.iloc[:5]  # choose the first five observations
some_labels = housing_labels.iloc[:5]
some_data
some_data_prepared = full_pipeline.transform(some_data)
some_data_prepared
print('Actual Prices:', list(some_labels))  # actual prices
# print predicted prices
print('Predicted Prices:', lin_reg.predict(some_data_prepared))<jupyter_output>Predicted Prices: [210644.60474898 317768.8075158  210956.43310747  59218.98838141
 189747.558595  ]<jupyter_text>As you can see, the sample predictions are not very good, i.e., the first one is off by (286600-210644)/286600 ~= 26.5%! Let's calculate the RMSE (root-mean-square error) on the whole training set.<jupyter_code>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
housing_labels.describe()<jupyter_output><empty_output><jupyter_text>25% and 75% quantile are `$120,000`  and  `$264,000` respectively, which means 50% of the house prices are between those two values. Therefore, $6,8628 error is not very good, which is a typical example of **underfitting**.Three ways you can potentially improve the results:1. try a different model2. try better features, i.e., feature engineering3. reduce the constraints on the model is any (for example, if the model is regularized)<jupyter_code># Try Decision Tree
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>No error at all! This is a typical example of **overfitting**. You cannot use the same set of data for both training and validation. Cross-Validation (CV) can help with model validation.Scikit-Learn CV features expect a utility function (greater is better) than a cost function (lower is better), which is the reason for having `-scores`:<jupyter_code># 10-fold cross validation
from sklearn.model_selection import cross_val_score

# for decision tree
tree_scores = cross_val_score(tree_reg, housing_prepared, housing_labels, scoring='neg_mean_squared_error', cv=10)
tree_rmse_scores = np.sqrt(-tree_scores)

# for linear regression
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)
print('Scores:', tree_rmse_scores)
print('Mean:', tree_rmse_scores.mean())
print('Standard Deviation:', tree_rmse_scores.std())
print('Scores:', lin_rmse_scores)
print('Mean:', lin_rmse_scores.mean())
print('Standard Deviation:', lin_rmse_scores.std())<jupyter_output>Standard Deviation: 2731.674004792122<jupyter_text>Now, Decision Tree Model's performance is actually worse than the Linear Regression Model: mean rmse 71493 vs. 69051 (the numbers differ everytime you run the models)<jupyter_code># Try Random Forest, which is an Ensemble Learning model
from sklearn.ensemble import RandomForestRegressor
forest_reg = RandomForestRegressor()
forest_reg.fit(housing_prepared, housing_labels)

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)
print('Mean:', forest_rmse_scores.mean())<jupyter_output>Mean: 52879.811612746016<jupyter_text>Random Forest is much better that the previous two models. You can try other models, such as SVM, NN, etc. We assume you choose Random Forest as the model and discuss how you can fine tune the model for better performance.   Fine Tune the ModelThere are different ways you can fine tune your model:1. try different combinations of hyperparameters of a model:     a. the following example trys 18 different combinations of hyperparameters and find the best one    b. you can also use Randomized Search to try more combinations when the search space is very large2. combine the best performing models<jupyter_code># use GridSearch to find best hyperparameter combinations
from sklearn.model_selection import GridSearchCV

param_grid = [
    {'n_estimators':[3, 10, 30], 'max_features': [2, 4, 6, 8]},  # try 3x4=12 combinations
    {'bootstrap': [False], 'n_estimators':[3, 10], 'max_features': [2, 3, 4]},  # try 2x3=6 combinations
]

forest_reg = RandomForestRegressor()
grid_search = GridSearchCV(forest_reg, param_grid, cv=5, scoring='neg_mean_squared_error')  # each model is trained 5 times, so (12+6)*5 = 80 rounds of training in total
grid_search.fit(housing_prepared, housing_labels)
grid_search.best_params_  # best parameters
grid_search.best_estimator_  # best estimators
# The importance of the features
feature_importances = grid_search.best_estimator_.feature_importances_
feature_importances
extra_attribs = ['rooms_per_hhold', 'pop_per_hhold', 'bedrooms_per_room']
cat_one_hot_attribs = list(encoder.classes_)
attributes = num_attribs = num_attribs + extra_attribs + cat_one_hot_attribs
sorted(zip(feature_importances, attributes), reverse=True)<jupyter_output><empty_output><jupyter_text>Based on the feature importance, you can choose to drop some features such as the last four ocean proximity features to simplify the model. In the following example, the performanc on the test set is actually better than the validation set.  Evaluation via the Test SetThis step is to see how the model performs on unknow data. As long as the result is not way off from the validation result, you should go ahead lauch the model.<jupyter_code>final_model = grid_search.best_estimator_  # best model

# see the best rmse on the validation set
best_valiation_score = grid_search.best_score_
best_validation_rmse = np.sqrt(-best_valiation_score)
best_validation_rmse
# see the final rmse on the test set
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)  # note DO NOT USE fit_transform!! not need to fit anymore
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.08090906530738685, 'fraction_numerical': 0.02354332704384804, 'mean_word_length': 4.496904024767802, 'pattern_counts': {'":': 2, '<': 100, '<?xml version=': 0, '>': 97, 'https://': 6, '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': '25', 'dir': 'jupyter-structured-clean-dedup', 'id': '648613', 'license': 'MIT', 'n_tokens_mistral': 10870, 'n_tokens_neox': 10015, 'n_words': 3659, 'path': 'house-price-prediction/house_price_prediction.ipynb', 'repo_name': 'BUAD625-Team4-2021/tutorial-buffet', 'size': '35509'}