Top > Python

実例:Pythonでデータ分析『相席ラウンジの女性人数を予測せよ!』- 重回帰分析 -



先日、『Pythonでデータ分析をする工程と主力ライブラリ』をご紹介しました。
今回は実例ということで、Data取得からData Cleaning→Data Analysis→Data Visualizeを分析テーマに沿って行っていきます。

そのテーマは、『相席ラウンジの女性人数を予測せよ』です。
某マッチングラウンジ相席トラッカーとして、Webスクレイピングをherokuから定期実行して取得しているデータがありますので、この蓄積データを使って予測モデルを作成してみます。
※相席ラウンジとは、エレガントな相席屋のことで、男女の出会いの場になっていて激混みしている人気スポットです。

相席ラウンジ 相席ラウンジ


Data取得

Data自体は1時間毎にスクレイピングを行って、Google Drive内のcsvに書き足し続けています。スクレイピングのスクリプトは以下のように書くことが出来ます。


from urllib import request as urllib2
from bs4 import BeautifulSoup
import pandas as pd

url = "某相席ラウンジのURL"
soup = BeautifulSoup(urllib2.urlopen(url), "lxml")
top_con = soup.find_all("a")

match_df_latest = pd.DataFrame()
for box in top_con:#店舗毎に処理を回す
    if "icon_mens.png" in str(box):#ソースの構造からポイントを指定
        try:
            temp_df = pd.DataFrame({"Box":[str(box.attrs['id']).split("_")[-1].capitalize()],#ラウンジの店舗
                                    "Date_Time":[str(datetime.datetime.now()).split(":")[0]+":00"],#1時間の代表値として格納
                                    "1_Male":[box.text.split("  ")[0]],#男性
                                    "2_Female":[box.text.split("  ")[-1]]})#女性
            match_df_latest = match_df_latest.append(temp_df)
        except:
            pass


サイトによっては明確にスクレイピングを禁止していますので、行う場合はサイトの規定を読んでから行いましょう。


Data Cleaning

継続的に行う分析の場合は、Google Drive上にあるデータの場合はGoogle Drive APIを使って直接読み込むのが良いのですが、今回はダウンロードしたcsvを読み込むこととします。
データ分析に進むために、データ内容の確認とデータクレンジング・項目追加をしていきます。


import pandas as pd
from pivottablejs import pivot_ui
import pandas_profiling

import datetime
import jpholiday

import matplotlib.pyplot as plt
import seaborn as sns
from pandas_highcharts.display import display_charts

df = pd.read_csv("latest_20180505.csv", header=0)
df.head()


Date_Time Box 1_Male 2_Female
0 2018/2/19 19:00 Sapporo 0 13
1 2018/2/19 19:00 Sendai 1 8
2 2018/2/19 19:00 Shinjuku 19 20
3 2018/2/19 19:00 Shibuya 3 30
4 2018/2/19 19:00 Machida 0 6


データはシンプルです。Timestamp、店舗、男性の人数、女性の人数です。では、続けてクリーニングを行います。


# Cleaning
df = df.dropna()

# 項目追加:DateとTimeを分割 (後で分析の切り口に使う)
df["date"] = df["Date_Time"].apply(lambda x: datetime.datetime.strptime(x.split(" ")[0],'%Y/%m/%d'))
df["date"] = df["date"].apply(lambda x: datetime.date(x.year, x.month, x.day))
df["time"] = df["Date_Time"].apply(lambda x: int(x.split(" ")[1].split(":")[0]))
df["time"] = df["time"].apply(lambda x: [x+24 if x < 12 else x][0])#新宿渋谷の1日は昼の12時起点なのだ
df["date"] = df.apply(lambda row: [row["date"]-datetime.timedelta(days=1) if row["time"]>=24 else row["date"]][0],axis=1)

# 項目追加:曜日、休日(土日+祝日)、休前日が分かるようにする (後で分析の切り口に使う)
national_holiday = list(dict(jpholiday.holidays(df["date"].min(), df["date"].max())).keys())#祝日APIで期間内の祝日を取得
df["day"] = df["date"].apply(lambda x:x.isoweekday())#曜日を付与
df["holiday"] = df["date"].apply(lambda x:[1 if (x.isoweekday()>5)|(x in national_holiday) else 0][0])#土日祝を休日とする
all_holiday = list(set(df[df["holiday"]==1]["date"]))#全休日のリストを生成
df["daybeforeholiday"] = df["date"].apply(lambda x:[1 if x  + datetime.timedelta(days=1) in all_holiday else 0][0])#翌日が休日の場合のフラグ

# Filtering:
df = df[df["date"] >= datetime.date(2018,4,1)]#最新の状況から予測をしたいので、4月以降に絞りました

# Filtering:店舗が多いので、主要どころの新宿と渋谷に絞込んでみる
box_filter = ["Shibuya","Shinjuku"]
df = df[df["Box"].str.contains("|".join(box_filter),na=False)]

# ダミー変数の追加(drop_first=Trueとして変数が従属関係になるのを防ぐことが多いのですが、今回は結果的に精度が高まったのでFalseにしています)
df_dummy_d = pd.get_dummies(df["day"], drop_first=False).rename(columns={1:"Mon",2:"Tue",3:"Wed",4:"Thu",5:"Fri",6:"Sat",7:"Sun"})
df_dummy_t = pd.get_dummies(df["time"], drop_first=False)
df = pd.merge(df, df_dummy_d, left_index=True, right_index=True)
df = pd.merge(df, df_dummy_t, left_index=True, right_index=True)

# ここまでで、pivot_ui(df) で可視化するとデータの状況が分かるのですが、
# ここでは散布図とヒストグラムを新宿と渋谷それぞれで可視化したいので、matplotlibを使ってみます

# histgram重ね書き
def hist(df):
    df = pd.DataFrame(df)
    a=0
    color_list=['b','r']
    for i in df.columns:
        df[i].hist(bins=50, range=(0, 150), normed=False, alpha=.3, color=color_list[a], label=i)
        a+=1

    plt.xlabel('Number of participants')
    plt.ylabel('The number of sample')
    plt.legend()
    plt.show()
    return  None

print("新宿 相席ラウンジ")#営業時間に絞る
hist(df[(df["time"] > 17)&(df["time"] < 29)&(df["Box"]=="Shinjuku")][["1_Male","2_Female"]])

print("渋谷 相席ラウンジ")#営業時間に絞る
hist(df[(df["time"] > 17)&(df["time"] < 29)&(df["Box"]=="Shibuya")][["1_Male","2_Female"]])


新宿 相席ラウンジ


渋谷 相席ラウンジ


上記のヒストグラムで人数の分布を見ると、新宿は想像していた分布に近かったのですが、渋谷は60人辺りに不自然な崖があります。
実地観察を行っていなくて恐縮ですが、おそらく店舗サイズによる人数制限があるものと思われます。ということで、値の幅が大きい新宿で予測を進めます!


# 散布図でもデータを確認しておく
sns.jointplot(x="1_Male", y="2_Female", data=df[(df["time"] > 17)&(df["time"] < 29)&(df["Box"]=="Shinjuku")])
sns.plt.show()




当然ですが、男女で相関関係は強いですね。そして女性が多い・・・まあ無料ですもんね・・・。
Highchartで時系列トレンドも確認してみます。


display_charts(df[df["Box"]=="Shinjuku"][["1_Male","2_Female"]],title="新宿相席ラウンジ")




上記を見ると、先程切り口として追加しておいた休日・休前日・曜日、あと当然ですが時間帯で傾向がありそうですので、作成した変数を元に重回帰分析をしてみます。


Data Analytics

ここでは、お題である『相席ラウンジの女性人数を予測せよ』を解決する予測値を得るために、重回帰分析を行います。
『は?時系列データを重回帰分析するの?』って思った方は正解だと思います。
時系列データの場合は、ARIMAモデルや状態空間モデルを使って予測するのが通常かと思うのですが、 ここでは時間の情報自体を説明変数として、さくっと予測値への寄与度も可視化するために無理やり重回帰分析を選択しました。

目的変数を女性の来客数とし、説明変数に休日区分・休前日区分・曜日のダミー変数・時間帯のダミー変数を使います。
なので、同じ曜日の同じ時間帯であれば、当日もしくは翌日が祝日でない限りは同じ値が出力されるわけです。
それ以上精度を高めようとすると、違う予測モデル使うか、追加変数としてイベントや天候などの情報を追加収集する必要がありそうです。
という前提のもと、statsmodelsを使用して進めてみます。ちなみにstatsmodelsはsklearnの回帰分析結果と比較すると出力変数が多いのが良いところです。


import statsmodels.api as sm
import math
import numpy as np

target_box = "Shinjuku"
target_gender = "2_Female"

# statsmodels
Y = df[df["Box"]==target_box][[target_gender]]
X = df[df["Box"]==target_box][['holiday', 'daybeforeholiday',
                               'Mon', 'Tue', 'Wed', 'Thu', #'Fri', 'Sat', 'Sun'は'holiday', 'daybeforeholiday'とかぶるので除外した
                               18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28]]

ols_model = sm.OLS(Y,sm.add_constant(X)).fit()
ols_model.summary()

OLS Regression Results
Dep. Variable: 2_Female R-squared: 0.850
Model: OLS Adj. R-squared: 0.847
Method: Least Squares F-statistic: 264.7
Date: Mon, 06 May 2018 Prob (F-statistic): 3.24e-313
Time: 21:04:17 Log-Likelihood: -3156.5
No. Observations: 813 AIC: 6349.
Df Residuals: 795 BIC: 6434.
Df Model: 17
Covariance Type: nonrobust
coef std err t P>|t| [0.025 0.975]
const -6.3527 1.424 -4.462 0.000 -9.147 -3.558
holiday 5.7171 1.161 4.925 0.000 3.438 7.996
daybeforeholiday 13.3847 1.149 11.653 0.000 11.130 15.639
Mon -0.5352 1.664 -0.322 0.748 -3.801 2.731
Tue -0.1483 1.752 -0.085 0.933 -3.587 3.291
Wed -0.6742 1.648 -0.409 0.683 -3.909 2.560
Thu -0.4485 1.528 -0.293 0.769 -3.449 2.552
18 14.6439 2.404 6.092 0.000 9.926 19.362
19 44.9921 2.061 21.834 0.000 40.947 49.037
20 65.6706 2.116 31.043 0.000 61.518 69.823
21 74.9358 2.061 36.367 0.000 70.891 78.981
22 73.7929 2.010 36.716 0.000 69.848 77.738
23 58.6875 2.061 28.473 0.000 54.642 62.734
24 37.7862 2.061 18.332 0.000 33.740 41.832
25 41.9353 2.116 19.823 0.000 37.783 46.088
26 36.2883 2.116 17.154 0.000 32.136 40.441
27 28.0143 2.088 13.417 0.000 23.916 32.113
28 17.9059 2.116 8.464 0.000 13.753 22.059
Omnibus: 7.903 Durbin-Watson: 0.465
Prob(Omnibus): 0.019 Jarque-Bera (JB): 8.434
Skew: 0.179 Prob(JB): 0.0147
Kurtosis: 3.348 Cond. No. 9.02

R-squaredは0.85ということでそこそこ高い値になっています。変数別に見ると、休日や休前日はプラスの寄与がある等、納得感があります。
ステップワイズ変数選定は省略します。上記から得られたモデルと元に実際の値と並べてみます。
重回帰分析の性質上、営業時間外の時間帯でも0人とならないので、ここは感と経験と真心で値の付け替えを行うこととします。

もっとちゃんと勉強をしたいという方は『Python関連のオススメの本』をご参照ください。


# Actualの値とSimulatedの値をJoinしていきます
reg_df = pd.DataFrame()
reg_df["actual"] = df[df["Box"]==target_box][target_gender]
reg_df["time"] = df[df["Box"]==target_box]["time"]
reg_df["simulated"] = ols_model.predict(sm.add_constant(X))

# マイナスはおかしいので0に付け替えと、営業時間外は0人に付け替えします
reg_df["simulated"] = reg_df["simulated"].apply(lambda x:[0 if x<0 else x][0])
reg_df["simulated"] = reg_df.apply(lambda row:[0 if (row["time"]>28)|(row["time"]<18) else row["simulated"]][0],axis=1)
reg_df.index = df[df["Box"]=="Shinjuku"]["Date_Time"]

# 直近の数字を様子見
display_charts(reg_df.iloc[500:,:][["actual","simulated"]], kind="line", title="新宿相席ラウンジの女性客数予測")



それっぽい感じになっていますが、4月27日以降の20〜22時のスパイクで過小評価になっていますね。
GW中の盛況っぷりを予測できなかったためと思われます。

と、いうことでお題の『相席ラウンジの女性人数を予測せよ』を一応達成しました。
実際の予測値はWebアプリとして近日中に『某マッチングラウンジ相席トラッカー』に機能追加する予定です。

さらに精度を高めるには、予測モデルを変えるか、イベント情報・天候などを考慮すると良さそうです。
今回の結果は特定の店舗だけでなく、同業も類似の動きになっているように思えますので、出会いを求める肉食獣の方々は参考にしていただければと思います。

最後に、全実値と予測値の差異の分布をみるために、散布図で表してみましょう。


plt.subplots(figsize=(7, 6))
x_ = np.arange(0, 120, 1)[:,np.newaxis]
plt.plot(x_,x_,"b-", linewidth=1)
plt.scatter(reg_df["actual"], reg_df["simulated"], color="black")
plt.xlim(0,120);plt.ylim(0,120);plt.xlabel("Actual");plt.ylabel("Simulated")
plt.show()



はやりバラツキがありますね(=あまり精度がよく見えませんね)。予測精度の計算方法はいくつかありますが、Root Mean Squared Error (RMSE)か、
[ 予測精度 = 100% - (|予測値 - 実値|) / 実値 ] というような計算式を使うことが多いように思います。(※別の予測モデルも試してから、実際の予測精度を追記予定)


Data Visualize

今回は上記の通り分析しながら可視化とコメントを記載してしまったので、Appendix的に相席ラウンジさん側のためにも行くと良い時間帯を紹介します。

print("来客女性の平均余剰人数  *1:月〜7:日")
df["Chance"] = df["2_Female"] - df["1_Male"]
df[(df["time"]>=18)&(df["time"]<=28)&(df["Box"]=="Shinjuku")].pivot_table("Chance",columns="day",index="time",aggfunc="mean").astype(int)


来客女性の平均余剰人数
day Mon Tue Wed Thu Fri Sat Sun
time
18 7 11 11 12 9 17 18
19 19 19 19 28 35 44 39
20 19 23 28 22 44 37 34
21 7 14 17 18 36 15 28
22 12 6 6 3 11 6 19
23 -1 -3 -1 -8 4 8 7
24 -3 -10 -6 -9 -18 -12 -2
25 3 -4 -3 0 -10 1 5
26 1 0 -1 0 -3 -3 6
27 0 0 2 0 2 7 3
28 2 2 2 5 0 17 3

どの曜日も19-20時は特に女性が多いようです。金曜は21時台でも同様の傾向が見られまし、土日は特に持て余している女性の方が多いようです。
逆に金曜の深夜帯は男性が過剰気味ですので、行くなら早めの時間が良さそうですよ。土曜の深夜というか日曜の朝方に女性増えているのはなんでしょうね。


ということで、『相席ラウンジの女性人数を予測せよ』をテーマに分析の実例をご紹介しました。次回はもう少しちゃんとレポーティングっぽくしようかと思っています。ご期待ください。




Ads by Google