<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Prediction model
Step1:  Import data
Might still need to clean up the files some after import
Step2:  Make fuel price a ratio of the coal price to the natural gas price
Step3:  One-hot encoding of the cluster variable
I'm trying to make this easy for using with different numbers of clusters
Step4:  Add free capacity of every group for that hour.
Turns out that this doesn't help prediction. It actually makes prediction much worse...
Step5:  Drop unnecessary columns and replace nan's with 0
Step6:  Split into training, validation, testing
Step7:  Need scaled versions of the X data for some of the models
Step8:  Check size of all arrays
Step9:  Linear Regression (OLS)
Step10:  XGBoost
Step11:  Validation curve for n_estimators
Step12:  Validation curve for n_estimators
Step13:  Validation curve for reg_alpha
Step14:  Learning curve for n_estimators=250 and max_depth=3
Step15:  Try XGBoost on non-scaled data
Turns out this works better
Step16:  Try ratio of fuel prices
Step17:  Linear Regression (OLS) with ratio of fuel prices
Slight improvement from first OLS try, but still nowhere near as good as gradient boosting
<ASSISTANT_TASK:>
Python Code:
%matplotlib inline
import matplotlib.pyplot as plt
import seaborn as sns
import pandas as pd
import numpy as np
import os
from xgboost import XGBRegressor
from sklearn.linear_model import LinearRegression
from sklearn.svm import SVR, LinearSVR
from sklearn.model_selection import GridSearchCV
from sklearn.preprocessing import StandardScaler
from sklearn.ensemble import GradientBoostingRegressor
from sklearn.model_selection import validation_curve, learning_curve
Explanation: Prediction model
End of explanation
path = '../Final report'
X_fn = 'x.csv'
y_fn = 'y.csv'
X_path = os.path.join(path, X_fn)
y_path = os.path.join(path, y_fn)
X = pd.read_csv(X_path)
y = pd.read_csv(y_path)
X.head()
Explanation: Import data
Might still need to clean up the files some after import
End of explanation
for fuel in ['All coal', 'Lignite', 'Subbituminous']:
    X.loc[:,fuel] = X.loc[:,fuel].values/X.loc[:,'NG Price ($/mcf)'].values
    
X.drop('NG Price ($/mcf)', axis=1, inplace=True)
Explanation: Make fuel price a ratio of the coal price to the natural gas price
End of explanation
cluster_ids = X['cluster'].unique()
for cluster in cluster_ids:
    X['cluster_{}'.format(cluster)] = np.eye(len(cluster_ids))[X['cluster'],cluster]
X.head()
X.tail()
Explanation: One-hot encoding of the cluster variable
I'm trying to make this easy for using with different numbers of clusters
End of explanation
free_cap_dict = {}
for cluster in range(6):
    free_cap_dict[cluster] = X.loc[X['cluster'] == cluster, ['DATETIME', 'nameplate_capacity', 'GROSS LOAD (MW)']]
    col_name = 'cluster_' +  str(cluster) + ' free capacity'
    free_cap_dict[cluster].loc[:,col_name] = (free_cap_dict[cluster].loc[:,'nameplate_capacity'].values - 
                                              free_cap_dict[cluster].loc[:,'GROSS LOAD (MW)'].values)
free_cap_dict[0].head()
for cluster in range(6):
    col_name = 'cluster_' +  str(cluster) + ' free capacity'
    X = pd.merge(X, free_cap_dict[cluster].loc[:,['DATETIME', col_name]], on='DATETIME')
X.head(n=10)
for idx in X.index:
datetime = X.loc[idx, 'DATETIME']
    for cluster in range(6):
        col_name = 'cluster_' +  cluster + ' free capacity'
        X.loc[idx, col_name] = 
y.tail()
Explanation: Add free capacity of every group for that hour.
Turns out that this doesn't help prediction. It actually makes prediction much worse...
End of explanation
X_cols = ['nameplate_capacity', 'GROSS LOAD (MW)', 'ERCOT Load, MW',
          'Total Wind Installed, MW', 'Total Wind Output, MW', 'Net Load Change (MW)',
          'All coal', 'Lignite', 'Subbituminous']
X_cluster_cols = ['cluster_{}'.format(cluster) for cluster in cluster_ids]
# X_cluster_free_cols  = ['cluster_{} free capacity'.format(cluster) for cluster in cluster_ids]
X_clean = X.loc[:,X_cols+X_cluster_cols]#+X_cluster_free_cols]
X_clean.fillna(0, inplace=True)
y_clean = y.loc[:,'Gen Change (MW)']
y_clean.fillna(0, inplace=True)
print X_clean.shape
print y_clean.shape
X_clean.head()
Explanation: Drop unnecessary columns and replace nan's with 0
End of explanation
X_train = X_clean.loc[(X['Year']<2012),:]
y_train = y_clean.loc[(X['Year']<2012)]
X_va = X_clean.loc[X['Year'].isin([2012, 2013]),:]
y_va = y_clean.loc[X['Year'].isin([2012, 2013])]
X_test = X_clean.loc[X['Year']>2013,:]
y_test = y_clean.loc[X['Year']>2013]
Explanation: Split into training, validation, testing
End of explanation
X_train_scaled = StandardScaler().fit_transform(X_train)
X_va_scaled = StandardScaler().fit_transform(X_va)
X_test_scaled = StandardScaler().fit_transform(X_test)
Explanation: Need scaled versions of the X data for some of the models
End of explanation
print X_train_scaled.shape, y_train.shape
print X_va_scaled.shape, y_va.shape
print X_test_scaled.shape, y_test.shape
Explanation: Check size of all arrays
End of explanation
lm = LinearRegression()
lm.fit(X_train_scaled, y_train)
lm.score(X_va_scaled, y_va)
y_pr = lm.predict(X_va_scaled)
y_va.values.shape, y_pr.shape, X.loc[X['Year'].isin([2012, 2013]),'cluster'].values.shape
y_lm_resids = pd.DataFrame(dict(zip(['Gen Change (MW)', 'y_pr', 'cluster'],
                               [y_va.values, y_pr, X.loc[X['Year'].isin([2012, 2013]),'cluster'].values])))
# y_lm_resids['y_pr'] = y_pr
# y_lm_resids['cluster'] = X.loc[:,'cluster']
y_lm_resids.head()
y_lm_resids.loc[:,'residuals'] = y_lm_resids.loc[:,'y_pr'] - y_lm_resids.loc[:,'Gen Change (MW)']
g = sns.FacetGrid(y_lm_resids, hue='cluster', col='cluster',
                  col_wrap=3)
g.map(plt.scatter, 'Gen Change (MW)', 'residuals')
g.add_legend()
Explanation: Linear Regression (OLS)
End of explanation
from xgboost import XGBRegressor
Explanation: XGBoost
End of explanation
param_values = [25, 100, 250, 350]
train_scores, valid_scores = validation_curve(XGBRegressor(), X_train, y_train, "n_estimators", param_values,
                                              n_jobs=-1, verbose=3)
train_scores_mean = np.mean(train_scores, axis=1)
train_scores_std = np.std(train_scores, axis=1)
valid_scores_mean = np.mean(valid_scores, axis=1)
valid_scores_std = np.std(valid_scores, axis=1)
plt.title("Validation Curve with XGBoost", size=15)
plt.xlabel("n_estimators", size=15)
plt.ylabel("Score", size=15)
plt.ylim(0.0, 1.1)
lw = 2
plt.plot(param_values, train_scores_mean, label="Training score",
             color="darkorange", lw=lw)
plt.fill_between(param_values, train_scores_mean - train_scores_std,
                 train_scores_mean + train_scores_std, alpha=0.2,
                 color="darkorange", lw=lw)
plt.plot(param_values, valid_scores_mean, label="Cross-validation score",
             color="navy", lw=lw)
plt.fill_between(param_values, valid_scores_mean - valid_scores_std,
                 valid_scores_mean + valid_scores_std, alpha=0.2,
                 color="navy", lw=lw)
plt.legend(loc="best")
plt.savefig('XGBoost n_estimators validation curve.pdf', bbox_inches='tight')
Explanation: Validation curve for n_estimators
End of explanation
param_values = [1,3,5,9,15]
train_scores, valid_scores = validation_curve(XGBRegressor(n_estimators=250), X_train, y_train, "max_depth", param_values,
                                              n_jobs=-1, verbose=3)
train_scores_mean = np.mean(train_scores, axis=1)
train_scores_std = np.std(train_scores, axis=1)
valid_scores_mean = np.mean(valid_scores, axis=1)
valid_scores_std = np.std(valid_scores, axis=1)
plt.title("Validation Curve with XGBoost", size=15)
plt.xlabel("max_depth", size=15)
plt.ylabel("Score", size=15)
plt.ylim(0.0, 1.1)
lw = 2
plt.plot(param_values, train_scores_mean, label="Training score",
             color="darkorange", lw=lw)
plt.fill_between(param_values, train_scores_mean - train_scores_std,
                 train_scores_mean + train_scores_std, alpha=0.2,
                 color="darkorange", lw=lw)
plt.plot(param_values, valid_scores_mean, label="Cross-validation score",
             color="navy", lw=lw)
plt.fill_between(param_values, valid_scores_mean - valid_scores_std,
                 valid_scores_mean + valid_scores_std, alpha=0.2,
                 color="navy", lw=lw)
plt.legend(loc="best")
plt.savefig('XGBoost max_depth validation curve.pdf', bbox_inches='tight')
Explanation: Validation curve for n_estimators
End of explanation
param_values = np.logspace(-5, 1, 7)
train_scores, valid_scores = validation_curve(XGBRegressor(n_estimators=250), X_train, y_train, "reg_alpha", param_values,
                                              n_jobs=-1, verbose=3)
train_scores_mean = np.mean(train_scores, axis=1)
train_scores_std = np.std(train_scores, axis=1)
valid_scores_mean = np.mean(valid_scores, axis=1)
valid_scores_std = np.std(valid_scores, axis=1)
plt.title("Validation Curve with XGBoost")
plt.xlabel("reg_alpha")
plt.ylabel("Score")
plt.ylim(0.0, 1.1)
lw = 2
plt.semilogx(param_values, train_scores_mean, label="Training score",
             color="darkorange", lw=lw)
plt.fill_between(param_values, train_scores_mean - train_scores_std,
                 train_scores_mean + train_scores_std, alpha=0.2,
                 color="darkorange", lw=lw)
plt.semilogx(param_values, valid_scores_mean, label="Cross-validation score",
             color="navy", lw=lw)
plt.fill_between(param_values, valid_scores_mean - valid_scores_std,
                 valid_scores_mean + valid_scores_std, alpha=0.2,
                 color="navy", lw=lw)
plt.legend(loc="best")
Explanation: Validation curve for reg_alpha
End of explanation
param_values = [1,3,5,9,15]
train_sizes, train_scores, valid_scores = learning_curve(XGBRegressor(n_estimators=250), X_train, y_train,
                                              n_jobs=-1, verbose=3)
train_scores_mean = np.mean(train_scores, axis=1)
train_scores_std = np.std(train_scores, axis=1)
valid_scores_mean = np.mean(valid_scores, axis=1)
valid_scores_std = np.std(valid_scores, axis=1)
plt.title("Learning Curve with XGBoost", size=15)
plt.xlabel("Sample size", size=15)
plt.ylabel("Score", size=15)
plt.ylim(0.0, 1.1)
lw = 2
plt.plot(train_sizes, train_scores_mean, label="Training score",
             color="darkorange", lw=lw)
plt.fill_between(train_sizes, train_scores_mean - train_scores_std,
                 train_scores_mean + train_scores_std, alpha=0.2,
                 color="darkorange", lw=lw)
plt.plot(train_sizes, valid_scores_mean, label="Cross-validation score",
             color="navy", lw=lw)
plt.fill_between(train_sizes, valid_scores_mean - valid_scores_std,
                 valid_scores_mean + valid_scores_std, alpha=0.2,
                 color="navy", lw=lw)
plt.legend(loc="best")
plt.savefig('XGBoost learning curve.pdf', bbox_inches='tight')
xgbr = XGBRegressor(n_estimators=250)
xgbr.fit(X_train, y_train)
y_pr = xgbr.predict(X_va)
y_xgbr_resids = pd.DataFrame(dict(zip(['Gen Change (MW)', 'y_pr', 'cluster'],
                               [y_va.values, y_pr, X.loc[X['Year'].isin([2012, 2013]),'cluster'].values])))
y_xgbr_resids.loc[:,'residuals'] = y_xgbr_resids.loc[:,'y_pr'] - y_xgbr_resids.loc[:,'Gen Change (MW)']
plt.scatter()
with sns.axes_style('whitegrid'):
    g = sns.FacetGrid(y_xgbr_resids, hue='cluster', col='cluster',
                      col_wrap=3)
    g.map(plt.scatter, 'y_pr', 'residuals', s=5, alpha=.3)
    g.set_xlabels(size=15)
    g.set_ylabels(size=15)
    plt.savefig('XGBR residuals.pdf')
model = XGBRegressor()
subsample = [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 1.0]
param_grid = dict(subsample=subsample)
grid_search = GridSearchCV(model, param_grid, n_jobs=-1, verbose=3)
result = grid_search.fit(X_train_scaled, y_train)
result.cv_results_
model = XGBRegressor()
colsample_bytree = [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 1.0]
param_grid = dict(colsample_bytree=colsample_bytree)
grid_search = GridSearchCV(model, param_grid, n_jobs=-1, verbose=3)
result = grid_search.fit(X_train_scaled, y_train)
result.cv_results_
model = XGBRegressor()
colsample_bylevel = [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 1.0]
param_grid = dict(colsample_bylevel=colsample_bylevel)
grid_search = GridSearchCV(model, param_grid, n_jobs=-1, verbose=3)
result = grid_search.fit(X_train_scaled, y_train)
result.cv_results_
model = XGBRegressor()
max_depth = [3, 6, 9]
n_estimators = [100, 250, 500]
reg_alpha = [1e-5, 1e-3, 0.1]
reg_lambda = [1e-3, 0.1, 1]
param_grid = dict(max_depth=max_depth, n_estimators=n_estimators,
                  reg_alpha=reg_alpha, reg_lambda=reg_lambda)
grid_search = GridSearchCV(model, param_grid, n_jobs=-1, verbose=2)
result = grid_search.fit(X_train_scaled, y_train)
import cPickle as pickle
pickle.dump((grid_search, result), open( "xgb gridsearch and results.pkl", "wb" ) )
result.cv_results_
result.best_estimator_
grid_search.score(X_va_scaled, y_va)
Explanation: Learning curve for n_estimators=250 and max_depth=3
End of explanation
xgb = XGBRegressor(n_estimators=250, reg_alpha=0.1)
xgb.fit(X_train, y_train)
xgb.score(X_va, y_va)
y_pr = xgb.predict(X_va)
y_xgb_resids = pd.DataFrame(dict(zip(['Gen Change (MW)', 'y_pr', 'cluster'],
                               [y_va.values, y_pr, X.loc[X['Year'].isin([2012, 2013]),'cluster'].values])))
y_xgb_resids.loc[:,'residuals'] = y_xgb_resids.loc[:,'y_pr'] - y_xgb_resids.loc[:,'Gen Change (MW)']
g = sns.FacetGrid(y_xgb_resids, hue='cluster', col='cluster',
                  col_wrap=3)
g.map(plt.scatter, 'y_pr', 'residuals')
g.add_legend()
y_xgb_resids.describe()
Explanation: Try XGBoost on non-scaled data
Turns out this works better
End of explanation
X_train.columns
X_train_ratio = X_train.copy()
X_va_ratio = X_va.copy()
for fuel in ['All coal', 'Lignite', 'Subbituminous']:
    X_train_ratio.loc[:,fuel] = X_train_ratio.loc[:,fuel].values/X_train_ratio.loc[:,'NG Price ($/mcf)'].values
    X_va_ratio.loc[:,fuel] = X_va.loc[:,fuel]/X_va.loc[:,'NG Price ($/mcf)']
    
X_train_ratio.drop('NG Price ($/mcf)', axis=1, inplace=True)
X_va_ratio.drop('NG Price ($/mcf)', axis=1, inplace=True)
X_train.head()
X_train_ratio.head()
xgb_ratio = XGBRegressor(n_estimators=250, reg_alpha=0.1)
xgb_ratio.fit(X_train_ratio, y_train)
xgb_ratio.score(X_va_ratio, y_va)
y_pr = xgb_ratio.predict(X_va_ratio)
y_xgb_resids = pd.DataFrame(dict(zip(['Gen Change (MW)', 'y_pr', 'cluster'],
                               [y_va.values, y_pr, X.loc[X['Year'].isin([2012, 2013]),'cluster'].values])))
y_xgb_resids.loc[:,'residuals'] = y_xgb_resids.loc[:,'y_pr'] - y_xgb_resids.loc[:,'Gen Change (MW)']
g = sns.FacetGrid(y_xgb_resids, hue='cluster', col='cluster',
                  col_wrap=3)
g.map(plt.scatter, 'y_pr', 'residuals')
g.add_legend()
y_xgb_resids.describe()
from xgboost import plot_importance
plot_importance(xgb_ratio)
Explanation: Try ratio of fuel prices
End of explanation
lm = LinearRegression(normalize=True)
lm.fit(X_train_ratio, y_train)
lm.score(X_va_ratio, y_va)
y_pr = lm.predict(X_va_scaled)
y_va.values.shape, y_pr.shape, X.loc[X['Year'].isin([2012, 2013]),'cluster'].values.shape
y_lm_resids = pd.DataFrame(dict(zip(['Gen Change (MW)', 'y_pr', 'cluster'],
                               [y_va.values, y_pr, X.loc[X['Year'].isin([2012, 2013]),'cluster'].values])))
# y_lm_resids['y_pr'] = y_pr
# y_lm_resids['cluster'] = X.loc[:,'cluster']
y_lm_resids.head()
y_lm_resids.loc[:,'residuals'] = y_lm_resids.loc[:,'y_pr'] - y_lm_resids.loc[:,'Gen Change (MW)']
g = sns.FacetGrid(y_lm_resids, hue='cluster', col='cluster',
                  col_wrap=3)
g.map(plt.scatter, 'Gen Change (MW)', 'residuals')
g.add_legend()
Explanation: Linear Regression (OLS) with ratio of fuel prices
Slight improvement from first OLS try, but still nowhere near as good as gradient boosting
End of explanation
<END_TASK>