Vova's Blog Homepage
Hi, I’m Vladimir Eremin, I am a data scientist during the day and sometimes at night as well. Welcome to my blog, where I share stories hidden behind data. Here, I explore the realms of data science, quantitative finance, and problem-solving. This blog has started as my offline journal and and some point I decided to start publishing in the hope it may be helpful for others as well. Please look around if something catches your eye. Feel free to reach out to me if you have questions or suggestions. Additionally, for those of you in NYC, I co-host a data science meetup - if you’re around, I’ll be happy to meet you at our next event!
Bike Sharing Demand Analysis
by Vladimir
Description: Kaggle playground competition
In [ ]:
import pandas as pd
from pprint import pprint
import numpy as np
import random
from tqdm import tqdm
import gc
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.preprocessing import LabelEncoder, OneHotEncoder
from sklearn.model_selection import train_test_split
#from scipy.sparse import vstack, csr_matrix, save_npz, load_npz
import warnings
import xgboost as xgb
from collections import defaultdict
import time
import seaborn as sns
import statsmodels.api as sm
import scipy
sns.set(style="ticks", color_codes=True)
from XGBoost.XGBoostOptimizer import XGBoostOptimizer
from DataProcessing import DataStats
dtypes = {'datetime': 'O',
'season': 'category',
'holiday': 'category',
'workingday': 'category',
'weather': 'category',
'temp': 'float64',
'atemp': 'float64',
'humidity': 'int64',
'windspeed': 'float64',
'casual': 'int64',
'registered': 'int64',
'count': 'int64'}
plt.style.use('ggplot')
warnings.filterwarnings('ignore')
In [ ]:
traindf = pd.read_csv('/home/vlad/csv/Social/Bike-sharing-demand/train.csv'
,dtype=dtypes,parse_dates=['datetime'])
print('Df size: ',traindf.shape,'\nMem used: ',round(sum(traindf.memory_usage()/1024/1024),2),'Mb')
Df size: (10886, 12) Mem used: 0.71 Mb
In [ ]:
traindf.head()
Out[ ]:
| datetime | season | holiday | workingday | weather | temp | atemp | humidity | windspeed | casual | registered | count | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 2011-01-01 00:00:00 | 1 | 0 | 0 | 1 | 9.84 | 14.395 | 81 | 0.0 | 3 | 13 | 16 |
| 1 | 2011-01-01 01:00:00 | 1 | 0 | 0 | 1 | 9.02 | 13.635 | 80 | 0.0 | 8 | 32 | 40 |
| 2 | 2011-01-01 02:00:00 | 1 | 0 | 0 | 1 | 9.02 | 13.635 | 80 | 0.0 | 5 | 27 | 32 |
| 3 | 2011-01-01 03:00:00 | 1 | 0 | 0 | 1 | 9.84 | 14.395 | 75 | 0.0 | 3 | 10 | 13 |
| 4 | 2011-01-01 04:00:00 | 1 | 0 | 0 | 1 | 9.84 | 14.395 | 75 | 0.0 | 0 | 1 | 1 |
In [ ]:
DataStats.get_df_stats(traindf)
Out[ ]:
| Feature | Unique_values | Percentage of missing values | Percentage of values in the biggest category | dtype | |
|---|---|---|---|---|---|
| 0 | datetime | 10886 | 0.0 | 0.009186 | datetime64[ns] |
| 1 | season | 4 | 0.0 | 25.114826 | category |
| 2 | holiday | 2 | 0.0 | 97.143120 | category |
| 3 | workingday | 2 | 0.0 | 68.087452 | category |
| 4 | weather | 4 | 0.0 | 66.066507 | category |
| 5 | temp | 49 | 0.0 | 4.289914 | float64 |
| 6 | atemp | 60 | 0.0 | 6.163880 | float64 |
| 7 | humidity | 89 | 0.0 | 3.380489 | int64 |
| 8 | windspeed | 28 | 0.0 | 12.061363 | float64 |
| 9 | casual | 309 | 0.0 | 9.057505 | int64 |
| 10 | registered | 731 | 0.0 | 1.791292 | int64 |
| 11 | count | 822 | 0.0 | 1.552453 | int64 |
In [ ]:
traindf.describe()
Out[ ]:
| temp | atemp | humidity | windspeed | casual | registered | count | |
|---|---|---|---|---|---|---|---|
| count | 10886.00000 | 10886.000000 | 10886.000000 | 10886.000000 | 10886.000000 | 10886.000000 | 10886.000000 |
| mean | 20.23086 | 23.655084 | 61.886460 | 12.799395 | 36.021955 | 155.552177 | 191.574132 |
| std | 7.79159 | 8.474601 | 19.245033 | 8.164537 | 49.960477 | 151.039033 | 181.144454 |
| min | 0.82000 | 0.760000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 1.000000 |
| 25% | 13.94000 | 16.665000 | 47.000000 | 7.001500 | 4.000000 | 36.000000 | 42.000000 |
| 50% | 20.50000 | 24.240000 | 62.000000 | 12.998000 | 17.000000 | 118.000000 | 145.000000 |
| 75% | 26.24000 | 31.060000 | 77.000000 | 16.997900 | 49.000000 | 222.000000 | 284.000000 |
| max | 41.00000 | 45.455000 | 100.000000 | 56.996900 | 367.000000 | 886.000000 | 977.000000 |
In [ ]:
#Define a null model - mean count value
null_pred = traindf.describe().loc['mean','count']
print('MAE for Null model: ',round(np.mean(np.abs(traindf['count'] - null_pred)),2))
#traindf['season', 'holiday', 'workingday' 'weather' 'temp' 'atemp' 'humidity' 'windspeed']
MAE for Null model: 142.71
In [ ]:
#Seasonality exploring
g = sns.FacetGrid(traindf, col='season', hue="weather", palette="Set1",size=5)
g = (g.map(sns.distplot, "count", hist=False, rug=True)
.add_legend()
.set_ylabels('density'))
In [ ]:
#Seasonality exploring
g = sns.FacetGrid(traindf, col='season'#, hue="weather"
, palette="Set1",size=5)
g = (g.map(sns.barplot,"weather", "count",palette="Blues_d")
.add_legend()
.set_ylabels('mean count')
)
plt.subplots_adjust(top=0.9)
g.fig.suptitle('Average Count by Season and Waether') # can also get the figure from plt.gcf()
Out[ ]:
Text(0.5,0.98,'Average Count by Season and Waether')
In [ ]:
#Tota count by season
plt.figure(figsize=(5,))
traindf.groupby('season')['count'].sum().plot.bar(ax=plt.gca())
plt.ylabel('total count')
plt.title('Total Count by Season')
#.loc[(traindf.weather=='4') &(traindf.season=='1')]['count'].mean()
Out[ ]:
Text(0.5,1,'Total Count by Season')
In [ ]:
#Holiday exploring
g = sns.FacetGrid(traindf, col='season'#, hue="weather"
, palette="Set1",size=5)
g = (g.map(sns.barplot,"holiday", "count",palette="Blues_d")
.add_legend()
.set_ylabels('mean count')
)
plt.subplots_adjust(top=0.9)
g.fig.suptitle('Average Count by Season and Holiday')
print('\n\n\n')
#Holiday exploring
g = sns.FacetGrid(traindf, col='season'#, hue="weather"
, palette="Set1",size=5)
g = (g.map(sns.barplot,"workingday", "count",palette="Blues_d")
.add_legend()
.set_ylabels('mean count')
)
plt.subplots_adjust(top=0.9)
g.fig.suptitle('Average Count by Workingday and Holiday')
Out[ ]:
Text(0.5,0.98,'Average Count by Workingday and Holiday')
In [ ]:
traindf_numerical = traindf[['humidity','windspeed','temp','count']]
g = sns.pairplot(traindf_numerical.sample(1000)
,y_vars=['count'], x_vars=['humidity','windspeed','temp']
,kind="reg"
,size=5
,palette="Blues_d")
plt.subplots_adjust(top=0.9)
g.fig.suptitle('Continues Variables vs. Count');
In [ ]:
corr = traindf_numerical.corr(); corr
Out[ ]:
| humidity | windspeed | temp | count | |
|---|---|---|---|---|
| humidity | 1.000000 | -0.318607 | -0.064949 | -0.317371 |
| windspeed | -0.318607 | 1.000000 | -0.017852 | 0.101369 |
| temp | -0.064949 | -0.017852 | 1.000000 | 0.394454 |
| count | -0.317371 | 0.101369 | 0.394454 | 1.000000 |
In [ ]:
sns.heatmap(corr,
xticklabels=corr.columns.values,
yticklabels=corr.columns.values)
#plt.subplots_adjust(top=0.9)
plt.title('Corelation Matrix')
Out[ ]:
Text(0.5,1,'Corelation Matrix')
In [ ]:
#Outlier analysis
fig ,axs = plt.subplots(1,2)
fig.set_size_inches(12,5)
sns.boxplot(data=traindf,y="count",orient="v",ax=axs[0],
)
sns.boxplot(data=traindf,y='count',x='workingday',ax=axs[1])
axs[0].set(title='Box Plot on Responce Variable')
axs[1].set(title='Box Plot on RV by Day Type')
Out[ ]:
[Text(0.5,1,'Box Plot on RV by Day Type')]
In [ ]:
fig,axs = plt.subplots(1,2)
fig.set_size_inches(12, 5)
sns.distplot(traindf["count"],ax=axs[0])
sm.qqplot(traindf["count"],ax=axs[1])
axs[0].set(title='Responce Variable Distribution',ylabel='Probability')
axs[1].set(title='QQ Plot on Variable Distribution')
Out[ ]:
[Text(0.5,1,'QQ Plot on Variable Distribution')]
In [ ]:
#Feature Engineering
#create a feature "TimeCat" -{0-3}
import time
from time import mktime
import datetime as dt
def getTimeCat(Datetime):
# extract time categoriesa
#str_time = datetime.strptime(Datetime, "%m/%j/%y %H:%M").time()
str_time = Datetime.time()
# --> Morning = 0400-1000
mornStart = dt.time(4, 0, 1)
mornEnd = dt.time(10, 0, 0)
# --> Midday = 1000-1600
midStart = dt.time(10, 0, 1)
midEnd = dt.time(16, 0, 0)
# --> Evening = 1600-2200
eveStart = dt.time(16, 0, 1)
eveEnd = dt.time(22, 0, 0)
# --> Late Night = 2200-0400
lateStart = dt.time(22, 0, 1)
Midnighy = dt.time(0, 0, 0)
lateEnd = dt.time(4, 0, 0)
if str_time >= mornStart and str_time <= mornEnd:
timecat = 0 #morning
elif str_time >= midStart and str_time <= midEnd:
timecat = 1 #midday
elif str_time >= eveStart and str_time <= eveEnd:
timecat = 2 #evening
elif str_time >= eveStart or (str_time >= Midnighy and str_time <= lateEnd):
timecat = 3 #late night
else:
timecat = -1
return timecat
#traindf['day_part'] =
In [ ]:
traindf['timecat_'] = traindf.datetime.map(getTimeCat)
In [ ]:
traindf.groupby('timecat_')['count'].mean().sort_values().plot.bar()
plt.ylabel('mean count')
plt.title('Count by Time Category')
Out[ ]:
Text(0.5,1,'Count by Time Category')
In [ ]:
traindf['humidity_temp_'] = (1/traindf.humidity) * traindf.temp**2
traindf['temp_'] = traindf['temp'].astype(int)
traindf['atemp_'] = traindf['atemp'].astype(int)
traindf = traindf.replace([np.inf, -np.inf], 1e6)
In [ ]:
traindf.head()
Out[ ]:
| datetime | season | holiday | workingday | weather | temp | atemp | humidity | windspeed | casual | registered | count | timecat_ | humidity_temp_ | temp_ | atemp_ | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 2011-01-01 00:00:00 | 1 | 0 | 0 | 1 | 9.84 | 14.395 | 81 | 0.0 | 3 | 13 | 16 | 3 | 1.195378 | 9 | 14 |
| 1 | 2011-01-01 01:00:00 | 1 | 0 | 0 | 1 | 9.02 | 13.635 | 80 | 0.0 | 8 | 32 | 40 | 3 | 1.017005 | 9 | 13 |
| 2 | 2011-01-01 02:00:00 | 1 | 0 | 0 | 1 | 9.02 | 13.635 | 80 | 0.0 | 5 | 27 | 32 | 3 | 1.017005 | 9 | 13 |
| 3 | 2011-01-01 03:00:00 | 1 | 0 | 0 | 1 | 9.84 | 14.395 | 75 | 0.0 | 3 | 10 | 13 | 3 | 1.291008 | 9 | 14 |
| 4 | 2011-01-01 04:00:00 | 1 | 0 | 0 | 1 | 9.84 | 14.395 | 75 | 0.0 | 0 | 1 | 1 | 3 | 1.291008 | 9 | 14 |
In [ ]:
In [ ]:
In [ ]:
#prepare a data
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
from sklearn.dummy import DummyRegressor
from sklearn.linear_model import SGDRegressor
from sklearn.metrics import mean_absolute_error
import xgboost as xgb
scaler = StandardScaler()
#traindf_float = traindf.copy()
#traindf_float['Date'] = traindf_float['Date'].values.astype(float)
'''
Index(['datetime', 'season', 'holiday', 'workingday', 'weather', 'temp',
'atemp', 'humidity', 'windspeed', 'casual', 'registered', 'count',
'timecat_', 'humidity_temp_', 'temp_', 'atemp_']
'''
columns_use = [
'weather'
# ,'temp'
# ,'humidity'
# ,'windspeed'
,'season'
,'workingday'
,'holiday'
,'timecat_'
,'humidity_temp_'
# ,'temp_'
,'atemp_'
,'count']
#remove outliers
#traindf_no_out = traindf.loc[abs(scipy.stats.zscore(traindf['count']))<2]
feature_matrix = scaler.fit_transform(traindf[columns_use])
# feature_matrix = traindf[columns_use].values
X_ss = feature_matrix[:,:-1]
y_ss = feature_matrix[:,-1]
sig = 3
# X_l = feature_matrix[abs(feature_matrix[:,-1])<sig,:-1]
X_l = np.array(traindf.loc[abs(scipy.stats.zscore(traindf['count']))<sig,columns_use[:-1]].astype(float))
y_l = np.log1p(traindf.loc[abs(scipy.stats.zscore(traindf['count']))<sig,'count'])
print('Original shape:',X_ss.shape,'\nOut. cleand shape: ',X_l.shape)
Original shape: (10886, 7) Out. cleand shape: (10739, 7)
In [ ]:
fig,axs = plt.subplots(3,2)
fig.set_size_inches(15,12)
sns.distplot(y,ax=axs[0][0])
sm.qqplot(y,ax=axs[0][1],line='s')
axs[0][0].set(title='RV Distribution (StandardScaler)',ylabel='Probability')
axs[0][1].set(title='QQ Plot on RV Distribution (StandardScaler)')
sns.distplot(np.log1p(traindf['count']),ax=axs[1][0])
sm.qqplot(np.log1p(traindf['count']),ax=axs[1][1],line='s')
axs[1][0].set(title='RV Distribution (log1p)',ylabel='Probability')
axs[1][1].set(title='QQ Plot on RV Distribution (log1p)')
sns.distplot(y_l,ax=axs[2][0])
sm.qqplot(y_l,ax=axs[2][1],line='s')
axs[2][0].set(title='RV Distribution (log1p+3sig rem)',ylabel='Probability')
axs[2][1].set(title='QQ Plot on RV Distribution (log1p+3sig rem)')
#scipy.stats.normaltest(traindf['count'])
print('normaltest p-values for RV:\nRaw:'
,scipy.stats.normaltest(y).pvalue
,'\nlog1p transformed:',scipy.stats.normaltest(np.log1p(traindf['count'])).pvalue
,'\nlog1p with 3sig removed:',scipy.stats.normaltest(y_l).pvalue)
normaltest p-values for RV: Raw: 0.0 log1p transformed: 9.264576012810486e-221 log1p with 3sig removed: 1.4561149236527596e-225
In [ ]:
X_train, X_test , y_train , y_test = train_test_split(X_ss,y_ss,test_size=0.35,random_state=10)
X_train_l, X_test_l , y_train_l , y_test_l = train_test_split(X_l
,y_l,test_size=0.35,random_state=10)
dtrain = xgb.DMatrix(X_train, label=y_train)
dtest = xgb.DMatrix(X_test, label=y_test)
dtrain_lg = xgb.DMatrix(X_train_l, label=y_train_l)
dtest_lg = xgb.DMatrix(X_test_l, label=y_test_l)
In [ ]:
#define a baseline model
baseline_est = DummyRegressor(strategy='median')
baseline_est.fit(X_train_l,y_train_l)
print("\t---Baseline Score---\n\t R^2 test = "+str(round(baseline_est.score(X_test_l,y_test_l),2)))
---Baseline Score--- R^2 test = -0.08
In [ ]:
reg_optimizer = SGDRegressor(penalty=None)
reg_optimizer_l = SGDRegressor(penalty=None)
reg_optimizer.fit(X_train,y_train)
reg_optimizer_l.fit(X_train_l,y_train_l)
y_pred = reg_optimizer.predict(X_test)
y_pred_l = reg_optimizer_l.predict(X_test_l)
print("\t---GB (no reg) Score---\n\t R^2 test = "+str(round(reg_optimizer.score(X_test,y_test),1))+
'\n\tMAE='+str(round(mean_absolute_error(y_test,y_pred),2))
+'\n\tMAE (l)='+str(round(mean_absolute_error(y_test_l,y_pred_l),2))
)
---GB (no reg) Score--- R^2 test = 0.2 MAE=0.67 MAE (l)=1.3662604833008612e+18
In [ ]:
#defould XGBoost
params = {'max_depth': 6,
'min_child_weight': 1,
'subsample': 1,
'colsample_bytree': 1,
'eta': 0.3,
'objective': 'reg:squarederror',
'eval_metric': 'mae'
}
num_boost_round = 999
early_stopping_round = 10
evals = [(dtest, "Test")]
xgb_optimizer = xgb.train(params=params,
dtrain=dtrain,
num_boost_round=num_boost_round,
evals=evals,
early_stopping_rounds=early_stopping_round
)
[0] Test-mae:0.797428
In [ ]:
xgb.plot_importance(xgb_optimizer)
pprint({i: columns_use[i] for i in range(len(columns_use))})
{0: 'weather',
1: 'season',
2: 'workingday',
3: 'holiday',
4: 'timecat_',
5: 'humidity_temp_',
6: 'atemp_',
7: 'count'}
In [ ]:
#XGB model
xgb_opt = XGBoostOptimizer(dtrain=dtrain_lg,dtest=dtest_lg,cv_metrics = "mae"
,params = params)
xgb_opt.run_optimizer(refit=True)
xgb_opt.plot_optimization()
fin_model = xgb_opt.get_best_model()
xgb.plot_importance(fin_model)
pprint({i: columns_use[i] for i in range(len(columns_use))})
{0: 'weather',
1: 'windspeed',
2: 'season',
3: 'workingday',
4: 'holiday',
5: 'timecat_',
6: 'humidity_temp_',
7: 'atemp_',
8: 'count'}
In [ ]:
#XGB model
xgb_opt = XGBoostOptimizer(dtrain=dtrain,dtest=dtest,cv_metrics = "mae"
,params = params)
xgb_opt.run_optimizer(refit=True)
xgb_opt.plot_optimization()
fin_model = xgb_opt.get_best_model()
xgb.plot_importance(fin_model)
pprint({i: columns_use[i] for i in range(len(columns_use))})
Level selected = complexity Solving best parameters ...
Best stage params: 0.01, -1, mae: 0.4396688
Finished optimization.
Best params: {'max_depth': 9, 'min_child_weight': 6, 'subsample': 1.0, 'colsample_bytree': 0.9, 'eta': 0.01, 'objective': 'reg:squarederror', 'eval_metric': 'mae'}
{0: 'weather',
1: 'windspeed',
2: 'season',
3: 'workingday',
4: 'holiday',
5: 'timecat_',
6: 'humidity_temp_',
7: 'atemp_',
8: 'count'}