In [121]:
import pandas as pd
src_data = pd.read_csv("joto_data_summary.csv")
src_data.head()
Out[121]:
車の目視台数とパケット数の散布図¶
回帰分析の前に、実測データ関係を見てみる
In [122]:
e_df = src_data.loc[:,['wifi_east','car_east']]
e_df.rename(columns={'wifi_east': 'Wi-Fiアドレス数', 'car_east': '台数(目視)'}, inplace=True)
e_df['direction'] = "東向き"
w_df = src_data.loc[:,['wifi_west','car_west']]
w_df.rename(columns={'wifi_west': 'Wi-Fiアドレス数', 'car_west': '台数(目視)'}, inplace=True)
w_df['direction'] = "西向き"
df4plot = pd.concat([e_df, w_df])
import seaborn as sns
sns.set(font_scale=1.5, font='IPAexGothic')
sns.pairplot(y_vars=['Wi-Fiアドレス数'], x_vars=['台数(目視)'], data=df4plot, hue='direction', height=6)
#src_data.plot(kind='scatter', x=['wifi_east','wifi_west'], y=['car_east','car_west'])
Out[122]:
- 東向きでアドレス数が多いのは11/29(金)18時. ダン⇒風月堂の歩行者数が多いのだと考えられる
- 車の台数以上にアドレス数の多い4点は、金、土の17,18時 （渋滞によるものと考えられる）
- 日曜日の15時に西向きアドレス数が多いのは不明（歩行者？）
- 全体としてダン⇔風月堂間の「Wi-Fi流量」は車によるものである。
Google 所要時間データの規格化¶
捕捉パケット数は「目視車台数ｘ渋滞率」に比例することを想定し、所要時間を最大値でわったデータを作成（「規格化」）
In [123]:
# 所要時間 (east_bound, west_bound)を最大値で規格化
src_data['norm_w_bound'] = src_data['west_bound'] / src_data['west_bound'].max()
src_data['norm_e_bound'] = src_data['east_bound'] / src_data['west_bound'].max()
# 規格化所要時間を車台数にかけて捕捉車台数（補正値）とする
src_data['scaled_car_east'] = src_data['car_east'] * src_data['norm_e_bound']
src_data['scaled_car_west'] = src_data['car_west'] * src_data['norm_w_bound']
src_data
Out[123]:
In [124]:
# 機械学習ライブラリsklearnでは、データはnumpyなので
# np.ndarray配列として作成
import numpy as np
# target （目的変数）Wi-Fi east & west
target_east = src_data['wifi_east'].values
target_west = src_data['wifi_west'].values
# type(target_east) # typeで確かめ
# (1) features (説明変数) pedes_total_* （歩行量調査データ）, car_* （車の目視データ）
features_east = src_data[['pedes_total_east', 'car_east']].values
features_west = src_data[['pedes_total_west', 'car_west']].values
# (2) features (説明変数) pedes_total_* （歩行量調査データ）, scaled_car_* （車の目視データを所要時間で補正したもの）
features_east_sc = src_data[['pedes_total_east', 'scaled_car_east']].values
features_west_sc = src_data[['pedes_total_west', 'scaled_car_west']].values
features_east[0:5,:] # 最初の5行をみてみる
Out[124]:
回帰分析¶
- 切片なし
- 渋滞係数なしの分析 (歩行量調査データと車通行量目視データ)
In [125]:
# ライブラリをロード
from sklearn.linear_model import LinearRegression
# 線形回帰機
regression = LinearRegression(fit_intercept=False)
# 渋滞係数なし
e_model = regression.fit(features_east, target_east)
print(e_model.coef_)
#print(e_model.intercept_)
w_model = regression.fit(features_west, target_west)
print(w_model.coef_)
#print(w_model.intercept_)
east_coef = e_model.coef_
west_coef = w_model.coef_
# 渋滞係数あり
s_e_model = regression.fit(features_east_sc, target_east)
print(s_e_model.coef_)
s_w_model = regression.fit(features_west_sc, target_west)
print(s_w_model.coef_)
s_east_coef = s_e_model.coef_
s_west_coef = s_w_model.coef_
s_e_model.predict(features_east_sc)
Out[125]:
Wi-Fiパケットデータと回帰係数からの推定値の比較¶
In [126]:
# 歩行量、車量に回帰係数をかけてWi-Fi推定値を算出
estimate_e = features_east[:,0]* east_coef[0] + features_east[:,1]* east_coef[1]
estimate_w = features_west[:,0]* west_coef[0] + features_west[:,1]* west_coef[1]
estimate_s_e = features_east_sc[:,0]* s_east_coef[0] + features_east_sc[:,1]* s_east_coef[1]
estimate_s_w = features_west_sc[:,0]* s_west_coef[0] + features_west_sc[:,1]* s_west_coef[1]
In [127]:
# 描画
from matplotlib import pyplot as plt
fig = plt.figure(figsize=(12,6))
ax1 = fig.add_subplot(1,2,1)
ax1.scatter(estimate_e, src_data['wifi_east'].values, c = 'red', label='東向き')
ax1.scatter(estimate_w, src_data['wifi_west'].values, c = 'blue', label='西向き')
ax1.legend(fontsize=14)
ax1.set_title('渋滞補正なし')
ax2 = fig.add_subplot(1,2,2)
ax2.scatter(estimate_s_e, src_data['wifi_east'].values, c = 'red', label='東向き')
ax2.scatter(estimate_s_w, src_data['wifi_west'].values, c = 'blue', label='西向き')
ax2.legend(fontsize=14)
ax2.set_title('渋滞補正あり')
fig.show()
In [128]:
# (1) features (説明変数) car_* （車の目視データ）
features_east = src_data[['car_east']].values
features_west = src_data[['car_west']].values
# (2) features (説明変数) scaled_car_* （車の目視データを所要時間で補正したもの）
features_east_sc = src_data[['scaled_car_east']].values
features_west_sc = src_data[['scaled_car_west']].values
In [129]:
# 線形回帰機
regression = LinearRegression(fit_intercept=False) # 切片の有無 (外れ値があるためか切片を入れると推定値がずれる)
# 渋滞係数なし
e_model = regression.fit(features_east, target_east)
print(e_model.coef_)
print(e_model.intercept_)
w_model = regression.fit(features_west, target_west)
print(w_model.coef_)
print(w_model.intercept_)
east_coef = e_model.coef_
west_coef = w_model.coef_
# 渋滞係数あり
s_e_model = regression.fit(features_east_sc, target_east)
print(s_e_model.coef_)
print(s_e_model.intercept_)
s_w_model = regression.fit(features_west_sc, target_west)
print(s_w_model.coef_)
print(s_w_model.intercept_)
s_east_coef = s_e_model.coef_
s_west_coef = s_w_model.coef_
In [130]:
s_e_model.predict(features_east_sc) # 下のように計算しなくても、predict()で推定値を得られる
Out[130]:
In [131]:
# 歩行量、車量に回帰係数をかけてWi-Fi推定値を算出
estimate_e = features_east[:,0]* east_coef[0] + e_model.intercept_
estimate_w = features_west[:,0]* west_coef[0] + w_model.intercept_
estimate_s_e = features_east_sc[:,0]* s_east_coef[0] + s_e_model.intercept_
estimate_s_w = features_west_sc[:,0]* s_west_coef[0] + s_w_model.intercept_
In [132]:
estimate_s_e
Out[132]:
In [133]:
# 描画
from matplotlib import pyplot as plt
fig = plt.figure(figsize=(14,6))
ax1 = fig.add_subplot(1,2,1)
ax1.scatter(estimate_e, target_east, c = 'red', label='東向き')
ax1.scatter(estimate_w, target_west, c = 'blue', label='西向き')
ax1.legend(fontsize=14)
ax1.set_title('渋滞補正なし')
ax1.set_xlabel('推定値')
ax1.set_ylabel('流動Wi-Fiアドレス数 / 時')
ax2 = fig.add_subplot(1,2,2)
ax2.scatter(estimate_s_e, target_east, c = 'red', label='東向き')
ax2.scatter(estimate_s_w, target_west, c = 'blue', label='西向き')
ax2.legend(fontsize=14)
ax2.set_title('渋滞補正あり')
ax2.set_xlabel('推定値')
ax2.set_ylabel('流動Wi-Fiアドレス数 / 時')
ax1.set_xlim([0,62])
ax1.set_ylim([0,62])
ax2.set_xlim([0,62])
ax2.set_ylim([0,62])
fig.suptitle("自動車目視数による単回帰 (切片なし)")
fig.show()
In [134]:
src_data['estimated_e'] = estimate_s_e
src_data['estimated_w'] = estimate_s_w
In [135]:
src_data
Out[135]:
In [139]:
fig, axes = plt.subplots(1, 3, figsize=(15, 6))
a = src_data.groupby(['date'])
a.get_group((20191129)).plot(x='hour', y='wifi_west', label='西向き', ax=axes[0])
a.get_group((20191129)).plot(x='hour', y='estimated_w', label='西向き(推定値)', ax=axes[0])
a.get_group((20191129)).plot(x='hour', y='wifi_east', label='東向き', ax=axes[0])
a.get_group((20191129)).plot(x='hour', y='estimated_e', label='東向き(推定値)', ax=axes[0])
a.get_group((20191130)).plot(x='hour', y='wifi_west', label='西向き', ax=axes[1])
a.get_group((20191130)).plot(x='hour', y='estimated_w', label='西向き(推定値)', ax=axes[1])
a.get_group((20191130)).plot(x='hour', y='wifi_east', label='東向き', ax=axes[1])
a.get_group((20191130)).plot(x='hour', y='estimated_e', label='東向き(推定値)', ax=axes[1])
a.get_group((20191201)).plot(x='hour', y='wifi_west', label='西向き', ax=axes[2])
a.get_group((20191201)).plot(x='hour', y='estimated_w', label='西向き(推定値)', ax=axes[2])
a.get_group((20191201)).plot(x='hour', y='wifi_east', label='東向き', ax=axes[2])
a.get_group((20191201)).plot(x='hour', y='estimated_e', label='東向き(推定値)', ax=axes[2])
axes[0].set_ylim([0,65])
axes[1].set_ylim([0,65])
axes[2].set_ylim([0,65])
axes[0].set_xticks(np.linspace(10, 19, 4))
axes[1].set_xticks(np.linspace(10, 19, 4))
axes[2].set_xticks(np.linspace(10, 19, 4))
axes[0].set_ylabel("アドレス数")
#axes[1].set_ylabel("アドレス数")
#axes[2].set_ylabel("アドレス数")
axes[0].set_title("2020/11/29")
axes[1].set_title("2020/11/30")
axes[2].set_title("2020/12/01")
axes[1].legend().set_visible(False)
axes[2].legend().set_visible(False)
fig.suptitle("自動車目視数によるWi-Fiパケット流量推定 (LM 切片なし)")
Out[139]:
サポートベクトル回帰 (SVM)¶
In [141]:
from sklearn import svm
# 学習を行う
svr = svm.SVR(kernel='rbf', gamma='scale')
# 渋滞係数あり
s_e_model = svr.fit(features_east_sc, target_east)
s_w_model = svr.fit(features_west_sc, target_west)
# 推定値
s_e_model.predict(features_east_sc)
s_w_model.predict(features_west_sc)
# 推定値をもとのデータフレームの列に追加
src_data['estimated_e'] = estimate_s_e
src_data['estimated_w'] = estimate_s_w
# 回帰曲線を描く (LMの時と同じ)
fig, axes = plt.subplots(1, 3, figsize=(15, 6))
a = src_data.groupby(['date'])
a.get_group((20191129)).plot(x='hour', y='wifi_west', label='西向き', ax=axes[0])
a.get_group((20191129)).plot(x='hour', y='estimated_w', label='西向き(推定値)', ax=axes[0])
a.get_group((20191129)).plot(x='hour', y='wifi_east', label='東向き', ax=axes[0])
a.get_group((20191129)).plot(x='hour', y='estimated_e', label='東向き(推定値)', ax=axes[0])
a.get_group((20191130)).plot(x='hour', y='wifi_west', label='西向き', ax=axes[1])
a.get_group((20191130)).plot(x='hour', y='estimated_w', label='西向き(推定値)', ax=axes[1])
a.get_group((20191130)).plot(x='hour', y='wifi_east', label='東向き', ax=axes[1])
a.get_group((20191130)).plot(x='hour', y='estimated_e', label='東向き(推定値)', ax=axes[1])
a.get_group((20191201)).plot(x='hour', y='wifi_west', label='西向き', ax=axes[2])
a.get_group((20191201)).plot(x='hour', y='estimated_w', label='西向き(推定値)', ax=axes[2])
a.get_group((20191201)).plot(x='hour', y='wifi_east', label='東向き', ax=axes[2])
a.get_group((20191201)).plot(x='hour', y='estimated_e', label='東向き(推定値)', ax=axes[2])
axes[0].set_ylim([0,65])
axes[1].set_ylim([0,65])
axes[2].set_ylim([0,65])
axes[0].set_xticks(np.linspace(10, 19, 4))
axes[1].set_xticks(np.linspace(10, 19, 4))
axes[2].set_xticks(np.linspace(10, 19, 4))
axes[0].set_ylabel("アドレス数")
#axes[1].set_ylabel("アドレス数")
#axes[2].set_ylabel("アドレス数")
axes[0].set_title("2020/11/29")
axes[1].set_title("2020/11/30")
axes[2].set_title("2020/12/01")
axes[1].legend().set_visible(False)
axes[2].legend().set_visible(False)
fig.suptitle("自動車目視数によるWi-Fiパケット流量推定 (SVM)")
Out[141]:
SVMと切片無しのLMとはほとんど同じ結果になる。
付録：練習¶
In [2]:
# 練習 python機械学習クックブック １３．２
# -*- coding: utf-8 -*-
# ライブラリをロード
from sklearn.linear_model import LinearRegression
from sklearn.datasets import load_boston
from sklearn.preprocessing import PolynomialFeatures
# データをロード with only two features
boston = load_boston()
features = boston.data[:,0:2]
target = boston.target
# boston {'data': [[]], 'target': [], 'feature_names': []}
# features は dataの最初の2列
# 交互作用の項を作成
interaction = PolynomialFeatures(
degree=3, include_bias=False, interaction_only=True)
features_interaction = interaction.fit_transform(features)
# 線形回帰器を作成
regression = LinearRegression()
# 線形回帰器を訓練
model = regression.fit(features_interaction, target)
##########
# 最初の観測値の特徴量の値を表示
features[0]
##########
# ライブラリをロード
import numpy as np
# 個々の観測値に対して、最初の特徴量値と2つ目の特徴量値を掛け合わせる
interaction_term = np.multiply(features[:, 0], features[:, 1])
##########
# 最初の観測値の交互作用項を表示
interaction_term[0]
##########
# 最初の観測値の値を表示
features_interaction[0]
Out[2]:
In [10]:
model.coef_
Out[10]:
In [11]:
model.intercept_
Out[11]:
Googleの道路状況がかなり即時的に表示されていることを示す面白い実験例
In [ ]: