データとかITとか統計とか好きなこと

データや統計に関すること&ITなどが好きで仕事にも生かせてます。音楽とかyoutubeも好き。

約束のネバーランドが連載完結したので、Googleトレンドを分析してみた。対鬼滅の刃

f:id:rautaku:20200615234051j:plain
main_never

本日のジャンプで約束のネバーランドが終了しましたね。日曜日の夜にYahooトップニュースで取り上げられていて知りました。 AmazonPrimeでアニメを見て好きになった作品です。アニメ第2期は2021年放送されるらしいので楽しみです。

鬼滅の刃に続いて人気作を4年程で終わらせたことも今回話題になった理由であると思いますが、ここで「鬼滅の刃と比べて約束のネバーランドってどのくらい人気だったんだろう」とふと気になりました。 個人的には約束のネバーランドの方が好みですが、鬼滅の刃の方が圧倒的に世間を賑わせているイメージでした。

そこで今回は完全に思いつきで、Googleトレンドを使って、「約束のネバーランド」と「鬼滅の刃」のキーワードのトレンド分析をしてみようと思います。思いつきなのでさっくりです。 分析した際のPythonコードも載せているので同じように分析をしたい人などは参考にしてください

GoogleトレンドAPIや時系列分析などの勉強にもなるかと思ったので)

作品基本情報

累計発行部数: 2020年6月時点 2100万部超
連載: 2016年8月1日 - 2020年6月15日
アニメ: 2019年1月 - 3月
累計発行部数: 2020年6月時点 6000万部超
連載: 2016年2月15日 - 2020年5月18日
アニメ: 2019年4月 - 9月

連載開始もアニメ化も割と近い印象の両作ですが、発行部数は鬼滅の刃がトリプルスコア

GoogleトレンドデータをAPIで簡単に取得

pytrendsというパッケージを使えば簡単にデータフレームとして取得できます!今回は「約束のネバーランド」と「鬼滅の刃」の2つのキーワード比較です。

日曜始まりの1週間のトレンド値を取得できます(注意点はGoogleトレンドで取得できるデータは絶対値算出ではなく、あくまで0-100で標準化された相対値です)

import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
%matplotlib inline
sns.set()
import warnings
warnings.filterwarnings('ignore')
from pytrends.request import TrendReq
from datetime import datetime

pytrends = TrendReq(hl='ja-JP', tz=360)
kw_list = ['約束のネバーランド', '鬼滅の刃']
pytrends.build_payload(kw_list, timeframe='2016-06-01 2018-12-31', geo='JP')
df = pytrends.interest_over_time().reset_index().rename(columns={'約束のネバーランド': 'neverland', '鬼滅の刃': 'kimetsu'})
df['date'] = pd.to_datetime(df['date'])
df
date neverland kimetsu isPartial
0 2016-06-05 0 3 False
1 2016-06-12 0 3 False
2 2016-06-19 0 3 False
3 2016-06-26 0 2 False
4 2016-07-03 0 1 False
... ... ... ... ...
130 2018-12-02 32 19 False
131 2018-12-09 28 24 False
132 2018-12-16 36 28 False
133 2018-12-23 56 27 False
134 2018-12-30 100 30 False

135 rows × 4 columns

連載開始 ~ アニメ化前

f:id:rautaku:20200615234241j:plain
never_start

では早速可視化などして分析してみましょう。まずは連載開始からアニメ化前の比較。お互いの連載が開始ており、アニメが始まる前の2016年8月1日 〜 2018年12月31日で見てみましょう

アニメ化前までは「約束のネバーランド」の方がトレンドだった

pytrends = TrendReq(hl='ja-JP', tz=360)
kw_list = ['約束のネバーランド', '鬼滅の刃']
pytrends.build_payload(kw_list, timeframe='2016-08-01 2018-12-31', geo='JP')
df = pytrends.interest_over_time().reset_index().rename(columns={'約束のネバーランド': 'neverland', '鬼滅の刃': 'kimetsu'})
df['date'] = pd.to_datetime(df['date'])
fig, ax = plt.subplots(figsize=(16, 4))
df.set_index('date')['neverland'].plot(ax=ax, label='neverland', color='orange')
df.set_index('date')['kimetsu'].plot(ax=ax, label='kimetsu', color='black')
plt.legend()
plt.ylabel('%')
plt.show()

f:id:rautaku:20200615233505p:plain
kimetsu_never_1

約束のネバーランドの方が連載半年ほどで目立ち始め、鬼滅の刃よりも高い人気だったことが読み取れます。

周期的なピークは単行本などの発売日でしょうが、約束のネバーランド初期はしっかり都度節目ごとに推されていてトレンドになるように仕組まれていたのではないか、また読者側もしっかりそれに反応するほど人気を得ていたのではないか、と思います。

一方で、鬼滅の刃はピークも少し弱い印象です。 しかし後半は追いついてきており、ここら辺から鬼滅の刃もアニメ化の話などが動き出したんですかね?(推測です)

約束のネバーランドアニメ開始

f:id:rautaku:20200615234314j:plain
never_anime
次は、約束のネバーランドのアニメが放送されていた2019年1月1日 〜 2019年4月1日で見てみましょう

アニメ化された「約束のネバーランド」は毎週コンスタントに話題に

pytrends = TrendReq(hl='ja-JP', tz=360)
kw_list = ['約束のネバーランド', '鬼滅の刃']
pytrends.build_payload(kw_list, timeframe='2019-01-01 2019-4-1', geo='JP')
df = pytrends.interest_over_time().reset_index().rename(columns={'約束のネバーランド': 'neverland', '鬼滅の刃': 'kimetsu'})
df['date'] = pd.to_datetime(df['date'])
fig, ax = plt.subplots(figsize=(16, 4))
df.set_index('date')['neverland'].plot(ax=ax, label='neverland', color='orange')
df.set_index('date')['kimetsu'].plot(ax=ax, label='kimetsu', color='black')
plt.axvline(datetime(2019, 1, 11), color='red')
plt.legend()
plt.ylabel('%')
plt.show()

f:id:rautaku:20200615233838p:plain
outpu8

約束のネバーランドが毎週トレンドをコンスタントに維持しています。初回話題になり、少し半ばでだれ始めるが、最終回に向けてはまた持ち直していった様子が見えます。このような波形を描くのが普通のアニメじゃないかなという印象です。

この頃は鬼滅の刃は潜伏期間と言えるでしょうか。アニメ放送に向けて準備している時期ですかね。

鬼滅の刃アニメ開始

f:id:rautaku:20200615234345j:plain
kimetsu_anime

次は、いよいよの鬼滅の刃のアニメが放送されていた期間を見てみましょう。比較のために、ネバーランドのアニメ化時期も含めてみます

鬼滅の刃がアニメ化によって一気に逆転!しかも伸び率が段違い

pytrends = TrendReq(hl='ja-JP', tz=360)
kw_list = ['約束のネバーランド', '鬼滅の刃']
pytrends.build_payload(kw_list, timeframe='2019-01-01 2019-10-1', geo='JP')
df = pytrends.interest_over_time().reset_index().rename(columns={'約束のネバーランド': 'neverland', '鬼滅の刃': 'kimetsu'})
df['date'] = pd.to_datetime(df['date'])
fig, ax = plt.subplots(figsize=(16, 4))
df.set_index('date')['neverland'].plot(ax=ax, label='neverland', color='orange')
df.set_index('date')['kimetsu'].plot(ax=ax, label='kimetsu', color='black')
plt.axvline(datetime(2019, 4, 6), color='red')
plt.legend()
plt.ylabel('%')
plt.show()

f:id:rautaku:20200615233923p:plain
output11

はい、鬼滅の刃のアニメ化による伸びが圧倒的です。こんなに圧倒的な伸びだったとは、、! しかもさっきの約束のネバーランドとは違って、後半最終回に向けて指数関数的にぐんぐん伸びていったんだなというのがすごいです。 自分はアニメを鬼滅の刃のアニメを見てないんですが、後半めちゃくちゃ面白くなっていったんじゃんないかな、と見たくなりました!

一方で約束のネバーランドはアニメ終了後、下降トレンドになってしまっていったようです。 もし第2期をすぐ開始していればもう少しこのトレンドを維持できていたのではないか、とか思ったりします。

アニメ終了から最近まで

そして最後、お互いのアニメも終わってからの直近までを見てみましょう

鬼滅の刃がアニメ化終了後もそのままの勢いでぐんぐん伸びていき、大大大人気作品に

pytrends = TrendReq(hl='ja-JP', tz=360)
kw_list = ['約束のネバーランド', '鬼滅の刃']
pytrends.build_payload(kw_list, timeframe='2019-01-01 2020-06-01', geo='JP')
df = pytrends.interest_over_time().reset_index().rename(columns={'約束のネバーランド': 'neverland', '鬼滅の刃': 'kimetsu'})
df['date'] = pd.to_datetime(df['date'])
fig, ax = plt.subplots(figsize=(16, 4))
df.set_index('date')['neverland'].plot(ax=ax, label='neverland', color='orange')
df.set_index('date')['kimetsu'].plot(ax=ax, label='kimetsu', color='black')
plt.axvline(datetime(2019, 4, 6), color='red')
plt.axvline(datetime(2019, 9, 28), color='red')
plt.legend()
plt.ylabel('%')
plt.show()

f:id:rautaku:20200615233952p:plain
output14

鬼滅の刃はもう完全に波に乗りましたね。アニメが終了してもほとんど影響がなく、人気が人気を呼ぶ感じ。こうなったらビジネスとしては大成功でしょう。 もうトレンドとしては約束のネバーランドはかなり開きがついてしまいました。

初期は、約束のネバーランドが人気だったのは確実だと思いますが、鬼滅の刃がアニメ化によってここまで圧倒的な人気作品になるとはびっくりです。 アニメ化が大成功するとこんなことになるんですね。

これを実現させたマーケティング手法とかとても興味が湧きます。

今回は以上です。あまり時系列分析はできなかった。。

おまけ

対ワンピース

f:id:rautaku:20200615234519j:plain
onepiece

自分はワンピナルト鰤世代なので、ワンピースとの比較はやはり気になったので調べてみた

pytrends = TrendReq(hl='ja-JP', tz=360)
kw_list = ['約束のネバーランド', '鬼滅の刃', 'ワンピース']
pytrends.build_payload(kw_list, timeframe='2019-01-01 2020-06-01', geo='JP')
df = pytrends.interest_over_time().reset_index().rename(columns={'約束のネバーランド': 'neverland', '鬼滅の刃': 'kimetsu', 'ワンピース': 'onepiece'})
df['date'] = pd.to_datetime(df['date'])
fig, ax = plt.subplots(figsize=(16, 4))
df.set_index('date')['neverland'].plot(ax=ax, label='neverland', color='orange')
df.set_index('date')['kimetsu'].plot(ax=ax, label='kimetsu', color='black')
df.set_index('date')['onepiece'].plot(ax=ax, label='onepiece', color='blue')
plt.axvline(datetime(2019, 4, 6), color='red')
plt.axvline(datetime(2019, 9, 28), color='red')
plt.legend()
plt.ylabel('%')
plt.show()

f:id:rautaku:20200615234011p:plain
output17

ワンピースは安定したハイトレンドですね。さすが看板。 でもこの短期感でも一気に抜ききった鬼滅の方がすごすぎるという見方も。

鬼滅のアニメ、後半までちゃんとみてみよ。。

Python 勾配ブースティングにおけるパラメータと調整方法について

また機械学習ネタです。 機械学習の醍醐味である予測モデル作製において勾配ブースティング(Gradient Boosting)について今回は勉強したいと思います。 以前のメルカリ価格予測(http://rautaku.hatenablog.com/entry/2017/12/22/195649/ )でもLightGBMという勾配ブースティングを使ったんですが使ってみたものの結局中身がどういう動きをしているのかだとかパラメーターなどもどういうものがあるのかよく理解できていなかったのと、あまり検索してもパラメータとその調整についてわかりやすく書いているものがあまり見つからなかったので。 このページ(https://www.analyticsvidhya.com/blog/2016/02/complete-guide-parameter-tuning-gradient-boosting-gbm-python/) がわかりやすかったので参考にして実際のデータも使いながら書いていきます。

勾配ブースティングとは

まず勾配ブースティングとは複数の弱学習器を組み合わせたアンサンブル学習の一種で、その中でも1つずつ順番に弱学習器を構築していく手法です。 一つずつ順番作っていくので、新しい弱学習器を作るときには、バギングと呼ばれるすべての弱学習器を独立に学習する方法と比べて並列計算ができず時間がかかります。時間はかかるのですが、過学習(高バリアンス)と学習不足(高バイアス)のバランスをうまく取ることができます。イメージとしては、ブースティングは最適化、バギングは平均化。  

各ステップごとに弱学習器を構築して、前のステップで間違って識別されたものへのウェイトを重くして,次のステップで間違ったものをうまく識別できるようにしていきます。視覚的に見るとこのような感じです。

左から見ていくと、最初の弱学習結果で2個の青丸と6個の赤丸をうまく分けることができていますが、一方で3個の青丸が間違えてしまっています。そこで間違えた3個のプラスを重み付けして次の弱学習器ではこの3個を間違わないように分類して見ます。そうすると次は、赤丸2個と青丸1個が間違えてしまいました。今度はこれらを重み付けして分類してみます。3回目では赤丸2個と青丸1個が間違えてしまっています。

このような作業を繰り返していき、最終的な境界線を描いていきます。今回の絵では3回試行した後での境界線でも青丸1個が間違えてしまっています。実際ぽくていいのではないかなと思います。この辺の最終モデルのちょうど良いバランス(適切なバリアンスとバイアス)を探すために分類回数や数値など適切なパラメータ調整が必要になります。

パラメータについて

GBMの調整できるパラメータはハイパーパラメータと呼ばれるそうですが、このハイパーパラメータは大まかにいうと3種類に分けることができます。 1. Tree-Specific Parameters: ツリー固有パラメータ
ツリーの形状や動作などを規定するパラメータ 2. Boosting Parameters: ブースティングパラメータ
勾配動作などを規定するパラメータ 3. Miscellaneous Parameters: 他パラメータ
上記2枠以外のパラメータ

ツリー固有パラメータ

ツリー固有パラメータは以下のようなものがあります (イメージ図を書いてみましたが間違いがあればご指摘下さい><)

ブースティングパラメータ

ブースティングパラメータは以下のようなものがあります

他パラメータ

他パラメータは以下のようなものがあります。今回はこれらのパラメータの調整は行いませんが参考に。

実データを用いたパラメータ調整

各パラメータを見たところで、実際のデータでのGBMのパラメータ調整の影響を確認していきたいと思います。
今回は参考サイトのと同じ、『インドでの貸付者情報に基づくローン支払い』に関するデータを使います。 生データはこちらから

インドのクレジットカード会社かどこかのデータぽいです。ローンを組んだ人の収入情報や住所、年齢、その他個人情報からその人がちゃんとローンを返したかという情報が入っています。クレジットカード会社側はこのような個人情報から顧客にどのくらい貸付できるかどのくらいで返せるだろうかなど試算しているんですね。怖い怖い。。 一つデータの中に聞き慣れない単語が入っていました。"EMI"という単語ですがこれは、Equated Monthly Installmentの略で、インドでのクレジットカードを利用した定額月分割払いの総称らしいです。言わば月払いのローンのことです。インドはクレジットカードの需要がかなり高いらしく、EMIという制度が大流行したそうです(ちょっと調べていると、借金が膨らみすぎる人が多発したらしく今は禁止されてるような内容もありました。)「EMI Liya Hai To Chukana Padega(EMIで借りた金はちゃんと返せよ)」という映画もあったそうです。ポスターがインパクトありすぎたので貼っておきますね。

データ前処理

今回、前処理については参考サイトですでに行われたものを用います。<a href="https://github.com/aarshayj/Analytics_Vidhya/tree/master/Articles/Parameter_Tuning_GBM_with_Example">こちら(github)からダウンロードします。
前処理前のデータで一つだけ遊んでみました。地域毎の支払率の可視化です。

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import folium
from folium import plugins
%matplotlib inline

train = pd.read_csv("train.csv", encoding="iso-8859-1")
num = 50
city_list = train["City"].value_counts().head(num).index
disbursed_df = pd.DataFrame(index=range(num), columns=["City", "ratio", "lat", "lon"])
disbursed_df["City"] = city_list
disbursed_ratio_list = []
for c in city_list:
    city_count = train.loc[train["City"]==c, "City"].count()
    disbursed_count = train.loc[train["City"]==c, "Disbursed"].sum()
    disbursed_ratio = disbursed_count/city_count
    disbursed_ratio_list.append(disbursed_ratio)
disbursed_df["ratio"] = disbursed_ratio_list

# http://openweathermap.org/help/city_list.txt: get each city's lattitude and longtitude
world_city_df = pd.read_table("http://openweathermap.org/help/city_list.txt", encoding="iso-8859-1")
india_city_df = world_city_df.loc[world_city_df["countryCode"]=="IN"]

def addLatLon(disbursed_df, india_city_df):
    disbursed_df_city = list(disbursed_df["City"])
    disbursed_df_lat = []
    disbursed_df_lon = []
    for d in disbursed_df_city:
        lat_v = india_city_df[india_city_df["nm"]==d]["lat"].values
        lon_v = india_city_df[india_city_df["nm"]==d]["lon"].values
        try:
            lat_v = lat_v[0]
            lon_v = lon_v[0]
        except:
            lat_v = np.nan
            lon_v = np.nan
        disbursed_df_lat.append(lat_v)
        disbursed_df_lon.append(lon_v)
    disbursed_df["lat"] = disbursed_df_lat
    disbursed_df["lon"] = disbursed_df_lon
    return disbursed_df

disbursed_df = addLatLon(disbursed_df, india_city_df)
latlon_df = disbursed_df.dropna().loc[:, ["lat","lon","ratio"]]

# For Folium Heatmap, make tuple list
latlon_list = []
for l in latlon_df.values:
    latlon_list.append(tuple(l))

# Make Heatmap by using Folium
map = folium.Map(location=[disbursed_df["lat"].median(), disbursed_df["lon"].median()], zoom_start=5)
map.add_child(plugins.HeatMap(latlon_list))
map

Foliumというパッケージを使って、地域別のローン支払い率のヒートマップを地図上に可視化してみました。 インドの知識はほとんどないのでこのマッピングを見て何も感じるものがないのですが、都市部が高いとかあるのではないでしょうか。また今度日本のデータ系で何か可視化してみたいです。

まぁ確かにインドのデータだとわかったところで、本題のパラメータ調整とその影響の調査に戻りたいと思います。

まずはデータのインポートやモデルフィット関数の準備です。 modelfit関数の中身は、分類アルゴリズム(今回はGradientBoostingClassifier)を使ってデータを分類して、精度評価として正確度評価(Accuracy)、AUC評価(Area under the ROC curve)、交差検定でスコア評価し、特徴量の影響度グラフと共に出力します。

import pandas as pd
import numpy as np
from sklearn.ensemble import GradientBoostingClassifier  #GBM algorithm
from sklearn import cross_validation, metrics   #Additional scklearn functions
from sklearn.grid_search import GridSearchCV   #Perforing grid search
import matplotlib.pylab as plt
%matplotlib inline
from matplotlib.pylab import rcParams
import seaborn as sns
rcParams['figure.figsize'] = 10, 4

train = pd.read_csv('train_modified.csv')
target = 'Disbursed'
IDcol = 'ID'

def modelfit(alg, dtrain, predictors, performCV=True, printFeatureImportance=True, cv_folds=5):
    #Fit the algorithm on the data
    alg.fit(dtrain[predictors], dtrain['Disbursed'])
        
    #Predict training set:
    dtrain_predictions = alg.predict(dtrain[predictors])
    dtrain_predprob = alg.predict_proba(dtrain[predictors])[:,1]
    
    #Perform cross-validation:
    if performCV:
        cv_score = cross_validation.cross_val_score(alg, dtrain[predictors], dtrain['Disbursed'], cv=cv_folds, scoring='roc_auc')
    
    #Print model report:
    print ("Model Report")
    print ("Accuracy : {:.4f}".format(metrics.accuracy_score(dtrain['Disbursed'].values, dtrain_predictions)))
    print ("AUC Score (Train): {:.4f}".format(metrics.roc_auc_score(dtrain['Disbursed'], dtrain_predprob)))
    
    if performCV:
        print ("CV Score : Mean - {:.6f} | Std - {:.6f} | Min - {:.6f} | Max - {:.6f}".format(np.mean(cv_score),np.std(cv_score),np.min(cv_score),np.max(cv_score)))
        
    #Print Feature Importance:
    if printFeatureImportance:
        feat_imp = pd.Series(alg.feature_importances_, predictors).sort_values(ascending=False)
        sns.set_palette("husl")
        sns.barplot(feat_imp.head(10).index, feat_imp.head(10).values)
        plt.title('Top10 Feature Importances')
        plt.ylabel('Feature Importance Score')
        plt.xticks(rotation=60)
        plt.show()

準備ができたところで、まずは標準値のまま勾配ブースティングをかけて見ます。random_stateは"10"で固定します。

#Choose all predictors except target & IDcols
predictors = [x for x in train.columns if x not in [target, IDcol]]
gbm0 = GradientBoostingClassifier(random_state=10)
modelfit(gbm0, train, predictors)
Model Report
Accuracy : 0.9856
AUC Score (Train): 0.8584
CV Score : Mean - 0.830063 | Std - 0.010523 | Min - 0.815374 | Max - 0.843712

交差検定スコアとして、0.830063を得ました。これをパラメータ調整により上回ることを目標とします。

パラメータ調整アプローチ

前述のように調整すべきパラメータは2種類あります。ツリー固有パラメータとブースティングパラメータです。 learning_rateには最適値というものはなくツリー数(n_estimatiors)が十分にあるのであれば小さければ小さいほど良いですが、あまりに小さすぎると計算負荷が大きく処理に時間がかかってしまいます。

以下の点に気をつけるべきです。
・比較的高めのlearnig_rateを設定する。標準値は0.1だが0.05~0.2あたりを試すべき
・設定したlearnig_rateに対してツリー数を設定する。40〜70くらい。いろんなパターンを試すためにも処理スピードにも気をつける ・learning_rateとツリー数を決めたらツリー固有パラメータを調整する。 ・learning_rateを低くしたなら、より堅牢なモデルを得る為、それに反比例してツリー数を増やす。

ツリー固有パラメータ調整のためのlearning_rateとツリー数(n_estimators)の検証

learning_rateを0.1に固定し、n_estimatorsを20~80まで10個間隔で検証して見ます。

他パラメータは以下で固定します。

min_samples_split = 500 : 総サンプル数の0.5~1%が一般的で、今回のようなクラスに属するサンプル数に偏りがある不均衡データの場合低い方(0.5%)を使います。
min_samples_leaf = 50 : 直感で選びます。過学習を防ぐために使うのと不均衡データなのであまりに小さい数は採用しません。
max_depth = 8 : サンプル数に応じて5~8ぐらいで選択。今回は87000×49と結構大きいので8。
max_features = ‘sqrt’: sqrtで始めるのが一般的である。
subsample = 0.8 : これも一般的に使われる値

GridSearchCVを使ってn_estimatorsを振った際の結果変化を検証します。

predictors = [x for x in train.columns if x not in [target, IDcol]]
param_test1 = {'n_estimators':list(range(20,81,10))}
gsearch1 = GridSearchCV(estimator = GradientBoostingClassifier(learning_rate=0.1, min_samples_split=500,min_samples_leaf=50,max_depth=8,max_features='sqrt',subsample=0.8,random_state=10), 
param_grid = param_test1, scoring='roc_auc',n_jobs=4,iid=False, cv=5)
gsearch1.fit(train[predictors],train[target])
mean_list = []
for s in gsearch1.grid_scores_:
    mean_list.append(s[1])
plt.figure(figsize=(4,2))
plt.plot(list(range(20,81,10)), mean_list)
plt.title("GridSearchCV Score")
plt.xlabel("n_estimators")
plt.show()
gsearch1.grid_scores_, gsearch1.best_params_

([mean: 0.83234, std: 0.01247, params: {'n_estimators': 20},
  mean: 0.83432, std: 0.01280, params: {'n_estimators': 30},
  mean: 0.83564, std: 0.01126, params: {'n_estimators': 40},
  mean: 0.83666, std: 0.01144, params: {'n_estimators': 50},
  mean: 0.83682, std: 0.01193, params: {'n_estimators': 60},
  mean: 0.83612, std: 0.01242, params: {'n_estimators': 70},
  mean: 0.83564, std: 0.01238, params: {'n_estimators': 80}],
 {'n_estimators': 60})

ツリー数60がベストスコアを出しました。この先はn_estimators=60として進めていきます。

from mpl_toolkits.mplot3d import Axes3D
import time
start_time = time.time()
param_test2 = {'max_depth':list(range(5,16,2)), 'min_samples_split':list(range(200,1001,200))}
gsearch2 = GridSearchCV(estimator = GradientBoostingClassifier(learning_rate=0.1, n_estimators=60, max_features='sqrt', subsample=0.8, random_state=10), 
param_grid = param_test2, scoring='roc_auc',n_jobs=4,iid=False, cv=5)
gsearch2.fit(train[predictors],train[target])
x=[]
y=[]
z=[]
for s in gsearch2.grid_scores_:
    x.append(s[0]["max_depth"])
    y.append(s[0]["min_samples_split"])
    z.append(s[1])
fig = plt.figure()
ax = Axes3D(fig)
ax.plot(x, y, z)
plt.show()
print('Showing results took {} secs.'.format(time.time() - start_time))
gsearch2.grid_scores_, gsearch2.best_params_, gsearch2.best_score_

Showing results took 593.9939701557159 secs.





([mean: 0.83094, std: 0.01075, params: {'max_depth': 5, 'min_samples_split': 200},
  mean: 0.83225, std: 0.01140, params: {'max_depth': 5, 'min_samples_split': 400},
  mean: 0.83265, std: 0.01255, params: {'max_depth': 5, 'min_samples_split': 600},
  mean: 0.83239, std: 0.01121, params: {'max_depth': 5, 'min_samples_split': 800},
  mean: 0.83203, std: 0.01215, params: {'max_depth': 5, 'min_samples_split': 1000},
  mean: 0.83260, std: 0.00662, params: {'max_depth': 7, 'min_samples_split': 200},
  mean: 0.83504, std: 0.01166, params: {'max_depth': 7, 'min_samples_split': 400},
  mean: 0.83517, std: 0.00963, params: {'max_depth': 7, 'min_samples_split': 600},
  mean: 0.83528, std: 0.01013, params: {'max_depth': 7, 'min_samples_split': 800},
  mean: 0.83471, std: 0.01230, params: {'max_depth': 7, 'min_samples_split': 1000},
  mean: 0.83094, std: 0.01107, params: {'max_depth': 9, 'min_samples_split': 200},
  mean: 0.83384, std: 0.01082, params: {'max_depth': 9, 'min_samples_split': 400},
  mean: 0.83522, std: 0.01034, params: {'max_depth': 9, 'min_samples_split': 600},
  mean: 0.83252, std: 0.00916, params: {'max_depth': 9, 'min_samples_split': 800},
  mean: 0.83356, std: 0.01169, params: {'max_depth': 9, 'min_samples_split': 1000},
  mean: 0.82438, std: 0.01229, params: {'max_depth': 11, 'min_samples_split': 200},
  mean: 0.82908, std: 0.01067, params: {'max_depth': 11, 'min_samples_split': 400},
  mean: 0.82865, std: 0.01038, params: {'max_depth': 11, 'min_samples_split': 600},
  mean: 0.83004, std: 0.01081, params: {'max_depth': 11, 'min_samples_split': 800},
  mean: 0.83232, std: 0.01079, params: {'max_depth': 11, 'min_samples_split': 1000},
  mean: 0.81861, std: 0.01130, params: {'max_depth': 13, 'min_samples_split': 200},
  mean: 0.82383, std: 0.00630, params: {'max_depth': 13, 'min_samples_split': 400},
  mean: 0.82776, std: 0.00844, params: {'max_depth': 13, 'min_samples_split': 600},
  mean: 0.83018, std: 0.01263, params: {'max_depth': 13, 'min_samples_split': 800},
  mean: 0.82855, std: 0.01021, params: {'max_depth': 13, 'min_samples_split': 1000},
  mean: 0.81289, std: 0.01212, params: {'max_depth': 15, 'min_samples_split': 200},
  mean: 0.82010, std: 0.01238, params: {'max_depth': 15, 'min_samples_split': 400},
  mean: 0.82385, std: 0.01026, params: {'max_depth': 15, 'min_samples_split': 600},
  mean: 0.82805, std: 0.01025, params: {'max_depth': 15, 'min_samples_split': 800},
  mean: 0.82918, std: 0.01099, params: {'max_depth': 15, 'min_samples_split': 1000}],
 {'max_depth': 7, 'min_samples_split': 800},
 0.8352763290658064)

'max_depth': 7, 'min_samples_split': 800がベストスコアでした。どんどん行きましょう、min_samples_leaf、max_featuresを順に確認して行きます。

param_test3 = {'min_samples_leaf':list(range(30,71,10))}
gsearch3 = GridSearchCV(estimator = GradientBoostingClassifier(learning_rate=0.1, n_estimators=60, min_samples_split=800, max_depth=7,max_features='sqrt', subsample=0.8, random_state=10), 
param_grid = param_test3, scoring='roc_auc',n_jobs=4,iid=False, cv=5)
gsearch3.fit(train[predictors],train[target])
mean_list = []
for s in gsearch3.grid_scores_:
    mean_list.append(s[1])
plt.figure(figsize=(4,2))
plt.plot(list(range(30,71,10)), mean_list)
plt.title("GridSearchCV Score")
plt.xlabel("min_samples_leaf")
plt.show()
gsearch3.grid_scores_, gsearch3.best_params_, gsearch3.best_score_

([mean: 0.83565, std: 0.01055, params: {'min_samples_leaf': 30},
  mean: 0.83599, std: 0.01158, params: {'min_samples_leaf': 40},
  mean: 0.83488, std: 0.00898, params: {'min_samples_leaf': 50},
  mean: 0.83634, std: 0.01103, params: {'min_samples_leaf': 60},
  mean: 0.83609, std: 0.00992, params: {'min_samples_leaf': 70}],
 {'min_samples_leaf': 60},
 0.8363437200261575)
param_test4 = {'max_features':list(range(7,20,2))}
gsearch4 = GridSearchCV(estimator = GradientBoostingClassifier(learning_rate=0.1, n_estimators=60,max_depth=7, min_samples_split=800, min_samples_leaf=60, subsample=0.8, random_state=10),
param_grid = param_test4, scoring='roc_auc',n_jobs=4,iid=False, cv=5)
gsearch4.fit(train[predictors],train[target])
mean_list = []
for s in gsearch4.grid_scores_:
    mean_list.append(s[1])
plt.figure(figsize=(4,2))
plt.plot(list(range(7,20,2)), mean_list)
plt.title("GridSearchCV Score")
plt.xlabel("max_features")
plt.show()
gsearch4.grid_scores_, gsearch4.best_params_, gsearch4.best_score_

([mean: 0.83442, std: 0.01179, params: {'max_features': 7},
  mean: 0.83606, std: 0.01068, params: {'max_features': 9},
  mean: 0.83704, std: 0.00725, params: {'max_features': 11},
  mean: 0.83873, std: 0.00966, params: {'max_features': 13},
  mean: 0.83671, std: 0.01033, params: {'max_features': 15},
  mean: 0.83629, std: 0.00860, params: {'max_features': 17},
  mean: 0.83582, std: 0.00988, params: {'max_features': 19}],
 {'max_features': 13},
 0.8387261096220344)

各パラメータのGridSearchCVによる調整した値が算出できました。
n_estimators = 60
max_depth = 7
min_samples_split = 800
min_samples_leaf = 60
max_features = 13
これらをGradientBoostingClassifierのパラメータに設定して再度評価してみましょう。標準値でのCVスコアは0.830063でしたのでこれを上回れば良いです。

#Choose all predictors except target & IDcols
predictors = [x for x in train.columns if x not in [target, IDcol]]
gbm1 = GradientBoostingClassifier(learning_rate=0.1, n_estimators=60,max_depth=7,min_samples_split=800, min_samples_leaf=60, subsample=0.8, random_state=10, max_features=13)
modelfit(gbm1, train, predictors)
Model Report
Accuracy : 0.9854
AUC Score (Train): 0.8871
CV Score : Mean - 0.838726 | Std - 0.009660 | Min - 0.829550 | Max - 0.851014

CVスコア=0.838726でした。無事上回りましたので良かったです!
という形で勾配ブースティングのパラメータ調整の影響を実感することができました。 より良いスコアを出そうとしたらsubsampleやlearning_rateとツリー数のバランス、warm_startなどを調整していくことで、どんどんスコアを上げていくことができると思います。 今回は以上にしておきます!

中古車の維持費を算出するサイト作ってみた

年末年始暇つぶしの2個目。 車って、結構大きな買い物ですよね。だいたいの人が買う前に「この車買ったらどのくらいの毎月負担になるんだろう」って計算すると思います。 計算をする際は、車の本体価格だけでなく、税金関係、保険代、ガソリン代などの要素もあります。 ガソリン代はハイブリッド車だったら安くつくでしょうから燃費は重要です。

またそれなりにいい車だったら手放す時に下取りでそれなりの金額で買い取ってもらえますから、購入金額との差額は自身の負担からは引いてしまって構わないです。

なので中古車などをカーセンサーなどのサイトで探していて、「この車は高いなー」って思っても意外と毎月の負担額は自分の収入でもいけるじゃんってことあるんじゃないかと思います。

そんな計算を、ざっくりですがしてくれるアプリです。

車の情報のサイトはいっぱいあるんですけど、結構バラバラになっていていろんなページを行き来しなくちゃいけなくて面倒だなーと思ったので、 カーセンサーのURLを貼るとその車に関する情報を3つのサイトから収集してきて、計算するようにしました。

以下のサイトを使ってます ・カーセンサー http://www.carsensor.netカーセンサーラボ http://www.carsensorlab.net (買取額を予想するために、その車の販売価格変動タイプを取ってきます。) ・自動車ランニングコスト https://car-life.adg7.com (燃費や重量税などを取得してきます。)

CALCボタンを押して、しばらく待ってください。10秒くらいかかるかもしれないです。。 そしたら以下のような結果が表示されます。

計算パラメータを変更してRECALCボタンを押せば、再計算してくれるので、調べたいパラメータに調整してやってください。

carlcur.herokuapp.com

まぁ車購入検討の際に参考程度に使ってやってください!

最寄り駅を教えてくれるLINEBotを作ってみた -Line Messaging API, Python Flask, Google Maps API-

年末年始の暇つぶしにLineBotを作ってみました。

最寄り駅を教えてくれるBotです。たまに慣れない場所などに行って帰ろうかーと思った時、「ここから一番近い駅ってどこだっけ?」みたいな時があるので、それを教えてくれるBotです。まぁ別にGoogleMapsとかあるしわざわざBotに聞く必要もない内容ですが、対話型で教えてくれた方が楽しいなーってのと、勉強のために作って見ました。

使ってるものは, LINEから提供されているMessagingAPI(https://developers.line.me/ja/docs/messaging-api/overview/), Python用のSDK(https://github.com/line/line-bot-sdk-python), スクリプト処理用のサーバにPython FlaskをHerokuに、現在値からの最寄駅を取得するためにSimpleAPI(http://map.simpleapi.net), 地図や位置情報を取得して地図などで受け取るためにGoogle Maps API(https://developers.google.com/maps/)を使ってみました。 早速ですが、こんな形になりました。

いざ作ってみるとこれだけ作るだけでもいろんな学びがありましたので、同じような初心者の人の助けになればといいかなと思います。 簡単な構成は以下の様になります。

MessagingAPIを介してHerokuへユーザーからのメッセージなどを送ります。Heroku上のFlaskで作ったWebアプリが送られてきたメッセージや位置情報に応じて処理を選択します。位置情報の場合、SimpleAPIへ現在地の緯度経度をPOSTして最寄り駅名を返してもらいます。受けた最寄り駅名をGoogleMapsAPIへ渡すことで、その駅までの距離や所要時間を教えてもらったり、現在地と最寄り駅の場所にピンを打った地図画像を作成してもらいます。 それらをユーザーへ返すことでやり取りします。 ざっくりやったことを自分の備忘録として残しておきます。

作りたいものの流れとしては ・FlaskでMessagingAPIを処理したい ・herokuに置きたい ・最寄り駅を知りたい ・最寄り駅までの徒歩時間と距離を知りたい です。

Flaskでオウム返しLINE botを実装し,herokuにdeployするまで

FlaskでLINE botを実装し,herokuにdeployするまでは、この記事がとても丁寧でわかりやすかったです。 https://qiita.com/suigin/items/0deb9451f45e351acf92
これでオウム返しをしてくれるBotが出来上がります。

メッセージによって返答を分岐する

次にオウム返しではなく、ユーザーからのメッセージに応じて返答をいくつか用意しておき少しだけ会話してるような感を持たせます 今回は特に単純な仕組みで、ユーザーからの返答があらかじめ用意したものと一致すればその返答を返すし、それ以外であればわからないとすぐ弱音を吐くものにします。

「帰るよー」 = 「お疲れ様です。位置情報を送ってください」
「ありがとう」 = 「どういたしまして」
「位置情報を教えて」 = 最寄り駅の位置情報を送る
その他 = 「まだその言葉はわからないです」

こんな単純な数パターンです。こんな感じです。

@handler.add(MessageEvent, message=TextMessage)
def handle_message(event):
    global near_station_name
    global near_station_address
    global near_station_geo_lat
    global near_station_geo_lon

    if event.type == "message":
        if (event.message.text == "帰るよー!") or (event.message.text == "帰るよ!") or (event.message.text == "帰る!") or (event.message.text == "帰るよ"):
            line_bot_api.reply_message(
                event.reply_token,
                [
                    TextSendMessage(text='お疲れ様です'+ chr(0x10002D)),
                    TextSendMessage(text='位置情報を送ってもらうと近くの駅を教えますよ'+ chr(0x10008D)),
                    TextSendMessage(text='line://nv/location'),
                ]
            )
        if (event.message.text == "ありがとう!") or (event.message.text == "ありがとう") or (event.message.text == "ありがと!") or (event.message.text == "ありがと"):
            line_bot_api.reply_message(
                event.reply_token,
                [
                    TextSendMessage(text="どういたしまして!気をつけて帰ってね" + chr(0x100033)),
                ]
            )
        if event.message.text == "位置情報教えて!":
            line_bot_api.reply_message(
                event.reply_token,
                [
                    LocationSendMessage(
                        title=near_station_name,
                        address=near_station_address,
                        latitude=near_station_geo_lat,
                        longitude=near_station_geo_lon
                    ),
                    TextSendMessage(text="タップした後右上のボタンからGoogleMapsなどで開けますよ"+ chr(0x100079)),
                    TextSendMessage(text="もし場所が間違えてたらもう一度地図画像をタップしてみたり位置情報を送り直してみてください"),    
                ]
            )
        else:
            line_bot_api.reply_message(
                event.reply_token,
                [
                    TextSendMessage(text="まだその言葉は教えてもらってないんです"+ chr(0x100029) + chr(0x100098)),
                ]
            )

LINEのSDKがとても簡単で使いやすいです。ユーザーから来たメッセージがテキストメッセージであれば、上記のハンドラー内のスクリプトを実行します。 ちなみに絵文字はhttps://devdocs.line.me/files/emoticon.pdf 内の絵文字が使えUnicodeで与えれているので、chr()メソッドで変更する必要があります。

位置情報を受けた時、SimpleAPIから最寄り駅名を受け取る

次に位置情報を受けた時に位置情報の緯度経度をSimpleAPIに投げます。返ってくるXML内から駅名を抽出します。 画像は東京駅付近の座標を渡した時の帰ってくるXMLです。

lat = event.message.latitude
lon = event.message.longitude
# SimpleAPIから最寄駅リストを取得
near_station_url = 'http://map.simpleapi.net/stationapi?x={}&y={}&output=xml'.format(lon, lat)
near_station_req = urllib.request.Request(near_station_url)
with urllib.request.urlopen(near_station_req) as response:
    near_station_XmlData = response.read()
near_station_root = ET.fromstring(near_station_XmlData)
near_station_list = near_station_root.findall(".//name")
near_station_n = len(near_station_list)

最寄駅名を得ることができたので、次はその駅までの徒歩時間と距離をgoogle maps APIに教えてもらいます。

# 最寄駅名から座標を取得
near_station_geo_url = 'https://maps.googleapis.com/maps/api/place/textsearch/xml?query={}&key={}'.format(urllib.parse.quote_plus(near_station_list[0].text, encoding='utf-8'), google_places_api_key);
near_station_geo_req = urllib.request.Request(near_station_geo_url) 
with urllib.request.urlopen(near_station_geo_req) as response:
    near_station_geo_XmlData = response.read() 
near_station_geo_root = ET.fromstring(near_station_geo_XmlData) 

#最寄駅情報(名前、住所、緯度経度)を取得
near_station_name = near_station_geo_root.findtext(".//name")
near_station_address = near_station_geo_root.findtext(".//formatted_address")
near_station_geo_lat = near_station_geo_root.findtext(".//lat") 
near_station_geo_lon = near_station_geo_root.findtext(".//lng")

#徒歩時間を取得
near_station_direction_url = 'https://maps.googleapis.com/maps/api/directions/xml?origin={},{}&destination={},{}&mode=walking&key={}'.format(lat, lon, near_station_geo_lat, near_station_geo_lon, google_directions_api_key);
near_station_direction_req = urllib.request.Request(near_station_direction_url) 
with urllib.request.urlopen(near_station_direction_req) as response:
    near_station_direction_XmlData = response.read() 
near_station_direction_root = ET.fromstring(near_station_direction_XmlData)
near_station_direction_time_second = int(near_station_direction_root.findtext(".//leg/duration/value"))
near_station_direction_distance_meter = int(near_station_direction_root.findtext(".//leg/distance/value"))
near_station_direction_time_min = near_station_direction_time_second//60
near_station_direction_distance_kilo = near_station_direction_distance_meter//1000 + ((near_station_direction_distance_meter//100)%10)*0.1

google Maps APIはいくつか種類があって、今回は駅名から座標を教えてくれるPlaces API、距離と時間を教えてくれるDirections APIを使います。 各APIはKEYを取得しなければいけないので、アカウントを適当に作って取得します。API KEYはHerokuの環境変数に登録しておきます。

heroku config:set GOOGLE_DIRECTIONS_API_KEY=AI*********************nfwbVTAdc

最後にGoogle Static Maps APIを使って地図に現在地と最寄り駅のマーカーを打って返してもらいます。
map_image_url = 'https://maps.googleapis.com/maps/api/staticmap?size=520x520&scale=2&maptype=roadmap&key={}'.format(google_staticmaps_api_key);
map_image_url += '&markers=color:{}|label:{}|{},{}'.format('red', '', near_station_geo_lat, near_station_geo_lon)
map_image_url += '&markers=color:{}|label:{}|{},{}'.format('blue', '', lat, lon)

こんな感じで現在地と目的地の地図を返してもらいます。
以上を組み合わせることで位置情報を受け取った時のハンドラーのコードは以下のようになります。

@handler.add(MessageEvent, message=LocationMessage)
def handle_location(event):
    global near_station_name
    global near_station_address
    global near_station_geo_lat
    global near_station_geo_lon
    
    lat = event.message.latitude
    lon = event.message.longitude

    zoomlevel = 18
    imagesize = 1040

    # SimpleAPIから最寄駅リストを取得
    near_station_url = 'http://map.simpleapi.net/stationapi?x={}&y={}&output=xml'.format(lon, lat)
    near_station_req = urllib.request.Request(near_station_url)
    with urllib.request.urlopen(near_station_req) as response:
        near_station_XmlData = response.read()
    near_station_root = ET.fromstring(near_station_XmlData)
    near_station_list = near_station_root.findall(".//name")
    near_station_n = len(near_station_list)

    # 最寄駅名から座標を取得
    near_station_geo_url = 'https://maps.googleapis.com/maps/api/place/textsearch/xml?query={}&key={}'.format(urllib.parse.quote_plus(near_station_list[0].text, encoding='utf-8'), google_places_api_key);
    near_station_geo_req = urllib.request.Request(near_station_geo_url) 
    with urllib.request.urlopen(near_station_geo_req) as response:
        near_station_geo_XmlData = response.read() 
    near_station_geo_root = ET.fromstring(near_station_geo_XmlData) 
    
    #最寄駅情報(名前、住所、緯度経度)を取得
    near_station_name = near_station_geo_root.findtext(".//name")
    near_station_address = near_station_geo_root.findtext(".//formatted_address")
    near_station_geo_lat = near_station_geo_root.findtext(".//lat") 
    near_station_geo_lon = near_station_geo_root.findtext(".//lng")

    #徒歩時間を取得
    near_station_direction_url = 'https://maps.googleapis.com/maps/api/directions/xml?origin={},{}&destination={},{}&mode=walking&key={}'.format(lat, lon, near_station_geo_lat, near_station_geo_lon, google_directions_api_key);
    near_station_direction_req = urllib.request.Request(near_station_direction_url) 
    with urllib.request.urlopen(near_station_direction_req) as response:
        near_station_direction_XmlData = response.read() 
    near_station_direction_root = ET.fromstring(near_station_direction_XmlData)
    near_station_direction_time_second = int(near_station_direction_root.findtext(".//leg/duration/value"))
    near_station_direction_distance_meter = int(near_station_direction_root.findtext(".//leg/distance/value"))
    near_station_direction_time_min = near_station_direction_time_second//60
    near_station_direction_distance_kilo = near_station_direction_distance_meter//1000 + ((near_station_direction_distance_meter//100)%10)*0.1


    map_image_url = 'https://maps.googleapis.com/maps/api/staticmap?size=520x520&scale=2&maptype=roadmap&key={}'.format(google_staticmaps_api_key);
    map_image_url += '&markers=color:{}|label:{}|{},{}'.format('red', '', near_station_geo_lat, near_station_geo_lon)
    map_image_url += '&markers=color:{}|label:{}|{},{}'.format('blue', '', lat, lon)

    actions = [
        MessageImagemapAction(
            text = "位置情報教えて!",
            area = ImagemapArea(
                x = 0,
                y = 0,
                width = 1040,
                height = 1040
        )
    )]

    line_bot_api.reply_message(
        event.reply_token,
        [
            ImagemapSendMessage(
                base_url = 'https://{}/imagemap/{}'.format(request.host, urllib.parse.quote_plus(map_image_url)),
                alt_text = '地図',
                base_size = BaseSize(height=imagesize, width=imagesize),
                actions = actions,
            ),
            TextSendMessage(text=near_station_list[0].text + 'が一番近いですね!'),
            TextSendMessage(text='歩いて約' + str(near_station_direction_time_min) + '分。距離は約'+ str(near_station_direction_distance_kilo) + 'kmです。'),
            TextSendMessage(text='画像をタップすれば位置情報を送ります'),
        ]
    )

書かなかったところとしては、地図に「位置情報教えて!」というテキスト情報を持たせることで、地図を押せばテキストを返すようにしていることや、取得した最寄り駅情報をグローバル変数とすることでテキストハンドラーでも代入された値を使用できるようにしています。 またリッチメニューを使って「帰るよ」のコマンド送信を簡易化してみました。リッチメニューはLineBusinessアカウントの設定画面から設定することができます。 かなりざっくりになりましたが、まとめとしては以上にしておきます。 また以下のQRコードを読み取って友達になってもらえれば誰でも試してもらうことができます。(やりとり内容や誰が友達かなど僕には全くわかりませんので特に気にすることはないと思います。)

最後に全体のコードを貼っておきます。またgithubにも置いておきます。 github.com

import os
import sys
from flask import Flask, request, abort, send_file

from linebot import (
    LineBotApi, WebhookHandler
)
from linebot.exceptions import (
    InvalidSignatureError
)
from linebot.models import (
    MessageEvent, TextMessage, TextSendMessage, LocationMessage, MessageImagemapAction, ImagemapArea, ImagemapSendMessage, BaseSize, LocationSendMessage
)
from PIL import Image
import requests
from io import BytesIO, StringIO
import urllib.parse
import urllib.request
import xml.etree.ElementTree as ET

app = Flask(__name__)

# 環境変数から各種KEYを取得
channel_secret = os.environ['LINE_CHANNEL_SECRET']
channel_access_token = os.environ['LINE_CHANNEL_ACCESS_TOKEN']
google_places_api_key = os.environ['GOOGLE_PLACES_API_KEY']
google_directions_api_key = os.environ['GOOGLE_DIRECTIONS_API_KEY']
google_staticmaps_api_key = os.environ['GOOGLE_STATICMAPS_API_KEY']

if channel_secret is None:
    print('Specify LINE_CHANNEL_SECRET as environment variable.')
    sys.exit(1)
if channel_access_token is None:
    print('Specify LINE_CHANNEL_ACCESS_TOKEN as environment variable.')
    sys.exit(1)

line_bot_api = LineBotApi(channel_access_token)
handler = WebhookHandler(channel_secret)

near_station_name = "東京駅"
near_station_address = "日本、〒100-0005 東京都千代田区丸の内1丁目"
near_station_geo_lat = 35.6811673
near_station_geo_lon = 139.7670516

@app.route("/")
def hello_world():
    return "hello world!"

@app.route("/callback", methods=['POST'])
def callback():
    # get X-Line-Signature header value
    signature = request.headers['X-Line-Signature']

    # get request body as text
    body = request.get_data(as_text=True)
    app.logger.info("Request body: " + body)

    # handle webhook body
    try:
        handler.handle(body, signature)
    except InvalidSignatureError:
        abort(400)

    return 'OK'

@app.route("/imagemap/<path:url>/<size>")
def imagemap(url, size):
    map_image_url = urllib.parse.unquote(url)
    response = requests.get(map_image_url)
    img = Image.open(BytesIO(response.content))
    img_resize = img.resize((int(size), int(size)))
    byte_io = BytesIO()
    img_resize.save(byte_io, 'PNG')
    byte_io.seek(0)
    return send_file(byte_io, mimetype='image/png')


@handler.add(MessageEvent, message=TextMessage)
def handle_message(event):
    global near_station_name
    global near_station_address
    global near_station_geo_lat
    global near_station_geo_lon

    if event.type == "message":
        if (event.message.text == "帰るよー!") or (event.message.text == "帰るよ!") or (event.message.text == "帰る!") or (event.message.text == "帰るよ"):
            line_bot_api.reply_message(
                event.reply_token,
                [
                    TextSendMessage(text='お疲れ様です'+ chr(0x10002D)),
                    TextSendMessage(text='位置情報を送ってもらうと近くの駅を教えますよ'+ chr(0x10008D)),
                    TextSendMessage(text='line://nv/location'),
                ]
            )
        if (event.message.text == "ありがとう!") or (event.message.text == "ありがとう") or (event.message.text == "ありがと!") or (event.message.text == "ありがと"):
            line_bot_api.reply_message(
                event.reply_token,
                [
                    TextSendMessage(text="どういたしまして!気をつけて帰ってね" + chr(0x100033)),
                ]
            )
        if event.message.text == "位置情報教えて!":
            line_bot_api.reply_message(
                event.reply_token,
                [
                    LocationSendMessage(
                        title=near_station_name,
                        address=near_station_address,
                        latitude=near_station_geo_lat,
                        longitude=near_station_geo_lon
                    ),
                    TextSendMessage(text="タップした後右上のボタンからGoogleMapsなどで開けますよ"+ chr(0x100079)),
                    TextSendMessage(text="もし場所が間違えてたらもう一度地図画像をタップしてみたり位置情報を送り直してみてください"),    
                ]
            )
        else:
            line_bot_api.reply_message(
                event.reply_token,
                [
                    TextSendMessage(text="まだその言葉は教えてもらってないんです"+ chr(0x100029) + chr(0x100098)),
                ]
            )

@handler.add(MessageEvent, message=LocationMessage)
def handle_location(event):
    global near_station_name
    global near_station_address
    global near_station_geo_lat
    global near_station_geo_lon
    
    lat = event.message.latitude
    lon = event.message.longitude

    zoomlevel = 18
    imagesize = 1040

    # SimpleAPIから最寄駅リストを取得
    near_station_url = 'http://map.simpleapi.net/stationapi?x={}&y={}&output=xml'.format(lon, lat)
    near_station_req = urllib.request.Request(near_station_url)
    with urllib.request.urlopen(near_station_req) as response:
        near_station_XmlData = response.read()
    near_station_root = ET.fromstring(near_station_XmlData)
    near_station_list = near_station_root.findall(".//name")
    near_station_n = len(near_station_list)

    # 最寄駅名から座標を取得
    near_station_geo_url = 'https://maps.googleapis.com/maps/api/place/textsearch/xml?query={}&key={}'.format(urllib.parse.quote_plus(near_station_list[0].text, encoding='utf-8'), google_places_api_key);
    near_station_geo_req = urllib.request.Request(near_station_geo_url) 
    with urllib.request.urlopen(near_station_geo_req) as response:
        near_station_geo_XmlData = response.read() 
    near_station_geo_root = ET.fromstring(near_station_geo_XmlData) 
    
    #最寄駅情報(名前、住所、緯度経度)を取得
    near_station_name = near_station_geo_root.findtext(".//name")
    near_station_address = near_station_geo_root.findtext(".//formatted_address")
    near_station_geo_lat = near_station_geo_root.findtext(".//lat") 
    near_station_geo_lon = near_station_geo_root.findtext(".//lng")

    #徒歩時間を取得
    near_station_direction_url = 'https://maps.googleapis.com/maps/api/directions/xml?origin={},{}&destination={},{}&mode=walking&key={}'.format(lat, lon, near_station_geo_lat, near_station_geo_lon, google_directions_api_key);
    near_station_direction_req = urllib.request.Request(near_station_direction_url) 
    with urllib.request.urlopen(near_station_direction_req) as response:
        near_station_direction_XmlData = response.read() 
    near_station_direction_root = ET.fromstring(near_station_direction_XmlData)
    near_station_direction_time_second = int(near_station_direction_root.findtext(".//leg/duration/value"))
    near_station_direction_distance_meter = int(near_station_direction_root.findtext(".//leg/distance/value"))
    near_station_direction_time_min = near_station_direction_time_second//60
    near_station_direction_distance_kilo = near_station_direction_distance_meter//1000 + ((near_station_direction_distance_meter//100)%10)*0.1


    map_image_url = 'https://maps.googleapis.com/maps/api/staticmap?size=520x520&scale=2&maptype=roadmap&key={}'.format(google_staticmaps_api_key);
    map_image_url += '&markers=color:{}|label:{}|{},{}'.format('red', '', near_station_geo_lat, near_station_geo_lon)
    map_image_url += '&markers=color:{}|label:{}|{},{}'.format('blue', '', lat, lon)

    actions = [
        MessageImagemapAction(
            text = "位置情報教えて!",
            area = ImagemapArea(
                x = 0,
                y = 0,
                width = 1040,
                height = 1040
        )
    )]

    line_bot_api.reply_message(
        event.reply_token,
        [
            ImagemapSendMessage(
                base_url = 'https://{}/imagemap/{}'.format(request.host, urllib.parse.quote_plus(map_image_url)),
                alt_text = '地図',
                base_size = BaseSize(height=imagesize, width=imagesize),
                actions = actions,
            ),
            TextSendMessage(text=near_station_list[0].text + 'が一番近いですね!'),
            TextSendMessage(text='歩いて約' + str(near_station_direction_time_min) + '分。距離は約'+ str(near_station_direction_distance_kilo) + 'kmです。'),
            TextSendMessage(text='画像をタップすれば位置情報を送ります'),
        ]
    )


if __name__ == "__main__":
    app.run()

Mercari Price Challenge -機械学習を使ったメルカリの価格予測 Ridge回帰 LightGBM

ポケモンデータ解析に続いて、またKaggleでのデータ解析ネタです。 今回の解析テーマはフリマアプリのメルカリです。提供されているのはアメリカで行われた商品毎の取引データです。データ内には出品者側が登録した商品名、商品コンディションランク、商品カテゴリ、ブランド名、価格、配送料負担サイド、商品概要が含まれています。自分もメルカリをたまに使用しますが、ユーザはカテゴリ欄や商品概要などは自由に選択したり書いているので、統率の取れているデータになっていない部分もあります。 現在賞金のかかったコンペが行われていて、一番精度の高い価格予測モデルを作成した人は10万ドル(約1,000万円)が貰えるそうです! まぁレベルが高すぎて、僕には関係ない話なのでコンペは無視してデータ解析の勉強をしていきたいと思います。 参加者の賢い人たちがコードを公開してくれているので、参考にしながら機械学習の方法を学んでいきます。 こちらの方を参考にします。(https://www.kaggle.com/iamprateek/submission-to-mercari-price-suggestion-challenge/notebook)
今回はけっこマニアックな内容になってしまいました。。出てくる主な学習ツールは、Ridge回帰、LightGBMです。

さてメルカリから提供されているデータは148万件もあります。今回148万件ものデータは処理するのが大変で、僕のPCでは処理計算に1時間ほどかかってしまうこともあり今回はサクサク進めていきたいので1万件のみを扱います。

まずはデータをざっくり俯瞰します。

import pandas as pd
train = pd.read_csv("train.csv", nrows=10000)
print ("train.shape: " + str(train.shape))
train.head()
train.shape: (10000, 8)

データにあるのはこんな情報ですね。以下写真の赤枠ですね。写真の取り方とか出品者のいいね数なども実際は影響しそうですが、今回はこれらの情報から価格が予想できるような学習モデルを導くのが目的です。

import matplotlib.pyplot as plt
%matplotlib inline
plt.figure(figsize=(10, 6))
plt.hist(train["price"])
plt.xlabel("price[$]")
plt.ylabel("count")
plt.title("Price Histogram")
plt.show()

ほとんどが200$以内の価格帯で取引されてますね。もうちょい拡大しますか。

plt.figure(figsize=(10, 6))
plt.hist(train["price"], bins=500)
plt.xlim(0, 200)
plt.xlabel("price[$]")
plt.ylabel("count")
plt.title("Price Histogram")
plt.show()

10~20$くらいが一番頻度の高い取引みたいです。1000~2000円程度ですので、実感的にも確かにそうだろうなって感じです。

次はカテゴリ別の価格分布を見てみたいです。カテゴリ名欄を分割して新カラムを作成し視覚化してみます。

import time
import seaborn as sns
def split_cat(text):
    try: return text.split("/")
    except: return ("No Label", "No Label", "No Label")
train["general_cat"], train["sub_cat1"], train["sub_cat2"] = \
    zip(*train['category_name'].apply(lambda x: split_cat(x)))
train.drop("category_name", axis=1)
start_time = time.time()
plt.figure(figsize=(16, 8))
ax = sns.violinplot(x="general_cat", y="price", data=train, inner=None)
ax = sns.swarmplot(x="general_cat", y="price", data=train, edgecolor="gray", hue="sub_cat1")
plt.xticks(rotation=30)
plt.ylim(0, 200)
plt.legend(loc=9, bbox_to_anchor=(0.5, -0.25), ncol=5)
plt.show()
print('Showing graph took {} secs.'.format(time.time() - start_time))

Showing graph took 1226.6749708652496 secs.

グラフを表示するだけに20分以上もかかってしまいました。。グラフをざっくりと考察してみると、キッズ用品や美容系は安値側にどっしりしたバイオリン図になっており比較的安価層に多く偏っていそうなことだったり、逆にElectronicsは安値から高値までバラついてそうなことがわかりまs。また安価帯の値段設定は, 10ドル、25ドル、40ドル、75ドルあたりなどに分布されてるものが多く、あまり連続的な分布ではなく階段状の分布がありそうですね。まぁ日本でも999円, 1999円, 9999円とか設定されやすい価格ってありますよね。そんな感じかなという印象です。

ちなみに一番高い価格でやり取りされた商品はどんな商品だったんでしょうか??見てみましょう。

train.sort_values(by="price", ascending=False).head(1)

Chanel Classic Flag Bag medium Caviar L 1506$
1500ドルですから約15万。まぁそこまで高くはなかったです。 シャネルのショルダーバッグみたいです。検索してみたところこれですかね

データ俯瞰はこれぐらいにして、本題である機械学習ツールを使った価格学習に入っていきます。 まずは機械学習のお決まり、学習器がちゃんと処理できるようにデータ整理です。

空欄になっているカラムについては"missing"という値を入れます。

def handle_missing_inplace(dataset):
    dataset['general_cat'].fillna(value='missing', inplace=True)
    dataset['sub_cat1'].fillna(value='missing', inplace=True)
    dataset['sub_cat2'].fillna(value='missing', inplace=True)
    dataset['brand_name'].fillna(value='missing', inplace=True)
    dataset['item_description'].fillna(value='missing', inplace=True)
handle_missing_inplace(train)
train['brand_name'].value_counts().head()
missing              4261
Nike                  350
PINK                  348
Victoria's Secret     322
LuLaRoe               201
Name: brand_name, dtype: int64

ブランド名のところは空白が多くて40%以上ものの欄にmissingが入りました。
次はブランド名やカテゴリ名なども出現頻度が極端に少ないものは学習しても仕方ないので、空欄と同じ扱いのmissing値に置換してしまいます。

def cutting(dataset):
    pop_brand = dataset['brand_name'].value_counts().loc[lambda x: x.index != 'missing'].index[:750]
    dataset.loc[~dataset['brand_name'].isin(pop_brand), 'brand_name'] = 'missing'
    pop_category1 = dataset['general_cat'].value_counts().loc[lambda x: x.index != 'missing'].index[:450]
    pop_category2 = dataset['sub_cat1'].value_counts().loc[lambda x: x.index != 'missing'].index[:450]
    pop_category3 = dataset['sub_cat2'].value_counts().loc[lambda x: x.index != 'missing'].index[:450]
    dataset.loc[~dataset['general_cat'].isin(pop_category1), 'general_cat'] = 'missing'
    dataset.loc[~dataset['sub_cat1'].isin(pop_category2), 'sub_cat1'] = 'missing'
    dataset.loc[~dataset['sub_cat2'].isin(pop_category3), 'sub_cat2'] = 'missing'
cutting(train)
train['brand_name'].value_counts().head()
missing              4272
Nike                  350
PINK                  348
Victoria's Secret     322
LuLaRoe               201
Name: brand_name, dtype: int64

ちょっとmissingが増えました。(参考コードではデータ量が多いのでこの作業の効力は大きいかもしれませんが、今回はデータ量はそれほど多くないので必要性はないかもしれません)
次はカテゴリ名と商品概要はテキストデータが入っており、次段階の解析で単語抽出を行いたいためデータ型をカテゴリ型に変換しておきます。

def to_categorical(dataset):
    dataset['general_cat'] = dataset['general_cat'].astype('category')
    dataset['sub_cat1'] = dataset['sub_cat1'].astype('category')
    dataset['sub_cat2'] = dataset['sub_cat2'].astype('category')
    dataset['item_condition_id'] = dataset['item_condition_id'].astype('category')
to_categorical(train)
train.dtypes
train_id                int64
name                   object
item_condition_id    category
category_name          object
brand_name             object
price                 float64
shipping                int64
item_description       object
general_cat          category
sub_cat1             category
sub_cat2             category
dtype: object

データの整理は終わりました。ここから少し難しくなります。
今回のデータ処理に対して回帰や分類をしようとする際に障壁となるのが、文字などの自然言語データです。これらを数学的に処理するのがひとつの壁になります。特に今回はユーザーがそれぞれの感覚で自由に書いているような文字列なので大変です。 このような自然言語データには単語の組み合わせやその傾向を見て、あぁこういうことを言ってるんだなとコンピュータに気づかせてあげるようにしなければなりません。そこでそういった処理に向いている便利な各種機械学習ライブラリを使っていきます。

まずは商品名に対してです。Scikit-learnという機械学習ライブラリにCountVectorizerというものがあり、これは単語の出現頻度を出すために使うようなクラスでこれを使ってみたいと思います。具体的には、fit関数に渡された文字列データに含まれる単語を単語毎にその出現回数をカウントして、それをベクトルとして表現してくれるクラスだそうです。 パラメータの min_dfは整数で設定し、学習データ全体で現れたサンプル数がmin_df未満の単語については特徴ベクトルからは除外する、というものです。

import numpy as np
from sklearn.feature_extraction.text import CountVectorizer
cv = CountVectorizer(min_df=5)
X_name = cv.fit_transform(train['name'])
print (X_name.shape)
occ = np.asarray(X_name.sum(axis=0)).ravel().tolist()
counts_cv = pd.DataFrame({'term': cv.get_feature_names(), 'occurrences': occ})
counts_cv.sort_values(by='occurrences', ascending=False).head(10)
(10000, 1508)

上記Outputは10,000個のテキストデータに対して1,508個のユニークな単語が含まれていたこと、出現数は"pink"という単語が563回で1番多く出てきたことを示します。5番目に多く出てきてる"lularoe"というのが聞きなれなかったので調べてみると女性用のファッションブランドぽいですね。このブランドの品が多く出品されているんだなとわかります。こういうものを使えば今流行っている服はなんだとかトレンドなんかも見えてきそうです。

同様にカテゴリ名に対しても同じ処理を行います。

cv = CountVectorizer(min_df=5)
combine_category = [train["general_cat"], train["sub_cat1"], train["sub_cat2"]]
X_category1 = cv.fit_transform(train['general_cat'])
print ("----general_cat----")
print (X_category1.shape)
occ = np.asarray(X_category1.sum(axis=0)).ravel().tolist()
counts_cv = pd.DataFrame({'term': cv.get_feature_names(), 'occurrences': occ})
print (counts_cv.sort_values(by='occurrences', ascending=False).head())
X_category2 = cv.fit_transform(train['sub_cat1'])
print ("----sub_cat1----")
print (X_category2.shape)
occ = np.asarray(X_category2.sum(axis=0)).ravel().tolist()
counts_cv = pd.DataFrame({'term': cv.get_feature_names(), 'occurrences': occ})
print (counts_cv.sort_values(by='occurrences', ascending=False).head())
X_category3 = cv.fit_transform(train['sub_cat2'])
print ("----sub_cat2----")
print (X_category3.shape)
occ = np.asarray(X_category3.sum(axis=0)).ravel().tolist()
counts_cv = pd.DataFrame({'term': cv.get_feature_names(), 'occurrences': occ})
print (counts_cv.sort_values(by='occurrences', ascending=False).head())
----general_cat----
(10000, 14)
    occurrences         term
13         4381        women
0          1508       beauty
5          1159         kids
2           816  electronics
7           642          men
----sub_cat1----
(10000, 109)
     occurrences         term
62           918       makeup
3            911  accessories
6            870      apparel
100          850         tops
9            849     athletic
----sub_cat2----
(10000, 317)
     occurrences         term
248          634       shirts
168          442     leggings
206          439        pants
2            438  accessories
288          396       tights

次は、商品概要が記入されているテキストデータの処理を行います。このデータは商品名やカテゴリ名とは違って多くの単語で構成される文章です。CountVectorizerは単語の出現数をカウントするのみでしたが、TfidfVectorizerというクラスを使えば各単語の出現頻度と希少性を計算して”重要そうな”単語を見つけ出すことができます。パラメータのmax_featuresは tfidf値を表示する単語数を指定します。tfidf値が高い順にとっていきます。ngramは複数単語の連続にも対応できます。

from sklearn.feature_extraction.text import TfidfVectorizer
tv = TfidfVectorizer(max_features=1000,
                         ngram_range=(1, 3),
                         stop_words='english')
X_description = tv.fit_transform(train['item_description'])
weights = np.asarray(X_description.mean(axis=0)).ravel().tolist()
weights_df = pd.DataFrame({'term': tv.get_feature_names(), 'weight': weights})
weights_df.sort_values(by='weight', ascending=False).head(10)

最もtfidf値が高いものが"description"となっていますが、これは"Not description yet"と"概要記載なし"のものが影響してそうです。あとは"新しい"だったり"中古"であったりという単語が重み付けされています。

次にブランド名の処理です。ブランド名はある程度限られたものしかないので各ブランドに該当するかどうかを0か1かでラベリングしていく方法が向いてそうです。例えばブランド名がNikeであれば1、Nikeでなければ0みたいな感じですね。このようなラベリングにはLabelBinarizerを使用します。

from sklearn.preprocessing import LabelBinarizer
lb = LabelBinarizer(sparse_output=True)
X_brand = lb.fit_transform(train['brand_name'])
occ = np.asarray(X_brand.sum(axis=0)).ravel().tolist()
counts_cv = pd.DataFrame({'term': lb.classes_, 'occurrences': occ})
counts_cv.sort_values(by='occurrences', ascending=False).head(10)

あとは商品コンディションと配送料負担ですがこれらはすでにラベリングがされているので、ダミー変数を使いましょう。ダミー変数については統計学用語ですのでわからない人は調べてください。またダミー変数を作成すると大量の疎行列(要素のほとんどが0の行列)が発生するためcsr_matrixという疎行列を効率よく処理できるクラスを使用します。(csr_matrixかcsc_matrixかどちらを使うべきかは、後工程で行に対する処理を行うか列に対する処理を行うかで使い分けるべきらしいです。今回は後工程で各行に対して回帰処理を行う予定なのでcsrを使います。)

from scipy.sparse import csr_matrix
X_dummies = csr_matrix(pd.get_dummies(train[['item_condition_id', 'shipping']],
                                          sparse=True).values)
X_dummies
<10000x6 sparse matrix of type '<class 'numpy.int64'>'
    with 14539 stored elements in Compressed Sparse Row format>

以上で10,000件のデータを行列群にすることができました。各10,000件の商品名、カテゴリ名を構成する単語出現数行列、商品概要を構成する単語の重み行列、ブランド名のラベリング、商品コンディション・配送料負担のダミー変数行列がありますのでこれらを学習器に入れていきましょう。それぞれ疎行列ですのでcsrで扱います。

from scipy.sparse import hstack
sparse_merge = hstack([X_dummies, X_description, X_brand, X_category1, X_category2, X_category3, X_name]).tocsr()
sparse_merge.shape
(10000, 3705)

3705個の因子を持ったデータフレームが出来ましたので、これらを回帰・分類学習器にかけることで各因子から予測して価格を算出してみます。ここでのミソとしては複数の学習器の結果を絶妙な割合で掛け合わせることで結果の精度を調整することです。今回使う学習器はRidge回帰とLightGBMという今流行りの勾配ブースティング学習器を使います。掛け合わせ割合は、Ridge: 0.3, LightGBM: 0.35 * 2個です。10,000件のデータは9,000個訓練用、1,000個検証用に分割します。

それぞれの学習器単体での対数平方平均二乗誤差(RMSLE: Root Mean Squared Logarithmic Error)を検証しながら進めていきます。

from sklearn.linear_model import Ridge
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_squared_error

X = sparse_merge
y = np.log1p(train["price"])
train_X, test_X, train_y, test_y = train_test_split(X, y, test_size = 0.1, random_state = 144) 

modelR = Ridge(alpha=.5, copy_X=True, fit_intercept=True, max_iter=100,
      normalize=False, random_state=101, solver='auto', tol=0.01)
modelR.fit(train_X, train_y)
predsR = modelR.predict(test_X)

def rmsle(y, y0):
     assert len(y) == len(y0)
     return np.sqrt(np.mean(np.power(np.log1p(y)-np.log1p(y0), 2)))

rmsleR = rmsle(predsR, test_y)
print ("Ridge Regression RMSLE = " + str(rmsleR))
Ridge Regression RMSLE = 0.150454850006
import lightgbm as lgb

train_XL1, valid_XL1, train_yL1, valid_yL1 = train_test_split(train_X, train_y, test_size = 0.1, random_state = 144) 
d_trainL1 = lgb.Dataset(train_XL1, label=train_yL1, max_bin=8192)
d_validL1 = lgb.Dataset(valid_XL1, label=valid_yL1, max_bin=8192)
watchlistL1 = [d_trainL1, d_validL1]
paramsL1 = {
        'learning_rate': 0.65,
        'application': 'regression',
        'max_depth': 3,
        'num_leaves': 60,
        'verbosity': -1,
        'metric': 'RMSE',
        'data_random_seed': 1,
        'bagging_fraction': 0.5,
        'nthread': 4
    }
modelL1 = lgb.train(paramsL1, train_set=d_trainL1, num_boost_round=8000, valid_sets=watchlistL1, \
early_stopping_rounds=5000, verbose_eval=500) 
predsL1 = modelL1.predict(test_X)
rmsleL1 = rmsle(predsL1, test_y)
print ("LightGBM1 RMSLE = " + str(rmsleL1))
Training until validation scores don't improve for 5000 rounds.
[500]   training's rmse: 0.411212   valid_1's rmse: 0.609769
[1000]  training's rmse: 0.35107    valid_1's rmse: 0.628343
[1500]  training's rmse: 0.31196    valid_1's rmse: 0.634885
[2000]  training's rmse: 0.283858   valid_1's rmse: 0.644041
[2500]  training's rmse: 0.260685   valid_1's rmse: 0.650578
[3000]  training's rmse: 0.243261   valid_1's rmse: 0.657922
[3500]  training's rmse: 0.229662   valid_1's rmse: 0.660836
[4000]  training's rmse: 0.217215   valid_1's rmse: 0.663572
[4500]  training's rmse: 0.206279   valid_1's rmse: 0.666011
[5000]  training's rmse: 0.195771   valid_1's rmse: 0.665814
Early stopping, best iteration is:
[96]    training's rmse: 0.528535   valid_1's rmse: 0.590428
LightGBM1 RMSLE = 0.153722873648
train_XL2, valid_XL2, train_yL2, valid_yL2 = train_test_split(train_X, train_y, test_size = 0.1, random_state = 101) 
d_trainL2 = lgb.Dataset(train_XL2, label=train_yL2, max_bin=8192)
d_validL2 = lgb.Dataset(valid_XL2, label=valid_yL2, max_bin=8192)
watchlistL2 = [d_trainL2, d_validL2]
paramsL2 = {
        'learning_rate': 0.85,
        'application': 'regression',
        'max_depth': 3,
        'num_leaves': 140,
        'verbosity': -1,
        'metric': 'RMSE',
        'data_random_seed': 2,
        'bagging_fraction': 1,
        'nthread': 4
    }
modelL2 = lgb.train(paramsL2, train_set=d_trainL2, num_boost_round=5500, valid_sets=watchlistL2, \
early_stopping_rounds=5000, verbose_eval=500) 
predsL2 = modelL2.predict(test_X)
rmsleL2 = rmsle(predsL2, test_y)
print ("LightGBM2 RMSLE = " + str(rmsleL2))
Training until validation scores don't improve for 5000 rounds.
[500]   training's rmse: 0.3948 valid_1's rmse: 0.660834
[1000]  training's rmse: 0.32773    valid_1's rmse: 0.677107
[1500]  training's rmse: 0.285268   valid_1's rmse: 0.687161
[2000]  training's rmse: 0.253586   valid_1's rmse: 0.693729
[2500]  training's rmse: 0.23149    valid_1's rmse: 0.702132
[3000]  training's rmse: 0.215441   valid_1's rmse: 0.707667
[3500]  training's rmse: 0.201827   valid_1's rmse: 0.710812
[4000]  training's rmse: 0.188848   valid_1's rmse: 0.713719
[4500]  training's rmse: 0.177358   valid_1's rmse: 0.715694
[5000]  training's rmse: 0.168306   valid_1's rmse: 0.71957
Early stopping, best iteration is:
[26]    training's rmse: 0.578606   valid_1's rmse: 0.62066
LightGBM2 RMSLE = 0.159046033321

各学習器単体での結果は、Ridge回帰: 0.1504、LightGBM1回目: 0.1524、LightGBM2回目: 0.1587でした。 これを組み合わせることで以下の様な結果が得られます。

preds = predsR*0.3 + predsL1*0.35 + predsL2*0.35
rmsle = rmsle(preds, test_y)
print ("Total RMSLE = " + str(rmsle))
Total RMSLE = 0.143722891253

掛け合わせることで、0.1435までスコアを上げることができました!ここからさらにスコアを上げていくには、各学習器のパラメータや学習器同士の割合をいろいろ調整していくことが必要になります。 では最後に実際の価格と予測価格の分布を見て終わります。

actual_price = np.expm1(test_y)
preds_price = np.expm1(preds)

plt.figure(figsize=(12,10))
cm = plt.cm.get_cmap('winter')
x_diff = np.clip(100 * ((preds_price - actual_price) / actual_price), -75, 75)
plt.scatter(x=actual_price, y=preds_price, c=x_diff, s=10, cmap=cm)
plt.colorbar()
plt.plot([0, 100], [0, 100], 'k--', lw=1)
plt.xlim(0, 100)
plt.ylim(0, 100)
plt.title('Actual vs. Predicted Prices')
plt.xlabel('Actual Prices [$]')
plt.ylabel('Predicted Prices [$]')
plt.show()

なんか今回はマニアックであまり面白くないものになりましたが、個人的には勉強になりました!以上です。

Pokemon Challenge - ポケモンデータセットを使った機械学習トライ

Pokemon Challenge - ポケモンデータセットを使った機械学習トライ

前回(http://rautaku.hatenablog.com/entry/2017/12/15/211915)の続きです。 前回はポケモン対戦データの解析を行いましたので、今回はそのデータから勝者予測を行ってみたいです。 変数やデータフレーム等は前回を引き継ぎますので、必要な方は前回を見てください。

勝者予測を行うにはデータの抜けなどがあると困りますので、まずはデータ補完を行なっていきます。 前回のポケモン個体ごとの勝率データを見返します。

  

battle_win.info()
<class 'pandas.core.frame.DataFrame'>
Index: 782 entries, Abomasnow to Zygarde Half Forme
Data columns (total 3 columns):
battle    782 non-null int64
win       781 non-null float64
ratio     781 non-null float64
dtypes: float64(2), int64(1)
memory usage: 24.4+ KB

上記情報を見ると、782体分しか情報がありません。またその内1体は勝星数、勝率情報 がありません。予測する側のデータにこのポケモン達がが出てきてしまうとデータ不足で予測できませんので補完していきます。まずは

battle_win[battle_win["win"].isnull()]

Shuckleこいつです。

ツボツボ。前回の勝星数算出処理の際に戦闘はしてるにも関わらず1勝もしていないため勝者リストに無く、正しく処理されなかったみたいです。ツボツボ悲しい。。win=0, ratio=0を追加してあげます。

battle_win.ix["Shuckle", ["win", "ratio"]] = 0
battle_win[battle_win.index=="Shuckle"]

次に戦闘をしていないためデータ抜けがある18体の補完を行いたいところですが、勝率などを補完するためには種族値などの情報も必要だと思われるので、図鑑に勝率を追加することで各ポケモンデータと勝率を結び付けます。

id_dict = dict(zip(pokemon['Name'], pokemon['#']))
battle_win["Name"] = battle_win.index
battle_win["#"] = battle_win["Name"].replace(id_dict)
ratio_dict = dict(zip(battle_win['#'], battle_win['ratio']))
pokemon["ratio"] = pokemon["#"].replace(ratio_dict)
pokemon.head()

図鑑に勝率情報が追加されました。図鑑ナンバーから勝率変換を行いましたので、データ抜けがあるポケモンたちの勝率欄には図鑑ナンバーが入ってしまっていると思います。調べます。

nobattle_pokemon = pokemon[pokemon["ratio"]>1]
print ("There are {} pokemons have NaN ratio.".format(len(nobattle_pokemon.index)))
nobattle_pokemon[["#", "Name", "ratio"]]
There are 18 pokemons have NaN ratio.

はい予想通り18体が該当しました。眺めて見ると1体名前が入ってないポケモンがいますね。#=63の格闘ポケモンです。普通に検索で調べてみるとPrimeapeオコリザルみたいです。追加しといてあげましょう。

pokemon.loc[62, "Name"] = "Primeape"
pokemon[pokemon["Name"]=="Primeape"][["#", "Name", "ratio"]]

さて勝率の補完ですが、前回の解析から種族値が勝率と相関がありそうでした。図示して見ます。

battle_pokemon = pokemon[pokemon["ratio"] <= 1]
sns.lmplot(x="stats_sum", y="ratio", data=battle_pokemon)
<seaborn.axisgrid.FacetGrid at 0x10eeaf048>

ばらつきはあるものの比較的相関はありそうですので、この相関を用いて、戦闘履歴のないポケモンの勝率を種族値から補完してみたいと思います。線形回帰を使います。

from sklearn.linear_model import LinearRegression
linreg = LinearRegression()
linreg.fit(battle_pokemon["stats_sum"].values.reshape(-1, 1), battle_pokemon["ratio"].values.reshape(-1, 1))
nobattle_pokemon["ratio"] = linreg.predict(nobattle_pokemon["stats_sum"].values.reshape(-1, 1))
nobattle_pokemon[["#", "Name", "ratio"]]

はい補完できました。これで予測作業のための補完は終了です。

ここからは勝者を予測するためのデータ整理を行っていきます。勝利に関与する因子として前回は種族値とタイプ相性を考察しましたが、他の情報としては、攻撃力、防御力、素早さなどもあるのでこれらも今回は入れ込みたいと思います。

戦闘結果combats.csvに今回因子としたい各種族値、タイプ相性情報を追加します。 【追加項目】 先攻側各種族値、勝率、後攻側各種族値、勝率、タイプ相性、先行側勝利判定

combats_add_data = combats.copy()
type_dict = dict(zip(pokemon['#'], pokemon['Type']))
hp_dict = dict(zip(pokemon['#'], pokemon['HP']))
attack_dict = dict(zip(pokemon['#'], pokemon['Attack']))
defense_dict = dict(zip(pokemon['#'], pokemon['Defense']))
spattack_dict = dict(zip(pokemon['#'], pokemon['Sp. Atk']))
spdefense_dict = dict(zip(pokemon['#'], pokemon['Sp. Def']))
speed_dict = dict(zip(pokemon['#'], pokemon['Speed']))
stats_sum_dict = dict(zip(pokemon['#'], pokemon['stats_sum']))
ratio_dict = dict(zip(pokemon['#'], pokemon['ratio']))
combats_add_data["First_pokemon_type"] = combats_add_data["First_pokemon"].replace(type_dict)
combats_add_data["First_pokemon_hp"] = combats_add_data["First_pokemon"].replace(hp_dict)
combats_add_data["First_pokemon_attack"] = combats_add_data["First_pokemon"].replace(attack_dict)
combats_add_data["First_pokemon_defense"] = combats_add_data["First_pokemon"].replace(defense_dict)
combats_add_data["First_pokemon_spattack"] = combats_add_data["First_pokemon"].replace(spattack_dict)
combats_add_data["First_pokemon_spdefense"] = combats_add_data["First_pokemon"].replace(spdefense_dict)
combats_add_data["First_pokemon_speed"] = combats_add_data["First_pokemon"].replace(speed_dict)
combats_add_data["First_pokemon_stats"] = combats_add_data["First_pokemon"].replace(stats_sum_dict)
combats_add_data["First_pokemon_ratio"] = combats_add_data["First_pokemon"].replace(ratio_dict)
combats_add_data["Second_pokemon_type"] = combats_add_data["Second_pokemon"].replace(type_dict)
combats_add_data["Second_pokemon_hp"] = combats_add_data["Second_pokemon"].replace(hp_dict)
combats_add_data["Second_pokemon_attack"] = combats_add_data["Second_pokemon"].replace(attack_dict)
combats_add_data["Second_pokemon_defense"] = combats_add_data["Second_pokemon"].replace(defense_dict)
combats_add_data["Second_pokemon_spattack"] = combats_add_data["Second_pokemon"].replace(spattack_dict)
combats_add_data["Second_pokemon_spdefense"] = combats_add_data["Second_pokemon"].replace(spdefense_dict)
combats_add_data["Second_pokemon_speed"] = combats_add_data["Second_pokemon"].replace(speed_dict)
combats_add_data["Second_pokemon_stats"] = combats_add_data["Second_pokemon"].replace(stats_sum_dict)
combats_add_data["Second_pokemon_ratio"] = combats_add_data["Second_pokemon"].replace(ratio_dict)


def calcTypeRelation(combats_add_data):
    r0 = 1
    first_type1 = combats_add_data["First_pokemon_type"].split("/")[0]
    first_type2 = combats_add_data["First_pokemon_type"].split("/")[1]
    second_type1 = combats_add_data["Second_pokemon_type"].split("/")[0]
    second_type2 = combats_add_data["Second_pokemon_type"].split("/")[1]
    if first_type2 != "None" and second_type2 != "None":
        r1 = df_type_relation[first_type1][second_type1]
        r2 = df_type_relation[first_type1][second_type2]
        r3 = df_type_relation[first_type2][second_type1]
        r4 = df_type_relation[first_type2][second_type2]
        r = r0 * r1 * r2 * r3 * r4
    elif first_type2 != "None" and second_type2 == "None":
        r1 = df_type_relation[first_type1][second_type1]
        r3 = df_type_relation[first_type2][second_type1]
        r = r0 * r1 * r3
    elif first_type2 == "None" and second_type2 != "None":
        r1 = df_type_relation[first_type1][second_type1]
        r2 = df_type_relation[first_type1][second_type2]
        r = r0 * r1 * r2
    elif first_type2 == "None" and second_type2 == "None":
        r1 = df_type_relation[first_type1][second_type1]
        r = r0 * r1
    return r

combats_add_data["Relation"] = combats_add_data.apply(lambda x: calcTypeRelation(x), axis = 1)
combats_add_data["First_win"] = combats_add_data.apply(lambda x: 1 if x["First_pokemon"]==x["Winner"] else 0, axis=1)
noneed_cols = ["First_pokemon", "Second_pokemon", "Winner", "First_pokemon_type", "Second_pokemon_type"]
combats_add_data = combats_add_data.drop(noneed_cols, axis=1)
combats_add_data.head()

先攻後攻それぞれのポケモンの情報、お互いのタイプ相性係数、その戦闘結果として先攻ポケモンが勝利したかどうかという情報が揃いました。 これを使って予測のためのモデル作成をしたいと思います。 モデルを作成した後、予測精度検証を行いたいので現在のcombats_add_dataデータフレームを訓練データ(全体の80%)とテストデータ(全体の20%)に分割します。

from sklearn.model_selection import train_test_split

X = combats_add_data.drop("First_win", axis=1)
y = combats_add_data["First_win"]
X_train, X_test, y_train, y_test = train_test_split(X, y, train_size=0.8, random_state=1)
print("X_train.shape = " + str(X_train.shape))
print("X_test.shape = " + str(X_test.shape))
print("y_train.shape = " + str(y_train.shape))
print("y_test.shape = " + str(y_test.shape))
X_train.shape = (40000, 17)
X_test.shape = (10000, 17)
y_train.shape = (40000,)
y_test.shape = (10000,)

準備が全て整いましたので学習アルゴリズムを用いて予測していきます。今回は勉強のため数種類の教師あり学習器を用います。どの学習器が一番精度が良いか検証します。 比較するアルゴリズムは、ロジスティック回帰、k近傍法、ガウシアンナイーブベイズパーセプトロン、決定木、ランダムフォレストを使ってみます。

# Logistic Regression
from sklearn.linear_model import LogisticRegression
logreg = LogisticRegression()
logreg.fit(X_train, y_train)
acc_log = round(logreg.score(X_test, y_test)*100, 2)
acc_log
88.659999999999997
from sklearn.neighbors import KNeighborsClassifier
knn = KNeighborsClassifier(n_neighbors = 3)
knn.fit(X_train, y_train)
acc_knn = round(knn.score(X_test, y_test) * 100, 2)
acc_knn
85.659999999999997
# Gaussian Naive Bayes
from sklearn.naive_bayes import GaussianNB
gaussian = GaussianNB()
gaussian.fit(X_train, y_train)
acc_gaussian = round(gaussian.score(X_test, y_test) * 100, 2)
acc_gaussian
77.810000000000002
# Perceptron
from sklearn.linear_model import Perceptron
perceptron = Perceptron()
perceptron.fit(X_train, y_train)
acc_perceptron = round(perceptron.score(X_test, y_test) * 100, 2)
acc_perceptron
80.969999999999999
# Decision Tree
from sklearn.tree import DecisionTreeClassifier
decision_tree = DecisionTreeClassifier()
decision_tree.fit(X_train, y_train)
acc_decision_tree = round(decision_tree.score(X_test, y_test) * 100, 2)
acc_decision_tree
93.159999999999997
# Random Forest
from sklearn.ensemble import RandomForestClassifier
random_forest = RandomForestClassifier(n_estimators=100)
random_forest.fit(X_train, y_train)
acc_random_forest = round(random_forest.score(X_test, y_test) * 100, 2)
acc_random_forest
94.959999999999994

各モデルの精度順に並べてみますと以下の様になります。

models = pd.DataFrame({
    'Model': ['Logistic Regression', 'KNN', 
              'Naive Bayes', 'Perceptron', 'Decision Tree', 'Random Forest'],
    'Score': [acc_log, acc_knn, acc_gaussian, acc_perceptron, 
              acc_decision_tree, acc_random_forest]})
models.sort_values(by='Score', ascending=False)

ランダムフォレストが一番精度が良く95%という成果を出すことができました。結構いい水準だと思います。 ちなみに各因子の影響度はどの程度だったんでしょうか?

effective = pd.DataFrame()
effective["feature_name"] = X.columns.tolist()
effective["feature_importance"] = random_forest.feature_importances_
effective.sort_values("feature_importance",ascending=False)

意外な結果が出ました。一番影響度の高い因子は、種族値や勝率ではなく素早さみたいですね。 種族値やタイプ相性が多少不利な状況でも素早さが高ければ勝つ確率は高いという結果が得られました。 ちなみに素早さ10傑は以下だそうです!勝率10傑とは少しメンツが変わりましたね。

pokemon.sort_values("Speed",ascending=False).head(10)

今回作成したモデルをtests.csvに反映させれば予測完了です。モデル反映は非常に簡単な手順ですので今回は省略させてもらいます。 以上、ポケモンデータセットを使った機械学習のトライでした!

Pokemon Challenge - ポケモンデータセットを使った解析トライ

Pokemon Challenge - ポケモンデータセットを使った解析トライ

Kaggleを見ていると面白そうなデータセットを見つけたので、暇つぶしに考察してみます。素人ですので本格的な解析を求める方はKaggleに行ってください。 https://www.kaggle.com/terminus7/pokemon-challenge/kernels ポケモンの対戦データセットです。データセットの中身は,combats.csv, pokemon.csv, tests.csvの3ファイル。 combats.csvはバトル結果が記録されています。お互いの図鑑番号と勝者側番号です。先に書いてある方が先制側です。50,000戦の結果があります。 pokemon.csvポケモン図鑑ですね。番号、名前、タイプ1、タイプ2、各種族値(攻撃、防御、特攻、特防、素早さ)、世代、伝説系か。800体分あります。 tests.csvはテストデータです。対戦組合せが提示されているのでその勝敗を予想しなさいということです。 とりあえず今回はデータからどの様な傾向があるのかなーっていう解析をしてみたいと思います。

ポケモンルビサファまでしかやってない軽度ユーザーの自分からしても、努力値のデータも技構成のデータなどもない種族値だけのデータでは何を予想しても仕方ない気がするのですが、とりあえず触ってみます。この対戦結果はどの様な環境から取ってきたものなんですかね。

データセットを軽く見ていきます。combats.csvから。

import pandas as pd
combats = pd.read_csv('combats.csv')
combats.head(3)

266番と298番が闘って298番が勝ってます。266番と298番ってどのポケモンでしょうか?pokemon.csvから調べます。

pokemon = pd.read_csv('pokemon.csv')
pokemon_266_298 = pokemon[pokemon['#'].isin([266, 298])]
pokemon_266_298

LavitarとNuzleafです。和名はLavitarヨーギラスとNuzleafコノハナです。こいつらです。

勝ったのはヨーギラスヨーギラスは岩地面で、コノハナは草悪。タイプ面でも普通にコノハナの方が有利ですね。ちなみにタイプ相性はこちら

そんな感じで対戦結果の上から3つはこんな感じ

names_dict = dict(zip(pokemon['#'], pokemon['Name']))
cols = ["First_pokemon","Second_pokemon","Winner"]
combats_name = combats[cols].replace(names_dict)
combats_name.head(3)

Virizionビリジオン(格闘草):Terrakionテラキオン(格闘岩)-> テラキオン勝利 Togeticトゲチック(フェアリー飛行):Beheeyemオーベムエスパー)-> オーベム勝利 そんな感じらしいです。うーんあまり傾向読めませんね。でも比較的同レベル同士が対戦してる感じがします。ポケモン種族値合計で700族や600族などのクラス分けされる傾向があるので調べてみましょう。 その前にまずはそもそもの全ポケモン種族値分布を。

pokemon["stats_sum"] = pokemon["HP"] + pokemon["Attack"] + pokemon["Defense"] + pokemon["Sp. Atk"] + pokemon["Sp. Def"] + pokemon["Speed"]
import matplotlib.pyplot as plt
import seaborn as sns
%matplotlib inline
print(pokemon["stats_sum"].describe())
sns.distplot(pokemon["stats_sum"])
plt.show()
count    800.00000
mean     435.10250
std      119.96304
min      180.00000
25%      330.00000
50%      450.00000
75%      515.00000
max      780.00000
Name: stats_sum, dtype: float64

300族と500族あたりに山がありますね。さてでは戦闘組み合わせの種族値差分布を見ます。

stats_sum_dict = dict(zip(pokemon['#'], pokemon['stats_sum']))
combats_stats_sum = combats[cols].replace(stats_sum_dict)
diff_stats_sum = abs(combats_stats_sum["First_pokemon"]-combats_stats_sum["Second_pokemon"])
print(diff_stats_sum.describe())
sns.distplot(diff_stats_sum)
plt.xlabel("diff_stats_sum")
plt.show()
count    50000.000000
mean       136.551440
std        101.221212
min          0.000000
25%         54.000000
50%        118.000000
75%        200.000000
max        590.000000
dtype: float64

うーん差が少ない方向に偏ってはいますが、種族値差が開いている組合せも結構あります。中間値が118です。100の差ってどうなんですかね、結構重要そうに思えます。次は種族値が高い方がちゃんと勝ってるのか調べましょう。

combats_stats_sum["Loser"] = combats_stats_sum.apply(lambda x: x["First_pokemon"] if x["First_pokemon"] !=  x["Winner"] else x["Second_pokemon"], axis = 1)
diff_win_lose_stats = combats_stats_sum["Winner"] - combats_stats_sum["Loser"]
print(diff_win_lose_stats.describe())
sns.distplot(diff_win_lose_stats)
plt.xlabel("diff_win_lose_stats")
plt.show()
count    50000.000000
mean        79.776320
std        150.093351
min       -562.000000
25%        -20.000000
50%         80.000000
75%        186.000000
max        590.000000
dtype: float64

まぁ中間値80と比較的種族値が高い方が勝つ傾向は見て取れます。種族値というものが大きな勝利要因になっていることは言えるでしょう。でも種族値200以上も差がついていながら低い方が勝利するという下剋上な結果もあります。 下剋上を起こした組み合わせはどの様な組合せか見てみます。ここでの下剋上の定義は種族値100以上差がありながら低い方が勝つことととします。

stats_sum_dict_re = dict(zip(pokemon['stats_sum'], pokemon['#']))
combats_stats_sum["diff"] = diff_win_lose_stats
surpassing_stats_sum = combats_stats_sum[combats_stats_sum["diff"] < -100]
print ("Surpassing one's superiors Battle number : " + str(len(surpassing_stats_sum)))
surpassing_id = surpassing_stats_sum[cols].replace(stats_sum_dict_re)
surpassing_name = surpassing_id[cols].replace(names_dict)
surpassing_name.join(combats_stats_sum["diff"]).head(8)
Surpassing one's superiors Battle number : 5716

下剋上の組合せは5716組あります。全体の約10%です。上から8個取り出すとこんな感じ。

結構面白いですね。自分でも見たことある映画伝説勢のボルケニオンがいかにも一般勢なゴチミルというポケモンに負けてます。最大の下剋上の組み合わせが気になってきたので調べてみます。

surpassing_name.join(combats_stats_sum["diff"]).sort_values(by="diff").head(4)

メガレックウザやられすぎですね。ピィがメガレックウザに勝つことなんてあるんでしょうか??562の種族値差です。

ではこれらの下剋上はなぜ起こったのかを探っていきましょう。思いつく理由としては、やはり先制攻撃ではないでしょうか。最近のポケモンは技も多様化しているので、種族値差があっても先攻を取ることで一気に落としてしまうことも見られる気がします。(まぁ上記のメガレックウザvsピィは後攻のピィが連勝していますが)

combats_stats_sum["First_Win"] =  combats_stats_sum.apply(lambda x: 1 if x["First_pokemon"] ==  x["Winner"] else 0, axis = 1)
surpassing_stats_sum = combats_stats_sum[combats_stats_sum["diff"] < -100]
sns.distplot(surpassing_stats_sum[surpassing_stats_sum["First_Win"]==0]["diff"], label="First_pokemon_win=0")
sns.distplot(surpassing_stats_sum[surpassing_stats_sum["First_Win"]==1]["diff"], label="First_pokemon_win=1")
plt.legend()
plt.show()

うーん全く差がないですね。。下剋上組合せだけでなく全体でも同じでしょうか

sns.distplot(combats_stats_sum[combats_stats_sum["First_Win"]==0]["diff"], label="First_Win=0")
sns.distplot(combats_stats_sum[combats_stats_sum["First_Win"]==1]["diff"], label="First_Win=1")
print ("-First_Win=0-")
print (combats_stats_sum[combats_stats_sum["First_Win"]==0]["diff"].describe())
print ("-First_Win=1-")
print (combats_stats_sum[combats_stats_sum["First_Win"]==1]["diff"].describe())
plt.legend()
plt.show()
-First_Win=0-
count    26065.000000
mean        76.155764
std        152.292707
min       -562.000000
25%        -27.000000
50%         79.000000
75%        185.000000
max        585.000000
Name: diff, dtype: float64
-First_Win=1-
count    23935.000000
mean        83.719072
std        147.563182
min       -544.000000
25%        -13.000000
50%         82.000000
75%        190.000000
max        590.000000
Name: diff, dtype: float64

全体でもほとんど差がない上に、むしろ先攻の方が勝ち数少ないですw なんなんでしょうかね1ターン目はステータスアップ系の技でも使う傾向が高いんですかね。とりあえず先制後攻はほとんど関係ないみたいです。 じゃあ次はタイプ相性を見ます。ポケモンといえばやはりタイプ相性でしょう。まずは存在するタイプ一覧を回収。

print ("There are {} Types.".format(len(pokemon["Type 1"].drop_duplicates())))
list(pokemon["Type 1"].drop_duplicates())
There are 18 Types

['Grass',
 'Fire',
 'Water',
 'Bug',
 'Normal',
 'Poison',
 'Electric',
 'Ground',
 'Fairy',
 'Fighting',
 'Psychic',
 'Rock',
 'Ghost',
 'Ice',
 'Dragon',
 'Dark',
 'Steel',
 'Flying']

18種類あるみたいですね。今はタイプが2つ持てますからその組み合わせまで考えると何通りあるんでしょうか。

type_cols = ["Type 1", "Type 2"]
print ("There are {} type-combinations.".format(len(pokemon[type_cols].drop_duplicates())))
There are 154 type-combinations.

154種類です。考えられる組合せは18*18/2=162ですからほぼ全ての組合せがいますね。これらの分布を可視化してみます。

pokemon["Type 2"] = pokemon["Type 2"].fillna("None")
type_cross = pd.crosstab(pokemon["Type 1"], pokemon["Type 2"])
type_cross.plot.bar(stacked=True, figsize=(8,6))
plt.legend(bbox_to_anchor=(0.01, 0.99), loc='upper left', ncol=3, fontsize=8, title="Type 2")
plt.show()

これだけの組合せがありますから比較検証するのは大変そうですね。。 ひとまずタイプ相関も数値化します。基本的にはダメージ変動は効果バツグンなら2倍、効果はいまひとつで0.5倍なので、それに応じた対応表を作ります。対応表見ながら手づくりします。

Normal = {"Normal": 1, "Fighting": 1, "Poison": 1, "Ground": 1, "Flying": 1, "Bug": 1, "Rock": 0.5, "Ghost": 0, "Steel": 0.5, "Fire": 1, "Water": 1, "Electric": 1, "Grass": 1, "Ice": 1, "Psychic": 1, "Dragon": 1, "Dark": 1, "Fairy": 1}
Fighting = {"Normal": 2, "Fighting": 1, "Poison": 0.5, "Ground": 1, "Flying": 0.5, "Bug": 0.5, "Rock": 2, "Ghost": 0, "Steel": 2, "Fire": 1, "Water": 1, "Electric": 1, "Grass": 1, "Ice": 2, "Psychic": 0.5, "Dragon": 1, "Dark": 2, "Fairy": 0.5}
Poison = {"Normal": 1, "Fighting": 1, "Poison": 0.5, "Ground": 0.5, "Flying": 1, "Bug": 1, "Rock": 0.5, "Ghost": 0.5, "Steel": 0, "Fire": 1, "Water": 1, "Electric": 1, "Grass": 2, "Ice": 1, "Psychic": 1, "Dragon": 1, "Dark": 1, "Fairy": 2}
Ground = {"Normal": 1, "Fighting": 1, "Poison": 2, "Ground": 1, "Flying": 0, "Bug": 0.5, "Rock": 2, "Ghost": 1, "Steel": 2, "Fire": 2, "Water": 1, "Electric": 2, "Grass": 0.5, "Ice": 1, "Psychic": 1, "Dragon": 1, "Dark": 1, "Fairy": 1}
Flying = {"Normal": 1, "Fighting": 2, "Poison": 1, "Ground": 1, "Flying": 1, "Bug": 2, "Rock": 0.5, "Ghost": 1, "Steel": 0.5, "Fire": 1, "Water": 1, "Electric": 0.5, "Grass": 2, "Ice": 1, "Psychic": 1, "Dragon": 1, "Dark": 1, "Fairy": 1}
Bug = {"Normal": 1, "Fighting": 0.5, "Poison": 0.5, "Ground": 1, "Flying": 0.5, "Bug": 1, "Rock": 1, "Ghost": 0.5, "Steel": 0.5, "Fire": 0.5, "Water": 1, "Electric": 1, "Grass": 2, "Ice": 1, "Psychic": 2, "Dragon": 1, "Dark": 2, "Fairy": 0.5}
Rock = {"Normal": 1, "Fighting": 0.5, "Poison": 1, "Ground": 0.5, "Flying": 2, "Bug": 2, "Rock": 1, "Ghost": 1, "Steel": 0.5, "Fire": 2, "Water": 1, "Electric": 1, "Grass": 1, "Ice": 2, "Psychic": 1, "Dragon": 1, "Dark": 1, "Fairy": 1}
Ghost = {"Normal": 0, "Fighting": 1, "Poison": 1, "Ground": 1, "Flying": 1, "Bug": 1, "Rock": 1, "Ghost": 2, "Steel": 1, "Fire": 1, "Water": 1, "Electric": 1, "Grass": 1, "Ice": 1, "Psychic": 2, "Dragon": 1, "Dark": 0.5, "Fairy": 1}
Steel = {"Normal": 1, "Fighting": 1, "Poison": 1, "Ground": 1, "Flying": 1, "Bug": 1, "Rock": 2, "Ghost": 1, "Steel": 0.5, "Fire": 0.5, "Water": 0.5, "Electric": 0.5, "Grass": 1, "Ice": 2, "Psychic": 1, "Dragon": 1, "Dark": 1, "Fairy": 0.5}
Fire = {"Normal": 1, "Fighting": 1, "Poison": 1, "Ground": 1, "Flying": 1, "Bug": 2, "Rock": 0.5, "Ghost": 1, "Steel": 2, "Fire": 0.5, "Water": 0.5, "Electric": 1, "Grass": 2, "Ice": 2, "Psychic": 1, "Dragon": 0.5, "Dark": 1, "Fairy": 1}
Water = {"Normal": 1, "Fighting": 1, "Poison": 1, "Ground": 2, "Flying": 1, "Bug": 1, "Rock": 2, "Ghost": 1, "Steel": 1, "Fire": 2, "Water": 0.5, "Electric": 1, "Grass": 0.5, "Ice": 1, "Psychic": 1, "Dragon": 0.5, "Dark": 1, "Fairy": 1}
Electric = {"Normal": 1, "Fighting": 1, "Poison": 1, "Ground": 0, "Flying": 2, "Bug": 1, "Rock": 1, "Ghost": 1, "Steel": 1, "Fire": 1, "Water": 2, "Electric": 0.5, "Grass": 0.5, "Ice": 1, "Psychic": 1, "Dragon": 0.5, "Dark": 1, "Fairy": 1}
Grass = {"Normal": 1, "Fighting": 1, "Poison": 0.5, "Ground": 2, "Flying": 0.5, "Bug": 0.5, "Rock": 2, "Ghost": 1, "Steel": 0.5, "Fire": 0.5, "Water": 2, "Electric": 1, "Grass": 0.5, "Ice": 1, "Psychic": 1, "Dragon": 0.5, "Dark": 1, "Fairy": 1}
Ice = {"Normal": 1, "Fighting": 1, "Poison": 1, "Ground": 2, "Flying": 2, "Bug": 1, "Rock": 1, "Ghost": 1, "Steel": 0.5, "Fire": 0.5, "Water": 0.5, "Electric": 1, "Grass": 2, "Ice": 0.5, "Psychic": 1, "Dragon": 2, "Dark": 1, "Fairy": 1}
Psychic = {"Normal": 1, "Fighting": 1, "Poison": 2, "Ground": 2, "Flying": 1, "Bug": 1, "Rock": 1, "Ghost": 1, "Steel": 0.5, "Fire": 1, "Water": 1, "Electric": 1, "Grass": 1, "Ice": 1, "Psychic": 0.5, "Dragon": 1, "Dark": 0, "Fairy": 1}
Dragon = {"Normal": 1, "Fighting": 1, "Poison": 1, "Ground": 1, "Flying": 1, "Bug": 1, "Rock": 1, "Ghost": 1, "Steel": 0.5, "Fire": 1, "Water": 1, "Electric": 1, "Grass": 1, "Ice": 1, "Psychic": 1, "Dragon": 2, "Dark": 1, "Fairy": 0}
Dark = {"Normal": 1, "Fighting": 0.5, "Poison": 1, "Ground": 1, "Flying": 1, "Bug": 1, "Rock": 1, "Ghost": 2, "Steel": 1, "Fire": 1, "Water": 1, "Electric": 1, "Grass": 1, "Ice": 1, "Psychic": 2, "Dragon": 1, "Dark": 0.5, "Fairy": 0.5}
Fairy = {"Normal": 1, "Fighting": 2, "Poison": 0.5, "Ground": 1, "Flying": 1, "Bug": 1, "Rock": 1, "Ghost": 1, "Steel": 0.5, "Fire": 0.5, "Water": 1, "Electric": 1, "Grass": 1, "Ice": 1, "Psychic": 1, "Dragon": 2, "Dark": 2, "Fairy": 1}

type_relation = {"Normal": Normal, "Fighting": Fighting, "Poison": Poison, "Ground": Ground, "Flying": Flying, "Bug": Bug, "Rock": Rock, "Ghost": Ghost, "Steel": Steel, "Fire": Fire, "Water": Water, "Electric": Electric, "Grass": Grass, "Ice": Ice, "Psychic": Psychic, "Dragon": Dragon, "Dark": Dark, "Fairy": Fairy}
df_type_relation = pd.DataFrame(type_relation)
print ("Row is Diffender, Column is Attacker")
df_type_relation
Row is Diffender, Column is Attacker

対応表ができたのでまずはタイプの勝敗結果への影響を見てみます。勝利側を攻撃側、敗退側を防御側として各タイプに対する数値を掛け合わせたもので比較してみます。無効の関係性にあるものは数値0だと掛け合わせ時に困るので0.25に変更します。

pokemon["Type"] = pokemon.apply(lambda x: x["Type 1"]+"/"+x["Type 2"], axis=1)
type_dict = dict(zip(pokemon['#'], pokemon['Type']))
combats_type = combats[cols].replace(type_dict)
combats_type["Loser"] = combats_type.apply(lambda x: x["First_pokemon"] if x["First_pokemon"] !=  x["Winner"] else x["Second_pokemon"], axis = 1)

zero_dict = {0: 0.25}
df_type_relation = df_type_relation[:].replace(zero_dict)

def calcRelation(combats_type):
    r0 = 1
    win_type1 = combats_type["Winner"].split("/")[0]
    win_type2 = combats_type["Winner"].split("/")[1]
    lose_type1 = combats_type["Loser"].split("/")[0]
    lose_type2 = combats_type["Loser"].split("/")[1]
    if win_type2 != "None" and lose_type2 != "None":
        r1 = df_type_relation[win_type1][lose_type1]
        r2 = df_type_relation[win_type1][lose_type2]
        r3 = df_type_relation[win_type2][lose_type1]
        r4 = df_type_relation[win_type2][lose_type2]
        r = r0 * r1 * r2 * r3 * r4
    elif win_type2 != "None" and lose_type2 == "None":
        r1 = df_type_relation[win_type1][lose_type1]
        r3 = df_type_relation[win_type2][lose_type1]
        r = r0 * r1 * r3
    elif win_type2 == "None" and lose_type2 != "None":
        r1 = df_type_relation[win_type1][lose_type1]
        r2 = df_type_relation[win_type1][lose_type2]
        r = r0 * r1 * r2
    elif win_type2 == "None" and lose_type2 == "None":
        r1 = df_type_relation[win_type1][lose_type1]
        r = r0 * r1
    return r

combats_type["Relation"] = combats_type.apply(lambda x: calcRelation(x), axis = 1)
print (combats_type["Relation"].describe())
sns.distplot(combats_type["Relation"])
plt.show()
count    50000.000000
mean         1.146062
std          0.869392
min          0.031250
25%          0.500000
50%          1.000000
75%          1.000000
max         16.000000
Name: Relation, dtype: float64

あまり特徴的な点はない様ですね。むしろ不利なタイプ相手の方が若干勝星多そうな感じがしますね。まぁ今回の係数の組合せ程度の比較ではあまり見えてこないということですね。この辺はまた時間ある時にでも考察できればと思います。

もうなんだか光が見えないですし、技構成とか戦闘環境とかわからなくて情報が少なすぎる!ってことで、最後に戦闘勝率が高いポケモンが何か調べて終わります! まずは組み合わせに出てくる頻度の高いポケモンは何でしょうか。

import numpy as np
from wordcloud import WordCloud, ImageColorGenerator
from PIL import Image

combats_names = combats[cols].replace(names_dict)
print (combats_names["Winner"].value_counts()[:10])
winners = list(combats_names["Winner"])
winners_str = [str(i) for i in winners]
winners_text = (",").join(winners_str)
maskPicture = np.array(Image.open("Pikachu.png", "r"))
imageColor = ImageColorGenerator(maskPicture)
wc = WordCloud(background_color= "black", random_state=1, margin=3, mask=maskPicture).generate(winners_text)
plt.figure(figsize=(10,10))
plt.axis("off")
plt.imshow(wc.recolor(color_func=imageColor))
plt.show()
Mewtwo                152
Aerodactyl            136
Infernape             136
Jirachi               134
Deoxys Speed Forme    133
Slaking               133
Murkrow               130
Mega Absol            130
Mega Houndoom         128
Mega Aerodactyl       127
Name: Winner, dtype: int64

WordCloud上で目立っているTherian FormeとIncarnate Formeは霊獣フォルムと化身フォルムで、ランドロスボルトロストルネロスという奴らの共通変化形態であり重複してる様です。

勝星数10傑はMewtwoミュウツー, Aerodactylプテラ, Infernapeゴウカザル, Jirachiジラーチ, Deoxys Speed Formeデオキシススピードフォルム, Slakingケッキング, Murkrowムクロー, Mega Absolメガアブソル, Mega Houndoomメガヘルガー, Mega Aerodactylメガプテラです。 ただこれらは単純に勝星数でありたくさん戦えばいいものなので求めるものとは違う気がするので勝率観点で調べてみましょう。

first_num = combats_names["First_pokemon"].value_counts()
second_num = combats_names["Second_pokemon"].value_counts()
battle_num = first_num + second_num
battle_win = pd.DataFrame({"battle": battle_num, "win": combats_names["Winner"].value_counts()}, columns=["battle", "win"])
battle_win["ratio"] = battle_win["win"]/battle_win["battle"]
battle_win.sort_values(by=["ratio"], ascending=False).head(10)

出ました。これが勝率10傑です!とりあえず今回のデータセット内ではこいつらを使えば勝ちやすくなるみたいです!

今回は以上にしておきます。次回は機械学習を使ってtests.csvを簡単に予測などしてみたいなと思っています。