【炫技】 用python对人们使用自行车情况分析与预测

这篇博客中，主要用到了pandas的数据清洗和分析工作，同时也用到了sklearn中回归预测的知识，非常的简单，但是产生了较好的预测效果。所有的数据都是可以下载的，重复这些代码也是能够完全重现以上的这些结果的，如果你有疑问，那么可以参考英文原博客[blog1] [blog2]，和原作者的github可以下载完整的代码和数据，

原文：

https://jakevdp.github.io/blog/2014/06/10/is-seattle-really-seeing-an-uptick-in-cycling/

https://jakevdp.github.io/blog/2015/07/23/learning-seattles-work-habits-from-bicycle-counts/

Part 1: 研究问题:

在美国西雅图市，好像人们对自行车越来越喜欢了，从越来越多的自行车俱乐部可以看出端倪。在我们的传统印象中，似乎骑自行车只是作为业余爱好，那么在西雅图是不是也是这种情况呢，自行车的使用情况随着周一到周末会有怎么样具体的变化呢，天气又对人们使用自行车的决定有多大的影响呢，下面我将尝试着回答这些问题。

Part 2:研究工具

本文使用的是python3.4+ipython notebook + pandas + numpy +sklearn，，其实以上的这些只用装一个Anaconda就可以完全解决了，数据是开源的，所有的结果是完全可重现的。

Part 3: 让数据说话

%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

hourly = pd.read_csv(r”E:研究生阶段课程作业python好玩的数据分析SeattleBike-masterFremontHourly.csv”,parse_dates=True,index_col=”Date”)#读入数据，同时默认“Date作为index”
hourly.sample(n = 10)  #随机抽取10行查看情况

随机抽样10行

上面的几行代码，首先是我们读入数据，这些数据可以在作者的github直接下载，https://github.com/jakevdp/SeattleBike

在这里我们读入数据的时候，做了一些小的处理，把csv文件中的”Date”字段当做日期处理。

hourly.columns = [‘northbound’, ‘southbound’]       #把列名改的简单一些，两个列名代表自行车经过时的方向。
hourly[“total”] = hourly.northbound + hourly.southbound    #新加一列，计算每个时刻自行车的总数
daily = hourly.resample(“d”,how=’sum’)  #对数据框的日期按照天进行重采样，属于同一天的加在一起
weekly = daily.resample(‘w’, how = ‘sum’) #对数据框的日期按照天进行重采样，属于同周的加在一起
weekly.plot()

从2013年到2014年西雅图市自行车行驶情况

从可视化结果上看，我们的第一眼直觉告诉我们，在2014年5月左右，西雅图街上的自行车最多，超过了32000/每周,上图，我们可以隐约看出来和2013年相比，2014年的自行车数量有所增加.
另外我们发现，夏天好像人们使用自行车的数量会异常的多，难道是夏天人们更愿意骑自行车吗，还是因为夏天，白昼时间会更长一些导致的呢。

def hours_of_daylight(date, axis=23.44, latitude=47.61):
    “””Compute the hours of daylight for the given date
       date：  exmple:2012-10-07
       axis:   地区的经度
       latitude: 地区的纬度
return :how many hours a day
    “””
    diff = date – pd.datetime(2000, 12, 21)
    day = diff.total_seconds() / 24. / 3600
    day %= 365.25
    m = 1. – np.tan(np.radians(latitude)) * np.tan(np.radians(axis) * np.cos(day * np.pi / 182.625))
    m = max(0, min(m, 2))
    return 24. * np.degrees(np.arccos(1 – m)) / 180.

“hours_of_daylight”只是用来计算某一天的白天有多少个小时的，这和天文学相关了，不用管具体原理了。

daily["day_of_hours"] = daily.index.map(hours_of_daylight) weekly["day_of_hours"] = weekly.index.map(hours_of_daylight) weekly.day_of_hours.plot()   #绘出白天时长随着时间的变化图

白天时长随着时间变化图

果然美帝也是夏天有16个小时白天，冬天只有将近8个小时的白天啊。我们再研究一下2013-2014年，西雅图市每周自行车的数量变化和白天的长度是什么关系。

plt.scatter(weekly.day_of_hours,weekly.total)
plt.xlabel('daylight hours')
plt.ylabel('weekly bicycle traffic');

我们可以很明显的看出，随着白天的时间增加，活跃的自行车数量也在增加，而且感觉上还存在线性关系，所以先前我们说夏天人们骑自行车更多，可能只是因为夏天的白天时间长，并不是因为天气炎热之类的，当然这需要进一步的研究。
既然，白昼时间的自行车总数可能存在线性关系，那么我们来验证一下。
得益于scikit-learn，我们日常需要的机器学习方法在里面全部可以找到。我强烈建议翻看它的文档，我也在按步骤把scikit-learn翻译成中文。

from sklearn.linear_model import LinearRegression X = weekly[['day_of_hours']]   
#白昼时间作为自变量

y = weekly['total']                   
#自行车数量作为应变量
clf = LinearRegression(fit_intercept=True).fit(X, y)  
  #对上面两个变量进行线性拟合print(clf.coef_)    out: 2056.44964989            
  #按照这个线性模型，白昼时间每增加一个小时，那么一周内观察到的自行车数量增加2056辆

print(clf.score(X,y)) out:0.74184335208140473 
   #计算线性模型对数据的拟合情况，结果越接近1，说明线性模型越精确

#为了更加直观看出我们建立的线性模型的预测效果，我们对预测结果和真实结果进行可视化对比

weekly["day_of_hours_trend"] = clf.predict(X)   
#利用模型去预测
"total"weekly[["day_of_hours_trend","total"]].plot()  
    #利用模型预测值和真实值进行对比

预测结果和真实结果可视化对比

我们可以认为我们建立的线性模型较为可靠，可以反应出街上的自行车数量随着白昼时间的变化。

回到我们当初的问题，一般在我们眼里，美帝人民使用自行车，那肯定用来休闲的，或者锻炼身体，那事情情况真是这样吗，我们继续挖掘。

days = ['Mon', 'Tues', 'Wed', 'Thurs', 'Fri', 'Sat', 'Sun'] def get_day_of_week(index):    return days[index.dayofweek] daily["day_of_week"] = daily.index.map(get_day_of_week) daily.head(10) day_count = daily.groupby("day_of_week").total.sum().plot(kind = 'bar')

上面的代码就是新添一个列，对应每一个日期的星期几。，然后对星期1-7进行分租，按组求出自行车出行的总次数。结果如下：

按星期几查看自行车总数情况

这个结果出乎我们的想象，竟然周一到周五人们使用自行车的情况最多，周六和周末最少，看来美帝人民大部分是用自行车来上班的，并不是我们想象的只是骑着自行车来休闲的。到目前为止，影响西雅图人民选择自行车出行的方式已经有白昼的时长、星期几这两个因素了，下面我们再结合天气因素进行深一步的了解。

weather = pd.read_csv(r"E:pythonSeattleBike-masterSeaTacWeather.csv",index_col="DATE",parse_dates = True,usecols=[2,3,6,7])  #读入天气文件

关于西雅图市的天气情况，我们有现成的weather数据。https://github.com/jakevdp/SeattleBikeweather.sample(n = 10)  #随机查看10行数据

西雅图的天气情况

 #对温度和降雨量进行单位转化，转化到我们日常生活中的单位
weather['TMIN'] = 0.18 * weather['TMIN'] + 32
weather['TMAX'] = 0.18 * weather['TMAX'] + 32
weather['PRCP'] /= 254
weather.TMIN.plot()      #绘出每日最高温随着时间的变化
weather.TMAX.plot()     #绘出每日最低温随着时间的变化

气温随着时间变化图

这种图简直就是密集恐惧症患者的灾难，为了更直观显示温度变化，我们对每周重采样，绘出每周的最低温和最高温随着时间的变化。

#按照每周进行时间规划，选出每周温度最大的和温度最小的，可视化
weather.TMIN.resample(‘w’,how=’min’).plot()
weather.TMAX.resample(‘w’,how=’max’).plot()
plt.xlabel(“日期”)
plt.ylabel(“每周的最大气温和最小气温”)

重采样后图

和上面的图相比，我觉得整个世界都清爽了呢。

我们现在拥有了天气、白昼时长、星期几这三个变量，下面的工作就是围绕这三个变量去预测自行车的数目变化，首先我们把天气信息放在一旁，只考虑每天的白昼时长和星期几，以这两个因素为自变量去预测街上自行车数量的变化.

daily_1 = dailyfor i in range(7):     daily_1[days[i]] = (daily_1.index.dayofweek == i).astype(float) daily_1.head(10)

这几行代码比较精简，也不容易理解，我来解释一下，首先是使用daily_1去备份daily这个数据框，然后为数据框增加7个列，分别是[‘Mon’, ‘Tues’, ‘Wed’, ‘Thurs’, ‘Fri’, ‘Sat’, ‘Sun’]，为什么要这么做呢，因为在我们的线性预测模型中，我们使用的自变量必须是数值，但是我们的周几却是字符串变量，因此我们使用了一个小技巧对其进行了量化。我们可以直接看一下量化的结果:

离散变量转化成连续变量

我们观察变化后的结果，可以发现：“day_of_week”转化成了7列，而这7列的值都是数值，因此我们利用 “day_of_hours”和[‘Mon’, ‘Tues’, ‘Wed’, ‘Thurs’, ‘Fri’, ‘Sat’, ‘Sun’]总共八个变量去预测自行车总数的变化。

X = daily_1[days + ['day_of_hours']] y = daily_1['total'] clf = LinearRegression().fit(X, y) daily_1['dayofweek_trend'] = clf.predict(X)print(clf.score(X,y)) 
   #检验线性拟合的效果out： 0.68 (np.abs(daily_1.total - daily_1.dayofweek_trend)).mean()  
   #计算模型预测值和真实值之间的平均误差out: 532.11

从计算的结果可以看出，我们利用两个给出的信息，分别是星期几和白昼时长，我们就可以判断，这个某一天西雅图市街上的自行车总数，平均误差在532辆。

因为我们手里也有西雅图市的天气资料，我们决定加上天气这个选项，联合上面的两个因素，继续做预测。理想的情况是，我们手里的信息越多，预测的结果会越来越准确，那么结果真的是这样吗，我们继续向下探索。

如果你是沿着这篇博客从开始看到这里，你就应该知道，我们现在有两个数据框:
daily和weather，现在我们需要合并这两个数据框，使用降水量、温度、白昼时长、和星期几去预测街上自行车的数目。

daily_new = daily_1.join(weather,how = 'inner')
##对daily 和 weather按照index，进行合并

#使用合并后的数据进行计算：包括周几，每天白昼有几个小时，温度的最大和最小值，
降雨量作为自变量去预测每天的自行车总数

columns = days + ['day_of_hours', 'TMIN', 'TMAX', 'PRCP'] X = daily_new[columns] y = daily_new['total'] clf = LinearRegression().fit(X, y) daily_new['overall_trend'] = clf.predict(X)

(np.abs(daily_new.total - daily_new.overall_trend)).mean()
out:  365

这次我们使用了五个变量去预测某一天街上自行车的数量，和前面只用两个因素得到的532辆平均误差相比，现在的平均误差只有365，我们的预测结果越来越好了哇。

我们来看降雨量这个因素的影响，首先我们需要控制其他变量不变的情况

ind = columns.index('PRCP') slope = clf.coef_[ind]print(slope)out: -824.64

表明在保持其他变量不变的情况下，当降雨量增加一英尺，街上的自行车数量差不多减少824辆.

#因为我们想知道人们在不同日期使用自行车的习惯，我们对数据进行透视
data_pivoted = data.pivot_table(values = [“West”,”East”], index = data.index.date,columns = data.index.hour,fill_value = 0)
data_pivoted.sample(n = 10)

经过透视表操作后的数据

经过处理后的数据列有48维，为了可视化需要，我们对48列进行降维

经过处理后的数据列有48维，为了可视化需要，我们对48列进行降维

from sklearn.decomposition import PCA          #使用自带的pca进行降维
xpca = PCA(n_components = 2).fit_transform(x)    #只保留2维

total_trips = data_pivoted.sum(1)       #计算每一天的经过自行车总数
plt.scatter(xpca[:, 0], xpca[:, 1], c=total_trips,
            cmap=’cubehelix’)
plt.colorbar(label=’total trips’);

很明显，我们可以看出这些数据明显可以分成两类，下面要做的就很简单了，聚类开始登场。首先我们使用kmeans

#首先使用kmeans，聚类的数目为2：
from sklearn.cluster import KMeans
kmeans_model = KMeans(n_clusters=2, random_state=1)
kmeans_model.fit(xpca)
labels = kmeans_model.labels_
#对聚类结果进行可视化
plt.scatter(xpca[:,0],xpca[:,1],c = labels)
plt.colorbar(label=’total trips’);

kmeans

果然kmeans表现的还是那么菜，聚类的效果并不完美，下面我们使用高斯混合模型聚类

#使用高斯混合模型，进行聚类
from sklearn.mixture import GMM
gmm = GMM(2, covariance_type=’full’, random_state=0)
gmm.fit(xpca)
cluster_label = gmm.predict(xpca)
plt.scatter(xpca[:, 0], xpca[:, 1], c=cluster_label);

还是高斯混合模型靠谱啊，多试几种方法，找出最适合你的方法吧

data_pivoted[“Cluster”] = cluster_label
data_new = data.join(data_pivoted[“Cluster”],on = data.index.date)  
data_new.sample(10)

增加了一个所属类别

#对两个聚类按照时间分别作图
#data_new_0        #包含cluster为0的所有数据
#data_new_1        #包含cluser为1的所有数据
data_new_0 = data_new[data_new.Cluster == 0]
data_new_1 = data_new[data_new.Cluster == 1]

by_hour_0 = data_new_0.groupby(data_new_0.index.time).mean()
by_hour_1 = data_new_1.groupby(data_new_1.index.time).mean()

by_hour_0[[“West”,”East”,”total”]].plot()
plt.xlabel(“time”)
plt.ylabel(“mean of bikes”)

第一个类中，自行车数量变化图

这个结果相当耐人寻味，我们可以看到，差不多在早上八点左右，自行车数量到达第一个巅峰，然后下午五点左右，自行车数量到达第二个巅峰。是不是这些日子都是工作日呢，大家骑着自行车去上班，换句话说，第一个聚类难道都是工作日。这个问题有待继续解答,我们再看第二个聚类的情况.

by_hour_1[["West","East","total"]].plot() plt.xlabel("time") plt.ylabel("mean of bikes")

12.png

和第一个聚类相比，这个结果没有出现明显的波峰，大约在下午两年的时候，自行车数量最多，好吧，你肯定在猜测，这些日子是不是周末呢。大家吃过中饭，骑车出来浪呢。为了搞清楚这些疑问，我们计算每个日期对应的星期几。

plt.scatter(xpca[:,0],xpca[:,1],
c = data_pivoted.index.dayofweek,cmap=plt.cm.get_cmap('jet', 7)) cb = plt.colorbar(ticks=range(7)) cb.set_ticklabels(['Mon', 'Tues', 'Wed', 'Thurs', 'Fri', 'Sat', 'Sun'])

我们可以得出这样的结论，周六和周末，人们对自行车的使用有着很大的相似，而周一到周五人们对自行车的使用也很相似，结合前面的聚类结果
 但是我们很奇怪的发现一个现象：有一些工作日的人们表现的和周末很相似，这些特别的日子具体是神马日子的，是不是节假日，另外和其他的工作日相比,周五表现的和周末很暧昧不清，这我们需要思考
 另外在工作日的聚类中，我们发现竟然没有一个非工作日的(至少从图中没有发现特例)，结果真是这样吗，我们需要进一步的使用数据进行分析

days = [‘Mon’, ‘Tues’, ‘Wed’, ‘Thurs’, ‘Fri’, ‘Sat’, ‘Sun’]     #分别对应“dayofweek”:[0,1,2,3,4,5,6]
def get_weekday(index):
    return days[index.dayofweek]
data_new_0[“weekday”] = data_new_0.index.map(get_weekday)
data_new_0_exception = data_new_0[data_new_0.weekday.isin([“Sat”,”Sun”])]
#在第一个聚类中，找特例,换句话说，就是找出这样的周六周末，人们对自行车的使用像工作日一样
len(data_new_0_exception)  

 #结果和我们在上图可视化的结果一样，没有一个周六周末，人们使用自行车像工作日一样

out:0

没有一个周末，人们使用自行车和工作日一样，这也能从侧面看出，看来美帝真心不加班啊，不像天朝，加班累成狗。

data_new_1[‘weekday’] = data_new_1.index.map(get_weekday)
data_new_1_exception = data_new_1[data_new_1.weekday.isin([‘Mon’, ‘Tues’, ‘Wed’, ‘Thurs’, ‘Fri’])]#在第2个聚类中，找特例
len(data_new_1_exception):
out:600

倒是有不少天，人们在工作日的时候和周六周末使用自行车的习惯差不多，我们猜测这些工作日很可能是假期，真的是这样吗，我们来验证一下。

date = set(data_new_1_exception.index.date) #列出从2012-2016年，美国的所有假期 from pandas.tseries.holiday import USFederalHolidayCalendar cal = USFederalHolidayCalendar() holidays = cal.holidays('2012', '2016', return_name=True) holidays_all = pd.concat([holidays,                        
 "Day Before " + holidays.shift(-1, 'D'),                       
   "Day After " + holidays.shift(1, 'D')]) holidays_all = holidays_all.sort_index() holidays_all.ix[(date)]

不出意外，这些表现反常的工作日，全部都在假期中。大家都放假了，当然开始骑车去浪了

最后一个问题，为什么周五的数据，可视化的时候，有几个点表现的特别反常，这几天究竟发生了什么

fri_day = (data_pivoted.index.dayofweek == 4)   #周5为true，其他为false
plt.scatter(xpca[:, 0], xpca[:, 1], c=’gray’, alpha=0.2)
plt.scatter(xpca[fri_day, 0], xpca[fri_day, 1], c=’yellow’);

weird_fridays = pivoted.index[fridays & (Xpca[:, 0] < -600)]weird_fridays
weird_fridays
out: Index([2013-05-17, 2014-05-16, 2015-05-15], dtype=’object’)

果然是，周五的这几个奇异点，果然有情况，我查阅了一下资料，这三天是一年一度的
自行车日。。。。。果然没有无缘无故的爱

总结

关于西雅图市人们使用自行车习惯的数据分析到此就结束了，数据蕴含着很多信息等待我们去挖掘。英文过得去的话，建议直接看原文章，如果有不懂的话，再结合这中文译文进行参照，如果你对python和机器学习，数据挖掘等感兴趣，我们也可以一起学习，一起分享好玩的文章。

来源：数据挖掘入门与实战  公众号： datadw

End

原文始发于微信公众号（PPV课数据科学社区）：【炫技】 用python对人们使用自行车情况分析与预测