Python業(yè)務(wù)分析實(shí)戰(zhàn)｜共享單車數(shù)據(jù)挖掘

數(shù)據(jù)集簡介

關(guān)于共享單車數(shù)據(jù)集

自行車共享系統(tǒng)是傳統(tǒng)自行車租賃的新一代，從注冊會員、租賃到歸還的整個過程都是自動化的。通過這些系統(tǒng)，用戶可以很容易地從一個特定的位置租用自行車，并在另一個位置歸還。目前，全球大約有500多個共享單車項(xiàng)目，這些項(xiàng)目由50多萬輛自行車組成。今天，由于它們在交通、環(huán)境和健康問題上的重要作用，人們對這些系統(tǒng)產(chǎn)生了極大的興趣。

除了自行車共享系統(tǒng)在現(xiàn)實(shí)世界的有趣應(yīng)用之外，眾多研究者們對這些系統(tǒng)所產(chǎn)生的數(shù)據(jù)產(chǎn)生濃厚的興趣。與其他運(yùn)輸服務(wù)（如公共汽車或地鐵）不同，共享自行車使用的持續(xù)時間、出發(fā)時間和到達(dá)位置都明確地記錄在系統(tǒng)中。這一功能將自行車共享系統(tǒng)變成了一個虛擬傳感器網(wǎng)絡(luò)，可用于感知城市中的流動性。因此，通過監(jiān)測這些數(shù)據(jù)，預(yù)計(jì)可以檢測到城市中的大多數(shù)重要事件。

今天我們就運(yùn)用這些數(shù)據(jù)集，挖掘出蘊(yùn)含在其中的有效的信息。接下來從探索數(shù)據(jù)屬性，清洗數(shù)據(jù)，到模型開發(fā)，一起來學(xué)習(xí)，共同進(jìn)步。

注意，該數(shù)據(jù)集是國外共享單車數(shù)據(jù)集，并非國內(nèi)的共享單車數(shù)據(jù)集。但不影響我們學(xué)習(xí)數(shù)據(jù)挖掘相關(guān)知識和技術(shù)。數(shù)據(jù)集獲取可以聯(lián)系原文作者云朵君（Mr_cloud_data）獲取。

屬性信息

hour.csv和?day.csv都有以下字段，day.csv中沒有?hr?字段

• instant：記錄索引
• dteday?：日期
• season?：季節(jié) （1：春天, 2：夏天, 3：秋天, 4：冬天）
• yr?：年份 （0：2011, 1：2012）
• mnth：月份 （ 1 to 12）
• hr：小時 （0 to 23）
• holiday?：是否是假期
• weekday?：星期幾
• workingday?：工作日，如果日既不是周末也不是假日，則為1，否則為0。
• weathersit：
• 1：晴，少云，部分云，無云
• 2：薄霧+多云，薄霧+碎云，薄霧+少量云，薄霧
• 3：小雪，小雨+雷暴+散云，小雨+散云
• 4：大雨+冰板+雷暴+霧，雪+霧
• temp：標(biāo)準(zhǔn)化溫度數(shù)據(jù)，單位為攝氏度。這些值是通過(t-t_min)/(t_max-t_min, t_min=-8, t_max=+39（僅在小時范圍內(nèi)）得到的
• atemp：以攝氏度為單位的正常體感溫度。這些值是通過(t-t_min)/(t_max-t_min), t_min=-16, t_max=+50（僅在小時范圍內(nèi)）得到的
• hum：標(biāo)準(zhǔn)化濕度。這些值被分割到100（最大值）
• windspeed：歸一化的風(fēng)速數(shù)據(jù)。這些值被分割到 67 （最大值）
• casual：注銷用戶數(shù)量
• registered：已注冊用戶數(shù)量
• cnt：出租自行車總數(shù)，包括注銷和注冊自行車

前期準(zhǔn)備

導(dǎo)入模塊

import?seaborn?as?sns
import?matplotlib.pyplot?as?plt
from?prettytable?import?PrettyTable
import?numpy?as?np  
import?pandas?as?pd  
from?sklearn.model_selection?import?RandomizedSearchCV
from?sklearn.metrics?import?mean_squared_error,?mean_absolute_error,?mean_squared_log_error
from?sklearn.linear_model?import?Lasso,?ElasticNet,?Ridge,?SGDRegressor
from?sklearn.svm?import?SVR,?NuSVR
from?sklearn.ensemble?import?BaggingRegressor,?RandomForestRegressor
from?sklearn.neighbors?import?KNeighborsClassifier
from?sklearn.cluster?import?KMeans
from?sklearn.ensemble?import?RandomForestClassifier
from?sklearn.ensemble?import?GradientBoostingClassifier
from?sklearn.linear_model?import?LinearRegression
import?random
%matplotlib?inline

random.seed(100)

定義數(shù)據(jù)獲取函數(shù)

class?Dataloader():
????'''
????自行車共享數(shù)據(jù)集數(shù)據(jù)加載器。
????'''
????def?__init__(self,?csv_path):
????????'''?
????????初始化自行車共享數(shù)據(jù)集數(shù)據(jù)加載器。
????????param：?csv_path?{str}?
????????--?自行車共享數(shù)據(jù)集CSV文件的路徑。
????????'''
????????self.csv_path?=?csv_path
????????self.data?=?pd.read_csv(self.csv_path)
????????#?Shuffle
????????self.data.sample(frac=1.0,?replace=True,?random_state=1)

????def?getHeader(self):
????????'''?
????????獲取共享單車CSV文件的列名。
????????return：?[list?of?str]
????????--CSV文件的列名
????????'''
????????return?list(self.data.columns.values)
????????
????def?getData(self):
????????'''???
????????劃分訓(xùn)練、驗(yàn)證和測試集
????????返回：?pandas?DataFrames
????????--?劃分后的不同數(shù)據(jù)集pandas?DataFrames
????????'''
????????#?將數(shù)據(jù)按60：20：20的比例劃分為訓(xùn)練、驗(yàn)證和測試集
??????
????????split_train?=?int(60?/?100?*?len(self.data))?
????????split_val?=?int(80?/?100?*?len(self.data))?
????????train?=?self.data[：split_train]
????????val?=?self.data[split_train：split_val]
????????test?=?self.data[split_val：]
????????return?train,?val,??test
????????
????def?getFullData(self):
????????'''?
????????在一個DataFrames中獲取所有數(shù)據(jù)。
????????return：?pandas?DataFrames
????????--?完整的共享數(shù)據(jù)集數(shù)據(jù)
????????'''
????????return?self.data

描述性分析

劃分訓(xùn)練、驗(yàn)證和測試數(shù)據(jù)集

dataloader?=?Dataloader('../data/bike/hour.csv')
train,?val,?test?=?dataloader.getData()
fullData?=?dataloader.getFullData()

category_features?=?['season',?'holiday',?'mnth',?'hr',?  
???????????????????? 'weekday',?'workingday',?'weathersit']
number_features?=?['temp',?'atemp',?'hum',?'windspeed']

features=?category_features?+?number_features
target?=?['cnt']

features

['season',
 'holiday',
 'mnth',
 'hr',
 'weekday',
 'workingday',
 'weathersit',
 'temp',
 'atemp',
 'hum',
 'windspeed']

獲取DataFrame的列名：

print(list(fullData.columns))

['instant', 'dteday', 'season', 'yr', 'mnth',
'hr', 'holiday', 'weekday', 'workingday', 
'weathersit', 'temp', 'atemp', 'hum', 
'windspeed', 'casual', 'registered', 'cnt']

打印數(shù)據(jù)集的前五個示例來探索數(shù)據(jù)：

fullData.head(5)

獲取每列的數(shù)據(jù)統(tǒng)計(jì)信息：

fullData[number_features].describe()

for?col?in?category_features：
????fullData[col]?=?fullData[col].astype('category')
fullData[category_features].describe()

缺失值分析

缺失值分析可參見往期文章：缺失值處理，你真的會了嗎？

檢查數(shù)據(jù)中的NULL值：

print(fullData.isnull().any())

instant       False
dteday        False
season        False
yr            False
mnth          False
hr            False
holiday       False
weekday       False
workingday    False
weathersit    False
temp          False
atemp         False
hum           False
windspeed     False
casual        False
registered    False
cnt           False
dtype：bool

異常值分析

箱形圖

sns.set(font_scale=1.0)
fig,?axes?=?plt.subplots(nrows=3,ncols=2)
fig.set_size_inches(15,?15)
sns.boxplot(data=train,y="cnt",orient="v",ax=axes[0][0])
sns.boxplot(data=train,y="cnt",x="mnth",orient="v",ax=axes[0][1])
sns.boxplot(data=train,y="cnt",x="weathersit",orient="v",ax=axes[1][0])
sns.boxplot(data=train,y="cnt",x="workingday",orient="v",ax=axes[1][1])
sns.boxplot(data=train,y="cnt",x="hr",orient="v",ax=axes[2][0])
sns.boxplot(data=train,y="cnt",x="temp",orient="v",ax=axes[2][1])

axes[0][0].set(ylabel='Count',title="Box?Plot?On?Count")
axes[0][1].set(xlabel='Month',?ylabel='Count',title="Box?Plot?On?Count?Across?Months")
axes[1][0].set(xlabel='Weather?Situation',?ylabel='Count',title="Box?Plot?On?Count?Across?Weather?Situations")
axes[1][1].set(xlabel='Working?Day',?ylabel='Count',title="Box?Plot?On?Count?Across?Working?Day")
axes[2][0].set(xlabel='Hour?Of?The?Day',?ylabel='Count',title="Box?Plot?On?Count?Across?Hour?Of?The?Day")
axes[2][1].set(xlabel='Temperature',?ylabel='Count',title="Box?Plot?On?Count?Across?Temperature")

for?tick?in?axes[2][1].get_xticklabels()：
????tick.set_rotation(90)

解析：?工作日和節(jié)假日箱形圖表明，正常工作日出租的自行車比周末或節(jié)假日多。每小時的箱形圖顯示當(dāng)?shù)卦缟?點(diǎn)最大，下午5點(diǎn)最大，這表明大多數(shù)自行車租賃服務(wù)的用戶使用自行車上班或上學(xué)。另一個重要因素似乎是溫度：較高的溫度導(dǎo)致自行車租賃數(shù)量增加，而較低的溫度不僅降低了平均租賃數(shù)量，而且在數(shù)據(jù)中顯示出更多的異常值。

從數(shù)據(jù)中去除異常值

sns.distplot(train[target[-1]]);

計(jì)數(shù)值的分布圖顯示，計(jì)數(shù)值不符合正態(tài)分布。我們將使用中位數(shù)和四分位區(qū)間（IQR）來識別和去除數(shù)據(jù)中的異常值。（另一種方法是將目標(biāo)值轉(zhuǎn)換為正態(tài)分布，并使用平均值和標(biāo)準(zhǔn)偏差。）

print("帶有異常值的列車集合中的樣本：?{}".format(len(train)))
q1?=?train.cnt.quantile(0.25)
q3?=?train.cnt.quantile(0.75)
iqr?=?q3?-?q1
lower_bound?=?q1?-(1.5?*?iqr)?
upper_bound?=?q3?+(1.5?*?iqr)?
train_preprocessed?=?train.loc[(train.cnt?>=?lower_bound)?&?(train.cnt?<=?upper_bound)]
print("沒有異常值的訓(xùn)練樣本集：?{}".format(len(train_preprocessed)))
sns.distplot(train_preprocessed.cnt);

帶有異常值的列車集合中的樣本：10427
沒有異常值的訓(xùn)練樣本集：10151

相關(guān)分析

matrix?=?train[number_features?+?target].corr()
heat?=?np.array(matrix)
heat[np.tril_indices_from(heat)]?=?False
fig,ax=?plt.subplots()
fig.set_size_inches(15,8)
sns.set(font_scale=1.0)
sns.heatmap(matrix,?mask=heat,vmax=1.0,?vmin=0.0,?square=True,annot=True,?cmap="Reds")

結(jié)論：?在描述性分析總結(jié)如下幾點(diǎn)：

• 變量"Casual"和"registered"包含關(guān)于共享自行車計(jì)數(shù)直接信息，而如果將這些信息用于預(yù)測（數(shù)據(jù)泄漏）。因此，它們不在特征集中考慮。
• 變量"temp"和"atemp"是高度相關(guān)的。為了降低預(yù)測模型的維數(shù)，可以刪除特征"atemp"。
• 變量"hr"和"temp"似乎是預(yù)測自行車共享數(shù)量的貢獻(xiàn)較大的特征。

features.remove('atemp')

模型選擇

所呈現(xiàn)的問題的特點(diǎn)是：

• 回歸：目標(biāo)變量是一個連續(xù)型數(shù)值。
• 數(shù)據(jù)集小：小于100K的樣本量。
• 少數(shù)特征應(yīng)該是重要的：相關(guān)矩陣表明少數(shù)特征包含預(yù)測目標(biāo)變量的信息。

這些特點(diǎn)給予了嶺回歸、支持向量回歸、集成回歸、隨機(jī)森林回歸等方法大展身手的好機(jī)會。有關(guān)回歸模型可參見往期文章??。

線性回歸中的多重共線性與嶺回歸
機(jī)器學(xué)習(xí) | 簡單而強(qiáng)大的線性回歸詳解
機(jī)器學(xué)習(xí) | 深度理解Lasso回歸分析
一文掌握sklearn中的支持向量機(jī)
集成算法 ｜ 隨機(jī)森林回歸模型
萬字長文，演繹八種線性回歸算法最強(qiáng)總結(jié)！
原理+代碼，總結(jié)了 11 種回歸模型
機(jī)器學(xué)習(xí)中回歸算法的基本數(shù)學(xué)原理

我們將對這些模型的性能進(jìn)行如下評估：

x_train?=?train_preprocessed[features].values
y_train?=?train_preprocessed[target].values.ravel()
#?Sort?validation?set?for?plots
val?=?val.sort_values(by=target)
x_val?=?val[features].values
y_val?=?val[target].values.ravel()
x_test?=?test[features].values

table?=?PrettyTable()
table.field_names?=?["Model",?
?????????????????????"Mean?Squared?Error",?
?????????????????????"R2?score"]

models?=?[
????SGDRegressor(max_iter=1000,?tol=1e-3),
????Lasso(alpha=0.1),
????ElasticNet(random_state=0),
????Ridge(alpha=.5),
????SVR(gamma='auto',?kernel='linear'),
????SVR(gamma='auto',?kernel='rbf'),
????BaggingRegressor(),
????BaggingRegressor(KNeighborsClassifier(),?
?????????????????????max_samples=0.5,?
?????????????????????max_features=0.5),
????NuSVR(gamma='auto'),
????RandomForestRegressor(random_state=0,?
??????????????????????????n_estimators=300)
]

for?model?in?models：
????model.fit(x_train,?y_train)?
????y_res?=?model.predict(x_val)

????mse?=?mean_squared_error(y_val,?y_res)
????score?=?model.score(x_val,?y_val)????

????table.add_row([type(model).__name__,?
??????????????????format(mse,?'.2f'),?
??????????????????format(score,?'.2f')])

print(table)

+-----------------------+--------------------+----------+
|?????????Model?????????|?Mean?Squared?Error?|?R2?score?|
+-----------------------+--------------------+----------+
|??????SGDRegressor?????|??????43654.10??????|???0.06???|
|?????????Lasso?????????|??????43103.36??????|???0.07???|
|???????ElasticNet??????|??????54155.92??????|??-0.17???|
|?????????Ridge?????????|??????42963.88??????|???0.07???|
|??????????SVR??????????|??????50794.91??????|??-0.09???|
|??????????SVR??????????|??????41659.68??????|???0.10???|
|????BaggingRegressor???|??????18959.70??????|???0.59???|
|????BaggingRegressor???|??????52475.34??????|??-0.13???|
|?????????NuSVR?????????|??????41517.67??????|???0.11???|
|?RandomForestRegressor?|??????18949.93??????|???0.59???|
+-----------------------+--------------------+----------+

#?Table?setup
table?=?PrettyTable()
table.field_names?=?["Model",?"Dataset",?"MSE",?"MAE",?'RMSLE',?"R2?score"]
#?模型訓(xùn)練
model?=?RandomForestRegressor(bootstrap=True,?criterion='mse',?max_depth=None,
???????????max_features='auto',?max_leaf_nodes=None,
???????????min_impurity_decrease=0.0,?min_impurity_split=None,
???????????min_samples_leaf=1,?min_samples_split=4,
???????????min_weight_fraction_leaf=0.0,?n_estimators=200,?n_jobs=None,
???????????oob_score=False,?random_state=None,?verbose=0,?warm_start=False)
model.fit(x_train,?y_train)?

def?evaluate(x,?y,?dataset)：
????pred?=?model.predict(x)
????mse?=?mean_squared_error(y,?pred)
????mae?=?mean_absolute_error(y,?pred)
????score?=?model.score(x,?y)????
????rmsle?=?np.sqrt(mean_squared_log_error(y,?pred))
????table.add_row([type(model).__name__,?dataset,
??????????????????format(mse,?'.2f'),?
??????????????????format(mae,?'.2f'),?
??????????????????format(rmsle,?'.2f'),?
??????????????????format(score,?'.2f')])
????
evaluate(x_train,?y_train,?'training')
evaluate(x_val,?y_val,?'validation')

print(table)

+-----------------------+------------+----------+-------+-------+----------+
|?????????Model?????????|??Dataset???|???MSE????|??MAE??|?RMSLE?|?R2?score?|
+-----------------------+------------+----------+-------+-------+----------+
|?RandomForestRegressor?|??training??|??298.21??|?10.93?|??0.21?|???0.98???|
|?RandomForestRegressor?|?validation?|?18967.46?|?96.43?|??0.47?|???0.59???|
+-----------------------+------------+----------+-------+-------+----------+

特征重要性

importances?=?model.feature_importances_
std?=?np.std([tree.feature_importances_?for?tree?in?model.estimators_],?axis=0)
indices = np.argsort(importances)[：：-1]

#?打印出特征排序
print("Feature ranking：")

for?f?in?range(x_val.shape[1])：
????print("%d.?feature?%s?(%f)"?%?(f?+?1,?features[indices[f]],?importances[indices[f]]))

Feature ranking：
1. feature hr (0.633655)
2. feature temp (0.160043)
3. feature hum (0.049495)
4. feature workingday (0.046797)
5. feature weathersit (0.027105)
6. feature windspeed (0.026471)
7. feature weekday (0.020293)
8. feature mnth (0.019940)
9. feature season (0.012995)
10. feature holiday (0.003206)

#?繪制出隨機(jī)森林的特征重要性
plt.figure(figsize=(14,5))
plt.title("Feature?importances")
plt.bar(range(x_val.shape[1]),?importances[indices],?
????????color="cornflowerblue",?yerr=std[indices],?align="center")
plt.xticks(range(x_val.shape[1]),?[features[i]?for?i?in?indices])
plt.xlim([-1,?x_val.shape[1]])
plt.show()）

解析：?結(jié)果對應(yīng)于特征相關(guān)矩陣中變量"hour"和變量"temperature"與自行車共享計(jì)數(shù)的高度相關(guān)。

寫在最后

以下是進(jìn)一步提高數(shù)據(jù)模型性能的一些思路：

• 目標(biāo)變量的分布調(diào)整：有些預(yù)測模型假設(shè)目標(biāo)變量的分布為正態(tài)分布，在數(shù)據(jù)預(yù)處理中進(jìn)行轉(zhuǎn)換可以提高這些方法的性能。
• 大規(guī)模數(shù)據(jù)集隨機(jī)森林的實(shí)現(xiàn)。對于大規(guī)模數(shù)據(jù)集（>10 Mio. 樣本），如果不能在工作內(nèi)存中保存所有的樣本，或者會遇到嚴(yán)重的內(nèi)存問題，那么使用python實(shí)現(xiàn)sklearn中的隨機(jī)森林將會非常慢。一個解決方案可以是woody實(shí)現(xiàn)，其中包含用于預(yù)分類的頂樹，以及在頂樹的葉子處用C語言實(shí)現(xiàn)的平坦隨機(jī)森林。