O conjunto de variáveis (features) inicialmente disponíveis na base de dados são:
1) datetime - hourly date (data por hora) + timestamp (registro de data e hora)
2) season , tal que: 1 = spring, 2 = summer, 3 = fall, 4 = winter
3) holiday - se o dia é considerado feriado
4) workingday - se o dia não é considerado nem final de semana nem feriado
weather -
1: Clear, Few clouds, Partly cloudy, Partly cloudy (claro ou pouco nublado)
2: Mist + Cloudy, Mist + Broken clouds, Mist + Few clouds, Mist (nevoeiro, neblina, nublado)
3: Light Snow, Light Rain + Thunderstorm + Scattered clouds, Light Rain + Scattered clouds (Neve fraca, Chuva fraca + Trovoada + Nuvens dispersas, Chuva fraca + Nuvens dispersas)
4: Heavy Rain + Ice Pallets + Thunderstorm + Mist, Snow + Fog (chuva torrencial, Blocos de neve, trovoadas + neve + nevoeiro intenso)
temp - temperatura em Celsius
atemp - "sensação térmica" em Celsius
humidity - umidade relativa
windspeed - velocidade do vento
casual - número de aluguéis realizados por usuários não registrados
registered - número de aluguéis realizados por usuários registrados
count - número total de aluguéis
Estatísticas descritivas
Começamos o exercício com uma análise descritiva dos dados para, em sequência, avaliar quais são as características mais relevantes para o nosso modelo de previsão de demandas de bicicletas compartilhadas:
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from matplotlib import style
import seaborn as sns
train_df = pd.read_csv('train.csv', parse_dates=['datetime'])
test_df = pd.read_csv('test.csv', parse_dates=['datetime'])
df=train_df.copy()
train_df.info() # para saber se as caracteristicas são object, int ou float
train_df.sample(10)
train_df.head(10)
train_df.tail(10)
train_df.shape
Assim, geramos a primeira saída, conforme figura abaixo:
Explora-se, a seguir, cada uma das variáveis discretas. A ver, Season, Holiday, Working Day e Weather. Note que a demanda por bicicletas compartilhadas é sazonal, e bastante sensível ao clima observado no dia. Além disso, nota-se uma maior demanda nos dias comerciais, i.e., Working Day tanto para o ano de 2011 quanto para 2012.
# SEASON
print(df.season.value_counts())
sns.factorplot(x='season',data=df,kind='count',size=5,aspect=0.75)
#WEATHER
df.weather.value_counts()
sns.factorplot(x='weather',data=df,kind='count',size=5,aspect=0.75)
#HOLIDAY vs WORKING DAY
def plot_by_hour(data, year=None, agg='sum'):
dd = data.copy()
if year: dd = dd[ dd.datetime.dt.year == year ]
dd.loc[:, ('hour')] = dd.datetime.dt.hour
by_hour = dd.groupby(['hour', 'workingday'])['count'].agg(agg).unstack()
return by_hour.plot(kind='bar', ylim=(0, 80000), figsize=(8,5), width=0.75, title="Year = {0}".format(year))
plot_by_hour(df, year = 2011);
plot_by_hour(df, year = 2012);
# ALUGUEL DE BICILETAS POR MÊS EM 2011 E 2012
def plot_by_year(agg_attr, title):
dd = df.copy()
dd['year'] = df.datetime.dt.year
dd['month'] = df.datetime.dt.month
dd['hour'] = df.datetime.dt.hour
by_year = dd.groupby([agg_attr, 'year'])['count'].agg('sum').unstack()
return by_year.plot(kind='bar', figsize=(8,5), width=0.9, title=title)
plot_by_year('month', "Aluguel de bikes por mês - 2011 e 2012");
plot_by_year('hour', "Aluguel de bikes por hora - 2011 e 2012");
# ALUGUEL DE BICILETAS POR HORA EM 2011 E 2012
def plot_hours(data, message = ''):
dd = data.copy()
dd['hour'] = data.datetime.dt.hour
hours = {}
for hour in range(24):
hours[hour] = dd[ dd.hour == hour ]['count'].values
plt.figure(figsize=(8,5))
plt.ylabel("Quantidade de aluguel (Count rent)")
plt.xlabel("Horas")
plt.title("count vs hours\n" + message)
plt.boxplot( [hours[hour] for hour in range(24)] )
axis = plt.gca()
axis.set_ylim([1, 1100])
plot_hours( df[df.workingday == 1], 'working day')
plot_hours( df[df.workingday == 0], 'non working day')
plot_hours( df[df.datetime.dt.year == 2011], 'year 2011');
plot_hours( df[df.datetime.dt.year == 2012], 'year 2012');
Agora, apresentamos o boxplot da quantidade de demanda de bicicletas compartilhadas por hora e ao longo dos anos de 2011 e 2012 para os working days e non working days:
Por sua vez, a distribuição das variáveis contínuas pode ser vista a seguir:
#EXPLORANDO A DISTRIBUIÇÃO DAS VARIÁVEIS CONTINUAS
df.describe()
# Apenas para enriquecer a análise descritiva
sns.boxplot(data=df[['temp',
'atemp', 'humidity', 'windspeed', 'casual', 'registered', 'count']])
fig=plt.gcf()
fig.set_size_inches(10,10)
Feito isso, podemos agora analisar as inter-relações entre as nossas variáveis:
#matriz de correlação
cor_mat= df[:].corr()
mask = np.array(cor_mat)
mask[np.tril_indices_from(mask)] = False
fig=plt.gcf()
fig.set_size_inches(30,12)
sns.heatmap(data=cor_mat,mask=mask,square=True,annot=True,cbar=True)
Pela figura acima, podemos inferir:
a) A relação entre uma característica e ela mesma é igual a 1 (esperado);
b) temp e atemp são altamente correlacionadas (esperado);
c) humidity e count são negativamente correlacionadas, o que reforça o fato
de as pessoas demandarem menos bicicletas quando há umidade elevada;
d) casual e working day são fortemente correlacionadas de forma inversa
e) count e holiday são altamente correlacionadas de forma inversa
f) temp(ou atemp) tem alto impacto sobre count.
g) weather e count são altamente correlacionados de forma inversa. Isso é consequência de a característica weather ser medida de 1 até 4. Na medida em que as condições climáticas se tornam menos favoráveis, o n. de bicicletas alugadas diminui.
h) registered, casual e count são altamente correlacionadas (de forma positiva) isso sugere que as bicicletas que são alugadas também são registradas.
Feature Engineering
Aqui, inicia-se o processo de avaliação e filtragem das características que são mais relevantes. Desta forma, poderemos fazer um drop, isto é, descartar do modelo aquelas características não relevantes para se aprimorar o poder de previsão da demanda por bicicletas compartilhadas.
Iniciamos seperando a variável season de acordo com seus valores. Isso tornará possível o realce de outras características que afetam o modelo de previsão.
season=pd.get_dummies(df['season'],prefix='season')
df=pd.concat([df,season],axis=1)
df.head()
season=pd.get_dummies(test_df['season'],prefix='season')
test_df=pd.concat([test_df,season],axis=1)
test_df.head()
Fazemos o mesmo para a característica weather
weather=pd.get_dummies(df['weather'],prefix='weather')
df=pd.concat([df,weather],axis=1)
df.head()
weather=pd.get_dummies(test_df['weather'],prefix='weather')
test_df=pd.concat([test_df,weather],axis=1)
test_df.head()
Feito isso, uma vez que foram redefinidas, podemos eliminar season e weather das características do modelo.
df.drop(['season','weather'],inplace=True,axis=1)
df.head()
test_df.drop(['season','weather'],inplace=True,axis=1)
test_df.head()
O mesmo se aplica a registered e casual. IMPORTANTE: a divisão entre date e time deve ser bastante criteriosa, dado a importância em saber quais horas do dia tem-se maior demanda por bicicletas. Conforme visto anteriormente, as horas de maior demanda coincidem com o horário padrão de entrada e saída do trabalho (working days).
Agora, iremos avaliar como nossa variável a ser prevista, count, se relaciona com as novas características. É difícil enxergar claramente a correlação entre count e temp, por exemplo. O mais indicado pode ser discretizar temp em intervalos denominados bins e redefinir essa característica da seguinte forma:
new_df=df.copy()
new_df.temp.describe()
new_df['temp_bin']=np.floor(new_df['temp'])//5
new_df['temp_bin'].unique()
# então, podemos avaliar a partir da frequência
sns.factorplot(x="temp_bin",y="count",data=new_df,kind='bar')
A demanda é maior para os intervalos "bins" 6 e 7 que correspondem as temperaturas 30-35(bin 6) e 35-40 (bin 7). Agora adicionaremos colunas extras para casual, registered e count na base test_df:
test_df["casual"] = np.nan
test_df["registered"] = np.nan
test_df["count"] = np.nan
test_df.head(10)
train_df.dtypes
test_df.dtypes
train_df.describe() # resumo das estatisticas descritivas da amostra de treino
train_df.info() # verifica missing values
test_df.describe() # resumo das estatísticas descritivas da amostra de teste
train_df.info() # verifica missing values
Filtrando as características
year, month, hour e day of the week (DOW)" serão adicionadas ao data frame de treinamento como novas variáveis.
train_df['year'] = train_df['datetime'].dt.year
train_df['month'] = train_df['datetime'].dt.month
train_df['hour'] = train_df['datetime'].dt.hour
train_df['DOW'] = train_df['datetime'].dt.dayofweek
test_df['year'] = test_df['datetime'].dt.year
test_df['month'] = test_df['datetime'].dt.month
test_df['hour'] = test_df['datetime'].dt.hour
test_df['DOW'] = test_df['datetime'].dt.dayofweek
O método iloc será utilizado para avaliar essas novas variáveis:
train_df.iloc[:,0]
train_df.head(10)
Para encontrar o count do número de bikes alugadas no dia 4, por exemplo, temos que acessar as informações disponíveis na 4a linha e na 12a coluna:
train_df.iloc[3, 11]
Agora, vamos treinar os modelos com as variáveis workinday, temp, year, month, hour, DOW. A variável count é a quantidade que estamos tentando prever.
ind_variables_selected = ['workingday', 'temp', 'year', 'month', 'hour', 'DOW']
X_orig_train = train_df[ind_variables_selected]
y_orig_train = train_df['count']
X_test = test_df[ind_variables_selected]
y_test = test_df['count']
Criando os sets de treinamento e validação
#classifiaction.
from sklearn.linear_model import LogisticRegression
from sklearn.svm import LinearSVC,SVC
from sklearn.neighbors import KNeighborsClassifier
from sklearn.ensemble import RandomForestClassifier,GradientBoostingClassifier
from sklearn.tree import DecisionTreeClassifier
from sklearn.naive_bayes import GaussianNB
from sklearn.tree import DecisionTreeRegressor
from sklearn.ensemble import RandomForestRegressor
import xgboost as xgb
#regression
from sklearn.linear_model import LinearRegression,Ridge,Lasso,RidgeCV
from sklearn.ensemble import RandomForestRegressor,BaggingRegressor,GradientBoostingRegressor,AdaBoostRegressor
from sklearn.svm import SVR
from sklearn.neighbors import KNeighborsRegressor
#model selection
from sklearn.model_selection import train_test_split,cross_validate
from sklearn.model_selection import KFold
from sklearn.model_selection import GridSearchCV
#evaluation metrics
from sklearn.metrics import mean_squared_log_error,mean_squared_error, r2_score,mean_absolute_error # for regression
from sklearn.metrics import accuracy_score,precision_score,recall_score,f1_score # for classification
#configure
# sets matplotlib to inline and displays graphs below the corressponding cell.
style.use('fivethirtyeight')
sns.set(style='whitegrid',color_codes=True)
## dividimos o set original de treinamento em um set 75% (treinamento) e 25%
## validação. O seed foi escolhido de modo aleatório (201)
X_train, X_valid, y_train, y_valid = train_test_split (X_orig_train,
y_orig_train,
test_size = 0.25,
random_state = 201)
X_train.head()
X_valid.head()
Scoring rule
O kaggle usa a RMSLE para avaliar a capacidade preditiva dos modelos. Agora apresentamos uma função para calcular essa RMSLE. O método clip trunca valores abaixo de 0.
def RMSLE(predictions, realizations):
predictions_use = predictions.clip(0)
rmsle = np.sqrt(np.mean(np.array(np.log(predictions_use + 1) -
np.log(realizations + 1))**2))
return rmsle
Ajustando uma Regression Tree
rt = DecisionTreeRegressor(min_samples_split=25, random_state=201)
rt_model = rt.fit(X_train, y_train)
rt_pred = rt_model.predict(X_valid)
pd.DataFrame(rt_model.feature_importances_, index=ind_variables_selected)
Que nos dará como output a relevância das características para a previsão:
- workingday = 0.036680
- temp = 0.082574
- year = 0.096813
- month = 0.070877
- hour = 0.660788
- DOW = 0.052269
Para fazer a plotagem da predictions vs. actuals:
plt.figure(figsize=(6,6))
plt.scatter(rt_pred, y_valid, s = 0.2)
plt.xlim(-50,1200)
plt.ylim(-50, 1200)
plt.plot([-50, 1200], [-50, 1200], color='r', linestyle='-', linewidth=3)
plt.suptitle('Dados - Previsão x Observados', fontsize=18)
plt.xlabel('Previsão - Random Tree', fontsize = 16)
plt.ylabel('y_valid', fontsize=16)
plt.show()
RMSLE(rt_pred, y_valid)
RMSLE(rtpred, yvalid) = 0.3966703555172774
Random Forest
rf = RandomForestRegressor(n_estimators=500, max_features=4,
min_samples_leaf=5, random_state=201)
rf_model = rf.fit(X_train, y_train)
rf_pred = rf_model.predict(X_valid)
pd.DataFrame(rf_model.feature_importances_, index=ind_variables_selected)
plt.figure(figsize=(6,6))
plt.scatter(rf_pred, y_valid, s = 0.2)
plt.xlim(-50,1200)
plt.ylim(-50, 1200)
plt.plot([-50, 1200], [-50, 1200], color='r', linestyle='-', linewidth=3)
plt.suptitle('Dados - Previsão x Observados', fontsize=18)
plt.xlabel('Previsão - Random Forest', fontsize = 16)
plt.ylabel('y_valid', fontsize=16)
plt.show()
RMSLE(rf_pred, y_valid)
RMSLE(rfpred, yvalid) = 0.37406876043550474
Boosted Tree
xgb_train = xgb.DMatrix(X_train, label = y_train)
xgb_valid = xgb.DMatrix(X_valid)
num_round_for_cv = 500
param = {'max_depth':6, 'eta':0.1, 'seed':201, 'objective':'reg:linear'}
xgb.cv(param,
xgb_train,
num_round_for_cv,
nfold = 5,
show_stdv = False,
verbose_eval = True,
as_pandas = False)
num_round = 400
xgb_model = xgb.train(param, xgb_train, num_round)
xgb_pred = xgb_model.predict(xgb_valid)
xgb_model.get_fscore()
xgb.plot_importance(xgb_model)
plt.figure(figsize=(6,6))
plt.scatter(xgb_pred, y_valid, s = 0.2)
plt.xlim(-50,1200)
plt.ylim(-50, 1200)
plt.plot([-50, 1200], [-50, 1200], color='r', linestyle='-', linewidth=3)
plt.suptitle('Dados - Previsão x Observados', fontsize=18)
plt.xlabel('Previsão - Boosted Tree', fontsize = 16)
plt.ylabel('y_valid', fontsize=16)
plt.show()
RMSLE(xgb_pred, y_valid)
RMSLE(xgbpred, yvalid) = 0.4704333738143786
Logo, o modelo de previsão com o menor RMSLE foi o Random Forest.
Outros Modelos de Previsão
Nos modelos anteriores, utilizou-se uma relação entre amostra de treino e teste de previsão (80%, 20%). Redimensionando para (75%, 25%) e adicionando novos modelos para prever a variável count, temos:
df.head()
df.columns.to_series().groupby(df.dtypes).groups
x_train,x_test,y_train,y_test=train_test_split(df.drop('count',axis=1),df['count'],test_size=0.25,random_state=42)
models=[RandomForestRegressor(),AdaBoostRegressor(),BaggingRegressor(),SVR(),KNeighborsRegressor()]
model_names=['RandomForestRegressor','AdaBoostRegressor','BaggingRegressor','SVR','KNeighborsRegressor']
rmsle=[]
d={}
for model in range (len(models)):
clf=models[model]
clf.fit(x_train,y_train)
test_pred=clf.predict(x_test)
rmsle.append(np.sqrt(mean_squared_log_error(test_pred,y_test)))
d={'Modelling Algo':model_names,'RMSLE':rmsle}
print(d)
rmsle_frame=pd.DataFrame(d)
rmsle_frame
sns.factorplot(y='Modelling Algo',x='RMSLE',data=rmsle_frame,kind='bar',size=5,aspect=2)
sns.factorplot(x='Modelling Algo',y='RMSLE',data=rmsle_frame,kind='point',size=5,aspect=2)
Que nos fornece o seguinte ranking em termos do RMSLE:
- RandomForestRegressor (0.327785 )
- BaggingRegressor (0.337443)
- KNeighborsRegressor (0.861661)
- AdaBoostRegressor (0.932090)
- SVR (1.434169)
