前言：
共享单车项目中，我们需要根据天气，风速，湿度，温度，时段等信息来预测华盛顿的自行车租赁需求。
训练集提供了11,12两年中每月前19天的每小时自行车租赁数据以及当时的时段，天气，风速，温度等信息。
测试集则提供了11,12两年每月20号至月底的每小时时段，风速，天气，温度等信息，不包括租赁数据（需要预测）。
本次项目说明：
- 编程语言：Python
- 编译工具：jupyter notebook
- 涉及到的库：pandas,numpy,matplotlib,sklearn,seaborn,datetime,dateutil
1.0 数据概览
#设置jupyter可以多行输出
from IPython.core.interactiveshell import InteractiveShell
InteractiveShell.ast_node_interactivity = 'all' #默认为'last'
#让在jupyter绘图可以自由拖拽，实时显示坐标的魔术方法
%matplotlib notebook
# 导入相关库
import numpy as np
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
from dateutil.parser import parse
from datetime import datetime
from sklearn.ensemble import RandomForestRegressor
# 读取数据
train_data = pd.read_csv('train.csv')
test_data = pd.read_csv('test.csv')
train_data.shape
test_data.shape
output:
(10886, 12)
(6493, 9)
train_data[:5]
output:
image
部分字段解释：
- season：季度
- 1：春
- 2：夏
- 3：秋
- 4：东
- holiday：当天是否是节假日
- working：当天是否是工作日
- weather：天气等级
- 1：晴/多云
- 2：阴天
- 3：小雨/小雪
- 4：大雨/大雪
- temp：温度
- atemp：体感温度
- humidity：湿度
- windspeed：风速
- casual：未注册用户租赁数量
- registered：组测用户租赁数量
- count：总的租赁数量
2.0 缺失值/重复值处理
在处理数据之前，我们首先将datetime列进行拆分，方便后续的处理与分析。
#训练集datetime列处理
train_data['date']=train_data['datetime'].map(lambda x:parse(x).date())
train_data['year']=train_data['datetime'].map(lambda x:parse(x).year)
train_data['month']=train_data['datetime'].map(lambda x:parse(x).month)
train_data['hour']=train_data['datetime'].map(lambda x:parse(x).hour)
train_data['weekday']=train_data['datetime'].map(lambda x:parse(x).isoweekday())
#测试集datetime列处理
train_data['date']=train_data['datetime'].map(lambda x:parse(x).date())
train_data['year']=train_data['datetime'].map(lambda x:parse(x).year)
train_data['month']=train_data['datetime'].map(lambda x:parse(x).month)
train_data['hour']=train_data['datetime'].map(lambda x:parse(x).hour)
train_data['weekday']=train_data['datetime'].map(lambda x:parse(x).isoweekday())
#检查训练集集有无缺失数据
train_data.info()
output:
<class 'pandas.core.frame.DataFrame'>
RangeIndex: 10886 entries, 0 to 10885
Data columns (total 17 columns):
datetime 10886 non-null object
season 10886 non-null int64
holiday 10886 non-null int64
workingday 10886 non-null int64
weather 10886 non-null int64
temp 10886 non-null float64
atemp 10886 non-null float64
humidity 10886 non-null int64
windspeed 10886 non-null float64
casual 10886 non-null int64
registered 10886 non-null int64
count 10886 non-null int64
date 10886 non-null object
year 10886 non-null int64
month 10886 non-null int64
hour 10886 non-null int64
weekday 10886 non-null int64
dtypes: float64(3), int64(12), object(2)
memory usage: 1.4+ MB
#检查训练集有无重复数据
train_data.duplicated().value_counts()
output:
False 10886
dtype: int64
可以看到，训练集和测试集都没有缺失值和重复值。
3.0 异常值处理
没有缺失值和重复值不代表数据集就是没有异常，我们先来看下训练集中的数值型数据情况。
train_data[['temp','atemp','humidity','windspeed','casual','registered','count']].describe()
output:
image
从上面这个表格中，我们可以发现count列的数值差异很大，最大值和最小值相差很多。为了进一步观察，我们可以绘制一份散点图矩阵来看下。
sns.pairplot(train_data[['temp','atemp','humidity','windspeed','count']])
image
通过散点图矩阵我们可以很明显的看到一些异常：
- atemp和temp的散点图中，有小部分数据的atemp值相等，但却明显游离在atemp和temp的回归方程之外（通过jupyter的坐标显示，此时的temp-atemp>10）
- atemp部分值缺失，导致散点图出现裂隙
- humidity有小部分数据值为0
- windspeed有大量的数据值为0
- count数据偏斜很严重，导致直方图尾部特别长
对应的解决方案如下：
- atemp游离的异常数据量较少，直接删掉
- atemp部分值缺失导致散点图裂隙，这个暂时还想到怎么处理
- humidity值为0的据量较小，直接删掉
- windspeed值为0的数据量较大，使用随机森林进行回补
- count列可以先剔除掉一些异常值看看
3.1 atemp列处理
len(train_data[(train_data['temp']-train_data['atemp'])>10])
output:
24
atemp列异常数据的量只有24个，数据量不大，可以直接剔除。
train_data = train_data[(train_data['temp']-train_data['atemp'])<=10]
3.2 humidity列处理
len(train_data[train_data['humidity']==0])
humidity列值为0的数据量共22个，也是可以直接剔除。
train_data = train_data[train_data['humidity']!=0]
3.3 count列处理
count列的问题在于数据数据偏斜很严重，count列分布的尾部很长，波动很大，这样容易造成模型过度拟合，所以我们先将一些极端值剔除掉看看，一般我们认为3个标准差以外的数据属于粗大误差，应该被剔除掉。
#剔除极端值
train_data_withoutoutliers = train_data[np.abs(train_data['count'] - train_data['count'].mean()) < (3 * train_data['count'].std())]
去除异常值之后再看下count列的分布
fig = plt.figure()
ax1 = fig.add_subplot(121)
ax2 = fig.add_subplot(122)
sns.distplot(train_data_withoutoutliers['count'],ax=ax1)
ax1.set(title='Distribution of count without outliers')
ax1.set_xticks(range(0,1200,200))
sns.distplot(train_data['count'],ax=ax2)
ax2.set(title='Distribution of count')
ax1.set_xticks(range(0,1200,200))
image
剔除掉异常值之后数据偏斜依旧很严重，所以这里我们需要将count列进行对数转换，减小波动。
train_data = train_data_withoutoutliers
train_data['count_log'] = np.log(train_data['count'])
image
经过对数转换后数据波动更小了，而且大小差异也变小了
3.4 windspeed列处理
这里我们使用随机森林对训练集中windspeed值为0的数据进行回补。
train_data['windspeed_rfr'] = train_data['windspeed']
#将训练集数据分为风速值为0和不为0两部分
windspeed0_data = train_data[train_data['windspeed_rfr']==0]
windspeed1_data = train_data[train_data['windspeed_rfr']!=0]
#设定训练集数据
windspeed_columns = ['season','weather','temp','atemp','humidity','year','month','hour']
X_train = windspeed1_data[windspeed_columns].as_matrix()
y_train = windspeed1_data['windspeed_rfr'].as_matrix()
#生成回归器并进行拟合
windspeed_rfr = RandomForestRegressor(n_estimators=1000,random_state=42)
windspeed_rfr.fit(X_train,y_train)
#预测风速值
X_test = windspeed0_data[windspeed_columns].as_matrix()
windspeed_pred = windspeed_rfr.predict(X_test)
#将预测值填充到风速为0的数据中，并将两个df合并
windspeed0_data['windspeed_rfr'] = windspeed_pred
train_data = pd.concat([windspeed0_data,windspeed1_data]).sort_values('datetime')
3.5 测试集数据处理
对训练集数据进行处理之后，我们同样需要对测试集进行数据清洗，需要注意的一点就是测试集我们不能删除数据，不然提交结果时会出错。
先来看下测试集数据的散点图矩阵
#测试集中没有count列
sns.pairplot(test_data[['temp','atemp','humidity','windspeed']])
image
测试集中atemp列并没有出现明显的异常数据，但是仍有值的缺失导致散点图出现裂隙（这个暂时还不知道怎么处理），humidity列也没有出现值为0的异常数据。
但是windspeed列仍有大量值为0的异常数据，和训练集一样，我们使用随机森林对测试集windspeed值为0的数据进行回补。
test_data['windspeed_rfr'] = test_data['windspeed']
#将测试集数据分为风速值为0和不为0两部分
windspeed0_data = test_data[test_data['windspeed_rfr']==0]
windspeed1_data = test_data[test_data['windspeed_rfr']!=0]
#设定训练集数据
windspeed_columns = ['season','weather','temp','atemp','humidity','year','month','hour']
X_train = windspeed1_data[windspeed_columns].as_matrix()
y_train = windspeed1_data['windspeed_rfr'].as_matrix()
#生成回归器并进行拟合
windspeed_rfr = RandomForestRegressor(n_estimators=1000,random_state=42)
windspeed_rfr.fit(X_train,y_train)
#预测风速值
X_test = windspeed0_data[windspeed_columns].as_matrix()
windspeed_pred = windspeed_rfr.predict(X_test)
#将预测值填充到风速为0的数据中，并将两个df合并
windspeed0_data['windspeed_rfr'] = windspeed_pred
test_data = pd.concat([windspeed0_data,windspeed1_data]).sort_values('datetime')
4.0 数据分析与可视化
在数据分析这部分我们除了可以观察各个变量与租赁数量之间的关系之外，还可以观察各个变量随时间变换的变化趋势。
观察各个变量与租赁数量之间的关系，我们对训练集进行分析即可。观察各个变量随时间变换的变化趋势，这里需要我们将训练集和测试集进行合并然后再进行分析。
#将训练集和测试集合并
bike_data = pd.concat([train_data,test_data],ignore_index=True)
4.1 整体观察
sns.pairplot(data=bike_data,x_vars=['holiday','workingday','weekday','weather','season','hour','month','year','temp','atemp','humidity','windspeed_rfr'],y_vars=['registered','casual','count'],plot_kws={'alpha': 0.1})
image
其实从上面的散点图矩阵我们可以得到很多信息了：
- 非节假日的用车里量多于节假日的用车量，当然这和非节假日的天数要多于节假日的天数有关
- 注册用户在工作日的用车量要多于周末的用车量，而非注册用户则恰好相反
- 天气等级越高，用车量越少
- 注册用户一天的用车高峰时上下班期间，另外中午也有个小高峰
- 非注册用户一天的用车分布基本呈正态分布，用车高峰在13点左右
- 1月和2月时一年中用车量最少的两个月
- 12年和11年相比，不管是注册用户还是非注册用户，用车量都有所增长
- 温度和体感温度对非注册用户影响较大，对注册用户影响较小
- 风速越高，用车量越少
为了更进一步理解数据，同时也是起到学习的作用，下面我们再一个一个的分析各个变量与自行车租赁之间的关系。
首先我们可以简单的将数据集的相关系数矩阵画出来，看下哪些变量与租赁数量关系最密切。
sns.heatmap(train_data.corr(), annot=True, vmax=1, square=True, cmap="Blues",fmt='.2f')
image
可以看到hour与count的相关系数0.41，是与count相关系数最高的一个变量，其次是temp的0.38，那么我们就先来具体看下hour和count之间的关系。
另外temp和atemp相关系数高达0.99，这是不是意味着我们在建模是可以将atemp或者temp排除出特征值。
4.2 时段对自行车租赁的影响
查看时段对自行车租赁影响时，我们可以将数据分为工作日和非工作日两个数据集，看下工作日和非工作日每小时的平均租赁量有什么不同。
fig,ax=plt.subplots(1,2,sharey=True)
workingday_df = train_data[train_data['workingday']==1]
workingday_df.groupby(['hour']).agg({'count':'mean','registered':'mean','casual':'mean'}).plot(ax=ax[0])
nworkingday_df = train_data[train_data['workingday']==0]
nworkingday_df.groupby(['hour']).agg({'count':'mean','registered':'mean','casual':'mean'}).plot(ax=ax[1])
image
在工作日时，注册用户的用车主要集中在早8点和下午17点这两个上下班高峰期，另外中午12点时也是一个小高峰，可能时中午出去吃饭。而非注册用在工作日的用车则是呈现一个正态分布趋势，在下午14点左右用车量最多。
而在非工作日时，注册用户和非注册用户都是在下午13点左右时用车最多，在凌晨4点时用车最少。
4.3 温度对自行车租赁的影响
我们可以先来看下温度随时间的变化趋势,将每天以及每月的平均气温绘制出来。
#按天对温度进行汇总，取一天气温的均数
temp_daily = bike_data.groupby(['date'],as_index=False).agg({'temp':'mean'})
#按月对温度进行汇总，取一月气温的均值
temp_monthly = bike_data.groupby(['year','month'],as_index=False).agg({'weekday':'min','temp':'mean'})
temp_monthly.rename(columns={'weekday':'day'},inplace=True)
temp_monthly['date'] = pd.to_datetime(temp_monthly[['year','month','day']])
fig = plt.figure()
ax = fig.add_subplot(1,1,1)
plt.plot(temp_daily['date'],temp_daily['temp'],label='daily average')
plt.plot(temp_monthly['date'],temp_monthly['temp'],marker='o',label='monthly average')
ax.legend()
image
可以明显的看到，温度呈周期性变化，在每年的7月平均温度最高，而1月则是全年温度最低的月份。
那么温度和单车租赁数量之间是什么关系呢？我们将按温度进行聚合，得到每个温度下的平均租赁数量来进行观察。
count_temp = train_data.groupby(['temp']).agg({'count':'mean','registered':'mean','casual':'mean'})
count_temp.plot()
image
可以发现，在4摄氏度左右时，用车量是最少的，之后随着温度的上升，用车量逐渐增多，到36度左右时，用车量达到了巅峰，再往后走，温度升高，用车量就逐渐减少了。
由于atemp和temp的相关系数高达0.99，所有这里我们就不研究atemp的变化趋势已经它和count之间的关系了，基本都一样。
4.4 湿度对自行车租赁的影响
同样的，我们可以先来看下湿度的全年变化趋势。将每天以及每月的平均湿度绘制出来。
#按天对湿度进行汇总，取一天的湿度均数
humidity_daily = bike_data.groupby(['date'],as_index=False).agg({'humidity':'mean'})
#按月对温度进行汇总，取一月的湿度均值
humidity_monthly = bike_data.groupby(['year','month'],as_index=False).agg({'weekday':'min','humidity':'mean'})
humidity_monthly.rename(columns={'weekday':'day'},inplace=True)
humidity_monthly['date'] = pd.to_datetime(humidity_monthly[['year','month','day']])
fig = plt.figure()
ax = fig.add_subplot(1,1,1)
plt.plot(humidity_daily['date'],humidity_daily['humidity'],label='daily average')
plt.plot(humidity_monthly['date'],humidity_monthly['humidity'],marker='o',label='monthly average')
ax.legend()
image
好像湿度没有像温度一样呈一个周期性变化，但是湿度的波动很大，同一个月内即可能很干燥，也可能特别湿润。
下面来看下湿度对自行车租赁的影响。同样的，我们将按湿度进行聚合，得到每个湿度下的平均租赁数量来进行观察。
count_temp = train_data.groupby(['humidity']).agg({'count':'mean','registered':'mean','casual':'mean'})
count_temp.plot()
image
4.5 风速对自行车租赁的影响
#按天对风速进行汇总，取一天的均数
windspeed_daily = bike_data.groupby(['date'],as_index=False).agg({'windspeed':'mean'})
#按月对风速进行汇总，取一月的均值
windspeed_monthly = bike_data.groupby(['year','month'],as_index=False).agg({'weekday':'min','windspeed':'mean'})
windspeed_monthly.rename(columns={'weekday':'day'},inplace=True)
windspeed_monthly['date'] = pd.to_datetime(humidity_monthly[['year','month','day']])
fig = plt.figure()
ax = fig.add_subplot(1,1,1)
plt.plot(windspeed_daily['date'],windspeed_daily['windspeed'],label='daily average')
plt.plot(humidity_monthly['date'],windspeed_monthly['windspeed'],marker='o',label='monthly average')
ax.legend()
image
风速的波动也是非常大，全年的话并没有一个很明显的变化趋势，似乎1/2/3/4这四个月再全年是风速最高的月份，再往后走，风速都是低于这4月的。
考虑到风速特别大的时候很少，分析风速与用车之间关系时，如果取平均值可能会出现异常，所以按风速对租赁数量取最大值。
count_windspeed = train_data.groupby(['windspeed']).agg({'count':'max','registered':'max','casual':'max'})
count_windspeed.plot()
image
整体来说，随着风速的增加，自行车租赁数量时逐渐下降的，但是在风速45和55左右出现了两次大的反弹，猜测可能是异常数据导致。
4.6 年份/月份/季节对自行车租赁的影响
先来看下11年和12年的会员和非会员的自行车租赁总量情况.
#按年份查看11和12年的自行车租赁总量
train_data.groupby(['year']).agg({'registered':'sum','casual':'sum'}).plot(kind='bar',stacked=True)
image
可以看到，12年相比11年，总的自行车租赁数量和会员租赁数量都有明显的增长，非会员的租赁数量小幅增长。
下面看下11和12两年每月的自行车租赁总量情况。
#按月份查看11-12年每月的自行车租赁总量
count_month = train_data.groupby(['year','month'],as_index=False).agg({'count':'sum','registered':'sum','casual':'sum','weekday':'min'})
count_month.rename(columns={'weekday':'day'},inplace=True)
count_month['date'] = pd.to_datetime(count_month[['year','month','day']])
count_month.set_index('date',inplace=True)
count_month.drop(['year','month','day'],axis=1,inplace=True)
count_month.plot()
image
每年的5/6/7三个月份是用车的高峰期，而12/1/2三个月则时用车的低谷期。
下面我们来看下在不同季节，用户的出行会有什么不同。
#按季节查看各季节的自行车租赁总量
train_data.groupby(['season']).agg({'registered':'sum','casual':'sum'}).plot(kind='bar',stacked=True)
image
夏秋冬三个季节的自行车租赁总量相差不大，都是在一个较高的水平，春季则是一年中用车最少的季节。
4.7 天气对自行车租赁的影响
#按天气等级查看各等级天气的自行车租赁总量
train_data.groupby(['weather']).agg({'registered':'sum','casual':'sum'}).plot(kind='bar',stacked=True)
image
天气等级越高，自行车租赁数量越少，毕竟下雨下雪天，应该不会有太多人会想骑自行车。
4.8 日期对自行车租赁的影响
得到每天的自行车租赁数量
day_df = train_data.groupby(['date'],as_index=False).agg({'count':'sum','registered':'sum','casual':'sum','weekday':'min','workingday':'min','holiday':'min'})
day_df[:5]
output:
image
我们先来看下工作日和非工作日平均每天的自行车租赁数量。
#工作日和非工作日的平均每天自行车租赁数量
day_df.groupby(['workingday']).agg({'registered':'mean','casual':'mean'}).plot(kind='bar',stacked=True)
image
可以看到:
- 非工作日的自行车租赁数量要稍多与工作日
- 注册用户在非工作日的租赁数量要少于工作日
- 未注册用户在非工作日的租赁数量要多于工作日
我们先来看下节假日和非节假日平均每天的自行车租赁数量。
#节假日和非节假日的平均每天自行车租赁数量
day_df.groupby(['holiday']).agg({'registered':'mean','casual':'mean'}).plot(kind='bar',stacked=True)
image
这个图和上面的工作日和非工作日对自行车租赁影响的图基本差不多，所以两者对自行车的影响的结论也是一样，就不赘述了。
最后来看下，一周7天，哪些天用车最多，那些天用车最少？我们按天进行聚合，每天的平均租赁量
day_df.groupby(['weekday']).agg({'registered':'mean','casual':'mean'}).plot(kind='bar',stacked=True)
image
可以发现，在周末是，会员用车下降明显，而非会员用车则明显增长。
5.0 建模与调参
根据观察，我们决定将hour,temp,humidity,year,month,season,weather,windspeed_rfr,weekday,weekday,workingday,holiday作为特征值。
使用随机森林进行回归预测，先将多类别型数据进行分类。
dummies_month = pd.get_dummies(train_data['month'], prefix= 'month')
dummies_season=pd.get_dummies(train_data['season'],prefix='season')
dummies_weather=pd.get_dummies(train_data['weather'],prefix='weather')
dummies_year=pd.get_dummies(train_data['year'],prefix='year')
然后设定训练集，生成模型进行拟合
#把5个新的DF和原来的表连接起来，设定训练集
X_train = pd.concat([train_data,dummies_month,dummies_season,dummies_weather,dummies_year],axis=1).drop(['casual' , 'count' , 'datetime' , 'date' , 'registered',
'windspeed' , 'atemp' , 'month','season','weather', 'year','count_log'],axis=1)
y_train = train_data['count_log']
#生成模型并进行拟合
rfr = RandomForestRegressor(n_estimators=1000)
rfr.fit(X_train,y_train)
对测试集中的变量也进行分类，然后用模型对测试集进行预测
dummies_month = pd.get_dummies(test_data['month'], prefix= 'month')
dummies_season=pd.get_dummies(test_data['season'],prefix='season')
dummies_weather=pd.get_dummies(test_data['weather'],prefix='weather')
dummies_year=pd.get_dummies(test_data['year'],prefix='year')
#把5个新的DF和原来的表连接起来生成测试集
X_test = pd.concat([test_data,dummies_month,dummies_season,dummies_weather,dummies_year],axis=1).drop([ 'datetime' , 'date' ,'windspeed' , 'atemp' , 'month','season','weather', 'year'],axis=1)
#进行预测
y_pred = rfr.predict(X_test)
#生成结果过并保存
submission=pd.DataFrame({'datetime':test_data['datetime'] , 'count':[max(0,x) for x in np.exp(y_pred)]})
submission.to_csv('bike_predictions.csv',index=False)
提交到kaggle后成绩如下：（这里分数越低代表模型性能越好）
image
由于这个项目已经关闭个人排名，所以看不到名称，根据团队排名估计的话，应该是在10%以内。
这里我并没有使用网格搜索进行调参，大家可以试下，另外特征工程中变量的选择，分类，都可以多尝试下。
