Predicting DC's bikeshare demand

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.ensemble import RandomForestRegressor

pd.set_option('display.width', 500)
pd.set_option('display.max_columns', 30)

sns.set_style("darkgrid")
sns.set_context(rc={"figure.figsize": (6, 3)})

dfs = {}
for name in ['train', 'test']:
    df = pd.read_csv('../data/%s.csv' % name)
    df['_data'] = name
    dfs[name] = df

# combine train and test data into one df
df = dfs['train'].append(dfs['test'])

# lowercase column names
df.columns = map(str.lower, df.columns)

df.head()

# parse datetime colum & add new time related columns
dt = pd.DatetimeIndex(df['datetime'])

df['date'] = dt.date
df['day'] = dt.day
df['month'] = dt.month
df['year'] = dt.year
df['hour'] = dt.hour
df['dow'] = dt.dayofweek
df['woy'] = dt.weekofyear
df['weeks_since_start'] = (dt - min(dt)).days / 7

# logarithmic transformation of dependent cols
# (adding 1 first so that 0 values don't become -inf)
for col in ['casual', 'registered', 'count']:
    df['%s_log' % col] = np.log(df[col] + 1)

df.iloc[0]

# let's plot some stuff to get a feel for the data

# get training data
dg = df[df['_data'] == 'train'].copy()

# distribution of casual / registered rentals by day
by_day = dg.groupby('date')[['casual', 'registered']].agg(sum)

print by_day.describe()

by_day.hist(figsize=(10, 3));

# total rentals by month & year
dg.groupby(['year', 'month'])['count'].agg(sum).plot(kind='bar');

# casual / registered rentals by day of the week
by_dow = dg.groupby('dow')[['casual', 'registered']].agg(sum)

by_dow.plot(kind='bar', width=0.8);

# rentals by hour, split by working day (or not)
by_hour = dg.groupby(['hour', 'workingday'])['count'].agg('sum').unstack()

by_hour.plot(kind='bar', figsize=(8,4), width=0.8);

# avg hourly rentals by temperature
dg['temp_int'] = np.round(dg['temp'])

dg.groupby('temp_int')['count'].agg('median').plot();

# correlation matrix
cols = [
    'count', 'casual', 'registered',
    'weather', 'temp', 'atemp', 'humidity', 'windspeed',
    'holiday', 'workingday', 'season', 
    'day', 'month', 'year', 'hour', 'dow', 'woy',
]

dg[cols].corr()

# instead of randomly splitting our training data 
# for cross validation, let's construct a framework that's more
# in line with how the data is divvied up for this competition
# (given first 19 days of each month, what is demand for remaining days)
# so, let's split our training data into 2 time contiguous datasets
# for fitting and validating our model (days 1-14 vs. days 15-19).

# also, since submissions are evaluated based on the
# root mean squared logarithmic error (RMSLE), let's replicate
# that computation as we test and tune our model.

def get_rmsle(y_pred, y_actual):
    diff = np.log(y_pred + 1) - np.log(y_actual + 1)
    mean_error = np.square(diff).mean()
    return np.sqrt(mean_error)


def get_data():
    data = df[df['_data'] == 'train'].copy()
    return data


def custom_train_test_split(data, cutoff_day=15):
    train = data[data['day'] <= cutoff_day]
    test = data[data['day'] > cutoff_day]
    
    return train, test


def prep_data(data, input_cols):
    X = data[input_cols].as_matrix()
    y_r = data['registered_log'].as_matrix()
    y_c = data['casual_log'].as_matrix()
    
    return X, y_r, y_c

    
def predict(input_cols, model_params={}):
    data = get_data()
    
    train, test = custom_train_test_split(data)
    
    X_train, y_train_r, y_train_c = prep_data(train, input_cols)
    X_test, y_test_r, y_test_c = prep_data(test, input_cols)

    model_params.update({
        'n_jobs': -1,
        'random_state': 123,
    })
    model = RandomForestRegressor(**model_params)
    
    model_r = model.fit(X_train, y_train_r)
    y_pred_r = np.exp(model_r.predict(X_test)) - 1
    
    model_c = model.fit(X_train, y_train_c)
    y_pred_c = np.exp(model_c.predict(X_test)) - 1

    y_pred_comb = np.round(y_pred_r + y_pred_c)
    y_pred_comb[y_pred_comb < 0] = 0
    
    y_test_comb = np.exp(y_test_r) + np.exp(y_test_c) - 2

    score = get_rmsle(y_pred_comb, y_test_comb)
    return score

# now that we've set that up, we can quickly try out different 
# combinations of features and have a pretty good sense for how 
# it would perform if we made a submission with that model

# for example, here's the RMSLE using the original input variables 
# (minus the datetime column)
cols = ['weather', 'temp', 'atemp', 'humidity', 
        'windspeed', 'holiday', 'workingday', 'season']
print predict(cols)

# now let's add the hour variable
cols = cols + ['hour']
print predict(cols)

# big improvement!

# let's try out a few more possibilities

base_cols = [
    'weather', 'temp', 'atemp', 'humidity', 'windspeed', 
    'holiday', 'workingday', 'season'
]

time_cols = ['hour', 'dow', 'year', 'month', 'weeks_since_start', 'day']

for i in range(0, len(time_cols) + 1):
    new_cols = time_cols[:i]
    all_cols = base_cols + new_cols
    print 'cols: base_cols + {}\nrmse: {}\n'.format(
        new_cols, 
        predict(all_cols)
    )

# now, let's hold the feature set constant and tune the
# model parameters to further improve accuracy

# feature set
x_cols = [
    'weather', 'temp', 'atemp', 'humidity', 'windspeed', 
    'holiday', 'workingday', 'season', 
    'hour', 'dow', 'year',
]

# random forest param grid
n_estimators = [500, 1000, 1500]
min_samples_splits = [6, 8, 10, 12, 14]

best_score, best_params = inf, None

# loop through param grid & find top performer
for ne in n_estimators:
    for mss in min_samples_splits:
        params = {'n_estimators': ne, 'min_samples_split': mss}
        score = predict(x_cols, params)
        print 'trees: {}, mss: {}, rmse: {}'.format(ne, mss, score)
        
        if score < best_score:
            best_params = params
            best_score = score
            
print 'best params: {}, rmse: {}'.format(best_params, best_score)

# putting it all together

# features (performs better on leaderboard w/o month var)
x_cols = [
    'weather', 'temp', 'atemp', 'humidity', 'windspeed', 
    'holiday', 'workingday', 'season', 
    'hour', 'dow', 'year'
]

# prepare training set
df_train = df[df['_data'] == 'train'].copy()
X_train = df_train[x_cols].as_matrix()
y_train_cas = df_train['casual_log'].as_matrix()
y_train_reg = df_train['registered_log'].as_matrix()

# prepare test set
df_test = df[df['_data'] == 'test'].copy()
X_test = df_test[x_cols].as_matrix()

# fit on training set
model = RandomForestRegressor(
    n_estimators = 1000, 
    min_samples_split = 12, 
    n_jobs = -1,
    random_state = 123456, 
)

# predict on test set & transform output back from log scale
casual_model = model.fit(X_train, y_train_cas)
y_pred_cas = casual_model.predict(X_test)
y_pred_cas = np.exp(y_pred_cas) - 1

registered_model = model.fit(X_train, y_train_reg)
y_pred_reg = registered_model.predict(X_test)
y_pred_reg = np.exp(y_pred_reg) - 1

# add casual & registered predictions together
df_test['count'] = np.round(y_pred_cas + y_pred_reg)

# output predictions for submission
final_df = df_test[['datetime', 'count']].copy()
final_df.to_csv('../output/predicted.csv', index=False)
print 'boom.'