<jupyter_start><jupyter_text>Shell APD Learning eXchange DnA Notebook Predicting power usage in the eastern United StatesIn this notebook, we will go through the cycle of a typical data science project, starting with data collection, cleaning, feature engineering and finally predictive modeling. We will also explore the three different categories of machine learning  problems, and how they work. By the end of this exercise notebook you should have a basic understanding of how data science works.  PJM Hourly Energy Consumption DataPJM is a regional transmission organization in the United States. It is part of the Eastern Interconnection grid operating an electric transmission system serving all or parts of Delaware, Illinois, Indiana, Kentucky, Maryland, Michigan, New Jersey, North Carolina, Ohio, Pennsylvania, Tennessee, Virginia, West Virginia, and the District of Columbia.We will be using hourly power consumption data from PJM's East grid. This data set ranges from 2002-2018. The data came from PJM's website and are given in megawatts (MW).  Overview:0. [Introduction](intro)1. [Import packages](setup)2. [Get the data](data)3. [Explore the data](explore)4. [Split train and test sets](traintest)5. [Regression - Predict energy usage](model)    - [Linear regression](linear)    - [Gradient-boosted trees](xgbreg)6. [Engineer new features](featureengineer)7. [Clean the dataset](cleaning)8. [Retrain your model with new features](modelnew)REGRESSION | CLASSIFICATION | CLUSTERING- | - | - |  |   1. Import some python libraries and define functions we'll need<jupyter_code>%%capture

import numpy as np
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
from matplotlib import cm
import matplotlib.dates as mdates
import xgboost as xgb
from xgboost import plot_importance, plot_tree
from sklearn.metrics import mean_squared_error, mean_absolute_error, r2_score
from sklearn.linear_model import LinearRegression, SGDClassifier
from sklearn.ensemble import RandomForestClassifier
from sklearn.metrics import confusion_matrix
from sklearn.metrics import precision_score, recall_score
from sklearn.metrics import f1_score
from sklearn.cluster import KMeans
import seaborn as sns
import math
from datetime import datetime

def mean_absolute_percentage_error(y_true, y_pred): 
    """Calculates MAPE given y_true and y_pred"""
    y_true, y_pred = np.array(y_true), np.array(y_pred)
    return np.mean(np.abs((y_true - y_pred) / y_true)) * 100

def plot_confusion_matrix(y_true, y_pred, classes,
                          normalize=False,
                          title=None,
                          cmap=plt.cm.Blues):
    """
    This function prints and plots the confusion matrix.
    Normalization can be applied by setting `normalize=True`.
    """
    if not title:
        if normalize:
            title = 'Normalized confusion matrix'
        else:
            title = 'Confusion matrix, without normalization'

    # Compute confusion matrix
    cm = confusion_matrix(y_true, y_pred)
    # Only use the labels that appear in the data
#     classes = classes[unique_labels(y_true, y_pred)]
    if normalize:
        cm = cm.astype('float') / cm.sum(axis=1)[:, np.newaxis]
        print("Normalized confusion matrix")
    else:
        print('Confusion matrix, without normalization')

    print(cm)

    fig, ax = plt.subplots(figsize=(7,7)) #
    fig.patch.set_facecolor('white')
    im = ax.imshow(cm, interpolation='nearest', cmap=cmap)
    ax.figure.colorbar(im, ax=ax)
    plt.figure(figsize=(10, 10)) 
    
    # We want to show all ticks...
    ax.set(xticks=np.arange(cm.shape[1]),
           yticks=np.arange(cm.shape[0]),
           # ... and label them with the respective list entries
           xticklabels=classes, yticklabels=classes,
           title=title,
           ylabel='True label',
           xlabel='Predicted label')

    # Rotate the tick labels and set their alignment.
    plt.setp(ax.get_xticklabels(), rotation=45, ha="right",
             rotation_mode="anchor")

    # Loop over data dimensions and create text annotations.
    fmt = '.2f' if normalize else 'd'
    thresh = cm.max() / 2.
    for i in range(cm.shape[0]):
        for j in range(cm.shape[1]):
            ax.text(j, i, format(cm[i, j], fmt),
                    ha="center", va="center",
                    color="white" if cm[i, j] > thresh else "black")
    fig.tight_layout()
    return ax

color_pal = ["#F8766D", "#D39200", "#93AA00", "#00BA38", "#00C19F", "#00B9E3", "#619CFF", "#DB72FB"]

np.set_printoptions(precision=2)
pd.set_option('display.max_columns', None)
plt.style.use('fivethirtyeight')

plt.rcParams.update({'font.size': 22})<jupyter_output><empty_output><jupyter_text>2. Download the data to our notebook<jupyter_code>url = 'https://github.com/robmoratore/ShellDnA/raw/master/data/PJM_Load.csv'
df_full = pd.read_csv(url, sep = ',', date_parser=["date"])
df_full.index = pd.DatetimeIndex(df_full["date"])
df_full = df_full.drop(columns=['date'])
pjme = df_full[['load', 'dayofyear']]<jupyter_output><empty_output><jupyter_text>3. Explore the data   Let's look at the first few rows of the dataset<jupyter_code>pjme.head()<jupyter_output><empty_output><jupyter_text>Plot energy usage over time<jupyter_code>_ = pjme["load"].plot(style='.', figsize=(15,5), color=color_pal[0], title='Energy consumption')
_.set(xlabel="Time", ylabel="Energy (MW)")<jupyter_output><empty_output><jupyter_text>4. Split the train and testing sets We do that to ensure our model is able to generalize. That means, perform well on unseen data. If we train and test using the same data, the model will memorize that specific output and not learn it.We will use the data from 2015 on as our test set.![alt](https://github.com/robmoratore/ShellDnA/blob/master/data/images/train_test.png?raw=true)<jupyter_code>split_date = '2015-01-01'
pjme_train = pjme.loc[pjme.index <= split_date].copy()
pjme_test = pjme.loc[pjme.index > split_date].copy()
ax = pd.DataFrame(pjme_test["load"]) \
    .rename(columns={'load': 'TEST SET'}) \
    .join(pd.DataFrame(pjme_train["load"]).rename(columns={'load': 'TRAINING SET'}), how='outer') \
    .plot(figsize=(15,5), title='Energy consumption', style='.')
_ = ax.set(xlabel="Time", ylabel="Energy (MW)")<jupyter_output><empty_output><jupyter_text>5. Regression - Can we predict energy use for a given hour?   5.1 Start with a simple model - Linear Regression<jupyter_code>lin_reg = LinearRegression()
lin_reg.fit(pjme_train["dayofyear"].values.reshape(-1, 1), pjme_train["load"].values.reshape(-1, 1))
pjme_test['load_Prediction_lin'] = lin_reg.predict(pjme_test["dayofyear"].values.reshape(-1, 1))
_ = pjme_test[['load','load_Prediction_lin']].plot(figsize=(15, 5))
_.set(xlabel="Time", ylabel="Energy (MW)")<jupyter_output><empty_output><jupyter_text>Check how well our model performsWhich error metrics would you choose?<jupyter_code>mse = mean_squared_error(y_true=pjme_test['load'], y_pred=pjme_test['load_Prediction_lin'])
print ("The Mean Squared Error (MSE) is ", round(mse, 2) , "MW\u00b2")
print ("The Root Mean Squared Error (RMSE) is ", round(np.sqrt(mse), 2), "MW")
print ("The Mean Absolute Error (MAE) is ", round(mean_absolute_error(y_true=pjme_test['load'], y_pred=pjme_test['load_Prediction_lin']), 2), "MW")
print ("The Mean Absolute Percentage Error is ", round(mean_absolute_percentage_error(y_true=pjme_test['load'], y_pred=pjme_test['load_Prediction_lin']), 2), "%")
print ("The R\u00b2 (coefficient of determination) regression score is ", round(r2_score(y_true=pjme_test['load'],y_pred=pjme_test['load_Prediction_lin']), 2))<jupyter_output>The Mean Squared Error (MSE) is  42420922.87 MW²
The Root Mean Squared Error (RMSE) is  6513.13 MW
The Mean Absolute Error (MAE) is  5242.46 MW
The Mean Absolute Percentage Error is  17.6 %
The R² (coefficient of determination) regression score is  -0.02<jupyter_text>5.2 Get fancy with some Machine Learning - Gradient-boosted trees This is one of the cleverest ways of combining multiple decision trees. Here, instead of having one single tree responsible for all the decisions, you use many trees on sequence. That means, the next tree does not try to predict the final result but the error of the previous tree. When you add all trees together, the results are generally very good<jupyter_code>include_variables = "dayofyear"
reg = xgb.XGBRegressor(n_estimators=50)
reg.fit(pjme_train[include_variables].values.reshape(-1, 1), pjme_train["load"].values.reshape(-1, 1))
pjme_test['load_Prediction_xgb'] = reg.predict(pjme_test[include_variables].values.reshape(-1, 1))
_ = pjme_test[['load','load_Prediction_xgb']].plot(figsize=(15, 5))
_.set(xlabel="Time", ylabel="Energy (MW)")<jupyter_output><empty_output><jupyter_text>Check how well our model performs<jupyter_code>mse = mean_squared_error(y_true=pjme_test['load'], y_pred=pjme_test['load_Prediction_xgb'])
print ("The Mean Squared Error (MSE) is ", round(mse, 2), "MW\u00b2")
print ("The Root Mean Squared Error (RMSE) is ", round(np.sqrt(mse), 2), "MW")
print ("The Mean Absolute Error (MAE) is ", round(mean_absolute_error(y_true=pjme_test['load'], y_pred=pjme_test['load_Prediction_xgb']), 2), "MW")
print ("The Mean Absolute Percentage Error is ", round(mean_absolute_percentage_error(y_true=pjme_test['load'], y_pred=pjme_test['load_Prediction_xgb']), 2), "%")
print ("The R\u00b2 (coefficient of determination) regression score is ", round(r2_score(y_true=pjme_test['load'],y_pred=pjme_test['load_Prediction_xgb']), 2))<jupyter_output>The Mean Squared Error (MSE) is  31031470.06 MW²
The Root Mean Squared Error (RMSE) is  5570.59 MW
The Mean Absolute Error (MAE) is  4354.54 MW
The Mean Absolute Percentage Error is  14.6 %
The R² (coefficient of determination) regression score is  0.25<jupyter_text>6. Engineer New Features   Which features would you create to help predict load?<jupyter_code>df_full.head().round(2)
_ = df_full.hist(bins=61, figsize=(20,20))<jupyter_output><empty_output><jupyter_text>7. Clean the dataset   Outliers<jupyter_code>_ = df_full["temp"].plot(style='k.', figsize=(15,5))
_.set(xlabel="Time", ylabel="Temperature (C)")
df_full = df_full[(df_full['temp']>=-30) & (df_full['temp']<=55)]
_ = df_full["temp"].plot(style='k.', figsize=(15,5))
_.set(xlabel="Time", ylabel="Temperature (C)")<jupyter_output><empty_output><jupyter_text>Missing values<jupyter_code>len(df_full.index)-df_full.count()
df_full = df_full.dropna()
len(df_full.index)-df_full.count()<jupyter_output><empty_output><jupyter_text>8. Retrain Gradient-boosted trees model with new features We can now retrain the Gradient-boosted trees model using the new features to see if performance improves.Here is a list of available features that can be used in the model.<jupyter_code>df_full.drop(columns=['load']).columns<jupyter_output><empty_output><jupyter_text>Select which features you would like to include in the model, and add them to the `include_variables` list bellow.Then check how the model performs.<jupyter_code># SELECT A COMBINATION OF FEATURES HERE
include_variables = ['hour', 'dayofyear', 'quarter']
#

split_date = '2015-01-01'
df_full_train = df_full.loc[df_full.index <= split_date].copy()
df_full_test = df_full.loc[df_full.index > split_date].copy()

df_full_test.index = pd.DatetimeIndex(df_full_test.index)
df_full_train.index = pd.DatetimeIndex(df_full_train.index)

reg = xgb.XGBRegressor(n_estimators=50)
reg.fit(df_full_train[include_variables], df_full_train["load"].values.reshape(-1, 1))
df_full_test['load_Prediction_xgb'] = reg.predict(df_full_test[include_variables])
_ = df_full_test[['load','load_Prediction_xgb']].plot(figsize=(15, 5))
_.set(xlabel="Time", ylabel="Energy (MW)")<jupyter_output><empty_output><jupyter_text>Check how well our model performsWhich Feature gives you the lowest MAE and highest R$^{2}$?<jupyter_code>mse = mean_squared_error(y_true=df_full_test['load'], y_pred=df_full_test['load_Prediction_xgb'])
print ("The Mean Squared Error (MSE) is ", round(mse, 2), "MW\u00b2")
print ("The Root Mean Squared Error (RMSE) is ", round(np.sqrt(mse), 2), "MW")
print ("The Mean Absolute Error (MAE) is ", round(mean_absolute_error(y_true=df_full_test['load'], y_pred=df_full_test['load_Prediction_xgb']), 2), "MW")
print ("The Mean Absolute Percentage Error is ", round(mean_absolute_percentage_error(y_true=df_full_test['load'], y_pred=df_full_test['load_Prediction_xgb']), 2), "%")
print ("The R\u00b2 (coefficient of determination) regression score is ", round(r2_score(y_true=df_full_test['load'],y_pred=df_full_test['load_Prediction_xgb']), 2))<jupyter_output>The Mean Squared Error (MSE) is  16502377.28 MW²
The Root Mean Squared Error (RMSE) is  4062.31 MW
The Mean Absolute Error (MAE) is  3209.28 MW
The Mean Absolute Percentage Error is  10.45 %
The R² (coefficient of determination) regression score is  0.6<jupyter_text>Feature ImportancesFeature importance is a great way to get a general idea about which features the model is relying on most to make the prediction. This is a metric that simply sums up how many times each feature is split on.We can see that the day of year was most commonly used to split trees, while hour came in next. Quarter has low importance (never used in any splits even) due to the fact that it could be created by different dayofyear splits.<jupyter_code>_ = plot_importance(reg, height=0.9)<jupyter_output><empty_output><jupyter_text>8.1 Look at Worst and Best Predicted Days<jupyter_code>df_full_test['error'] = df_full_test['load'] - df_full_test['load_Prediction_xgb']
df_full_test['abs_error'] = df_full_test['error'].apply(np.abs)
error_by_day = df_full_test.groupby(['year','month','dayofmonth']) \
    .mean()[['load','load_Prediction_xgb','error','abs_error']]<jupyter_output><empty_output><jupyter_text>Worst days<jupyter_code>error_by_day.sort_values('abs_error', ascending=False).round(2).head(10)
f, ax = plt.subplots(1)
f.set_figheight(5)
f.set_figwidth(10)

_ = df_full_test[['load','load_Prediction_xgb']][(df_full_test[['load','load_Prediction_xgb']].index>'2015-02-20') & (df_full_test[['load','load_Prediction_xgb']].index<'2015-02-21')].plot(ax=ax, style=['-','.'])
plot = plt.suptitle('Worst Predicted Day')
_.set(xlabel="Time", ylabel="Energy (MW)")<jupyter_output><empty_output><jupyter_text>Temperature in worst predicted day<jupyter_code>df_full[df_full.index=='2015-02-20 00:00:00']['temp'].head()<jupyter_output><empty_output><jupyter_text>Best days<jupyter_code>error_by_day.sort_values('abs_error', ascending=True).round(2).head(10)
f, ax = plt.subplots(1)
f.set_figheight(5)
f.set_figwidth(10)
_ = df_full_test[['load','load_Prediction_xgb']][(df_full_test[['load','load_Prediction_xgb']].index>'2018-05-17') & (df_full_test[['load','load_Prediction_xgb']].index<'2018-05-18')].plot(ax=ax, style=['-','.'])
plot = plt.suptitle('Best Predicted Day')
_.set(xlabel="Time", ylabel="Energy (MW)")<jupyter_output><empty_output>
StarcoderdataJupyterStructuredCleanDedup