<jupyter_start><jupyter_text>Predicting the Sale Price of Bulldozers using Machine Learning> In this notebook, we're going to go through an example machine learning project with the goal of predicting the sale price of bulldozers.  1. Problem definition> How well can we predict the future sale price for a bulldozer, given its characteristics and previous examples of how much similar bulldozers have been sold for? 2. DataThe data is downloaded from the Kaggle Bluebook for Bulldozers compitition: https://www.kaggle.com/c/bluebook-for-bulldozers/dataThere are 3 main datasets:* Train.csv is the training set, which contains data through the end of 2011.* Valid.csv is the validation set, which contains data from January 1, 2012 - April 30, 2012 You make predictions on this set throughout the majority of the competition. Your score on this set is used to create the public leaderboard.* Test.csv is the test set, which won't be released until the last week of the competition. It contains data from May 1, 2012 - November 2012. Your score on the test set determines your final rank for the competition. 3. EvaluationThe evaluation metric for this competition is the RMSLE (root mean squared log error) between the actual and predicted auction prices.For more on the evaluation of this project check: https://www.kaggle.com/c/bluebook-for-bulldozers/overview/evaluation**Note:** The goal for most regression evaluation metrics is to minimize the error. For example, our goal for this project will be to build a machine learning model which minimises RMSLE. 4. FeaturesKaggle provides a data dictionary detailing all of the features of the dataset. You can view this data dictionary on Google Sheets: https://docs.google.com/spreadsheets/d/18ly-bLR8sbDJLITkWG7ozKm8l3RyieQ2Fpgix-beSYI/edit?usp=sharing<jupyter_code>import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import sklearn
# Import training and validation sets
df = pd.read_csv("data/TrainAndValid.csv",
                 low_memory=False)
df.info()
df.isna().sum()
df.columns<jupyter_output><empty_output><jupyter_text>Do you know how to plot a figure using `matlotlib.plt`?Let's plot a figure with `saledate` as its x-label and `SalePrice` as its y-label.Since there are too many data, we only plot the first 1000.<jupyter_code>fig, ax = plt.subplots()
ax.scatter(df["saledate"][:1000], df["SalePrice"][:1000])<jupyter_output><empty_output><jupyter_text>Why is the x-label messed up? Let's check it.<jupyter_code>df.saledate[:1000]<jupyter_output><empty_output><jupyter_text>It looks weird! Let's check the data type of `saledate`...<jupyter_code>df.saledate.dtype<jupyter_output><empty_output><jupyter_text>Oh! It's `Object`, maybe that's why x-label messed up...But we are not sure, let's go further!<jupyter_code>df.SalePrice.plot.hist()<jupyter_output><empty_output><jupyter_text>The lower the price, the more its sales.  Parsing datesWhen we work with **time series data**, we want to **enrich the time & date component** as much as possible.We can do that by telling pandas which of our columns has dates in it using the `parse_dates` parameter.<jupyter_code># Import data again but this time parse dates
df = pd.read_csv("data/TrainAndValid.csv",
                 low_memory=False,
                 parse_dates=["saledate"])
df.saledate.dtype<jupyter_output><empty_output><jupyter_text>Before: Object | After: <M8[ns]But what is `<M8[ns]`?- Well, `<M8[ns]` & `datetime64[ns]` are basically the same thing, at least on my laptop (whose byte order is `little endian`): https://stackoverflow.com/questions/29206612/difference-between-data-type-datetime64ns-and-m8ns- `datetime64[ns]` is data type in **Numpy** that represents a date or a time. It was named this way because `datetime` has been taken by **Python**.<jupyter_code>np.dtype('datetime64[ns]') == np.dtype('<M8[ns]')<jupyter_output><empty_output><jupyter_text>So, this means after parsing dates, the data type of the `saledate` column is a kind of datetime.Let's check the first 1000 rows of `saledate` column:<jupyter_code>df.saledate[:1000]<jupyter_output><empty_output><jupyter_text>Wow, this data type looks clearner and more readable than what we have seen before.The purpose of this date parsing is to regulate the x-label in the figure.<jupyter_code>fig, ax = plt.subplots()
ax.scatter(df["saledate"][:1000], df["SalePrice"][:1000])<jupyter_output><empty_output><jupyter_text>Looks matplotlib can understand this `<M8[ns]` data type. Great!Let's do more EDA(exploratory data analysis):<jupyter_code>df.head()<jupyter_output><empty_output><jupyter_text>Everytime we deal with a dataframe with long width, we can do **Transpose(转置矩阵)**. https://en.wikipedia.org/wiki/Transpose<jupyter_code>df.head().T
df.saledate.head(20)<jupyter_output><empty_output><jupyter_text>Sort DataFrame by saledateWhen working with time series data, it's a good idea to sort it by date.<jupyter_code># Sort DataFrame in date order
df.sort_values(by=["saledate"], inplace=True, ascending=True)
df.saledate.head(20)<jupyter_output><empty_output><jupyter_text>Make a copy of the original DataFrameWe make a copy of the original dataframe so when we manipulate the copy, we've still got our original data.<jupyter_code># Make a copy of the original DataFrame to perform edits on
df_tmp = df.copy()<jupyter_output><empty_output><jupyter_text>What we are about to do is called **Feature Engineering**.  Add datetime parameters for `saledate` column<jupyter_code>df_tmp["saleYear"] = df_tmp.saledate.dt.year
df_tmp["saleMonth"] = df_tmp.saledate.dt.month
df_tmp["saleDay"] = df_tmp.saledate.dt.day
df_tmp["saleDayOfWeek"] = df_tmp.saledate.dt.dayofweek
df_tmp["saleDayOfYear"] = df_tmp.saledate.dt.dayofyear
df_tmp.head().T
# Now we've enriched our DataFrame with data time features, we can remove the "saledate" column.
df_tmp.drop("saledate", axis=1, inplace=True)
df_tmp.saledate
# Check the values of different columns
df_tmp.state.value_counts()<jupyter_output><empty_output><jupyter_text>5. ModellingWe've done enough EDA (we could always do more), but let's start to do some model-driven EDA.<jupyter_code># Let's build a machine learning model. (Regression)
from sklearn.ensemble import RandomForestRegressor

model = RandomForestRegressor(n_jobs=-1,
                              random_state=42)

model.fit(df_tmp.drop("SalePrice", axis=1), df_tmp["SalePrice"])
df_tmp.info()
df_tmp.isna().sum()<jupyter_output><empty_output><jupyter_text>Convert strings to categoriesOne way we can turn all of our data into numbers is by converting them into pandas categories.We can check the different datatypes compatible with pandas here:https://pandas.pydata.org/pandas-docs/stable/reference/general_utility_functions.htmldata-types-related-functionality<jupyter_code>df_tmp.head().T
pd.api.types.is_string_dtype(df_tmp["UsageBand"])
df_tmp.items()
# Find the columns which contain strings
for label, content in df_tmp.items():
    if pd.api.types.is_string_dtype(content):
        print(label)
# If you're wondering what df.items() does, here's an example
random_dict = {"key1": "Hello,",
               "key2": "world!"}

for key, value in random_dict.items():
    print(f"this is a key: {key}",
          f"this is a value: {value}")
# This will turn all of the string values into category values
for label, content in df_tmp.items():
    if pd.api.types.is_string_dtype(content):
        df_tmp[label] = content.astype("category").cat.as_ordered()
df_tmp.info()
df_tmp.state.cat.categories
df_tmp.state.cat
df_tmp.state.cat.codes<jupyter_output><empty_output><jupyter_text>Thanks to pandas Categories, we now have a way to access all of our data in the form of numbers.But we still have a bunch of missing data...<jupyter_code># Check missing data
df_tmp.isnull().sum()/len(df_tmp)<jupyter_output><empty_output><jupyter_text>Save preprocessed data<jupyter_code># Export current tmp dataframe
df_tmp.to_csv("data/train_tmp.csv",
              index=False)
# Import preprocessed data
df_tmp = pd.read_csv("data/train_tmp.csv",
                     low_memory=False)
df_tmp.head().T
df_tmp.isnull().sum()<jupyter_output><empty_output><jupyter_text>Fill missing values Fill numerical missing values first<jupyter_code>for label, content in df_tmp.items():
    if pd.api.types.is_numeric_dtype(content):
        print(label)
# Check for which numeric colums have null values
for label, content in df_tmp.items():
    if pd.api.types.is_numeric_dtype(content):
        if pd.isnull(content).sum():
            print(label)<jupyter_output>auctioneerID
MachineHoursCurrentMeter<jupyter_text>Looks like these two numerical columns have null values. Let's fill them with the median.<jupyter_code># Fill numeric rows with the median
for label, content in df_tmp.items():
    if pd.api.types.is_numeric_dtype(content):
        if pd.isnull(content).sum():
            # Add a binary column which tells us if the data was missing
            df_tmp[label+"_is_missing"] = pd.isnull(content)
            # Fill missing numeric values with median
            df_tmp[label] = content.fillna(content.median())<jupyter_output><empty_output><jupyter_text>After this, the new column will be `True` if it was null, and `False` if it was filled.<jupyter_code># Why median is more robust than mean
hundreds = np.full((1000,), 100)
hundreds_billion = np.append(hundreds, 1000000000)
np.mean(hundreds), np.mean(hundreds_billion), np.median(hundreds), np.median(hundreds_billion)
pd.DataFrame(hundreds).plot.hist();
pd.DataFrame(hundreds_billion).plot.hist();
# Check if there's any null numeric values:
for label, content in df_tmp.items():
    if pd.api.types.is_numeric_dtype(content):
        if pd.isnull(content).sum():
            print(label)<jupyter_output><empty_output><jupyter_text>None? Perfect! That's what we've expected.<jupyter_code># Check to see how many examples were missing
df_tmp.auctioneerID_is_missing.value_counts()<jupyter_output><empty_output><jupyter_text>20136 True's means there were 20136 missing values.Remember before the `auctioneerID` and `MachineHoursCurrentMeter` had missing values?Now they should all be filled.<jupyter_code>df_tmp.isna().sum()<jupyter_output><empty_output><jupyter_text>Look! The `auctioneerID` and `MachineHoursCurrentMeter` now don't have missing values!<jupyter_code># Check for columns which aren't numeric
for label, content in df_tmp.items():
    if not pd.api.types.is_numeric_dtype(content):
        print(label)
# Turn categorical variables into numbers and fill missing values
for label, content in df_tmp.items():
    if not pd.api.types.is_numeric_dtype(content):
        # Add binary column to indicate whether sample had missing value
        df_tmp[label+"_is_missing"] = pd.isnull(content)
        # Turn categories into numbers and add +1
        df_tmp[label] = pd.Categorical(content).codes + 1
df_tmp.head().T
df_tmp.isna().sum()<jupyter_output><empty_output><jupyter_text>Now we have no missing values and all values are numerical...Based on what you know about machine learning, what can we do now?We should be able to build a machine learning model.<jupyter_code>df_tmp.head()<jupyter_output><empty_output><jupyter_text>`%%` is a jupyter magic function.`%%time` calculates the time the current cell takes.<jupyter_code>len(df_tmp)
%%time
# Instantiate model
from sklearn.ensemble import RandomForestRegressor
model = RandomForestRegressor(n_jobs=-1,
                              random_state=42)

# Fit the model
model.fit(df_tmp.drop("SalePrice", axis=1), df_tmp["SalePrice"])
# Score the model
model.score(df_tmp.drop("SalePrice", axis=1), df_tmp["SalePrice"])<jupyter_output><empty_output><jupyter_text>**Question:** Why doesn't the above metric hold water? (why isn't the metric reliable)  Splitting data into train/validation sets<jupyter_code>df_tmp.head()
df_tmp.saleYear
df_tmp.saleYear.value_counts()
# Split data into training and validation
df_val = df_tmp[df_tmp.saleYear == 2012]
df_train = df_tmp[df_tmp.saleYear != 2012]

len(df_val), len(df_train)
# Split data into X & y
X_train, y_train = df_train.drop("SalePrice", axis=1), df_train.SalePrice
X_valid, y_valid = df_val.drop("SalePrice", axis=1), df_val.SalePrice

X_train.shape, y_train.shape, X_valid.shape, y_valid.shape<jupyter_output><empty_output><jupyter_text>Building an evaluation function<jupyter_code># Create evaluation function (the competition uses RMSLE)
from sklearn.metrics import mean_squared_log_error, mean_absolute_error, r2_score

def rmsle(y_test, y_preds):
    """
    Caculates root mean squared log error between predictions and true labels.
    """
    return np.sqrt(mean_squared_log_error(y_test, y_preds))

# Create function to evaluate model on a few different levels
def show_scores(model):
    train_preds = model.predict(X_train)
    val_preds = model.predict(X_valid)
    scores = {"Train MAE": mean_absolute_error(y_train, train_preds),
              "Valid MAE": mean_absolute_error(y_valid, val_preds),
              "Training RMSLE": rmsle(y_train, train_preds),
              "Valid RMSLE": rmsle(y_valid, val_preds),
              "Training R^2": r2_score(y_train, train_preds), 
              "Valid R^2": r2_score(y_valid, val_preds)}
    return scores<jupyter_output><empty_output><jupyter_text>Testing our model on a subset (to tune the hyperparameters)<jupyter_code># # This takes far too long ... for experimenting
# %%time
# model = RandomForestRegressor(n_jobs=-1,
#                               random_state=42)

# model.fit(X_train, y_train)
# Change max_samples value
model = RandomForestRegressor(n_jobs=-1,
                              random_state=42,
                              max_samples=10000)

model
%%time
# Cutting down on the max number of samples each estimator can see improves training time
model.fit(X_train, y_train)
show_scores(model)<jupyter_output><empty_output><jupyter_text>Hyperparamter tuning with RandomizedSearchCV<jupyter_code>%%time
from sklearn.model_selection import RandomizedSearchCV

# Different RandomForestRegressor hyperparameters
rf_grid = {"n_estimators": np.arange(10, 100, 10),
           "max_depth": [None, 3, 5, 10],
           "min_samples_split": np.arange(2, 20, 2),
           "min_samples_leaf": np.arange(1, 20, 2),
           "max_features": [0.5, 1, "sqrt", "auto"],
           "max_samples": [10000]}

# Instantiate RadomizedSearchCV model
rs_model = RandomizedSearchCV(RandomForestRegressor(n_jobs=-1,
                                                    random_state=42),
                             param_distributions=rf_grid,
                             n_iter=2,
                             cv=5,
                             verbose=True)

# Fit the RandomizedSearchCV model
rs_model.fit(X_train, y_train)
# Find the best model hyperparameters
rs_model.best_params_
# Evaluate the RandomizedSearch model
show_scores(rs_model)<jupyter_output><empty_output><jupyter_text>Train a model with the best hyperparamters**Note:** These were found after 100 iterations of `RandomizedSearchCV`.<jupyter_code>%%time

# Most ideal hyperparameters
ideal_model = RandomForestRegressor(n_estimators=40,
                                    min_samples_leaf=1,
                                    min_samples_split=14,
                                    max_features=0.5,
                                    n_jobs=-1,
                                    max_samples=None,
                                    random_state=42) # random state

# Fit the ideal model
ideal_model.fit(X_train, y_train)
# Scores for ideal_model (trained on all the data)
show_scores(ideal_model)
# Scores on rs_model (only trained on -10,000 examples)
show_scores(rs_model)
# Import the test data
df_test = pd.read_csv("data/Test.csv",
                      low_memory=False,
                      parse_dates=["saledate"])

df_test.head()
# Make predictions on the test dataset
test_preds = ideal_model.predict(df_test)<jupyter_output><empty_output><jupyter_text>Preprocessing the data (getting the test data in the same format as our training dataset)<jupyter_code>def preprocess_data(df):
    """
    Performs transformations on df and returns transformed df.
    """
    df["saleYear"] = df.saledate.dt.year
    df["saleMonth"] = df.saledate.dt.month
    df["saleDay"] = df.saledate.dt.day
    df["saleDayOfWeek"] = df.saledate.dt.dayofweek
    df["saleDayOfYear"] = df.saledate.dt.dayofyear
    
    df.drop("saledate", axis=1, inplace=True)
    
    # Fill the numeric rows with median
    for label, content in df.items():
        if pd.api.types.is_numeric_dtype(content):
            if pd.isnull(content).sum():
                # Add a binary column which tells us if the data was missing
                df[label+"_is_missing"] = pd.isnull(content)
                # Fill missing numeric values with median
                df[label] = content.fillna(content.median())
                
            # Fill categorical missing data and turn categories into numbers    
        if not pd.api.types.is_numeric_dtype(content):
                df[label+"_is_missing"] = pd.isnull(content)
                # We add +1 to the category code because pandas encodes missing values as -1
                df[label] = pd.Categorical(content).codes+1
    
    return df
# Process the test data
df_test = preprocess_data(df_test)
df_test.head()
# Make predictions on updated test data
test_preds = ideal_model.predict(df_test)
X_train.head()
# We can find how the columns differ using sets
set(X_train.columns) - set(df_test.columns)
# Manually adjust df_test to have auctioneerID_is_missing column
df_test["auctionID_is_missing"] = False
df_test.head()<jupyter_output><empty_output><jupyter_text>Finally now our test dataframe has the same features as our training dataframe, we can make predictions!<jupyter_code># Make predictions on the test data
test_preds = ideal_model.predict(df_test)
len(test_preds)
# Format predictions into the sae format Kaggle is after
df_preds = pd.DataFrame()
df_preds["SalesID"] = df_test["SalesID"]
df_preds["SalesPrice"] = test_preds
df_preds
# Export prediction data
df_preds.to_csv("data/test_predictions.csv", index=False)<jupyter_output><empty_output><jupyter_text>Feature ImportanceFeature importance seeks to figure out which different attributes of the data were most importance when it comes to predicting the **target variable** (SalePrice).<jupyter_code># Find feature importance of our best model
ideal_model.feature_importances_
# Helper function for plotting feature importance
def plot_features(columns, importances, n=20):
    df = (pd.DataFrame({"features": columns,
                        "feature_importances": importances})
          .sort_values("feature_importances", ascending=False)
          .reset_index(drop=True))
    
    # Plot the dataframe
    fig, ax = plt.subplots()
    ax.barh(df["features"][:n], df["feature_importances"][:20])
    ax.set_ylabel("Features")
    ax.set_xlabel("Fature importance")
    ax.invert_yaxis()
plot_features(X_train.columns, ideal_model.feature_importances_)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.09880678249947666, 'fraction_numerical': 0.011513502198032238, 'mean_word_length': 4.316917084028937, 'pattern_counts': {'":': 16, '<': 141, '<?xml version=': 0, '>': 138, 'https://': 6, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 1, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '34', 'dir': 'jupyter-structured-clean-dedup', 'id': '571153', 'license': 'Apache-2.0', 'n_tokens_mistral': 6272, 'n_tokens_neox': 5944, 'n_words': 1696, 'path': '_notebooks/2020-05-10-end-to-end-bulldozer-price-regression.ipynb', 'repo_name': 'CongLiu-CN/AIisLove.github.io', 'size': '19108'}