Primeira vez aqui? Seja bem vindo e cheque o FAQ!
x

Um exercício de Machine Learning usando a base "Bike Sharing Demand" do Kaggle

+2 votos
89 visitas
perguntada Jun 21 em Aprendizagem de Máquinas por Douglas Sad Silveira (121 pontos)  

Os sistemas de compartilhamento de bicicletas são um meio de alugar bicicletas em que o processo de obtenção de associação, locação e devolução de bicicletas é automatizado por meio de uma rede de locais de quiosques em toda a cidade. Usando esses sistemas, as pessoas podem alugar uma bicicleta de um local e devolvê-la em um local diferente conforme necessário. Atualmente, existem mais de 500 programas de compartilhamento de bicicletas em todo o mundo.

Os dados gerados por esses sistemas os tornam atraentes para os pesquisadores porque a duração da viagem, o local da partida, a localização da chegada e o tempo decorrido são explicitamente registrados. Os sistemas de compartilhamento de bicicletas, portanto, funcionam como uma rede de sensores, que pode ser usada para estudar mobilidade em uma cidade.

Nesta competição, os participantes são convidados a combinar os padrões históricos de uso com os dados meteorológicos, a fim de prever a demanda de aluguel de bicicletas no programa Capital Bikeshare em Washington, D.C.

link do exercício no kaggle

Compartilhe

1 Resposta

+2 votos
respondida Jun 21 por Douglas Sad Silveira (121 pontos)  
editado Jun 24 por Douglas Sad Silveira

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:

A imagem será apresentada aqui.

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');

A imagem será apresentada aqui. A imagem será apresentada aqui.
A imagem será apresentada aqui.
A imagem será apresentada aqui.
A imagem será apresentada aqui.
A imagem será apresentada aqui.

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:

A imagem será apresentada aqui.
A imagem será apresentada aqui.
A imagem será apresentada aqui.
A imagem será apresentada aqui.

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)

A imagem será apresentada aqui.

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)

A imagem será apresentada aqui.

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 imagem será apresentada aqui.

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)

A imagem será apresentada aqui.
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)

A imagem será apresentada aqui.

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)

A imagem será apresentada aqui.

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:

  1. RandomForestRegressor (0.327785 )
  2. BaggingRegressor (0.337443)
  3. KNeighborsRegressor (0.861661)
  4. AdaBoostRegressor (0.932090)
  5. SVR (1.434169)

A imagem será apresentada aqui.

comentou Jun 24 por danielcajueiro (5,726 pontos)  
Douglas, ficou super legal a sua resposta. Mas para mim nao ficou claro o processo que vc fez na seção que chama de Feature Engineering.
comentou Jun 24 por danielcajueiro (5,726 pontos)  
Outra curiosidade. Alguns dos modelos que vc usou tem varios parametros. Como vc escolheu os parametros?
comentou Jun 24 por Douglas Sad Silveira (121 pontos)  
editado Jun 25 por Douglas Sad Silveira
Daniel, obrigado pelos comentários. O meu código ficou enorme e a resposta pode ter sido comprometida pela minha pouca experiência em ML. Tentarei ser mais didático:

Feature Engeeniring

A ideia da seção foi basicamente considerar novas características a partir daquelas já existentes no próprio DataFrame, conforme apresentado na parte descritiva da resposta.

Inicialmente, extraiu-se de datetime o ano (2011 e 2012). Esse mesmo procedimento foi feito para a característica "month",  "hour" e "day of the week". Note que, neste caso,
o código se refere a uma coluna por inteira. Para acessar o valor de uma célula em específico, sim, foi utilizado o método iloc, que nos permite acessar uma combinação específica de linha e coluna.

Sobre os Parâmetros das árvores de decisão:

eu segui o default do min_samples_split (20), que vi nos códigos do Kaggle. Em palavras, isso significa que 20 samples eram necessárias em cada nó para que esse nó fosse então subdivido. Assim, quanto menor for seu set up pro min_sample_split, mais "profundo" é o fit da árvore de decisão. Para os demais parâmetros, me orientei pelos exemplos sugeridos no capítulo 8 do seguinte livro (gratuito), que também estava indicado no fórum do Kaggle:  

http://www-bcf.usc.edu/~gareth/ISL/

Mesmo com os exemplos computacionais executados no R, não foi difícil trazer a intuição para o ambiente do python.

Na boosted tree (xgb model), pra fazer a cross validation, a escolha para
num_round_for_cv = 500

foi baseada no critério do "número de rounds (num_rounds) com o menor
test-rmse". O set dos outros parâmetros foi guiado pela seguinte referência:

https://xgboost.readthedocs.io/en/latest/

param = {'max_depth':6, 'eta' = 0.1, 'seed': 201, 'objective': 'reg:linear'}

desta forma, a profundidade de cada árvore é 6 ('max_depth"), o parâmetro 'eta' corresponde a taxa de aprendizado (1-'eta'), que neste caso corresponde a um aprendizado lento. Como o conjunto de características são continuas, faz-se uma regressão linear, cujo seed é dado por 201 (mas poderia ser outro valor).

Daí então, foi possível encontrar o melhor número de árvores a se utilizar no modelo ajustado (que nos retornou um num-round=400).

Espero ter clarificado um pouco mais a resposta!
comentou Jun 24 por danielcajueiro (5,726 pontos)  
By the way, gostei das figuras.
comentou Jun 24 por danielcajueiro (5,726 pontos)  
Essa base nao tem problema de missing data?
comentou Jun 24 por Douglas Sad Silveira (121 pontos)  
Não há problema de missing data! Eu testei isso, também. Não coube na resposta, pois o código ficou muito grande. Você acha que caberia postar um link com ele no github, por exemplo?
comentou Jun 24 por Douglas Sad Silveira (121 pontos)  
As figuras deram um certo trabalho, mas foi interessante ir explorando as bibliotecas, especialmente a seaborn, que não conhecia.
comentou Jun 25 por danielcajueiro (5,726 pontos)  
Se vc quiser colocar um link para o github nao tem problema. Mas vc pode dividir tb a sua resposta em varias sub-respostas.
comentou Jun 25 por Douglas Sad Silveira (121 pontos)  
Boa ideia, Professor! Não tinha pensado nisso...
...