这篇文章来自Kaggle上的一位数据科学家，为New York City Taxi Trip Duration（纽约的士路程所花时间预测）项目写的一篇从分析到代码全部囊括在内的解答文章。原文
An MVP Aproach
1. 简介
作为一名（公认犯懒）的数据科学家，我喜欢使用最简单的方式来处理数据。Google是我最好的朋友。除了极少数事之外，几乎没有什么是不能通过一次简单的搜索解决的。
这个kernel使用IPython Notebook编写，使用XGBoost开发，目标是实现让读者看完后就能立刻提交结果（将结果提交后Kaggle会打分并进行排名）。这意味着任何读者都运行代码块、提交结果，并且获得高分。当然，如果你的目标是前10%的话，那需要你付出更多的努力和设计更复杂的算法。在最后一部分我会写一些可以帮助提高分数的建议和一些可优化的参数。如果你能把这些建议实现的话，我相信你提交的结果一定可以在排行榜上获得一个很高的排名。
还请注意，这篇文章是作为解决这一个特定问题而写的一篇引导教程（只适用于此问题，而不是给出一个对于数据分析问题的最优解决方案），还有足够的空间给每位参赛者提出更多创新的解决办法。
简而言之，整个项目是关于学习和分享的，希望我的这些分享能让你在数据科学家的路上走得更远。
1.1 The Problem
纽约有很多单向的街道，并且几乎每条街上每个时间点都有无数的行人。更不要说那些堵塞街道的汽车、摩托车和自行车的数量了。所以，如果你要从A地到B地，不管你多么希望自己准时到，你都会迟到。
在纽约从A地到B地是简单的，Taxi、Uber、Lyft等等都可以。坐出租车的话就不用在意交通状况和行人了，你可以利用这段时间做点别的事，比如收发邮件。虽然想想很美好，但这不意味着你能按时到你的目的地。所以在你选择怎么去的时候，要选择最快的路线。“最快”指的是时间最短，而不是路程最短。
为了弄清楚那条路线最短，我们需要首先预测一段路线走完要多久的时间。因此，这个Playground的比赛是为了预测测试数据中提供了开始和结束位置的每一段路走完所需的时间。
1.2 库和库的作用
Python版本为3.6.1, 请导入下面列出的库（建议使用Anaconda）。请注意%matplotlib inline，这是为了在iPython Notebook允许图标内联。
文档
使用sklearn是为了使用里面已经实现好了的一些机器学习方法。Pandas用于数据存取操作。Numpy是Python中科学计算的基础包。Seaborn是基于matplotlib实现的很友好的数据可视化工具。
import代码如下:
%matplotlib inline
import pandas as pd
from datetime import datetime
import pandas as pd
from sklearn.model_selection import train_test_split
import xgboost as xgb
from sklearn.linear_model import LinearRegression, Ridge,BayesianRidge
from sklearn.cluster import MiniBatchKMeans
from sklearn.metrics import mean_squared_error
from math import radians, cos, sin, asin, sqrt
import seaborn as sns
import matplotlib
import numpy as np
import matplotlib.pyplot as plt
plt.rcParams['figure.figsize'] = [16, 10]
1.3 Loading the Data
使用Pandas的read_csv函数载入数据：
train = pd.read_csv('../input/nyc-taxi-trip-duration/train.csv')
test = pd.read_csv('../input/nyc-taxi-trip-duration/test.csv')
1.4 数据探索
1.4.1 文件结构
Let's start off by exploring the files we just imported, using Pandas' `head` and `describe` functions:
In [3]:
pd.set_option('display.float_format', lambda x: '%.3f' % x)
train.head()
| id | vendor_id | pickup_datetime | dropoff_datetime | passenger_count | …… | |
|---|---|---|---|---|---|---|
| 0 | id2875421 | 2 | 2016-03-14 17:24:55 | 2016-03-14 17:32:30 | 1 | …… |
| 1 | id2377394 | 1 | 2016-06-12 00:43:35 | 2016-06-12 00:54:38 | 1 | …… |
| 2 | id3858529 | 2 | 2016-01-19 11:35:24 | 2016-01-19 12:10:48 | 1 | …… |
| 3 | id3504673 | 2 | 2016-04-06 19:32:31 | 2016-04-06 19:39:40 | 1 | …… |
| 4 | id2181028 | 2 | 2016-03-26 13:30:55 | 2016-03-26 13:38:10 | 1 | …… |
pd.set_option('display.float_format', lambda x: '%.3f' % x)的作用是设置float型数据的显示方式，显示3位小数。想了解更多请参考：这里.
你也可以使用train.tail()显示两个文件的最后5行。从输出中我们能知道数据是有结构的。也可以在Input Files中打开文件查看数据。
1.4.2 数据统计
通过Pandas获得数据一些统计学参数非常容易。使用函数describe即可：
pd.set_option('display.float_format', lambda x: '%.3f' % x)
train.describe()
| vendor_id | passenger_count | pickup_longitude | pickup_latitude | dropoff_longitude | …… | |
|---|---|---|---|---|---|---|
| count | 1458644.000 | 1458644.000 | 1458644.000 | 1458644.000 | 1458644.000 | …… |
| mean | 1.535 | 1.665 | -73.973 | 40.751 | -73.973 | …… |
| std | 0.499 | 1.314 | 0.071 | 0.033 | 0.071 | …… |
| min | 1.000 | 0.000 | -121.933 | 34.360 | -121.933 | …… |
| 25% | 1.000 | 1.000 | -73.992 | 40.737 | -73.991 | …… |
| 50% | 2.000 | 1.000 | -73.982 | 40.754 | -73.980 | …… |
| 75% | 2.000 | 2.000 | -73.967 | 40.768 | -73.963 | …… |
| max | 2.000 | 9.000 | -61.336 | 51.881 | -61.336 | …… |
下面分析一下我们得到了些什么数据。 开始时什么也看不出来，直到看到了 trip_duration。最小值是1秒，最大值是3 526 282秒 (大概980小时)，没有人会坐这么长时间的出租车的，否则账单会爆了的。所以，数据中还有很多的异常值等着我们处理掉。
最后一个查看数据的方法是Pandas的info函数。这个函数的优点是能组合describe和head/tail，同时还有shape（描述文件的行列数）。更多关于info的知识，请参考 这里.
train.info()
<class 'pandas.core.frame.DataFrame'>
RangeIndex: 1458644 entries, 0 to 1458643
Data columns (total 11 columns):
id 1458644 non-null object
vendor_id 1458644 non-null int64
pickup_datetime 1458644 non-null object
dropoff_datetime 1458644 non-null object
passenger_count 1458644 non-null int64
pickup_longitude 1458644 non-null float64
pickup_latitude 1458644 non-null float64
dropoff_longitude 1458644 non-null float64
dropoff_latitude 1458644 non-null float64
store_and_fwd_flag 1458644 non-null object
trip_duration 1458644 non-null int64
dtypes: float64(4), int64(3), object(4)
memory usage: 122.4+ MB
所以现在我们知道了数据大概的样子和我们要关注的变量。接下来要干什么呢？请注意，有趣的地方开始了。
2. 数据预处理
我最初的想法是从pickup_datetime字段中提取出小时、周几和几号这些数据，因为从逻辑上讲不同的时候交通状况就不同，有高峰期（如早高峰和晚高峰），另外，周末交通状况和工作日时的交通状况也不一样。同样，交通状况也可以和月份有关，有一定的季节性。例如，冬天时纽约路面湿滑，车开得就会慢一些，不管你选哪条路，行驶时间都会变长。
一个有趣的变量是passenger_count，乘客人数。乘客人数的增加会造成一些特殊的场景。比如，出租车停车的次数会增加，因为中途有人要下车，这会增加从A地到B地的时间。同时，从物理学的角度考虑，人数太多车太重，走得也会慢点（笑）。当然，这点影响微不足道。
vendor和store_and_fwd_flag一样。我们可以发现，有两个出租车公司，1和2，并且它们中的任何一个都对寻找“最短路线信息”没有帮助，因为要出租车公司的雇员在纽约寻找最佳路线几乎不可能。当然，我也觉得不可能。不过，这个可以也作为一个候选参数（至少也得先看看，最后再排除）。至于store_and_fwd_flag—— 在个别路线中和服务器没有连接这一特征，也可以用来说明几件事。例如，如果我们检查后发现行驶时间和服务器断连有很强的关联性，那么这个特征就可以用来训练模型预测那些个别路线要花费的行驶时间。
接下来要处理的是_lattitude和_longitute（经度和纬度）特征，这是我认为最有价值的特征。可选项有将它们聚类分族，同时，可以发现两个地点的距离和行驶方向。
2.1 行驶时长数据清洗
我们之前提到有一些异常的数据，它们的表现在于trip_duration，行驶时长太短（比如1秒）和太长（比如980小时），这样的数据需要去掉。我的做法是去掉那些偏离平均数两个标准差的数据。当然，探索一下如果是4个标准差会对结果有什么影响是很有价值的，建议大家试试。
m = np.mean(train['trip_duration'])
s = np.std(train['trip_duration'])
train = train[train['trip_duration'] <= m + 2*s]
train = train[train['trip_duration'] >= m - 2*s]
2.2 经纬度数据清洗
纽约城的坐标如下：
city_long_border = (-74.03, -73.75)
city_lat_border = (40.63, 40.85)
和我们用train.describe()输出的数据一比较，就会发现一些地点的坐标超过边界了。所以下面我们需要进行数据清洗，将那些超出边界的数据清洗掉。
train = train[train['pickup_longitude'] <= -73.75]
train = train[train['pickup_longitude'] >= -74.03]
train = train[train['pickup_latitude'] <= 40.85]
train = train[train['pickup_latitude'] >= 40.63]
train = train[train['dropoff_longitude'] <= -73.75]
train = train[train['dropoff_longitude'] >= -74.03]
train = train[train['dropoff_latitude'] <= 40.85]
train = train[train['dropoff_latitude'] >= 40.63]
2.3 日期数据清洗
最后一步数据准备是改变我们的日期数据格式（pickup_datetime和dropoff_datetime）。这在我们进行数据提取时会帮很大的忙。
train['pickup_datetime'] = pd.to_datetime(train.pickup_datetime)
test['pickup_datetime'] = pd.to_datetime(test.pickup_datetime)
train.loc[:, 'pickup_date'] = train['pickup_datetime'].dt.date
test.loc[:, 'pickup_date'] = test['pickup_datetime'].dt.date
train['dropoff_datetime'] = pd.to_datetime(train.dropoff_datetime) #Not in Test
3. 数据可视化和数据分析
接下来的步骤是数据可视化。其实大家可以发现，在日常生活中，一幅准确丰富的图比起一张表格可以给我们提供更直观的信息，让我们更好地了解数据的特征。当然，我们也不要忽视在第1部分时做的工作（数据统计信息，最大最小值、平均数等）。
3.1 初步分析
下面我们来画一个简单的直方图来看一下行驶时间数据，直方图会将行驶时间划分为一百等份。大家可以通过修改划分的等份数来获得对数据的一个大致的感觉。简单地说，这种划分会用到数据的最大最小值，最大最小值的差除以等分数就是每个等份的时间长度。然后将所有的行驶时间长度放入相应的等份内，统计每个等份内的数目作为纵轴。下面是代码：
plt.hist(train['trip_duration'].values, bins=100)
plt.xlabel('trip_duration')
plt.ylabel('number of train records')
plt.show()
这是一个思考如何对数据进行转换的好机会，我们可以去看使用特定的数据转换是否可以使数据出现一个我们熟悉的可以提前假设的模式。比如，log（对数）变换。在这个案例中，对行驶时间使用对数变换是有意义的。如下代码：
train['log_trip_duration'] = np.log(train['trip_duration'].values + 1)
plt.hist(train['log_trip_duration'].values, bins=100)
plt.xlabel('log(trip_duration)')
plt.ylabel('number of train records')
plt.show()
sns.distplot(train["log_trip_duration"], bins =100)
显示了两幅图，第一幅图是将行驶时间做了对数变换后的直方图，而第二幅图上的曲线是密度曲线，从图形可以看出，做了对数变换后，行驶时间已经符合正态分布了。
接下来我们开始研究一下时间。如前文所述，上车时间（包括年月日）可能会揭示一些季节性的规律和一些特定的趋势，但是也有可能会出现一些离群值（如果不是已经清理过的数据集的话），并且会出现一些缺失值（同上，只有不是已经清理过的数据集）。
下面我们对测试数据和训练数据都画一个简单的时序线性图，不止是要看一下数据的趋势和季节性，同时也得看一下两个数据集是否遵循一样的规律。我们可以合理地假设两个数据集会是一样的图形，因为测试集应该是随机从先前的数据集中选取的。通过随机地选取测试数据点，每个数据点都有一样的概率被选为测试数据，从而确保测试样本的一致性。
plt.plot(train.groupby('pickup_date').count()[['id']], 'o-', label='train')
plt.plot(test.groupby('pickup_date').count()[['id']], 'o-', label='test')
plt.title('Trips over Time.')
plt.legend(loc=0)
plt.ylabel('Trips')
plt.show()
很明显，训练数据集图形和测试数据集图形确实是如预料般是同样的形状。
同时，从图形上有几点我们可以一眼看出来。在二月初和六月末的时候行程数量有大规模的下降，但六月的没有二月的那么剧烈。第一个大幅下降也许和这个原因有关：那时是纽约的冬季，所以出行较少（谁会想在结冰的路面开车？）。然而，这好像又不太可能，因为下降的情况好像只在一两天。在我看来更大的可能是罢工。。。（不过咱们这儿是纽约，又不是南非，可能性也不大）。但不管是什么原因啦，这些异常值都值得注意。
下面我们看看两家出租车公司在行驶时间上重要性（或不重要性）吧。
import warnings
warnings.filterwarnings("ignore")
plot_vendor = train.groupby('vendor_id')['trip_duration'].mean()
plt.subplots(1,1,figsize=(17,10))
plt.ylim(ymin=800)
plt.ylim(ymax=840)
sns.barplot(plot_vendor.index,plot_vendor.values)
plt.title('Time per Vendor')
plt.legend(loc=0)
plt.ylabel('Time in Seconds')
所以好像公司1和公司2没有什么区别。也许会有人假设有公司可以选择更快的路，但这更多是交易策略，而不是真正的位置辨识。但这看起来有点不太确定，所以我们可以来先看看另一个特征store_and_fwd_flag，来看它在平均行驶时间上是否有显著的差异。
snwflag = train.groupby('store_and_fwd_flag')['trip_duration'].mean()
plt.subplots(1,1,figsize=(17,10))
plt.ylim(ymin=0)
plt.ylim(ymax=1100)
plt.title('Time per store_and_fwd_flag')
plt.legend(loc=0)
plt.ylabel('Time in Seconds')
sns.barplot(snwflag.index,snwflag.values)
可以看出store_and_fwd_flag特征还是能让行程时间有所区别的。很明显对那些不准确记录行程时间的司机有所倾斜。
前面我们还挖掘出一个变量，就是乘客数量，我假设这个变量会对行程时间有所影响。这个想法是假设乘客越多司机停车的时间也就越多，因此行程时间也就越长（除非不同的下车被记录成了单一行程）。
所以，为了弄清楚这个变量对行程时间的影响，下面我们根据乘客数量对行程时间进行分组，然后计算每组的平均时间。
pc = train.groupby('passenger_count')['trip_duration'].mean()
plt.subplots(1,1,figsize=(17,10))
plt.ylim(ymin=0)
plt.ylim(ymax=1100)
plt.title('Time per store_and_fwd_flag')
plt.legend(loc=0)
plt.ylabel('Time in Seconds')
sns.barplot(pc.index,pc.values)
所以其实乘客数量其实对每个行程没有影响。有趣的是，平均时间是4分多钟时没有乘客（乘客数目为0的那一列，200多秒，就是4分钟左右）。这可能是错误的数据，除非是司机自己给自己钱。
接下来，我们又一次地要看训练数据和测试数据是否对于不同数量的乘客的数据分布是否一致：
In [15]:
train.groupby('passenger_count').size()
Out[15]:
passenger_count
0 52
1 1018715
2 206864
3 58989
4 27957
5 76912
6 47639
dtype: int64
In [16]:
test.groupby('passenger_count').size()
Out[16]:
passenger_count
0 23
1 443447
2 90027
3 25686
4 12017
5 33411
6 20521
9 2
dtype: int64
乍一看好像没有什么东西，但我们要看的不是乘客数目的顺序而是计数。注意在测试集中出现一次坐了9个乘客，这也是我们还没有处理到的异常值。当然，除此之外，两个数据集相对而言是相关的。
下一步是挖掘经纬度数据，大家可以想象一幅地图，上面有各种由经纬度构成的点，点与点之间有路线，大家可以从这种特殊的空间结构感受数据。
3.2 坐标映射
就像我们比较训练集和测试集的行程时间和乘客人数一样，我们也能尝试确认上车地点数据是否在两个数据集中一致。
3.2.1 上车位置
我们使用上面提到过的纽约城的边界位置来创建一幅图，在上面根据经纬度画点。这就可以通过一幅简单的图来看实际的坐标分布。
In [17]:
city_long_border = (-74.03, -73.75)
city_lat_border = (40.63, 40.85)
fig, ax = plt.subplots(ncols=2, sharex=True, sharey=True)
ax[0].scatter(train['pickup_longitude'].values[:100000], train['pickup_latitude'].values[:100000],
color='blue', s=1, label='train', alpha=0.1)
ax[1].scatter(test['pickup_longitude'].values[:100000], test['pickup_latitude'].values[:100000],
color='green', s=1, label='test', alpha=0.1)
fig.suptitle('Train and test area complete overlap.')
ax[0].legend(loc=0)
ax[0].set_ylabel('latitude')
ax[0].set_xlabel('longitude')
ax[1].set_xlabel('longitude')
ax[1].legend(loc=0)
plt.ylim(city_lat_border)
plt.xlim(city_long_border)
plt.show()
我们可以看到两幅图表示的上车位置是相似的，只有一个显著的区别就是训练集的数据更多（这当然是很有意义的）。
3.2.2 距离和方向
下一步相对的有趣。感谢Beluga的文章，我们能通过上车位置和下车位置计算出行程的距离的方向。下面是三个用于计算的函数：
def haversine_array(lat1, lng1, lat2, lng2):
lat1, lng1, lat2, lng2 = map(np.radians, (lat1, lng1, lat2, lng2))
AVG_EARTH_RADIUS = 6371 # in km
lat = lat2 - lat1
lng = lng2 - lng1
d = np.sin(lat * 0.5) ** 2 + np.cos(lat1) * np.cos(lat2) * np.sin(lng * 0.5) ** 2
h = 2 * AVG_EARTH_RADIUS * np.arcsin(np.sqrt(d))
return h
def dummy_manhattan_distance(lat1, lng1, lat2, lng2):
a = haversine_array(lat1, lng1, lat1, lng2)
b = haversine_array(lat1, lng1, lat2, lng1)
return a + b
def bearing_array(lat1, lng1, lat2, lng2):
AVG_EARTH_RADIUS = 6371 # in km
lng_delta_rad = np.radians(lng2 - lng1)
lat1, lng1, lat2, lng2 = map(np.radians, (lat1, lng1, lat2, lng2))
y = np.sin(lng_delta_rad) * np.cos(lat2)
x = np.cos(lat1) * np.sin(lat2) - np.sin(lat1) * np.cos(lat2) * np.cos(lng_delta_rad)
return np.degrees(np.arctan2(y, x))
将这些函数应用于测试集和训练集上，我们就可以通过给定的经纬度计算出在两个点之间的半正矢距离，就是一个大圆球上的两个点之间的距离。我们可以计算出在Manhattan行驶的大致距离，最后我们就可以计算（三角函数）这段距离（位移）的方向。这些计算出的结果将作为单独的变量存储在数据集中。
下一步我们将从数据和显示上创建邻域（位置聚集点），比如Soho或者上东区。
train.loc[:, 'distance_haversine'] = haversine_array(train['pickup_latitude'].values, train['pickup_longitude'].values, train['dropoff_latitude'].values, train['dropoff_longitude'].values)
test.loc[:, 'distance_haversine'] = haversine_array(test['pickup_latitude'].values, test['pickup_longitude'].values, test['dropoff_latitude'].values, test['dropoff_longitude'].values)
train.loc[:, 'distance_dummy_manhattan'] = dummy_manhattan_distance(train['pickup_latitude'].values, train['pickup_longitude'].values, train['dropoff_latitude'].values, train['dropoff_longitude'].values)
test.loc[:, 'distance_dummy_manhattan'] = dummy_manhattan_distance(test['pickup_latitude'].values, test['pickup_longitude'].values, test['dropoff_latitude'].values, test['dropoff_longitude'].values)
train.loc[:, 'direction'] = bearing_array(train['pickup_latitude'].values, train['pickup_longitude'].values, train['dropoff_latitude'].values, train['dropoff_longitude'].values)
test.loc[:, 'direction'] = bearing_array(test['pickup_latitude'].values, test['pickup_longitude'].values, test['dropoff_latitude'].values, test['dropoff_longitude'].values)
3.2.3 创建"Neighborhoods"
也许大家认为需要有地图才能给坐标划分好区域，但这是不必要的。我们可以使用KMeans算法给数据点进行聚类，从而划分好区域。这很简单直接，因为Numpy可以帮我们创建上车和下车地点的数组，然后我们可以使用sklearn的MinBatchKMeans模块，设置一些参数即可完成聚类。
一共分为三步：
- 创建位置堆
- 配置KMeans聚类参数
- 创建实际的聚类
coords = np.vstack((train[['pickup_latitude', 'pickup_longitude']].values,
train[['dropoff_latitude', 'dropoff_longitude']].values))
sample_ind = np.random.permutation(len(coords))[:500000]
kmeans = MiniBatchKMeans(n_clusters=100, batch_size=10000).fit(coords[sample_ind])
train.loc[:, 'pickup_cluster'] = kmeans.predict(train[['pickup_latitude', 'pickup_longitude']])
train.loc[:, 'dropoff_cluster'] = kmeans.predict(train[['dropoff_latitude', 'dropoff_longitude']])
test.loc[:, 'pickup_cluster'] = kmeans.predict(test[['pickup_latitude', 'pickup_longitude']])
test.loc[:, 'dropoff_cluster'] = kmeans.predict(test[['dropoff_latitude', 'dropoff_longitude']])
当你完成了上述步骤之后，就可以显示上车地点的聚类而不是上车位置的点了，如下：
fig, ax = plt.subplots(ncols=1, nrows=1)
ax.scatter(train.pickup_longitude.values[:500000], train.pickup_latitude.values[:500000], s=10, lw=0,
c=train.pickup_cluster[:500000].values, cmap='autumn', alpha=0.2)
ax.set_xlim(city_long_border)
ax.set_ylim(city_lat_border)
ax.set_xlabel('Longitude')
ax.set_ylabel('Latitude')
plt.show()
这是一个KMeans聚类算法很好的可视化表示案例（我们使用了100个聚类，但是可以修改这个值，建议大家试试，然后看看结果有没有改变）。聚类有效地把Manhattan地区划分好了各个区域。
下一步是开始从日期里提取信息，并且要开始考虑给数据编码。
3.3 日期信息提取
提取出日期的不同部分信息的部分原因是为了让我们能够做独热编码中文&英文。这涉及到将分类变量变为二进制变量。此变换会使得训练模型变得很容易，因为机器理解0和1比理解2月和4月要容易得多。为了确保我们能使用我们提取出来的信息来作为特征，我们需要先确定两个数据集中的日期数据是否有相同的大小（比如月份、天数等）。
In [24]:
#Extracting Month
train['Month'] = train['pickup_datetime'].dt.month
test['Month'] = test['pickup_datetime'].dt.month
In [25]:
train.groupby('Month').size(),test.groupby('Month').size()
Out[25]:
(Month
1 226444
2 235054
3 252443
4 247855
5 244591
6 230741
dtype: int64, Month
1 97676
2 102314
3 109697
4 107432
5 107570
6 100445
dtype: int64)
训练集和测试集都是6个月，所以月份可以作为哑变量（编码）。
In [26]:
train['DayofMonth'] = train['pickup_datetime'].dt.day
test['DayofMonth'] = test['pickup_datetime'].dt.day
len(train.groupby('DayofMonth').size()),len(test.groupby('DayofMonth').size())
Out[26]:
(31, 31)
训练集和测试集都有31天，所以天数也可以作为哑变量（编码）。
In [27]:
train['Hour'] = train['pickup_datetime'].dt.hour
test['Hour'] = test['pickup_datetime'].dt.hour
len(train.groupby('Hour').size()),len(test.groupby('Hour').size())
Out[27]:
(24, 24)
训练集和测试集都有24小时，所以小时也可以作为哑变量（编码）。
In [28]:
train['dayofweek'] = train['pickup_datetime'].dt.dayofweek
test['dayofweek'] = test['pickup_datetime'].dt.dayofweek
len(train.groupby('dayofweek').size()),len(test.groupby('dayofweek').size())
Out[28]:
(7, 7)
训练集和测试集每周都有7天，所以星期几也可以编码。
综上所述，我们可以安全地使用日期的不同部分作为模型的特征。下面我们来看一下平均速度如何影响时间，特别是要关注小时、周几和月份如何影响平均速度。当然，平均速度是距离是时间的函数，而测试集中是没有时间的，所以我们训练模型之前要把平均速度删除。
In [29]:
train.loc[:, 'avg_speed_h'] = 1000 * train['distance_haversine'] / train['trip_duration']
train.loc[:, 'avg_speed_m'] = 1000 * train['distance_dummy_manhattan'] / train['trip_duration']
fig, ax = plt.subplots(ncols=3, sharey=True)
ax[0].plot(train.groupby('Hour').mean()['avg_speed_h'], 'bo-', lw=2, alpha=0.7)
ax[1].plot(train.groupby('dayofweek').mean()['avg_speed_h'], 'go-', lw=2, alpha=0.7)
ax[2].plot(train.groupby('Month').mean()['avg_speed_h'], 'ro-', lw=2, alpha=0.7)
ax[0].set_xlabel('Hour of Day')
ax[1].set_xlabel('Day of Week')
ax[2].set_xlabel('Month of Year')
ax[0].set_ylabel('Average Speed')
fig.suptitle('Average Traffic Speed by Date-part')
plt.show()
显然，不同的区域，平均速度是有不同的变化的。城市中心更大程度上更忙（这是因为大部分活动集中在城市中心），郊区平均速度会有所增加。
所以我们可以预计模型在加入了我们的聚类数据后会表现得更好，因为平均速度与区域（即各个类）有关。另外，我们需要使用吴小林直线算法创建数据，因为可以显著提高XGBoost模型的性能。这涉及到图形区域像素化，并且需要记录各条从上车地点到下车地点直线的交叉点。分辨率越高越好，最好能看到交通路口、红绿灯、桥等。使用"has crossed coordinate X"特征可以在模型训练算法时额外创建10000个特征。
注意： XGBoost可以在Mac Book Pro上除以大约10000个特征、1百万行的数据。
4. 数据扩充和哑变量
这个时候其实已经可以训练我们的模型了，但是在数据科学中我们经常还需要分析其它不同的数据库，我们需要看一看有没有其它的数据来源以提高模型的精确度。并且我们还需要创建哑变量来提高模型精确度。
4.1 数据扩充
对于这个问题，我们可以增加OSRM(Open Source Routing Machine)这一特征。这个数据集包含了纽约的从特定点起始的最快路线。
In [31]:
fr1 = pd.read_csv('../input/new-york-city-taxi-with-osrm/fastest_routes_train_part_1.csv', usecols=['id', 'total_distance', 'total_travel_time', 'number_of_steps'])
fr2 = pd.read_csv('../input/new-york-city-taxi-with-osrm/fastest_routes_train_part_2.csv', usecols=['id', 'total_distance', 'total_travel_time', 'number_of_steps'])
test_street_info = pd.read_csv('../input/new-york-city-taxi-with-osrm/fastest_routes_test.csv',
usecols=['id', 'total_distance', 'total_travel_time', 'number_of_steps'])
train_street_info = pd.concat((fr1, fr2))
train = train.merge(train_street_info, how='left', on='id')
test = test.merge(test_street_info, how='left', on='id')
In [32]:
train.shape, test.shape
Out[32]:
((1437128, 29), (625134, 22))
严格的说，我们并不需要看数据的内容，但出于理智我们还是检查了一下训练集和测试集的行数和列数。输出结果测试集比训练集少了5列，即少了5个特征，分别是：dropoff_datetime（下车时间）, avg_speed_m, avg_speed_h, pickup_lat_bin, pickup_long_bin和trip_duration。
下面开始创建哑变量。
4.2 创建哑变量
这一步就是说的我们之前提到的独热编码。这一步依然是使用我们可爱的Pandas。一个简单的函数即可将一个分类的数据转换成哑变量（指示器变量）（参考）。程序很简单，如下：
In [33]:
vendor_train = pd.get_dummies(train['vendor_id'], prefix='vi', prefix_sep='_')
vendor_test = pd.get_dummies(test['vendor_id'], prefix='vi', prefix_sep='_')
passenger_count_train = pd.get_dummies(train['passenger_count'], prefix='pc', prefix_sep='_')
passenger_count_test = pd.get_dummies(test['passenger_count'], prefix='pc', prefix_sep='_')
store_and_fwd_flag_train = pd.get_dummies(train['store_and_fwd_flag'], prefix='sf', prefix_sep='_')
store_and_fwd_flag_test = pd.get_dummies(test['store_and_fwd_flag'], prefix='sf', prefix_sep='_')
cluster_pickup_train = pd.get_dummies(train['pickup_cluster'], prefix='p', prefix_sep='_')
cluster_pickup_test = pd.get_dummies(test['pickup_cluster'], prefix='p', prefix_sep='_')
cluster_dropoff_train = pd.get_dummies(train['dropoff_cluster'], prefix='d', prefix_sep='_')
cluster_dropoff_test = pd.get_dummies(test['dropoff_cluster'], prefix='d', prefix_sep='_')
month_train = pd.get_dummies(train['Month'], prefix='m', prefix_sep='_')
month_test = pd.get_dummies(test['Month'], prefix='m', prefix_sep='_')
dom_train = pd.get_dummies(train['DayofMonth'], prefix='dom', prefix_sep='_')
dom_test = pd.get_dummies(test['DayofMonth'], prefix='dom', prefix_sep='_')
hour_train = pd.get_dummies(train['Hour'], prefix='h', prefix_sep='_')
hour_test = pd.get_dummies(test['Hour'], prefix='h', prefix_sep='_')
dow_train = pd.get_dummies(train['dayofweek'], prefix='dow', prefix_sep='_')
dow_test = pd.get_dummies(test['dayofweek'], prefix='dow', prefix_sep='_')
看到多简单了吗？一个函数就搞定！
我们检查一下输出 (在测试成功前不要相信程序！)。
In [34]:
vendor_train.shape,vendor_test.shape
Out[34]:
((1437128, 2), (625134, 2))
In [35]:
passenger_count_train.shape,passenger_count_test.shape
Out[35]:
((1437128, 7), (625134, 8))
In [36]:
store_and_fwd_flag_train.shape,store_and_fwd_flag_test.shape
Out[36]:
((1437128, 2), (625134, 2))
In [37]:
cluster_pickup_train.shape,cluster_pickup_test.shape
Out[37]:
((1437128, 100), (625134, 100))
In [38]:
cluster_dropoff_train.shape,cluster_dropoff_test.shape
Out[38]:
((1437128, 100), (625134, 100))
In [39]:
month_train.shape,month_test.shape
Out[39]:
((1437128, 6), (625134, 6))
In [40]:
dom_train.shape,dom_test.shape
Out[40]:
((1437128, 31), (625134, 31))
In [41]:
hour_train.shape,hour_test.shape
Out[41]:
((1437128, 24), (625134, 24))
In [42]:
dow_train.shape,dow_test.shape
Out[42]:
((1437128, 7), (625134, 7))
观察一下输出，你可以发现一切都是那么美好，当然，除了乘客人数以外。在测试集中，乘客人数有出现9的情况，所以测试集会比训练集多出一个哑变量特征，就是乘客数量为9。这是一个离群点（异常值），所以我们要删掉它：
In [43]:
passenger_count_test = passenger_count_test.drop('pc_9', axis = 1)
下面我们来到了通常被认为是数据科学最难的地方：将数据清洗操作成可用的格式或结构。
这个在开始模型训练前的最后一步就是去掉所有的分类变量（这些变量在前面那步已经用哑变量代替了），并且准备训练数据和测试数据的最后一个版本。
In [44]:
train = train.drop(['id','vendor_id','passenger_count','store_and_fwd_flag','Month','DayofMonth','Hour','dayofweek',
'pickup_longitude','pickup_latitude','dropoff_longitude','dropoff_latitude'],axis = 1)
Test_id = test['id']
test = test.drop(['id','vendor_id','passenger_count','store_and_fwd_flag','Month','DayofMonth','Hour','dayofweek',
'pickup_longitude','pickup_latitude','dropoff_longitude','dropoff_latitude'], axis = 1)
train = train.drop(['dropoff_datetime','avg_speed_h','avg_speed_m','pickup_lat_bin','pickup_long_bin','trip_duration'], axis = 1)
下面看看去掉分类变量后剩下的变量：
In [45]:
train.shape,test.shape
Out[45]:
((1437128, 11), (625134, 10))
然后开始加入所有的哑变量。
In [46]:
Train_Master = pd.concat([train,
vendor_train,
passenger_count_train,
store_and_fwd_flag_train,
cluster_pickup_train,
cluster_dropoff_train,
month_train,
dom_train,
hour_test,
dow_train
], axis=1)
In [47]:
Test_master = pd.concat([test,
vendor_test,
passenger_count_test,
store_and_fwd_flag_test,
cluster_pickup_test,
cluster_dropoff_test,
month_test,
dom_test,
hour_test,
dow_test], axis=1)
In [48]:
Train_Master.shape,Test_master.shape
Out[48]:
((1437128, 290), (625134, 289))
我们还可以删掉pickup_datetime和pickup_date，因为这两个变量的信息已经包含在其它的变量中了。
In [49]:
Train_Master = Train_Master.drop(['pickup_datetime','pickup_date'],axis = 1)
Test_master = Test_master.drop(['pickup_datetime','pickup_date'],axis = 1)
再检查一次：
In [50]:
Train_Master.shape,Test_master.shape
Out[50]:
((1437128, 288), (625134, 287))
这个就是我们所期待的最终版数据了！训练集比测试集多一个特征，多出来的那个特征就是目标特征，也就是所需时间。下一步就是讲训练集切分成子训练集和子测试集。因为我们需要调整模型参数以提高准确度，所以需要子测试集进行测试（降低RSME）。更重要的是，我们需要一个验证集验证XGBoost。所以XGBoost算法需要三个数据集：训练集、测试集和验证集。
验证集被用于持续地对模型的准确度进行评估，最终再把模型用于对测试集进行预测。所以分隔训练集到一个训练集和一个测试集用于验证，就是为了给我们一个测试样本。很有道理吧？
本教程只使用了10000个数据点，这是为了训练起来更快并且会少一些麻烦。同时，我们将训练集切分为80%训练集，20%测试集用于验证。
所以，你有两件事情可以做了：第一件是改变训练集和验证集的比例（比如三七开），看看结果有什么改变；另一件是提高训练集数据量（比如100K），这应该会提升模型准确度。
代码如下：
In [51]:
Train, Test = train_test_split(Train_Master[0:100000], test_size = 0.2)
接下来还要做几件额外的工作。首先是我们需要去掉log_trip_duration特征（这仅仅是其中一个特征的对数变换）；同时，我们需要去掉索引并重建，以确保各行可以精确地被引用。
In [52]:
X_train = Train.drop(['log_trip_duration'], axis=1)
Y_train = Train["log_trip_duration"]
X_test = Test.drop(['log_trip_duration'], axis=1)
Y_test = Test["log_trip_duration"]
Y_test = Y_test.reset_index().drop('index',axis = 1)
Y_train = Y_train.reset_index().drop('index',axis = 1)
在开始训练前，我们还有一步要做：我们需要创建XGBoost的指标，这会在我们使用XGBoost训练模型时用到。
注意，我们是使用新创建的测试集和训练集用于我们模型的输入（为了训练和验证），并且最终使用模型对测试集进行预测，输出预测结果。
In [53]:
dtrain = xgb.DMatrix(X_train, label=Y_train)
dvalid = xgb.DMatrix(X_test, label=Y_test)
dtest = xgb.DMatrix(Test_master)
watchlist = [(dtrain, 'train'), (dvalid, 'valid')]
我们要开始XGBoost了！！！你可能已经注意到了这一部分实际上很短。但这仅仅只是代码。实际上我们会要通过使用不同的参数以及参数的不同的值来获得性能最好的模型，即输出结果最优，准确度最高。下面我们就进入XGBoost建模吧！
5. XGBoost - 模型训练和准确度测试
如前文所述你可以通过给XGBoost算法设置不同的参数来调整模型的输出。下面是一个虽然很短但是很有效地使用迭代模型参数来调整模型的方法。实现很简单，只需要去掉注释，并且运行这个kernel就是了，可以参考XGBoost的文档来理解每个参数做了什么，又是如何影响训练过程的。
In [54]:
#md = [6]
#lr = [0.1,0.3]
#mcw = [20,25,30]
#for m in md:
# for l in lr:
# for n in mcw:
# t0 = datetime.now()
# xgb_pars = {'min_child_weight': mcw, 'eta': lr, 'colsample_bytree': 0.9,
# 'max_depth': md,
# 'subsample': 0.9, 'lambda': 1., 'nthread': -1, 'booster' : 'gbtree', 'silent': 1,
# 'eval_metric': 'rmse', 'objective': 'reg:linear'}
# model = xgb.train(xgb_pars, dtrain, 50, watchlist, early_stopping_rounds=10,
# maximize=False, verbose_eval=1)
[0] train-rmse:5.70042 valid-rmse:5.69993
[0] train-rmse:5.70042 valid-rmse:5.69993 Multiple eval metrics have been passed: 'valid-rmse' will be used for early stopping.
Will train until valid-rmse hasn't improved in 10 rounds.
[11] train-rmse: 3.25699 valid-rmse: 3.25698
...
...
[89] train-rmse:0.335358 valid-rmse:0.345624
[90] train-rmse:0.334614 valid-rmse:0.344972
[91] train-rmse:0.333921 valid-rmse:0.344405
建议使用下面这些参数：
In [55]:
xgb_pars = {'min_child_weight': 1, 'eta': 0.5, 'colsample_bytree': 0.9,
'max_depth': 6,
'subsample': 0.9, 'lambda': 1., 'nthread': -1, 'booster' : 'gbtree', 'silent': 1,
'eval_metric': 'rmse', 'objective': 'reg:linear'}
model = xgb.train(xgb_pars, dtrain, 10, watchlist, early_stopping_rounds=2,
maximize=False, verbose_eval=1)
print('Modeling RMSLE %.5f' % model.best_score)
[0] train-rmse:3.02558 valid-rmse:3.01507
Multiple eval metrics have been passed: 'valid-rmse' will be used for early stopping.
Will train until valid-rmse hasn't improved in 2 rounds.
[1] train-rmse:1.55977 valid-rmse:1.55301
[2] train-rmse:0.862101 valid-rmse:0.858751
[3] train-rmse:0.56334 valid-rmse:0.564431
[4] train-rmse:0.456502 valid-rmse:0.461474
[5] train-rmse:0.421733 valid-rmse:0.430346
[6] train-rmse:0.410094 valid-rmse:0.421733
[7] train-rmse:0.404835 valid-rmse:0.418769
[8] train-rmse:0.40078 valid-rmse:0.417905
[9] train-rmse:0.398338 valid-rmse:0.417041
Modeling RMSLE 0.41704
当然，如果我没有说明如何如何调整模型的话该怎么做呢？
下面是一些建议：
- 迭代次数超过10次
- 更低的eta-value.
- 提高最大深度
这里一定要注意不要过拟合了，过拟合就是模型训练得很适合训练集，但是对没有用到过的数据，效果却非常不好的一种情况。模型过度地拟合训练数据。建议使用colsample_bytree和subsample控制过拟合。
有兴趣的话，我们可以研究每个特征的重要性来理解行程时间是怎么被影响的。代码：
In [56]:
xgb.plot_importance(model, max_num_features=28, height=0.7)
Out[56]:
<matplotlib.axes._subplots.AxesSubplot at 0x7fe81a457ef0>
从上到下我们可以看到哪个特征对行程时间最有影响。逻辑上是觉得距离是最有影响的。你要去的地方越远，那么花费的时间也就越多。剩下的特征也遵循同样的逻辑，大家可以对着排名思考一下。
在提交结果前的最后一步就是使用我们训练好的模型来生成我们的预测。
In [57]:
pred = model.predict(dtest)
pred = np.exp(pred) - 1
这很简单。我们现在已经成功训练好了模型并且在测试数据上应用了，生成了对测试数据从A点到B点的行程时间预测。因此，激动人心的最后一步就到来了！上传你的结果，然后看看你在排行榜上的排名！开始吧！
提交结果
In [58]:
submission = pd.concat([Test_id, pd.DataFrame(pred)], axis=1)
submission.columns = ['id','trip_duration']
submission['trip_duration'] = submission.apply(lambda x : 1 if (x['trip_duration'] <= 0) else x['trip_duration'], axis = 1)
submission.to_csv("submission.csv", index=False)
这就是我在NYC这个问题上提出的基本模型。我非常希望每件事我都写清楚了，如果有没看懂的地方，欢迎评论并提问。同样如果你发现了一些错误（我也是人嘛。。。），也请给我写信提出，我会立刻改正的。
谢谢阅读，如果觉得有用的话就请点个赞并分享一下吧！
欢迎关注我的微信公众号“金融与人工智能”：
