2.- Bike Sharing: Predicción de Demanda Horaria

a) Cargue los datos de entrenamiento y pruebas como dataframes de pandas. Describa las variables involucradas en el problema, explorando el tipo de datos de que se trata, el número de valores distintos y, si corresponde, un gráfico (e.g. un histograma) que resuma su comportamiento. Su primera operación de pre-procesamiento de datos será obtener la hora del día desde el campo fecha (que en este momento es de tipo string), creando una nueva columna denominada hour y de tipo int. Para hacer esta operación se concatenarán los dataframes de entrenamiento y pruebas y luego se volverán a separar manteniendo la separación original.


In [14]:
import pandas as pd
import numpy as np

dftrain = pd.read_csv('data/bike_sharing_train.csv')
dfval = pd.read_csv('data/bike_sharing_val.csv')
dftest = pd.read_csv('data/bike_sharing_test.csv')
ntrain = len(dftrain)
nval = len(dftrain) + len(dfval)
df = pd.concat([dftrain,dfval,dftest])

print('\nSummary - dataframe completo:\n')
print(df.describe())
df.head()


Summary - dataframe completo:

         Unnamed: 0        season       holiday    workingday       weather  \
count  10886.000000  10886.000000  10886.000000  10886.000000  10886.000000   
mean    5442.500000      2.506614      0.028569      0.680875      1.418427   
std     3142.661849      1.116174      0.166599      0.466159      0.633839   
min        0.000000      1.000000      0.000000      0.000000      1.000000   
25%     2721.250000      2.000000      0.000000      0.000000      1.000000   
50%     5442.500000      3.000000      0.000000      1.000000      1.000000   
75%     8163.750000      4.000000      0.000000      1.000000      2.000000   
max    10885.000000      4.000000      1.000000      1.000000      4.000000   

              temp         atemp      humidity     windspeed        casual  \
count  10886.00000  10886.000000  10886.000000  10886.000000  10886.000000   
mean      20.23086     23.655084     61.886460     12.799395     36.021955   
std        7.79159      8.474601     19.245033      8.164537     49.960477   
min        0.82000      0.760000      0.000000      0.000000      0.000000   
25%       13.94000     16.665000     47.000000      7.001500      4.000000   
50%       20.50000     24.240000     62.000000     12.998000     17.000000   
75%       26.24000     31.060000     77.000000     16.997900     49.000000   
max       41.00000     45.455000    100.000000     56.996900    367.000000   

         registered         count  
count  10886.000000  10886.000000  
mean     155.552177    191.574132  
std      151.039033    181.144454  
min        0.000000      1.000000  
25%       36.000000     42.000000  
50%      118.000000    145.000000  
75%      222.000000    284.000000  
max      886.000000    977.000000  
Out[14]:
Unnamed: 0 datetime season holiday workingday weather temp atemp humidity windspeed casual registered count
0 3 2011-01-01 03:00:00 1 0 0 1 9.84 14.395 75 0.0000 3 10 13
1 4 2011-01-01 04:00:00 1 0 0 1 9.84 14.395 75 0.0000 0 1 1
2 5 2011-01-01 05:00:00 1 0 0 2 9.84 12.880 75 6.0032 0 1 1
3 6 2011-01-01 06:00:00 1 0 0 1 9.02 13.635 80 0.0000 2 0 2
4 7 2011-01-01 07:00:00 1 0 0 1 8.20 12.880 86 0.0000 1 2 3

In [15]:
import matplotlib.pyplot as plt

for column in ['season','workingday','weather','temp', 'atemp', 'humidity', 'windspeed','casual','registered']:
    g = df.groupby(column)['count','registered','casual'].mean()
    if len(g) > 4:
        g.plot()
    else:
        g.plot.bar()
    plt.show()



In [16]:
df['hour'] = pd.to_datetime(df['datetime']).apply(lambda x: x.strftime('%H'))
df['hour'] = pd.to_numeric(df['hour'])

df['day'] = pd.to_datetime(df['datetime']).apply(lambda x: x.strftime('%d'))
df['day'] = pd.to_numeric(df['day'])

df['month'] = pd.to_datetime(df['datetime']).apply(lambda x: x.strftime('%m'))
df['month'] = pd.to_numeric(df['month'])

df['year'] = pd.to_datetime(df['datetime']).apply(lambda x: x.strftime('%Y'))
df['year'] = pd.to_numeric(df['year'])

for column in ['hour','day','month','year']:
    g = df.groupby(column)['count'].mean()
    g.plot('bar')
    plt.show()



In [17]:
Xdf=df.loc[:,['season','holiday','workingday','weather','temp','atemp','humidity','windspeed','hour']]
Ydf=df.loc[:,'count']
X_train = Xdf[0:ntrain].values
X_val = Xdf[ntrain:nval].values
X_test = Xdf[nval:].values
Y_train = Ydf[0:ntrain].values
Y_val = Ydf[ntrain:nval].values
Y_test = Ydf[nval:].values

b) Entrene un árbol de regresión para resolver el problema usando parámetros por defecto. Con este fin, construya una matriz $X_{train}$ de forma $n_{train} \times d_1$ que contenga los datos de entrenamiento en sus filas, seleccionando las columnas que desee/pueda utilizar para el entrenamiento. Implemente además, la función de evaluación que hemos definido anteriormente para este problema. Evalúe el árbol de regresión ajustado a los datos de entrenamiento sobre el conjunto de entrenamiento y pruebas. Construya un gráfico que compare las predicciones con los valores reales. En este punto usted debiese tener un modelo con puntaje del orden de 0.59, lo que lo dejará más o menos en la posición 2140 de la competencia.


In [18]:
from sklearn.tree import DecisionTreeRegressor
import matplotlib.pyplot as plt

def eval_bikemodel(y_predict,y_true,**kwargs):
    diff = np.log(y_predict+1.0) - np.log(y_true+1.0)
    return np.sqrt(np.sum(np.square(diff))/len(y_predict))

model = DecisionTreeRegressor(random_state=0)
model.fit(X_train,Y_train)
score_test = model.score(X_test,Y_test)
print("SCORE TEST=%f"%score_test)

Y_pred_train = model.predict(X_train)
Y_pred_val = model.predict(X_val)
Y_pred_test = model.predict(X_test)
kagg_train = eval_bikemodel(Y_pred_train,Y_train)
kagg_val = eval_bikemodel(Y_pred_val,Y_val)
kagg_test = eval_bikemodel(Y_pred_test,Y_test)

print("KAGG EVAL TRAIN =%f"%kagg_train)
print("KAGG EVAL TEST =%f"%kagg_test)
plt.plot(Y_test,Y_pred_test,'.')
plt.show()


SCORE TEST=0.703388
KAGG EVAL TRAIN =0.028516
KAGG EVAL TEST =0.574239

c) Mejore el árbol de regresión definido en el punto anterior haciendo modificaciones a los hiper-parámetros del modelo. Por ejemplo, como estos modelos tienden a sobre-ajustar, podría intentar limitar la profundidad del árbol (¿Por qué esto debiese ayudar?). Naturalmente, está absolutamente prohibido tomar este tipo de decisiones en función del resultado de pruebas. Debe realizar estas elecciones evaluando sobre el conjunto de validación. Si no desea utilizarlo, y prefiere implementar validación cruzada u otra técnica automática, tiene la ventaja de poder usar el conjunto de validación como parte del entrenamiento. Con estas modificaciones debiese poder mejorar su ranking en unas 300 posiciones.


In [19]:
from sklearn.model_selection import cross_val_score
from sklearn.metrics import make_scorer

def train_tree(X,Y,crit=["mse","mae"],max_depth=30,min_depth=5):
    scores = []
    count,current = len(crit)*(max_depth-min_depth+1),1
    for d in range(min_depth,max_depth+1):
        for c in crit:
            model = DecisionTreeRegressor(random_state=0,max_depth=d,criterion=c)
            score = 0#cross_val_score(model, X,Y).mean()
            score_bike = cross_val_score(model, X,Y, scoring=make_scorer(eval_bikemodel)).mean()
            scores.append([model, score, score_bike, c, d])
            # print("\t Score:%f \t BikeScore:%f \t Max Depth: %d \t Criterion: %s"%(score,score_bike,d,c))
            print ("%d/%d"%(current,count), end="\r")
            current+=1
    return scores

scores = train_tree(np.concatenate([X_train,X_val]),np.concatenate([Y_train,Y_val]),max_depth=20,min_depth=5)


32/32

In [20]:
def top_scores(scores):
    sorted_scores = sorted(scores, key=(lambda x: -x[2]), reverse=False)
    print("Top 10 models:")
    for s in sorted_scores[-10:]:
        print("  Score:%f \t BikeScore:%f \t Max Depth: %d \t Criterion: %s"%(s[1],s[2],s[4],s[3]))
    return sorted_scores[-1]

t=top_scores(scores)


Top 10 models:
  Score:0.000000 	 BikeScore:0.592093 	 Max Depth: 13 	 Criterion: mae
  Score:0.000000 	 BikeScore:0.588915 	 Max Depth: 12 	 Criterion: mse
  Score:0.000000 	 BikeScore:0.583421 	 Max Depth: 9 	 Criterion: mae
  Score:0.000000 	 BikeScore:0.582569 	 Max Depth: 8 	 Criterion: mae
  Score:0.000000 	 BikeScore:0.580230 	 Max Depth: 11 	 Criterion: mse
  Score:0.000000 	 BikeScore:0.579067 	 Max Depth: 11 	 Criterion: mae
  Score:0.000000 	 BikeScore:0.577337 	 Max Depth: 10 	 Criterion: mae
  Score:0.000000 	 BikeScore:0.573863 	 Max Depth: 8 	 Criterion: mse
  Score:0.000000 	 BikeScore:0.570939 	 Max Depth: 10 	 Criterion: mse
  Score:0.000000 	 BikeScore:0.566530 	 Max Depth: 9 	 Criterion: mse

In [21]:
t[0].fit(np.concatenate([X_train,X_val]),np.concatenate([Y_train,Y_val]))
print("Top BikeScore Test: %f"%eval_bikemodel(t[0].predict(X_test),Y_test))


Top BikeScore Test: 0.490544

Se utilizó validación cruzada con la finalidad de no guiarse por una sola evaluación, la que podría ser "buena" por cosa de suerte.

Además se varió la profundidad del árbol, teniendo en cuanta que un número muy alto causará overfitting del modelo. Junto s esto, también se varió el criterio de decisión utilizado.

Como se puede ver en esta lista, las mejores configuraciones fueron mse ó friedman_mse con una profundidad de 9.

d) Mejore el árbol de regresión definido en el punto anterior haciendo modificaciones sobre la representación utilizada para aprender desde los datos. Por ejemplo, los histogramas que construyó en el punto (a) así como la forma especial de la función de evaluación, sugieren una cierta transformación de la variable respuesta. Podría intentar también normalizando los datos o normalizando la respuesta. Otra opción es intentar rescatar algo más acerca de la fecha (anteriormente sólo se extrajo la hora), como por ejemplo el año o el día de la semana ('lunes','martes', etc) que corresponde. Sea creativo, este paso le debiese reportar un salto de calidad muy significativo. Una observación importante es que si hace una transformación a la variable respuesta (por ejemplo raíz cuadrada), debe invertir esta transformación antes de evaluar el desempeño con eval_bikemodel (por ejemplo, elevar al cuadrado si tomó raíz cuadrada). Con modificaciones de este tipo, podría mejorar su ranking en unas 1000 posiciones, entrando ya al top-1000 con un score del orden de $0.45$.


In [22]:
df['cday'] = pd.to_datetime(df['datetime']).dt.dayofweek #0:lunes,6:domingo
df['cday'] = pd.to_numeric(df['cday'])

df['cyear'] = pd.to_datetime(df['datetime']).apply(lambda x: x.strftime('%Y'))
df['cyear'] = pd.to_numeric(df['cyear'])

df['cmonth'] = pd.to_datetime(df['datetime']).apply(lambda x: x.strftime('%m'))
df['cmonth'] = pd.to_numeric(df['cmonth'])

df['ghour'] = df['hour'].apply(lambda x: int(x/3))
df['gatemp'] = df['atemp'].apply(lambda x: int(x/3))

df['holiday'] = pd.to_numeric(df['holiday'])
df['workingday'] = pd.to_numeric(df['workingday'])
df['tday'] = df['holiday'] + 2*df['workingday']

df['humidity'] = pd.to_numeric(df['humidity'])

Xdf=df.loc[:,['season','weather','humidity','cday','hour','cyear','tday','gatemp']]
Xdf.head()


Out[22]:
season weather humidity cday hour cyear tday gatemp
0 1 1 75 5 3 2011 0 4
1 1 1 75 5 4 2011 0 4
2 1 2 75 5 5 2011 0 4
3 1 1 80 5 6 2011 0 4
4 1 1 86 5 7 2011 0 4

In [23]:
X_train2 = Xdf[0:nval].values
Y_train2 = Ydf[0:nval].values

X_test2 = Xdf[nval:].values
Y_test2 = Ydf[nval:].values

scores2 = train_tree(X_train2,Y_train2,crit=["mse","mae"],max_depth=15,min_depth=5)
t=top_scores(scores2)


Top 10 models:
  Score:0.000000 	 BikeScore:0.468551 	 Max Depth: 12 	 Criterion: mae
  Score:0.000000 	 BikeScore:0.467679 	 Max Depth: 10 	 Criterion: mae
  Score:0.000000 	 BikeScore:0.464680 	 Max Depth: 11 	 Criterion: mae
  Score:0.000000 	 BikeScore:0.464342 	 Max Depth: 15 	 Criterion: mse
  Score:0.000000 	 BikeScore:0.459172 	 Max Depth: 13 	 Criterion: mse
  Score:0.000000 	 BikeScore:0.459068 	 Max Depth: 14 	 Criterion: mse
  Score:0.000000 	 BikeScore:0.455209 	 Max Depth: 9 	 Criterion: mse
  Score:0.000000 	 BikeScore:0.451104 	 Max Depth: 12 	 Criterion: mse
  Score:0.000000 	 BikeScore:0.443767 	 Max Depth: 10 	 Criterion: mse
  Score:0.000000 	 BikeScore:0.441475 	 Max Depth: 11 	 Criterion: mse

In [24]:
t[0].fit(X_train2,Y_train2)
print("Top BikeScore Test: %f"%eval_bikemodel(t[0].predict(X_test2),Y_test2))


Top BikeScore Test: 0.412036

e) Entrene una SVM no lineal para resolver el problema midiendo el efecto de las distintas representaciones que haya descubierto hasta este punto. Un detalle importante es que antes de entrenar la SVM sería aconsejable hacer dos tipos de pre-procesamiento adicional de los datos: (i) codificar las variables categóricas en un modo apropiado - por ejemplo como vector binario con un 1 en la posición del valor adoptado-, (ii) escalar los atributos de modo que queden centrados y con rangos comparables. Usando parámetros por defecto para la SVM debiese obtener un score del orden de $0.344$, quedando definitivamente en el top-10 de la competencia.


In [25]:
def eval_bikemodel_log(y_predict,y_true,**kwargs):
    y_predict = np.exp(y_predict)-1
    y_true    = np.exp(y_true)-1
    diff = np.log(y_predict+1.0) - np.log(y_true+1.0)
    return np.sqrt(np.sum(np.square(diff))/len(y_predict))

df['day_weekend'] = df['tday'].apply(lambda x: x==0)
df['day_holiday'] = df['tday'].apply(lambda x: x==1)
df['day_work'] = df['tday'].apply(lambda x: x==2)

df['season_spr'] = df['season'].apply(lambda x: x==1)
df['season_sum'] = df['season'].apply(lambda x: x==2)
df['season_fal'] = df['season'].apply(lambda x: x==3)
df['season_win'] = df['season'].apply(lambda x: x==4)

df['weather_clear'] = df['weather'].apply(lambda x: x==1)
df['weather_mist'] = df['weather'].apply(lambda x: x==2)
df['weather_snow'] = df['weather'].apply(lambda x: x==3)
df['weather_rain'] = df['weather'].apply(lambda x: x==4)


Xdf=pd.get_dummies(df[['humidity','windspeed','hour','atemp','tday','season','weather','cmonth','cyear']],
                   columns=['hour','season','weather','tday','cmonth','cyear'])

Ydf=np.log(df.loc[:,'count']+1)
X_train = Xdf[0:ntrain].values
X_val = Xdf[ntrain:nval].values
X_test = Xdf[nval:].values
Y_train = Ydf[0:ntrain].values
Y_val = Ydf[ntrain:nval].values
Y_test = Ydf[nval:].values

In [26]:
from sklearn.preprocessing import StandardScaler
scalerX = StandardScaler()
X_train = scalerX.fit_transform(X_train)
X_val = scalerX.fit_transform(X_val)
X_test = scalerX.transform(X_test)

from sklearn.svm import SVR
model = SVR()
model.fit(X_train,Y_train)
Y_pred_train = model.predict(X_train)
Y_pred_val = model.predict(X_val)
Y_pred_test = model.predict(X_test)

print("Top BikeScore Test: %f"%eval_bikemodel_log(Y_pred_test,Y_test))


Top BikeScore Test: 0.343511

f) Mejore la SVM definida en el punto anterior haciendo modificaciones a los hiper-parámetros de la máquina ($C$, $\epsilon$ o la misma función de kernel). Naturalmente, está absolutamente prohibido tomar este tipo de decisiones de diseño mirando el resultado de pruebas. Debe realizar estas elecciones evaluando sobre el conjunto de validación. Si no desea utilizarlo, y prefiere implementar validación cruzada u otra técnica automática, tiene la ventaja de poder usar el conjunto de validación como parte del entrenamiento.


In [27]:
def train_svm(X,Y,kernel=['linear','poly','rbf','sigmoid','precomputed'],max_c=20,min_c=1):
    scores = []
    c_range = range(min_c,max_c+1)
    e_range = [0.01,0.05,0.1,0.5,1,10]
    count,current = len(kernel)*len(c_range)*len(e_range),1
    for c in c_range:
        for k in kernel:
            for e in e_range:
                model = SVR(C=c,epsilon=e,kernel=k)
                score_bike = cross_val_score(model, X,Y, scoring=make_scorer(eval_bikemodel_log)).mean()
                scores.append([model, 0, score_bike, c, k, e])
                print ("%d/%d"%(current,count), end="\r")
                current+=1
    return scores

def top_scores_svm(scores):
    sorted_scores = sorted(scores, key=(lambda x: -x[2]), reverse=False)
    print("Top 10 models:")
    for s in sorted_scores[-10:]:
        print("  BikeScore:%f \t C: %d \t kernel: %s \t e: %f"%(s[2],s[3],s[4],s[5]))
    return sorted_scores[-1]


scores_svm = train_svm(np.concatenate([X_train,X_val]),np.concatenate([Y_train,Y_val]),kernel=['rbf'],max_c=5)
t=top_scores_svm(scores_svm)


Top 10 models:
  BikeScore:0.349680 	 C: 5 	 kernel: rbf 	 e: 0.100000
  BikeScore:0.348753 	 C: 1 	 kernel: rbf 	 e: 0.010000
  BikeScore:0.347904 	 C: 3 	 kernel: rbf 	 e: 0.050000
  BikeScore:0.347745 	 C: 2 	 kernel: rbf 	 e: 0.010000
  BikeScore:0.347035 	 C: 4 	 kernel: rbf 	 e: 0.100000
  BikeScore:0.346859 	 C: 1 	 kernel: rbf 	 e: 0.050000
  BikeScore:0.345501 	 C: 1 	 kernel: rbf 	 e: 0.100000
  BikeScore:0.344689 	 C: 2 	 kernel: rbf 	 e: 0.050000
  BikeScore:0.344506 	 C: 3 	 kernel: rbf 	 e: 0.100000
  BikeScore:0.342479 	 C: 2 	 kernel: rbf 	 e: 0.100000

In [28]:
model = SVR(C=t[3],kernel=t[4],epsilon=t[5])
model.fit(np.concatenate([X_train,X_val]),np.concatenate([Y_train,Y_val]))
Y_pred_test_svr = model.predict(X_test)
print("Top BikeScore Test: %f"%eval_bikemodel_log(Y_pred_test_svr,Y_test))


Top BikeScore Test: 0.331970

g) Evaúe el efecto de utilizar el dataset de validación para entrenamiento y seleccionar los parámetros estructurales del árbol de clasificación y la SVM usando validación cruzada. El código de ejemplo para esto ha sido proporcionado en las tareas 1 y 2, pero se adjunta de nuevo a continuación

Hecho en f)

h) Evalúe el efecto de utilizar un ensamblado de 2 máquinas de aprendizaje para predecir la demanda total de bicicletas. Un modelo se especializará en la predicción de la demanda de bicicletas de parte de usuarios registrados y otra en la predicción de la demanda de usuarios casuales. Hay razones claras para pensar que los patrones son distintos.


In [29]:
# Ydf=df.ix[:,'count'] #demanda total
Ydf_r=np.log(df['registered']+1) #demanda registrada
Yr_train = Ydf_r[0:ntrain].values
Yr_val = Ydf_r[ntrain:nval].values
Yr_test = Ydf_r[nval:].values

Ydf_c=np.log(df['casual']+1) #demanda casual
Yc_train = Ydf_c[0:ntrain].values
Yc_val = Ydf_c[ntrain:nval].values
Yc_test = Ydf_c[nval:].values

model_r = SVR(C=t[3],kernel=t[4],epsilon=t[5])
model_c = SVR(C=t[3],kernel=t[4],epsilon=t[5])

model_r.fit(np.concatenate([X_train,X_val]),np.concatenate([Yr_train,Yr_val]))
model_c.fit(np.concatenate([X_train,X_val]),np.concatenate([Yc_train,Yc_val]))

Yr_pred_test = model_r.predict(X_test)
Yc_pred_test = model_c.predict(X_test)

print("Top BikeScore Test registered: %f"%eval_bikemodel_log(Yr_pred_test,Yr_test))
print("Top BikeScore Test casual    : %f"%eval_bikemodel_log(Yc_pred_test,Yc_test))
print("Top BikeScore Test count     : %f"%eval_bikemodel_log(Yr_pred_test+Yc_pred_test,Y_test))


Top BikeScore Test registered: 0.329733
Top BikeScore Test casual    : 0.539797
Top BikeScore Test count     : 2.796318

i) Evalúe el efecto de utilizar un algoritmo genérico para ensamblar máquinas de aprendizaje para predecir la demanda total de bicicletas. Puede experimentar con una sola técnica (e.g. Random Forest), discuta la evolución a medida que aumenta el número de máquinas.


In [47]:
from sklearn.ensemble import RandomForestRegressor

def train_random_forest(Xtrain,Ytrain,Xeval,Yval, estimators, depth):
    model = RandomForestRegressor(n_estimators=estimators,max_depth=depth, random_state=0)
    model.fit(Xtrain, Ytrain)

    Yeval_pred = model.predict(Xeval)
    
    return (model,eval_bikemodel_log(Yeval_pred,Yval),estimators,depth)

def top_scores_rfr(scores):
    sorted_scores = sorted(scores, key=(lambda x: -x[1]), reverse=False)
    print("Top 10 models:")
    for s in sorted_scores[-10:]:
        print("  BikeScore:%f \t estimators: %d \t depth: %d"%(s[1],s[2],s[3]))
    return sorted_scores[-1]

scores = []
count,current=20*20,1
for e in range(1,20):
    for d in range(1,20):
        scores.append(train_random_forest(X_train,Y_train,X_val,Y_val, e, d))
        print ("%d/%d"%(current,count), end="\r")
        current+=1
        
t=top_scores_rfr(scores)


Top 10 models:
  BikeScore:0.374417 	 estimators: 10 	 depth: 19
  BikeScore:0.371400 	 estimators: 11 	 depth: 19
  BikeScore:0.371047 	 estimators: 12 	 depth: 19
  BikeScore:0.370860 	 estimators: 19 	 depth: 19
  BikeScore:0.370166 	 estimators: 16 	 depth: 19
  BikeScore:0.370093 	 estimators: 15 	 depth: 19
  BikeScore:0.370052 	 estimators: 13 	 depth: 19
  BikeScore:0.370002 	 estimators: 18 	 depth: 19
  BikeScore:0.369732 	 estimators: 17 	 depth: 19
  BikeScore:0.369007 	 estimators: 14 	 depth: 19

In [54]:
t[0].fit(X_train,Y_train)
Ytest_pred = t[0].predict(X_test)

print("Top BikeScore Test: %f"%eval_bikemodel_log(Ytest_pred,Y_test))


Top BikeScore Test: 0.412570

In [ ]: