通过自行车流量了解西雅图骑行上班族的工作习惯

This notebook originally appeared as a post on the blog Pythonic Perambulations. The content is BSD licensed.

2015年7月25日更新:根据我的同事Ariel Rokem的建议,增加了一些图片。


去年我写过一篇文章,研究西雅图自行车流量的变化趋势,以及与天气、日照时间、日期等因素的关系。

这次,我从另外一个角度出发重新研究这些自行车流量数据。上次是先做出假设,再建立模型描述数据的规律,这次不做任何假设,看看从这些数据本身能够得到什么信息。或者说,上次是使用"有监督机器学习方法"进行数据建模,这次使用“无监督学习方法”进行数据挖掘。

本文将展示使用Python导入数据,转换格式,以及数据可视化和数据分析的技巧。主要使用PandasMatplotlibScikit-learn等扩展包。本文还将展示如何使用无监督机器学习算法挖掘数据信息,如主成分分析(PCA)、高斯混合模型(GMM)等。

当然,我们最关心的问题是,通过分析自行车流量数据,能否了解西雅图骑行上班族的工作习惯。

数据

2012年,西雅图市政府在Fremont桥两侧的车道安装了自行车自动计数器。我们使用的数据就是这些计数器记录的Fremont桥每小时自行车流量。年度、月度数据可以从http://data.seattle.gov/ 下载。这里是每小时流量数据的直接下载链接。 可以将下面这一行命令的注释去掉并执行来下载所需要的数据。

In [1]:
# !curl -o FremontBridge.csv https://data.seattle.gov/api/views/65db-xm6k/rows.csv?accessType=DOWNLOAD

下载之后,使用Pandas库读入数据。

In [2]:
import pandas as pd
data = pd.read_csv('FremontBridge.csv', index_col='Date', parse_dates=True)
data.head()
Out[2]:
Fremont Bridge West Sidewalk Fremont Bridge East Sidewalk
Date
2012-10-03 00:00:00 4 9
2012-10-03 01:00:00 4 6
2012-10-03 02:00:00 1 1
2012-10-03 03:00:00 2 3
2012-10-03 04:00:00 6 1

译注:由于2015年9月的自行车流量发生了显著变化(原因未知),为了保持与该文撰写时所用的数据集一致,增加如下操作。

In [3]:
# 删除2015年7月份以后的数据。
data = data[data.index<'20150701']

首先对数据进行一些预处理。将各栏标题改短,改为“West”和“East”,缺失的数据用0代替,并且增加一栏总数,“Total”。

In [4]:
data.columns = ['West', 'East']
data.fillna(0, inplace=True)
data['Total'] = data.eval('East + West')

对数据集进行简单的可视化,以便更好地了解数据整体的特征。比如,对数据按周重新采样,看看这两年多的变化趋势。

In [5]:
# first some standard imports
%matplotlib inline
import matplotlib.pyplot as plt
import seaborn; seaborn.set()  # plot styling
import numpy as np

data.resample('W', how='sum').plot()
plt.ylabel('weekly trips');

数据表现出很强的季节性,并且有一些局部特征,可能是受温度、日期、降水等因素的影响。

从数据中挖掘信息

我们可以根据直觉猜测哪些因素会影响自行车的流量,再做可种可视化展示,比如我去年那篇文章中讨论了工作日、天气等因素的影响。但是我们也可以用无监督机器学习的方法从数据中发掘有价值的信息。

我们将数据集中的每一天作为一个独立的实体(或者称为样本,这是机器学习的术语)。每天24小时中每个小时都有两个观测数据(东、西车道各一个),每天总共48个观测数据。通过研究每天光照时间,不依赖于任何假设,就能从数据中得到有价值的数量关系。

调整数据格式

第一步先调整数据格式。要得到一个矩阵,每一行代表一天,矩阵的列代表48个观测值。在Pandas中,使用pivot_table()透视表函数调整数据的格式。将原数据表中“东”、“西”两列数据以日期为索引组织,并且按小时分开。所有缺失的数据用0代替。

In [6]:
pivoted = data.pivot_table(['East', 'West'],
                           index=data.index.date,
                           columns=data.index.hour,
                           fill_value=0)
pivoted.head()
Out[6]:
East ... West
0 1 2 3 4 5 6 7 8 9 ... 14 15 16 17 18 19 20 21 22 23
2012-10-03 9 6 1 3 1 10 50 95 146 104 ... 77 72 133 192 122 59 29 25 24 5
2012-10-04 11 0 6 3 1 11 51 89 134 94 ... 63 73 114 154 137 57 27 31 25 11
2012-10-05 7 4 3 2 2 7 37 101 119 81 ... 63 80 120 144 107 42 27 11 10 16
2012-10-06 7 5 2 2 1 2 15 16 47 55 ... 89 115 107 107 41 40 25 18 14 15
2012-10-07 5 5 1 2 2 3 8 12 26 36 ... 126 122 132 118 68 26 19 12 9 5

5 rows × 48 columns

然后将原始数据放入一个矩阵。

In [7]:
X = pivoted.values
X.shape
Out[7]:
(1001, 48)

矩阵中有1000多天的数据,每一天包括48个观测值。

数据可视化

目前的数据可以视为生活在48维空间中的独立个体。每一维数据是桥的某一侧在特定的一小时内的流量值。很难形象地展示48维数据,可以先使用标准的降维技术降低数据的维数。

我们使用的方法称为主成分分析(PCA),这是一种快速的线性投影变换,通过旋转数据,降低维数,并且保持方差最大。我们要求主成分保持90%的方差信息。

In [8]:
from sklearn.decomposition import PCA
Xpca = PCA(0.9).fit_transform(X)
Xpca.shape
Out[8]:
(1001, 2)

输出结果只有两维,意味着,使用这两个主成分就能描述原数据集90%的方差信息。虽然48维数据很难展示,但是二维数据显然很容易绘图。做一个散点图,并且点的颜色由每天的流量总数决定。

In [9]:
total_trips = X.sum(1)
#plt.scatter(Xpca[:, 0], Xpca[:, 1], c=total_trips, cmap='cubehelix')
plt.scatter(Xpca[:, 0], Xpca[:, 1], c=total_trips,
            cmap='cubehelix')
plt.colorbar(label='total trips');

从上图可以看到,所有日子分成两组,每天的总流量数沿着各组的长轴增加。并且每天的流量越小,两组间的界线越不明显。

现在得到一个有趣的结论,可以断定,对于西雅图骑行上班族而言,所有的日子有两种基本类型。下面继续挖掘信息,看能不能找出是哪两类。

无监督聚类

当我们对数据进行分类时,最好不需要事先标记数据的类型,直接自动分类,这种方法称为聚类算法。有很多聚类算法,但是对于上面这种大致呈椭圆分布的数据,使用高斯混合模型(GMM)比较合适。我们仍然使用scikit-learn计算GMM聚类结果,并且标识出每一点所属的类别。

In [10]:
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);

聚类算法将所有样本分为两类。我们把分类信息代回原数据集,看看有什么发现。

In [11]:
pivoted['Cluster'] = cluster_label
data = data.join(pivoted['Cluster'], on=data.index.date)
data.head()
Out[11]:
West East Total Cluster
Date
2012-10-03 00:00:00 4 9 13 0
2012-10-03 01:00:00 4 6 10 0
2012-10-03 02:00:00 1 1 2 0
2012-10-03 03:00:00 2 3 5 0
2012-10-03 04:00:00 6 1 7 0

使用groupby函数得到每一类不同时间的平均趋势。

In [12]:
by_hour = data.groupby(['Cluster', data.index.time]).mean()
by_hour.head()
Out[12]:
West East Total
Cluster
0 00:00:00 5.312139 6.213873 11.526012
01:00:00 2.713873 2.969653 5.683526
02:00:00 2.294798 1.732659 4.027457
03:00:00 1.570809 1.426301 2.997110
04:00:00 4.179191 2.650289 6.829480

最后,把每一类每天各小时的平均流量数据画出来。

In [13]:
fig, ax = plt.subplots(1, 2, figsize=(14, 5))
hourly_ticks = 4 * 60 * 60 * np.arange(6)

for i in range(2):
    by_hour.ix[i].plot(ax=ax[i], xticks=hourly_ticks)
    ax[i].set_title('Cluster {0}'.format(i))
    ax[i].set_ylabel('average hourly trips')

上图直观地显示了两类数据的差异。第一类是较窄的双峰交通流量模式,而第二类是扁平单峰模式。

在双峰这一类中,上午8点左右的流量峰值主要是西侧车道的自行车流量,而下午5点的流量峰值主要是东侧车道的流量。这是典型的工作日模式,早上骑车去西雅图市区上班,下午返回。

而在单峰这一类中,两个方向的流量从早到晚一直比较平稳,峰值出现在下午2点。这是典型的休闲模式,一整天都有人骑车出行。

仅仅使用简单的无监督降维和聚类技术,就能发现数据分为截然不同的两类,并且每一类都有非常直观的解释。

西雅图市民的工作习惯

下面,我们进一步深入研究这些数据,看看能不能发现关于西雅图骑行上班族工作习惯的更多信息。容易想到,第一类很有可能是工作日,而第二类可能是休息日。重新画图,标记出每一个样本点是星期几就能检验这一猜想。

In [14]:
dayofweek = pd.to_datetime(pivoted.index).dayofweek
plt.scatter(Xpca[:, 0], Xpca[:, 1], c=dayofweek,
            cmap=plt.cm.get_cmap('jet', 7))
cb = plt.colorbar(ticks=range(7))
cb.set_ticklabels(['Mon', 'Tues', 'Wed', 'Thurs', 'Fri', 'Sat', 'Sun'])
plt.clim(-0.5, 6.5);

显然,关于工作日和周末的猜测基本正确,但是有一点例外,少数工作日服从周末的车流量模式。另外还有一个有趣的现象,在这张图上周五虽然属于工作日这一类,但是更靠近周末这一类。

让我们仔细看看这些“站错队”的工作日特殊在哪里。先将每一天的类别和星期几列出来。

In [15]:
results = pd.DataFrame({'cluster': cluster_label,
                        'is_weekend': (dayofweek > 4),
                        'weekday': pivoted.index.map(lambda x: x.strftime('%a'))},
                       index=pivoted.index)
results.head()
Out[15]:
cluster is_weekend weekday
2012-10-03 0 False Wed
2012-10-04 0 False Thu
2012-10-05 0 False Fri
2012-10-06 1 True Sat
2012-10-07 1 True Sun

查一下有多少周末被归入第一类,即工作日模式。

In [16]:
weekend_workdays = results.query('cluster == 0 and is_weekend')
len(weekend_workdays)
Out[16]:
0

没有!显然,西雅图骑行上班族作为一个整体在这几年的周末都没有去上班。

类似地,看看有多少工作日被归入第二类,即休息日模式。

In [17]:
midweek_holidays = results.query('cluster == 1 and not is_weekend')
len(midweek_holidays)
Out[17]:
23

过去几年中有23个工作日,西雅图骑行上班族没去上班。为了标记这些日子,先用Pandas读取美国的法定节假日。

In [18]:
from pandas.tseries.holiday import USFederalHolidayCalendar
cal = USFederalHolidayCalendar()
holidays = cal.holidays('2012', '2016', return_name=True)
holidays.head()
Out[18]:
2012-01-02                 New Years Day
2012-01-16    Dr. Martin Luther King Jr.
2012-02-20                Presidents Day
2012-05-28                   MemorialDay
2012-07-04                      July 4th
dtype: object

完整起见,将这些节假日之前和之后一天也标记出来。

In [19]:
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.head()
Out[19]:
2012-01-01                 Day Before New Years Day
2012-01-02                            New Years Day
2012-01-03                  Day After New Years Day
2012-01-15    Day Before Dr. Martin Luther King Jr.
2012-01-16               Dr. Martin Luther King Jr.
dtype: object

注意,这些假期是调休后的日期(observed holiday),比如2012年的元旦标注在1月2日。现在,将西雅图骑行上班族没去上班的所有非周末的日子列出来。

In [20]:
holidays_all.name = 'name'  # required for join
joined = midweek_holidays.join(holidays_all)
set(joined['name'])
Out[20]:
{'Christmas',
 'Day After Christmas',
 'Day After Thanksgiving',
 'Day Before Christmas',
 'July 4th',
 'Labor Day',
 'MemorialDay',
 'New Years Day',
 'Thanksgiving'}

同时也能得到,在下面这些法定节假日,西雅图骑行上班族还是去上班了。

In [21]:
set(holidays) - set(joined.name)
Out[21]:
{'Columbus Day',
 'Dr. Martin Luther King Jr.',
 'Presidents Day',
 'Veterans Day'}

更新:周五怎么了?

我的同事Ariel Rokem看了本文第一版之后,发现一些有趣的现象。在参数平面上,几乎所有的周五都位于工作日这一类的上边缘,更接近周末这一类。但是有三个异常的周五,位于这一类的另一侧。

把周五标出来看清楚一些。

In [22]:
fridays = (dayofweek == 4)
plt.scatter(Xpca[:, 0], Xpca[:, 1], c='gray', alpha=0.2)
plt.scatter(Xpca[fridays, 0], Xpca[fridays, 1], c='yellow');

上图左下方(译注:上图可能会左右颠倒)的黄点非常特殊,它们和其它周五很不一样,甚至和其他工作日也明显不同。找出来看看是哪几天。

In [23]:
weird_fridays = pivoted.index[fridays & (Xpca[:, 0] < -600)]
weird_fridays
Out[23]:
Index([], dtype='object')
In [24]:
# 译注:如果上图中异常的三个黄点在右下方,则使用以下命令
weird_fridays = pivoted.index[fridays & (Xpca[:, 0] > 600)]
weird_fridays
Out[24]:
Index([2013-05-17, 2014-05-16, 2015-05-15], dtype='object')

这三天都是5月中旬。有点意思!

画个图看看这三天的流量与所有日期平均值的对比情况。用透视表做:

In [25]:
all_days = data.pivot_table('Total', index=data.index.time, columns=data.index.date)
all_days.loc[:, weird_fridays].plot();
all_days.mean(1).plot(color='gray', lw=5, alpha=0.3,
                      xticks=hourly_ticks);

显然这三个特殊周五的自行车流量惊人,比平时大得多。为什么会这样?

在网上查了查,这三天是西雅图一年一度的骑车上班日活动。

总结

我们使用基本的可视化方法和非监督机器学习技术对Fremont桥的自行车流量数据进行了研究,得到一些关于经过该桥骑车上下班人群的工作习惯的结论。概括起来,主要有以下几点:

  • 西雅图的骑行上班族整体上在元旦、阵亡将士纪念日、独立日、美国劳动节,以及感恩节和圣诞节前后这些节日不去上班。
  • 但在一些小的节日,他们一般会去单位,如哥伦布日、马丁·路德·金日、总统日,以及美国退伍军人节等。
  • 西雅图骑行上班族整体上从来不在周末去上班。

本文完全用IPython notebook撰写。可以下载Notebook源文件,或者在线阅读

译注:中文版Notebook

Comments

Comments powered by Disqus