In typical environments our data would be available in a relational database (or some other common data store) and spread across multiple tables/documents/files. To access it, we would first need to get our credentials and access authorizations and familiarize ourselves with the data schema.

In thos notebook, however, things are much simpler: we will just download a single compressed file, housing.tgz, which contains a comma-separated values (CSV) file called housing.csv with all the data.

It is preferable to create a small function to do that. Having a function that downloads the data is useful in particular if the data changes regularly: we can write a small script that uses the function to fetch the latest data (or we can set up a scheduled job to do that automatically at regular inter‐ vals). Automating the process of fetching the data is also useful if we need to install the dataset on multiple machines.

The next cell is the function to fetch the data:

Now let's load the data using pandas

The info() method is useful to get a quick description of the data, in particular the total number of rows, each attribute’s type, and the number of nonnull values

All attributes are numerical, except the ocean_proximity field. Its type is object, so it could hold any kind of Python object. But since we loaded this data from a CSV file, we know that it must be a text attribute. When we looked at the top five rows, we probably noticed that the values in the ocean_proximity column were repetitive, which means that it is probably a categorical attribute. We can find out what categories exist and how many districts belong to each category by using the value_counts() method

Let’s look at the other fields. The describe() method shows a summary of the numerical attributes

Another quick way to get a feel of the type of data we are dealing with is to plot a histogram for each numerical attribute. A histogram shows the number of instances (on the vertical axis) that have a given value range (on the horizontal axis). We can either plot this one attribute at a time, or we can call the hist() method on the whole dataset (as shown in the following code example), and it will plot a histogram for each numerical attribute

It should be noted that the median income attribute does not look like it is expressed in US dollars (USD). In fact, the data has been scaled and capped at 15 (actually, 15.0001) for higher median incomes, and at 0.5 (actually, 0.4999) for lower median incomes. The numbers represent roughly tens of thousands of dollars (e.g., 3 actually means about $30,000). Working with preprocessed attributes is common in Machine Learning, and it is not necessarily a problem, but we should try to understand how the data was computed.

What's indeed a problem is that housing median age and the median house value were also capped. The latter may be a serious problem since it is our target attribute (our labels). Our Machine Learning algorithms may learn that prices never go beyond that limit. We need to check with our business team to see if this is a problem or not. If they tell us that they need precise predictions even beyond $500,000, then we have two options:

1. Collect proper labels for the districts whose labels were capped.
2. Remove those districts from the training set (and also from the test set, since our system should not be evaluated poorly if it predicts values beyond $500,000).

Finally, many histograms are tail-heavy: they extend much farther to the right of the median than to the left. This may make it a bit harder for some Machine Learning algorithms to detect patterns. We will try transforming these attributes later on to have more bell-shaped distributions.

Now before we look at the data any further, we need to create a test set, put it aside, and never look at it.

Creating a Test Set¶

It may sound strange to voluntarily set aside part of the data at this stage. After all, we have only taken a quick glance at the data, and surely we should learn a whole lot more about it before we decide what algorithms to use, right? This is true, but our brain is an amazing pattern detection system, which means that it is highly prone to overfitting: if we look at the test set, we may stumble upon some seemingly interesting pattern in the test data that leads us to select a particular kind of Machine Learning model. When we estimate the generalization error using the test set, our estimate will be too optimistic, and we will launch a system that will not perform as well as expected. This is called data snooping bias.

Creating a test set is theoretically simple: pick some instances randomly, typically 20% of the dataset (or less if your dataset is very large), and set them aside:

import numpy as np

def split_training_and_test_data(data, test_ratio):
    shuffled_indices = np.random.permutation(len(data))
    test_set_size = int(len(data) * test_ratio)
    test_indices = shuffled_indices[:test_set_size]
    training_indices = shuffled_indices[test_set_size:]
    return data.iloc[training_indices], data.iloc(test_indices)

training_set, test_set = split_training_and_test_data(housing, 0.2)

Well, this works, but it is not perfect: if we run the program again, it will generate a different test set! Over time, we (or our Machine Learning algorithms) will get to see the whole dataset, which is what we want to avoid.

To have a stable train/test split even after updating the dataset, a common solution is to use each instance’s identifier to decide whether or not it should go in the test set (assuming instances have a unique and immutable identifier). For example, we could compute a hash of each instance’s identifier and put that instance in the test set if the hash is lower than or equal to 20% of the maximum hash value. This ensures that the test set will remain consistent across multiple runs, even if we refresh the dataset. The new test set will contain 20% of the new instances, but it will not contain any instance that was previously in the training set.

Here is a possible implementation:

from zlib import crc32

def test_set_check(identifier, test_ratio):
    return crc32(np.int64(identifier)) & 0xffffffff < test_ratio * 2**32

def split_train_test_by_id(data, test_ratio, id_column):
    ids = data[id_column]
    in_test_set = ids.apply(lambda id_: test_set_check(id_, test_ratio)) return data.loc[~in_test_set], data.loc[in_test_set]

Unfortunately, the housing dataset does not have an identifier column. The simplest solution is to use the row index as the ID:

housing_with_id = housing.reset_index() # adds an `index` column
train_set, test_set = split_train_test_by_id(housing_with_id, 0.2, "index")

If we use the row index as a unique identifier, we need to make sure that new data gets appended to the end of the dataset and that no row ever gets deleted. If this is not possible, then we can try to use the most stable features to build a unique identifier. For example, a district’s latitude and longitude are guaranteed to be stable for a few million years, so we could combine them into an ID like so

housing_with_id["id"] = housing["longitude"] * 1000 + housing["latitude"]
train_set, test_set = split_train_test_by_id(housing_with_id, 0.2, "id")

Scikit-Learn provides a few functions to split datasets into multiple subsets in various ways. The simplest function is train_test_split(), which does pretty much the same thing as the function split_training_and_test_data(), with a couple of additional features. First, there is a random_state parameter that allows us to set the random generator seed. Second, we can pass it multiple datasets with an identical number of rows, and it will split them on the same indices (this is very useful, for example, if we have a separate DataFrame for labels):

from sklearn.model_selection import train_test_split
train_set, test_set = train_test_split(housing, test_size=0.2, random_state=42)

So far we have considered purely random sampling methods. This is generally fine if our dataset is large enough (especially relative to the number of attributes), but if it is not, we run the risk of introducing a significant sampling bias. We want a sampling that is representative of the whole dataset. For this, we need to use some domain knowledge.

Suppose we chatted with experts who told us that the median income is a very important attribute to predict median housing prices. We may want to ensure that the test set is representative of the various categories of incomes in the whole dataset. Since the median income is a continuous numerical attribute, we first need to create an income category attribute. Let’s look at the median income histogram above more closely: most median income values are clustered around 1.5 to 6 (i.e., $15,000–$60,000), but some median incomes go far beyond 6. It is important to have a sufficient number of instances in our dataset for each stratum, or else the estimate of a stratum’s importance may be biased. This means that we should not have too many strata, and each stratum should be large enough. The following code uses the pd.cut() function to create an income category attribute with five categories (labeled from 1 to 5): category 1 ranges from 0 to 1.5 (i.e., less than $15,000), category 2 from 1.5 to 3, and so on:

We use this to construct a stratified sampling, i.e. the number of samples of each median income category follows the same distribution shown above. For this, we can use Scikit-Learn’s StratifiedShuffleSplit class:

Let’s see if this worked as expected. We can start by looking at the income category proportions in the test set:

Now we should remove the income_cat attribute so the data is back to its original state:

Visualizing Data¶

First, make sure we have put the test set aside and we are only exploring the training set. Also, if the training set is very large, we may want to sample an exploration set, to make manipulations easy and fast. In our case, the set is quite small, so we can just work directly on the full set. Let’s create a copy so that we can play with it without harming the training set:

Now that’s much better: we can clearly see the high-density areas, namely the Bay Area and around Los Angeles and San Diego, plus a long line of fairly high density in the Central Valley, in particular around Sacramento and Fresno.

Now let’s look at the housing prices shown in the figure above. The radius of each circle represents the district’s population (option s), and the color represents the price (option c). We will use a predefined color map (option cmap) called jet, which ranges from blue (low values) to red (high prices)

Looking for Correlations¶

Since the dataset is not too large, we can easily compute the standard correlation coefficient (also called Pearson’s r) between every pair of attributes using the corr() method:

The correlation coefficient ranges from –1 to 1

• When the coefficient is close to 1, it means that there is a strong positive correlation
• When the coefficient is close to –1, it means that there is a strong negative correlation
• Coefficients close to 0 mean that there is no linear correlation

The correlation coefficient only measures linear correlations (“if x goes up, then y generally goes up/down”). It may completely miss out on nonlinear relationships (e.g., “if x is close to 0, then y gener‐ ally goes up”).

Another way to check for correlation between attributes is to use the pandas scatter_matrix() function, which plots every numerical attribute against every other numerical attribute. Since there are now 11 numerical attributes, we would get 11 x 11 = 121 plots, which would not fit on a page; so let’s just focus on a few promising attributes that seem most correlated with the median housing value:

Note that the main diagonal (top left to bottom right) would be full of straight lines if pandas plotted each variable against itself, which would not be very useful. So instead pandas displays a histogram of each attribute (other options are available; see the pandas documentation for more details).

Experimenting with Attribute Combinations¶

We see that the new bedrooms_per_room attribute is much more correlated with the median house value than the total number of rooms or bedrooms shown previsouly. Trying out various attribute combinations could also give us useful insights that will help us get a good prototype in the next training iteration.

Preparing the Data for Machine Learning Algorithms¶

After we have gained enough empirical ideas on what our data looks like, it's time to feed the data into a data preprocessing pipeline and use the preprosessed data to train a model.

What we mean by "preprocessing" is to make sure our training data doesn't contain any value that could potentially trigger a runtime error during model training. For our housing data, we will execute the following 4 preprocessings in order:

1. Handling the missing feautre values in dataset
2. Handling text and categorical attributes
3. Composing attributes
4. Feature scaling

We will show the code for each of these 4 steps separately below

But first let’s separate the predictors and the labels, since we don’t necessarily want to apply the same transformations to the predictors and the target values :

Handling Missing Values¶

Most Machine Learning algorithms cannot work with missing features. We saw earlier that the total_bedrooms attribute has some missing values. We have 3 options:

1. Get rid of the corresponding districts.
2. Get rid of the whole attribute.
3. Set the values to some value (zero, the mean, the median, etc.).

Scikit-Learn provides a handy class to take care of missing values: SimpleImputer:

Since the median can only be computed on numberical attributes, we need to create a copy of the data without the text attribute ocean_proximity and then fit the data. The result is a plain NumPy array containing the transformed features that are put back into a pandas DataFrame

Handling Text and Categorical Attributes¶

We see the only text attribute, ocean_proximity, is not random text, but a limited number of possible values, each of which represents a category. Moast machine-learning algorithms perfer working with numbers, so let's convert these categories to numbers. For this, we use scikit-learn's OrdinalEncoder:

One issue with this representation is that ML algorithms will assume that two nearby values are more similar than two distant values. This may be fine in some cases (e.g., for ordered categories such as “bad,” “average,” “good,” and “excellent”), but it is obviously not the case for the ocean_proximity column (for example, categories 0 and 4 are clearly more similar than categories 0 and 1). To fix this issue, a common solution is to create one binary attribute per category:

• one attribute equal to 1 when the category is <1H OCEAN and 0, otherwise
• one attribute equal to 1 when the category is 'INLAND' and 0, otherwise
• so on and so forth...

This is called one-hot encoding, because only one attribute will be equal to 1 (hot) while others will all be 0 (cold) for one row from the housing data. The new attribute are sometimes called dummy attribute. scikit-learn offers a OneHotEncoder class to convert categorical values into one-hot vectors:

If a categorical attribute has a large number of possible categories (e.g., country code, profession, species), then one-hot encoding will result in a large number of input features. This may slow down training and degrade performance. If this happens, we may want to replace the categorical input with useful numerical features related to the categories: for example, we could replace the ocean_proximity feature with the distance to the ocean (similarly, a country code could be replaced with the country’s population and GDP per capita). Alternatively, we could replace each category with a learnable, low-dimensional vector called an embedding. Each category’s representation would be learned during training. This is an example of representation learning

Composing Attributes (Custom Transformers)¶

Before moving on, let's pause for a second for one question: "How would we bundle up preprocessing steps together?" This is an important question because at the end of the day the increasing degree of automation is required so that

• we could reproduce these transformations easily on any dataset
• we will be able to gradually build transformation library with increasing complexity
• we are able to perform typerparameter tunning which greatly enhances the quality of Machine Learning outcome (we will briefly touch this in this section)

sciki-learn, for this purpose, introduced transformers. Let's pick up the attribute composition we have done early and re-do it using such transformer

From the housing.head() we used earlier, we need and know the indexes of following 4 attributes:

1. total_rooms = 3
2. total_bedrooms = 4
3. population = 5
4. households = 6

In this example the transformer has one hyperparameter, add_bedrooms_per_room, set to True by default (it is often helpful to provide sensible defaults). This hyperparameter will allow us to easily find out whether adding this attribute helps the Machine Learning algorithms or not. More generally, we can add a hyperparameter to gate any data preparation step that we are not 100% sure about. The more we automate these data preparation steps, the more combinations we can automatically try out, making it much more likely that we will find a great combination (and saving us a lot of time).

Feature Scaling¶

With few exceptions, Machine Learning algorithms don’t perform well when the input numerical attributes have very different scales. There are two common ways to get all attributes to have the same scale:

1. min-max scaling
2. standardization.

Min-max scaling (many people call this normalization) is the simplest: values are shifted and rescaled so that they end up ranging from 0 to 1. We do this by subtracting the min value and dividing by the max minus the min. Scikit-Learn provides a transformer called MinMaxScaler for this. It has a feature_range hyperparameter that lets us change the range if, for some reason, we don’t want 0 ~ 1.

Standardization is different: first it subtracts the mean value (so standardized values always have a zero mean), and then it divides by the standard deviation so that the resulting distribution has unit variance. Unlike min-max scaling, standardization does not bound values to a specific range, which may be a problem for some algorithms (e.g., neural networks often expect an input value ranging from 0 to 1). However, standardization is much less affected by outliers. For example, suppose a district had a median income equal to 100 (by mistake). Min-max scaling would then crush all the other values from 0 ~ 15 down to 0 ~ 0.15, whereas standardization would not be much affected. Scikit-Learn provides a transformer called StandardScaler for standardization.

Puting Things Together - Transformation Pipelines¶

Selecting and Training a Model¶

Thanks to all previous steps, we are ready to train a Linear Regression model

Done! We now have a working Linear Regression model. Let’s try it out on a few instances from the training set:

Let’s measure this regression model’s RMSE on the whole train‐ ing set using Scikit-Learn’s mean_squared_error() function:

Most districts’ median_housing_values range between \$120,000 and \$265,000, so a typical prediction error of \$68,628 is not very satisfying.

This is an example of a model underfitting the training data. When this happens it can mean that the features do not provide enough information to make good predictions, or that the model is not powerful enough.

The main approaches to fix underfitting are to

• select a more powerful model
• feed the training algorithm with better features, or
• reduce the constraints on the model.

This model is not regularized, which rules out the last option. We could try to add more features (e.g., the log of the population), but first let’s try a more complex model to see how it does.

The model has badly overfit the data. Note that we haven't touched test dataset yet. One way to evaluate the Decision Tree model would be to use such validation set.

A great alternative is to use scikit-learn’s K-fold cross-validation feature. The following code randomly splits the training set into 10 distinct subsets called folds, then it trains and evaluates the Decision Tree model 10 times, picking a different fold for evaluation every time and training on the other 9 folds. The result is an array containing the 10 evaluation scores:

scikit-learn’s cross-validation features expect a utility function (greater is better) rather than a cost function (lower is better), so the scoring function is actually the opposite of the MSE (i.e., a neg‐ ative value), which is why the preceding code computes -scores before calculating the square root.

Now the Decision Tree doesn’t look as good as it did earlier. In fact, it seems to perform worse than the Linear Regression model! Notice that cross-validation allows us to get not only an estimate of the performance of our model, but also a measure of how precise this estimate is (i.e., its standard deviation).

Let’s compute the same scores for the Linear Regression model just to be sure:

That’s right: the Decision Tree model is overfitting so badly that it performs worse than the Linear Regression model.

Building a model on top of many other models is called Ensemble Learning, and it is often a great way to push ML algorithms even further. So let's try Random Forests which trains many Decision Trees on random subsets of the features, then averaging out their predictions:

Fine-Tuning Models¶

Grid Search¶

One option would be to fiddle with the hyperparameters manually, until we find a great combination of hyperparameter values. This would be very tedious work, and we may not have time to explore many combinations.

Instead, we should get scikit-learn’s GridSearchCV to search for us. All we need to do is tell it which hyperparameters we want it to experiment with and what values to try out, and it will use cross-validation to evaluate all the possible combinations of hyperparameter values. For example, the following code searches for the best combination of hyperparameter values for the RandomForestRegressor:

This param_grid above tells scikit-learn to first evaluate all 3 × 4 = 12 combinations of n_estimators and max_features hyperparameter values, then try all 2 × 3 = 6 combinations of hyperparameter values in the second dict, but this time with the bootstrap hyperparameter set to False

As a result, the grid search will explore 12 + 6 = 18 combinations of RandomForestRegressor hyperparameter values, train each model cv=5 times. There will be 18 x 5 = 90 rounds of traninig.

When training is done, we get the besst parameter combination and best estimator (best model) by

And of course the evalueation scores are also available:

Randomized Search¶

The grid search approach is fine when we are exploring relatively few combinations, like in the previous setup, but when the hyperparameter search space is large, it is often preferable to use RandomizedSearchCV instead.

This class can be used in much the same way as the GridSearchCV class, but instead of trying out all possible combinations, it evaluates a given number of random combinations by selecting a random value for each hyperparameter at every iteration. This approach has two main benefits

1. If we let the randomized search run for, say, 1,000 iterations, this approach will explore 1,000 different values for each hyperparameter (instead of just a few values per hyperparameter with the grid search approach)
2. Simply by setting the number of iterations, we have more control over the computing budget we want to allocate to hyperparameter search.

Ensemble Methods¶

Another way to fine-tune our system is to try to combine the models that perform best. The group (or "ensemble") will often perform better than the best individual model (just like Random Forests perform better than the individual Decision Trees they rely on), especially if the individual models make very different types of errors.

Analyze the Best Models and Their Errors¶

You will often gain good insights on the problem by inspecting the best models. For example, the RandomForestRegressor can indicate the relative importance of each attribute for making accurate predictions:

With this information, you may want to try dropping some of the less useful features

Evaluate Model on the Test Set¶

After tweaking our models for a while, we eventually have a system that performs sufficiently well. Now is the time to evaluate the final model on the test set.