Uplift Modeling用のパッケージtools4upliftを使ってみた

はじめに

今回は、今後仕事で使いたいという思いもあり、RでUplift Modelingに関して便利なパッケージがないか探した結果、2019年に登場したばかりのtools4upliftの存在を知りました。アップリフトモデリングのモチベーションに関しても簡単に説明しながら、サンプルデータで実践してみようと思います。

・Uplift Modelingとはなにか
・Uplift Modelingの卑近な例え話
・Uplift Modelingのサンプルデータ
・tools4upliftについて
・tools4upliftでCriteoデータを試してみる
・『仕事ではじめる機械学習』の9章のコードをCriteoデータに試してみる
・おわりに
・参考文献

Uplift Modelingとはなにか

きちんとした説明は、あまりにも今更感があるので説明は端折りたいと思います。既出の文献がありますので、そちらを熟読ください。

Uplift Modelingの卑近な例え話

自分が吉野家のマーケティング担当だとしましょう。吉野家のアプリで割引クーポンを顧客にばらまくことができるとします。
マーケターとして重要なのは、割引クーポンを渡したことをきっかけとして吉野家に足を運び購入する顧客を増やせるかどうかになります。

マーケターの手元にあるのは、割引クーポンをばらまいた顧客とばらまかなかった顧客、そして吉野家で牛丼を食べたかどうかのデータです。
以前のマーケティング担当者がランダムにクーポンをばらまいていたことが重要なポイントです。

このデータから、顧客は以下の4分類に分かれます。

無関心：割引クーポンをばらまこうが我関せず。そもそも吉野家に行く気はない。
説得可能：普段、牛丼が安いすき屋にばかり行っているが、割高に感じている吉野家に負い目を感じている。割引クーポンで揺さぶられ来店する。
天の邪鬼：吉野家コピペのように、割引クーポンを握りしめた家族連れに遭遇したくないので、割引クーポンをばらまかれたら来店しないような客。
鉄板：毎日決まった時間に吉野家に行くことを心に決めている客。

マーケターは割引クーポンをばらまいた顧客と割引クーポンをばらまいていない顧客にデータを二分し、それぞれ機械学習のための訓練用データとテスト用データを用意します。

つまり、「割引クーポンをばらまいた顧客」の訓練用データとテスト用データと「割引クーポンをばらまいていない顧客」の訓練用データとテスト用データの計4つのデータセットを用意します。

まず、牛丼の購入の有無を教師とした訓練用データでロジスティック回帰モデルなどを推定します。
その結果、「割引クーポンをばらまいた顧客」から推定したモデルと、「割引クーポンをばらまいていない顧客」から推定したモデルが手元に残ります。

2つのテスト用データを1つにまとめて、先程推定したモデルを用いて、牛丼の購入確率を求めます。モデルは２つあるので、予測結果がテスト用データ1つに対して2つあることになります。

その予測結果の比（「割引クーポンをばらまいた顧客」モデルベースの予測値÷「割引クーポンをばらまいていない顧客」モデルベースの予測値）をアップリフトとみなします。

以下の図はこれまでの説明を図にしたものです。

アップリフトがどの程度の水準であれば、説得可能なユーザーが多いのかを探っていくことで、吉野家のアプリにおいて、どのユーザーに割引クーポンを発行するべきかがわかることになります。

Uplift Modelingのサンプルデータ

残念なことに吉野家のアプリのデータはありません。そこで今回は公開データを利用します。
以前より、The MineThatData E-Mail Analytics And Data Mining ChallengeのメールのデータがUplift Modelingで非常にしばしば取り上げられるデータでしたが、Twitterで他にデータないのかとぼやいたところ、2名の方にCriteo Uplift Prediction Datasetを紹介していただきました。

余談ですが、Criteo社と言えばディスプレイ広告のキング的な存在で、少し商品のリンクを踏んだだけであっという間に広告がレコメンドされますよね。自社で出稿用バナーを作っていましたが、CVRが高くなる良いクリエイティブを作ってきたのか、単にCriteo社のアルゴリズムが優秀なだけなのか非常に気になるところでしたね。

Criteo社が提供してくれている今回のデータは、2500万行に及ぶユーザーのデータで、プライバシー保護の観点から特徴量は復元できないような形式で提供されています。バイナリーのラベルとしては訪問やコンバージョンなどがあり、データ全体に占める処置群の割合は84.6%となっています。要は、吉野家で言う割引クーポンをばらまいた顧客が全体の84.6%に及ぶということです。

tools4upliftについて

2019年1月に公開されたRのUplift Modeling用のパッケージです。

特徴量における連続値をカテゴリ変数にする際に、最適な階級値を求めてくれる関数
アップリフトモデリングの可視化する関数
アップリフトモデリングにおける特徴量選択ができる関数
アップリフトモデリングにおけるモデルのバリデーションを行う関数

などが提供されており、ちょいとRを触れるマーケターにとって、アップリフトモデリングにおける試行錯誤がかなりしやすくなる便利なパッケージだと思いました。
なお、このパッケージで扱っているモデルはロジスティック回帰になります。介入データをもとに推定したモデルの条件付き確率と非介入データをもとに推定したモデルの条件付き確率の差をアップリフトとして推定しています。

このパッケージの解説論文においては、アップリフトモデリングの評価指標としてQini曲線というものが提案されていました。Qini曲線はローレンツ曲線のようなもので、Qini曲線とランダムに割り当てた際のアップリフト量の差分の合計をQini係数と定義しています。

tools4upliftでCriteoデータを試してみる

library(tools4uplift)
library(tidyverse)
library(data.table)

x = fread("gunzip -c criteo-uplift.csv.gz")

head(x,5)

# 1000万件でデータをサンプリング
sampleNum <- sample(nrow(x),10000000)
x_sub <- x[sampleNum,]

x_sub <- x_sub %>% select(-exposure,-visit)
summary(x_sub)

# Baseline models ---------------------------------------------------------
set.seed(123);
split.data1 <- SplitUplift(data = x_sub,
                           p = 0.7,
                           group = c("treatment", "conversion"))

train <- split.data1[[1]]
valid <- split.data1[[2]]

base.tm <- DualUplift(data = train,
                      treat = "treatment",
                      outcome = "conversion",
                      predictors = colnames(train[,1:12]))

# baseline model for control group
base.tm[[1]]

# baseline model for treatment group
base.tm[[2]]

# predict the uplift on the validation set
base.tm.valid <- DualPredict(data = valid,
                             treat = "treatment",
                             outcome = "conversion", 
                             model = base.tm, nb.group = 5)[[1]]

# evaluate the model performance
base.tm.perf <- QiniTable(data = base.tm.valid,
                          treat = "treatment",
                          outcome = "conversion",
                          prediction = "uplift_prediction",
                          nb.group = 5)

# Qini曲線の描画
QiniCurve(base.tm.perf, title = "")
# アップリフト量の棒グラフの描画
QiniBarPlot(base.tm.perf, title = "")
# Qini係数の算出
QiniArea(base.tm.perf)

library(tools4uplift)

library(tidyverse)

library(data.table)

x = fread("gunzip -c criteo-uplift.csv.gz")

head(x,5)

# 1000万件でデータをサンプリング

sampleNum <- sample(nrow(x),10000000)

x_sub <- x[sampleNum,]

x_sub <- x_sub %>% select(-exposure,-visit)

summary(x_sub)

# Baseline models ---------------------------------------------------------

set.seed(123);

split.data1 <- SplitUplift(data = x_sub,

p = 0.7,

group = c("treatment", "conversion"))

train <- split.data1[[1]]

valid <- split.data1[[2]]

base.tm <- DualUplift(data = train,

treat = "treatment",

outcome = "conversion",

predictors = colnames(train[,1:12]))

# baseline model for control group

base.tm[[1]]

# baseline model for treatment group

base.tm[[2]]

# predict the uplift on the validation set

base.tm.valid <- DualPredict(data = valid,

treat = "treatment",

outcome = "conversion",

model = base.tm, nb.group = 5)[[1]]

# evaluate the model performance

base.tm.perf <- QiniTable(data = base.tm.valid,

treat = "treatment",

outcome = "conversion",

prediction = "uplift_prediction",

nb.group = 5)

# Qini曲線の描画

QiniCurve(base.tm.perf, title = "")

# アップリフト量の棒グラフの描画

QiniBarPlot(base.tm.perf, title = "")

# Qini係数の算出

QiniArea(base.tm.perf)

こちらはアップリフト値の予測値の上位から右に並べた際のアップリフトの増大のグラフになります。20%あたりでピークになるようです。

こちらはアップリフト量の棒グラフです。20%の階級値を超えたらガクンと下がるのがわかります。

なお、Qini係数は0.03233551でした。

『仕事ではじめる機械学習』の9章のコードをCriteoデータに試してみる

tools4upliftの結果を鵜呑みにするのもあれなので、『仕事ではじめる機械学習』の9章のコードを使ってアップリフトモデリングを実践してみます。コードは丸パクリですが、謹んで掲載させていただきます。

%matplotlib inline

import pandas as pd
import matplotlib.pyplot as plt
from operator import  itemgetter
plt.style.use("ggplot")

from sklearn.model_selection import  train_test_split
from sklearn.linear_model import LogisticRegression

# https://ailab.criteo.com/criteo-uplift-prediction-dataset/
source_df = pd.read_csv("criteo-uplift.csv.gz")
source_df.head(10)

source_df.describe()

feature_vector_df = source_df.drop(["treatment","conversion","visit","exposure"],axis=1)

is_treat_list = list(source_df["treatment"] == 1)
is_cv_list = list(source_df["conversion"] == 1)

train_is_cv_list, test_is_cv_list, train_is_treat_list, \
        test_is_treat_list, train_feature_vector_df,\
        test_feature_vector_df = train_test_split(is_cv_list, is_treat_list,
                                                                                feature_vector_df,
                                                                                train_size=0.5,
                                                                                test_size=0.5,
                                                                                random_state=42)

treat_model = LogisticRegression(C=0.01)
control_model = LogisticRegression(C=0.01)

train_sample_num = len(train_is_cv_list)

treat_is_cv_list = [train_is_cv_list[i] for i in range(train_sample_num) if train_is_treat_list[i] == True]
treat_feature_vector_list = train_feature_vector_df[train_is_treat_list]

control_is_cv_list = [train_is_cv_list[i] for i in range(train_sample_num) if train_is_treat_list[i] == False]
control_feature_vector_list = train_feature_vector_df[list(map(lambda a:a == False ,train_is_treat_list))]

treat_model.fit(treat_feature_vector_list, treat_is_cv_list)
control_model.fit(control_feature_vector_list, control_is_cv_list)

treat_score = treat_model.predict_proba(test_feature_vector_df)
control_score = control_model.predict_proba(test_feature_vector_df)
score_list = treat_score[:,1] / control_score[:,1]

result = list(zip(test_is_cv_list, test_is_treat_list, score_list))
result.sort(key=itemgetter(2),reverse=True)
    

treat_uu = 0
control_uu = 0
treat_cv = 0
control_cv = 0
treat_cvr = 0.0
control_cvr = 0.0
lift = 0.0

stat_data = []

for is_cv, is_treat, score in result:
    if is_treat:
        treat_uu += 1
        if is_cv:
            treat_cv += 1
        treat_cvr = treat_cv / treat_uu
    else:
        control_uu += 1
        if is_cv:
            control_cv += 1
        control_cvr = control_cv / control_uu

    # コンバージョンレートの差に実験群の人数を掛けることでliftを算出
    lift = (treat_cvr - control_cvr) * treat_uu

    stat_data.append([is_cv, is_treat, score, treat_uu, control_uu, treat_cv, control_cv, treat_cvr, control_cvr, lift])


qdf = pd.DataFrame(columns=('treat_cvr', 'control_cvr'))

quantile_data = []
for n in range(10):
    start = int(n * len(result) / 10)
    end = int((n + 1) * len(result) / 10) - 1
    quantiled_result = result[start:end]
    
    treat_uu = list(map(lambda item:item[1], quantiled_result)).count(True)
    control_uu = list(map(lambda item:item[1], quantiled_result)).count(False)
    
    treat_cv = [item[0] for item in quantiled_result if item[1] == True].count(True)
    control_cv = [item[0] for item in quantiled_result if item[1] == False].count(True)
    
    treat_cvr = treat_cv / treat_uu
    control_cvr = control_cv / control_uu
    
    
    quantile_data.append([treat_uu, control_uu, treat_cv, control_cv, treat_cvr, control_cvr])
        
    label = "{}%~{}%".format(n*10, (n+1)*10)
    qdf.loc[label] = [treat_cvr, control_cvr]

qdf.plot.bar()
plt.xlabel("percentile")
plt.ylabel("conversion rate")

df = pd.DataFrame(stat_data)
df.columns = ["is_cv", "is_treat", "score", "treat_uu", "control_uu", "treat_cv", "control_cv", "treat_cvr", "control_cvr", "lift"]


# ベースラインを書き加える
df["base_line"] = df.index * df["lift"][len(df.index) - 1] / len(df.index)

df.plot(y=["treat_cv", "control_cv"])
plt.xlabel("uplift score rank")
plt.ylabel("conversion count")

df.plot(y=["treat_cvr", "control_cvr"], ylim=[0, 0.04])
plt.xlabel("uplift score rank")
plt.ylabel("conversion rate")

df.plot(y=["lift", "base_line"])
plt.xlabel("uplift score rank")
plt.ylabel("lift count")

df.plot(y=["treat_cv", "control_cv"], x="score", title="conversion count")
df.plot(y=["treat_cvr", "control_cvr"], ylim=[0, 0.04], x="score", title="conversion rate")
df.plot(y=["lift", "base_line"], x="score", title="lift")

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

%matplotlib inline

import pandas as pd

import matplotlib.pyplot as plt

from operator import itemgetter

plt.style.use("ggplot")

from sklearn.model_selection import train_test_split

from sklearn.linear_model import LogisticRegression

# https://ailab.criteo.com/criteo-uplift-prediction-dataset/

source_df = pd.read_csv("criteo-uplift.csv.gz")

source_df.head(10)

source_df.describe()

feature_vector_df = source_df.drop(["treatment","conversion","visit","exposure"],axis=1)

is_treat_list = list(source_df["treatment"] == 1)

is_cv_list = list(source_df["conversion"] == 1)

train_is_cv_list, test_is_cv_list, train_is_treat_list, \

test_is_treat_list, train_feature_vector_df,\

test_feature_vector_df = train_test_split(is_cv_list, is_treat_list,

feature_vector_df,

train_size=0.5,

test_size=0.5,

random_state=42)

treat_model = LogisticRegression(C=0.01)

control_model = LogisticRegression(C=0.01)

train_sample_num = len(train_is_cv_list)

treat_is_cv_list = [train_is_cv_list[i] for i in range(train_sample_num) if train_is_treat_list[i] == True]

treat_feature_vector_list = train_feature_vector_df[train_is_treat_list]

control_is_cv_list = [train_is_cv_list[i] for i in range(train_sample_num) if train_is_treat_list[i] == False]

control_feature_vector_list = train_feature_vector_df[list(map(lambda a:a == False ,train_is_treat_list))]

treat_model.fit(treat_feature_vector_list, treat_is_cv_list)

control_model.fit(control_feature_vector_list, control_is_cv_list)

treat_score = treat_model.predict_proba(test_feature_vector_df)

control_score = control_model.predict_proba(test_feature_vector_df)

score_list = treat_score[:,1] / control_score[:,1]

result = list(zip(test_is_cv_list, test_is_treat_list, score_list))

result.sort(key=itemgetter(2),reverse=True)

treat_uu = 0

control_uu = 0

treat_cv = 0

control_cv = 0

treat_cvr = 0.0

control_cvr = 0.0

lift = 0.0

stat_data = []

for is_cv, is_treat, score in result:

if is_treat:

treat_uu += 1

if is_cv:

treat_cv += 1

treat_cvr = treat_cv / treat_uu

else:

control_uu += 1

if is_cv:

control_cv += 1

control_cvr = control_cv / control_uu

# コンバージョンレートの差に実験群の人数を掛けることでliftを算出

lift = (treat_cvr - control_cvr) * treat_uu

stat_data.append([is_cv, is_treat, score, treat_uu, control_uu, treat_cv, control_cv, treat_cvr, control_cvr, lift])

qdf = pd.DataFrame(columns=('treat_cvr', 'control_cvr'))

quantile_data = []

for n in range(10):

start = int(n * len(result) / 10)

end = int((n + 1) * len(result) / 10) - 1

quantiled_result = result[start:end]

treat_uu = list(map(lambda item:item[1], quantiled_result)).count(True)

control_uu = list(map(lambda item:item[1], quantiled_result)).count(False)

treat_cv = [item[0] for item in quantiled_result if item[1] == True].count(True)

control_cv = [item[0] for item in quantiled_result if item[1] == False].count(True)

treat_cvr = treat_cv / treat_uu

control_cvr = control_cv / control_uu

quantile_data.append([treat_uu, control_uu, treat_cv, control_cv, treat_cvr, control_cvr])

label = "{}%~{}%".format(n*10, (n+1)*10)

qdf.loc[label] = [treat_cvr, control_cvr]

qdf.plot.bar()

plt.xlabel("percentile")

plt.ylabel("conversion rate")

df = pd.DataFrame(stat_data)

df.columns = ["is_cv", "is_treat", "score", "treat_uu", "control_uu", "treat_cv", "control_cv", "treat_cvr", "control_cvr", "lift"]

# ベースラインを書き加える

df["base_line"] = df.index * df["lift"][len(df.index) - 1] / len(df.index)

df.plot(y=["treat_cv", "control_cv"])

plt.xlabel("uplift score rank")

plt.ylabel("conversion count")

df.plot(y=["treat_cvr", "control_cvr"], ylim=[0, 0.04])

plt.xlabel("uplift score rank")

plt.ylabel("conversion rate")

df.plot(y=["lift", "base_line"])

plt.xlabel("uplift score rank")

plt.ylabel("lift count")

df.plot(y=["treat_cv", "control_cv"], x="score", title="conversion count")

df.plot(y=["treat_cvr", "control_cvr"], ylim=[0, 0.04], x="score", title="conversion rate")

df.plot(y=["lift", "base_line"], x="score", title="lift")

こちらの図はアップリフト値の階級値ごとのCVRです。最上位のアップリフト値はCVRの差が大きいですが、上位40~50%程度のアップリフト値のときにCVRの差が最も大きいようです。

アップリフト値の順位とCVRの図です。順位が低くても処置群のほうがCVRがわずかに高いようです。

アップリフトのスコアとCVRの関係です。2未満であればCVRは処置群が上回っていますが、一様な傾向はなさそうです。

コンバージョンレートの差に対象群の人数を掛けることでliftを算出したものです。アップリフトスコアが1~2点であれば儲かるようです。

tools4upliftと出している指標が違うので比較ができないのが難点に思いました。tools4upliftはオートマチックな感じで便利なのですが、『仕事ではじめる機械学習』の9章を正義として進めたいので、どうにか揃えれるようにしていきたいと思います。

おわりに

tools4upliftというマーケターにとって銀の弾丸になりそうなパッケージの存在を知ることができ、実際に非常に便利そうな関数が用意されているのがわかりました。ただ、開発されたばかりのパッケージなのでそこまで結果を信じていません。『仕事ではじめる機械学習』本の結果と揃えたいなと思いました。その点がはっきりすれば業務で使ってみるのも良いですし、任意のマーケターに安心して共有できると思います。

参考文献

[1] 有賀康顕・中山心太・西林孝 (2018) 『仕事ではじめる機械学習』オライリージャパン
[2] Mouloud Belbahri, Alejandro Murua, Olivier Gandouet, Vahid Partovi Nia (2019). “Uplift Regression: The R Package tools4uplift”, arXiv:1901.10867 [stat.AP]
[3] ohke (2019) 「Uplift modelingで施策が効く人を見極める」け日記
[4] usaito (2018) 「Uplift Modelingで介入効果を最適化する」 Qiita

ABEJA SIX 2019の1日目に行ってきましたレポート

今日は午後から有給をいただいて、ABEJA SIXの1日目に行ってきました。印象的だなと感じたものに関して、つらつらと雑記を載せておきたいと思います。

こちらは品川グランドプリンスホテルの庭園です。

こちらは会場の雰囲気です。

ブースの様子1です。

ブースの様子2です。

ABEJA SIX 2019

「食事画像認識モデル開発プロジェクトでの10個5個の教訓」株式会社FiNC Technologies 南野充則氏

FiNCは450万ダウンロードされているヘルスケア系のアプリを開発している会社。
ユーザーの継続率を高めるための施策として、機械学習を用いている。
今回の紹介事例ではユーザーの食事に関する情報を入力する手間を機械学習で短縮させ、短縮させることで継続率を高めることを狙っている。
食事の画像は1日に数万枚がアプリに投稿される。
食事の画像から栄養価などを計算することを目指している。
食事レシピ認識モデルでは、画像からレシピを識別し、メニューの量（グラム数）なども推定し、カテゴリ単位で決まっている栄養価から推定している。レシピ本の情報を入力したり、レシピサイトをクローリングし、レシピを一人あたりの栄養価になるように標準化などもしている。きれいな画像と栄養価（材料何グラムか）の伴ったクリーンなデータセットを用意するために自社のキッチンに料理人を呼び2000レシピ分の料理を作ったとのこと。
食材認識モデルでは食材一つ一つ（トマト一つとか、キャベツ一枚とか）を識別して、栄養価を素材単位で計算している。
学習の結果、管理栄養士よりも3%程度の誤差でメニューの栄養価を推定可能になった。
開発期間は6ヶ月間。
東大の松尾研にアドバイスをもらっているらしい。

5つの教訓

１．DL/ML人材をソフトウェアエンジニアから輩出すべき：
インフラ、サーバー、DB、パフォーマンスなどに明るいソフトウェアエンジニアが機械学習や深層学習を学ぶと、分析も実装もできる頼もしいメンバーになるので、ソフトウェアエンジニアのデータサイエンティスト化に注力しているらしい。目指すは論文のリプリケーションができるレベルとのこと。
2.データ取得から学習までのPDCAを最速にする：
ユーザーが画像を出したあとのフローをしっかりしていなかった。予期せぬデータが入ってくるので、そこへの対応も必要。アノテーションした項目を再学習するような仕組みを作り、そばの画像が苦手であれば、そばの画像を集中的に集めて学習させる。
3.オペレーションは自社で構築せよ：
泥臭い仕事と思い、丸投げしてはいけない。データセットの質が最も大事。データセットの質を担保するには評価手法を理解し細かいオペレーションを作る必要がある。アルバイトも自社で雇用、マネジャーもエンジニアとすることで当事者意識も芽生えやすい。
4.評価方法の決定からプロジェクトを始めよう：
AIを使えば、想像を超える何かが出てくると期待していまうフシがある。評価の仕方を決めたほうが、メンバーのゴールが見えるし。やりやすい。10%以内の誤差の難易度がどの程度なのかわからなかったりするし、解釈の多様性が生まれてしまうこともある。
5.プロジェクトはアジャイルで進めるべき：
作ったことのないモデルを作る際にスケジューリングを引くことは難しい。SOTAくらいいけますよと言ってしまい、自らを苦しめることになりかねない。

「機械学習におけるクラウド活用のポイント」アマゾンウェブサービスジャパン株式会社針原佳貴氏 & 宇都宮聖子氏

SageMakerいいぞというお話。
ビジネスにおいて、機械学習を進めるに際して重要なポイントは、
「ビジネス価値に落とし込む」
「データの流れを理解する」
「自分の力で頑張らない」
の3点が挙げられていた。
必要ではあるが、付加価値にはつながりにくい作業のことをUndifferentiated heavy liftingと呼ぶらしい。
機械学習プロジェクトを回す際に重要なこととして、
データ取得
↓
データ前処理
↓
モデルの開発・学習
↓
モデルの評価
↓
モデルの変換（エッジデバイスに送るにはデータを小さくする必要がある。）
↓
本番環境のデプロイ
↓
監視・評価データ変換
のループを繰り返すことが挙げられている。
S3(Simple Storage Service)に蓄積しているデータがあったとして、そのデータに対して、SageMakerで前処理やら機械学習を行い、学習済みの結果をS3にためれば、それを用いてエンドポイントの推論としてカジュアルに活用することができる。S3→SageMaker→S3のコンボが良いとのこと。
ここ1年間で200個くらいAWSのサービスやら機能が増えているので、それを知るだけでも大変そう。でもうまく使えば、Undifferentiated heavy liftingを避けることができる。
わからないことがあれば、ソリューションアーキテクトに質問したり、SageMakerのSlackで聞いたりすると良いらしい。
SageMakerでの学習の進め方としては3種類ある。１つ目は、TensorFlowなどでゴリゴリとアルゴリズムを書く。2つ目はAWS Marketplaceで販売されているアルゴリズムを時間単位で課金して使う。3つ目はAWSのビルトインのアルゴリズム（Object Detection、Semantic Segmentation、Factorization Machineなど）を使う。

「少数データからの学習法の展開とABEJAの取り組み」株式会社ABEJA 藤本敬介氏

データの質がモデルの結果を左右するが、きれいなデータを大量に集めるためにアノテーションをやるのは大変。少ないデータでも性能を出したい。
アプローチとしては、Data Augmentation、Transfer Learning、Meta learningの3つがある。

Data Augmentation（データ拡張）

データを擬似的に増やす手法
mixup:Beyond Empirical Risk Minimizationの研究では画像とラベルを合成してデータを水増しが行われている。
AutoAugment:LearningAugmentation Policies from Dataの研究では、強化学習で精度が上がりやすい水増し方法を学習している。
使わない手はない。

Transfer Learning（転移学習）

異なるデータセットで学習したものを再利用する。
Fine-tuning：別のデータで学習済みのモデルに対して、タスクに対してのデータに適用する。
Domain Adaptation：学習済みのモデルやデータの知識を再利用する。
Fine-tuningは有効な手段。

Meta learning

タスクの学習のしかたを学習する
少数のデータでのうまい学習方法を訓練しておいて、それを使い回す。

ABEJAの取り組み

データが少ない場合はFine-tuningで高精度を出しやすい。
External Network：中間層の情報を利用して、例外的な処理（ネットワークにバイパスみたいなものを通す）をすることで、Fine-tuningした際に精度が落ちないようにしている。不均衡データやクラス追加に対して強い手法とされている。データ数に応じてExternal Networkのサイズを調整でき、クラス1に大量のデータがある場合、1だけネットワークを深くして、2やら3はネットワークを浅くするなどの柔軟な対応が可能。これでもって不均衡データに対応できるとのこと。また、クラス追加に関しては、追加したクラスの分だけ学習すればいいようにネットワークの学習ができるらしい。ただし、学習に時間がかかるとのこと。
（よくわからないが）Model-Agnostic Meta-Learning（MAML（マムル））を応用したら精度が高まるらしい。

うーん、DNNは全然追いかけれていないので断片的にしかわからなかった。悔しいものです。

「Deep Learningの都市伝説と現実」株式会社ABEJA 白川達也氏

リサーチャーをする上で大事なこととしては、
1.先に見つけること
2.シンプルに解くこと
3.先に失敗する（大きな失敗は会社としてしないために）
の3つがある。
クリーンなデータで学習したほうが精度が高くなりやすく、過学習しにくい。ラベルの精度が高ければ、高いほどよい。Big Clean Data + DLで勝つる？
アノテーションは簡単ではない。アノテーターごとにわかりやすい情報がバラバラで、ブレるのが本質的。どこまでやるのか、どこが基準なのかというフレーミングとアンカーリングが重要。人間とかタスクを理解してすすめるのが良い。
半教師あり学習（アノテーションされていないデータを使って精度向上させる取り組み）も魅力的だが、教師データを増やしたほうが効率的。アノテーションできるならば、アノテーションしてしまおう。事前学習も意味があるので行う。
次にどんな技術がくるのか？ Graph Convolution、Annotation、Poincare Embeddings、ML in Hyperbolic Space
Taskonomyという研究が今後熱くなるかも。見たこともないタスクも解けるという柔軟性を持つモデルが構築できる？

感想

機械学習で精度を出すためにそこまで頑張るのか！という事例を聞けたり、知識として不足していたAWS系のサービスの話を聞けたり、自分の足りていない知識を補えた良いイベントだと思いました。

RのContextualパッケージをいじってみた際のメモ書き

はじめに

このブログの私の中での位置づけは、今後仕事で使いそうなものを調べて書き溜めるというところにあります。仕事で使っているものはブログに載せないというスタンスでもあるのですが、出来るだけ先回りしておきたいところです。今回は、昨年のJapan.RやTokyo.Rで紹介されていたcontextualパッケージを触ってみたというゆるふわな内容となっています。

・バンディット問題とは
・マーケティング関連でバンディット問題が役に立つ場面
・バンディット問題で出てくる数学的な知識と方策
・Contextual Bandit問題とは
・Contextualパッケージでできること
・サンプル実行
・おわりに
・参考情報

バンディット問題とは

「選択肢の集合から1つの要素を選択して、その選択肢に対する報酬を得るものの、他の選択肢の報酬情報は得られないというプロセスを繰り返す設定において、報酬の合計値を最大化することを目指す逐次決定問題」とされています。バンディットは昔ながらのスロットマシンが客からお金をむしり取ること（盗賊）にちなんでいるそうです。胴元は盗賊ということなんでしょうか？

大学時代の知人は毎日パチンコ屋に行ってから講義に行っていましたが、出そうな台・出そうな店を転々としていましたが、あれはバンディット問題を彼なりに解いていたのでしょう。当時はサクラの台というのがあったらしく、3000円ほど投資すれば大当たりになるのだとか。そしてその大当たりに釣られて他の客が頑張るという意味で、サクラの台だそうです。

マーケティング関連でバンディット問題が役に立つ場面

私はマーケティング×データ分析を生業としているので、マーケティング方面にしか関心がないのですが、バンディット問題は役立つ可能性が十分にあるというか既に一部の企業ではバリューを出しています。

・インターネット広告配信：オレシカナイトでSpeeeの方がトンプソン抽出で精度を増していた。
・推薦システムにおけるコールドスタート問題：ネットフリックスが情報推薦の際にContextual Banditを適用

バンディット問題とは異なるものの、最適腕識別問題においては、クックパッドのクリエイティブ出し分けやGoogleのウェブテスト（旧Webサイトオプティマイザー）などで使われています。ちなみに、バンディット問題と最適腕識別問題は似て非なるものであるということを『バンディット問題の理論とアルゴリズム』で知りました。

また、マーケティングとは違いますが、株価のトレーディングの際にバンディットアルゴリズムを使っているという事例（Bandits and Stocks）が当然ながらあるようです。

バンディット問題で出てくる数学的な知識と方策

バンディット問題の書籍を読もうとすると、数理統計学の知識が必要です。

あるスロットを何回引くべきかという意思決定の際に、「神のみぞ知る真の報酬」と「あるスロットの報酬」がどれくらい外れているか、そしてそのハズレ具合は許容できるのかということが重要になります。
「神のみぞ知る真の報酬と、あるスロットの報酬がΔだけ外れている確率」の推論の精度に関心があるということです。

バンディット問題において、「その時のベストのスロットを引いた際のリターン」と「その時実際に選んだスロットのリターン」の差の期間合計値をリグレットとして、そのリグレットを小さくするようにスロットを選びます。
そのリグレットに対して理論的な下限を求める際に、数理統計学の知識が必要になります。

具体的には、ヘフディングの不等式、その前提となるマルコフの不等式やチェビシェフの不等式やチェルノフ限界、積率母関数やイェンセンの不等式などです。
それらを駆使しながら、様々な施策の中で、理論的な下限がより小さくなるようなものを探そうという流れのようです。

『バンディット問題の理論とアルゴリズム』を読む上で前提となっていそうな知識として、スタンフォード大学の講義資料（CS229 Supplemental Lecture notes Hoeffding’s inequality）を運良く見つけることが出来たので、これをもとに学ぶと理解が捗ると思います。

リグレットの下限を低めることを目指して、様々なアプローチが議論されます。

ε-貪欲法

概要：スロットを回す回数のうち、一定割合（ε）をスロットの探索に当て、残りの期間を良いとされるスロットを回し続ける。
メリット：実装が容易でシステムに組み込み易い
デメリット：期待値が悪いスロットも良いスロットも同じ回数引いてしまうので性能が悪くなる。スロットの種類が多い際はより一層悪くなりやすい。

UCB(Upper Confidence Bound)方策

概要：標本平均に補正項を足した、UCBスコアを各時点ごとに計算し、最もスコアが高いスロットを回す。なお、補正項は選択回数の少ないスロットに対して大きくなります。
メリット：ε-貪欲法と異なり、リグレットの上限がεなどの水準に左右されない。ハイパーパラメータが少ない。
デメリット：真の期待値についての信頼区間を求めることは本質的ではない。

KL-UCB

概要：KLダイバージェンスを用いてUCBスコアを計算し、最もスコアが高いスロットを回す。
メリット：KLダイバージェンスを様々なモデルに応じて置き換えることができるなど、柔軟性がある。
デメリット：KLダイバージェンスの逆関数を計算する必要があり、毎回ニュートン法などを適用する必要がある。

MED(Minimum Empirical Divergence)方策

概要：期待値最大である際の尤度が一定以上のスロットを回すという方策。
メリット：KLダイバージェンスの逆関数を計算する必要がない。
デメリット：KL-UCBよりも性能が悪い。IMEDという方策であればその弱点を克服している。

トンプソン抽出

概要：期待値最大でないスロットの選択数の期待値を近似的に最小化するという取り組みを、ベイズ統計の枠組みで行ったもの。
メリット：経験的に高い性能となりやすい。
デメリット：？

Contextual Bandit問題とは

ある時点のあるスロットの報酬が、ユーザーの特徴量と誤差項により線形で表すことができるものを、線形バンディットと呼びます。
ユーザーの各行動の特徴量が時刻により異なる値を取ることを許すという設定を、文脈付きバンディット（Contextual Bandit）と呼びます。
つまり、Contextual Banditは時刻により異なるユーザーの特徴量が与えられたもとでの、利得の期待値の最大化問題となります。

具体的には、パチンコ店における期待値最大化の行動を考えるとすると、パチンコ台の大当たり確率は、午前か午後か、大当たりが既に他の台で出たか、その台がどれくらい回されているかなどの時間による文脈に左右されるという状況となります。

このContextual Banditにおいても、先程あげたようなリグレットを最小にするような様々な方策があります。LinUCB方策や、線形モデルのトンプソン抽出、ロジスティック回帰モデルのバンディットなどです。

Contextualパッケージでできること

こちらの資料にある通り、バンディットアルゴリズムのシミュレーションとオフライン評価が行えるパッケージです。
多様なバンディットアルゴリズムを試すことができます。
要となるデータですが、シミュレーションにより生成することもできれば、過去にランダムに出し分けたログなどのデータがあればそのデータをもとにアルゴリズムの検証をすることができます。

サンプル実行

さて、今回は完全に手抜きです。GitHubにあったサンプルコードを3つほど回すだけです。ただ、特徴量の突っ込み方などをサンプルコードから学べるので、ぜひ開発者のGitHubをご覧ください。

サンプル1：ABテストによる最適腕選択

パッケージのGitHub
にコードがありました。Bandit Algorithms for Website Optimizationという書籍に登場してきている例をRで実行できるサンプルです。
・ε-貪欲法を様々なεでシミュレーションして最適なスロットを見つける
・ソフトマックスによる方策に関しても様々なτに応じたシミュレーションをして最適なスロットを見つける
・UCB方策によりシミュレーションを行い、最適なスロットを見つける。ε-貪欲法やソフトマックスとの比較を行う
という実験ができます。シミュレーションの設定として、スロットごとの当たりの出る確率をベクトルで指定しています。

実行するのに10分くらいはかかるかもしれません。

library(contextual)

# Bandit algorithms for website optimization -----------------------------------------------------------------

## Simulation of the multi-armed Bandit examples in
## of "Bandit algorithms for website optimization"
## by John Miles White.

# The code from the book chooses the arm with the first index when all arms are equal.
# Contextuals policies correctly picks one of the max arms.
# That's why the plots below are slightly different from the book - they are correct, though.

# Chapter 4 - Debugging and epsilon greedy -------------------------------------------------------------------

prob_per_arm       <- c(0.1, 0.1, 0.1, 0.1, 0.9)
horizon            <- 250
simulations        <- 5000

bandit             <- BasicBernoulliBandit$new(prob_per_arm)

agents             <- list(Agent$new(EpsilonGreedyPolicy$new(0.1), bandit, "Epsilon = 0.1"),
                           Agent$new(EpsilonGreedyPolicy$new(0.2), bandit, "Epsilon = 0.2"),
                           Agent$new(EpsilonGreedyPolicy$new(0.3), bandit, "Epsilon = 0.3"),
                           Agent$new(EpsilonGreedyPolicy$new(0.4), bandit, "Epsilon = 0.4"),
                           Agent$new(EpsilonGreedyPolicy$new(0.5), bandit, "Epsilon = 0.5"))

simulation         <- Simulator$new(agents, horizon, simulations)
history            <- simulation$run()

# Figure 4-2. How often does the epsilon greedy algorithm select the best arm?

plot(history, type = "optimal", legend_position = "bottomright", ylim = c(0,1))

# Figure 4-3. How much reward does the epsilon greedy algorithm earn on average?

plot(history, type = "average", regret = FALSE, legend_position = "bottomright", ylim = c(0,1))

# Figure 4-4. How much reward has the epsilon greedy algorithm earned by trial t?

plot(history, type = "cumulative", regret = FALSE)

# Chapter 5 - Softmax ----------------------------------------------------------------------------------------

agents             <- list(Agent$new(SoftmaxPolicy$new(0.1), bandit, "Tau = 0.1"),
                           Agent$new(SoftmaxPolicy$new(0.2), bandit, "Tau = 0.2"),
                           Agent$new(SoftmaxPolicy$new(0.3), bandit, "Tau = 0.3"),
                           Agent$new(SoftmaxPolicy$new(0.4), bandit, "Tau = 0.4"),
                           Agent$new(SoftmaxPolicy$new(0.5), bandit, "Tau = 0.5"))

simulation         <- Simulator$new(agents, horizon, simulations)
history            <- simulation$run()

# Figure 5-2. How often does the softmax algorithm select the best arm?

plot(history, type = "optimal", legend_position = "bottomright", ylim = c(0,1))

# Figure 5-3. How much reward does the softmax algorithm earn on average?

plot(history, type = "average", regret = FALSE, legend_position = "bottomright", ylim = c(0,1))

# Figure 5-4. How much reward has the softmax algorithm earned by trial t?

plot(history, type = "cumulative", regret = FALSE)


# Chapter 6 - UCB --------------------------------------------------------------------------------------------

agents             <- list(Agent$new(SoftmaxPolicy$new(0.1), bandit, "Softmax"),
                           Agent$new(EpsilonGreedyPolicy$new(0.1), bandit, "EpsilonGreedy"),
                           Agent$new(UCB1Policy$new(), bandit, "UCB1"))

simulation         <- Simulator$new(agents, horizon, simulations)
history            <- simulation$run()

# Figure 6-3. How often does the UCB algorithm select the best arm?

plot(history, type = "optimal", legend_position = "bottomright", ylim = c(0,1))

# Figure 6-4. How much reward does the UCB algorithm earn on average?

plot(history, type = "average", regret = FALSE, legend_position = "bottomright", ylim = c(0,1))

# Figure 6-5. How much reward has the UCB algorithm earned by trial t?

plot(history, type = "cumulative", regret = FALSE)

library(contextual)

# Bandit algorithms for website optimization -----------------------------------------------------------------

## Simulation of the multi-armed Bandit examples in

## of "Bandit algorithms for website optimization"

## by John Miles White.

# The code from the book chooses the arm with the first index when all arms are equal.

# Contextuals policies correctly picks one of the max arms.

# That's why the plots below are slightly different from the book - they are correct, though.

# Chapter 4 - Debugging and epsilon greedy -------------------------------------------------------------------

prob_per_arm <- c(0.1, 0.1, 0.1, 0.1, 0.9)

horizon <- 250

simulations <- 5000

bandit <- BasicBernoulliBandit$new(prob_per_arm)

agents <- list(Agent$new(EpsilonGreedyPolicy$new(0.1), bandit, "Epsilon = 0.1"),

Agent$new(EpsilonGreedyPolicy$new(0.2), bandit, "Epsilon = 0.2"),

Agent$new(EpsilonGreedyPolicy$new(0.3), bandit, "Epsilon = 0.3"),

Agent$new(EpsilonGreedyPolicy$new(0.4), bandit, "Epsilon = 0.4"),

Agent$new(EpsilonGreedyPolicy$new(0.5), bandit, "Epsilon = 0.5"))

simulation <- Simulator$new(agents, horizon, simulations)

history <- simulation$run()

# Figure 4-2. How often does the epsilon greedy algorithm select the best arm?

plot(history, type = "optimal", legend_position = "bottomright", ylim = c(0,1))

# Figure 4-3. How much reward does the epsilon greedy algorithm earn on average?

plot(history, type = "average", regret = FALSE, legend_position = "bottomright", ylim = c(0,1))

# Figure 4-4. How much reward has the epsilon greedy algorithm earned by trial t?

plot(history, type = "cumulative", regret = FALSE)

# Chapter 5 - Softmax ----------------------------------------------------------------------------------------

agents <- list(Agent$new(SoftmaxPolicy$new(0.1), bandit, "Tau = 0.1"),

Agent$new(SoftmaxPolicy$new(0.2), bandit, "Tau = 0.2"),

Agent$new(SoftmaxPolicy$new(0.3), bandit, "Tau = 0.3"),

Agent$new(SoftmaxPolicy$new(0.4), bandit, "Tau = 0.4"),

Agent$new(SoftmaxPolicy$new(0.5), bandit, "Tau = 0.5"))

simulation <- Simulator$new(agents, horizon, simulations)

history <- simulation$run()

# Figure 5-2. How often does the softmax algorithm select the best arm?

plot(history, type = "optimal", legend_position = "bottomright", ylim = c(0,1))

# Figure 5-3. How much reward does the softmax algorithm earn on average?

plot(history, type = "average", regret = FALSE, legend_position = "bottomright", ylim = c(0,1))

# Figure 5-4. How much reward has the softmax algorithm earned by trial t?

plot(history, type = "cumulative", regret = FALSE)

# Chapter 6 - UCB --------------------------------------------------------------------------------------------

agents <- list(Agent$new(SoftmaxPolicy$new(0.1), bandit, "Softmax"),

Agent$new(EpsilonGreedyPolicy$new(0.1), bandit, "EpsilonGreedy"),

Agent$new(UCB1Policy$new(), bandit, "UCB1"))

simulation <- Simulator$new(agents, horizon, simulations)

history <- simulation$run()

# Figure 6-3. How often does the UCB algorithm select the best arm?

plot(history, type = "optimal", legend_position = "bottomright", ylim = c(0,1))

# Figure 6-4. How much reward does the UCB algorithm earn on average?

plot(history, type = "average", regret = FALSE, legend_position = "bottomright", ylim = c(0,1))

# Figure 6-5. How much reward has the UCB algorithm earned by trial t?

plot(history, type = "cumulative", regret = FALSE)

ε-貪欲法

・最適なスロットを選んだ確率

・平均報酬額

　・累積報酬額

ソフトマックスによる方策

・最適なスロットを選んだ確率

・平均報酬額

・累積報酬額

UCB方策

・最適なスロットを選んだ確率

・平均報酬額

・累積報酬額

サンプル2：文脈付きバンディット問題で映画のレーティングの最適化

同じGitHubにあるこちらのコードは、映画のデータセットに対して、文脈付きバンディット問題でオフラインテストをするためのコードです。映画のレーティングが4以上なら1そうでないなら0のデータを作り、特徴量として映画館で見たか家で見たか、一人で見たか家族と見たか、週末に見たかどうかなどの変数を7個ほど作成しています。方策としては、ランダムなもの、ε-貪欲法、トンプソン抽出、LinUCBをシミュレーションしています。

実行してから処理が止まるまで1時間程度はかかりましたが、LinUCBが累積の報酬が大きいようです。

library(contextual)
library(data.table)

# Import personalization data-set

# Info: https://d1ie9wlkzugsxr.cloudfront.net/data_irecsys_CARSKit/Movie_DePaulMovie/README.txt

url         <- "http://d1ie9wlkzugsxr.cloudfront.net/data_irecsys_CARSKit/Movie_DePaulMovie/ratings.csv"
data        <- fread(url, stringsAsFactors=TRUE)

# Convert data

data        <- contextual::one_hot(data, cols = c("Time","Location","Companion"), sparsifyNAs = TRUE)
data[, itemid := as.numeric(itemid)]
data[, rating := ifelse(rating <= 3, 0, 1)]

# Set simulation parameters.
simulations <- 10  # here, "simulations" represents the number of boostrap samples
horizon     <- nrow(data)

# Initiate Replay bandit with 10 arms and 100 context dimensions
log_S       <- data
formula     <- formula("rating ~ itemid | Time_Weekday + Time_Weekend + Location_Cinema + Location_Home +
                                          Companion_Alone + Companion_Family + Companion_Partner")
bandit      <- OfflineBootstrappedReplayBandit$new(formula = formula, data = data)

# Define agents.
agents      <-
  list(Agent$new(RandomPolicy$new(), bandit, "Random"),
       Agent$new(EpsilonGreedyPolicy$new(0.03), bandit, "EGreedy 0.05"),
       Agent$new(ThompsonSamplingPolicy$new(), bandit, "ThompsonSampling"),
       Agent$new(LinUCBDisjointOptimizedPolicy$new(0.37), bandit, "LinUCB 0.37"))

# Initialize the simulation.
simulation  <-
                Simulator$new(
                  agents           = agents,
                  simulations      = simulations,
                  horizon          = horizon
                )

# Run the simulation.
# Takes about 5 minutes: bootstrapbandit loops for arms x horizon x simulations (times nr of agents).
sim  <- simulation$run()

# plot the results
plot(sim, type = "cumulative", regret = FALSE, rate = TRUE,
     legend_position = "topleft", ylim=c(0.48,0.87))

library(contextual)

library(data.table)

# Import personalization data-set

# Info: https://d1ie9wlkzugsxr.cloudfront.net/data_irecsys_CARSKit/Movie_DePaulMovie/README.txt

url <- "http://d1ie9wlkzugsxr.cloudfront.net/data_irecsys_CARSKit/Movie_DePaulMovie/ratings.csv"

data <- fread(url, stringsAsFactors=TRUE)

# Convert data

data <- contextual::one_hot(data, cols = c("Time","Location","Companion"), sparsifyNAs = TRUE)

data[, itemid := as.numeric(itemid)]

data[, rating := ifelse(rating <= 3, 0, 1)]

# Set simulation parameters.

simulations <- 10 # here, "simulations" represents the number of boostrap samples

horizon <- nrow(data)

# Initiate Replay bandit with 10 arms and 100 context dimensions

log_S <- data

formula <- formula("rating ~ itemid | Time_Weekday + Time_Weekend + Location_Cinema + Location_Home +

Companion_Alone + Companion_Family + Companion_Partner")

bandit <- OfflineBootstrappedReplayBandit$new(formula = formula, data = data)

# Define agents.

agents <-

list(Agent$new(RandomPolicy$new(), bandit, "Random"),

Agent$new(EpsilonGreedyPolicy$new(0.03), bandit, "EGreedy 0.05"),

Agent$new(ThompsonSamplingPolicy$new(), bandit, "ThompsonSampling"),

Agent$new(LinUCBDisjointOptimizedPolicy$new(0.37), bandit, "LinUCB 0.37"))

# Initialize the simulation.

simulation <-

Simulator$new(

agents = agents,

simulations = simulations,

horizon = horizon

)

# Run the simulation.

# Takes about 5 minutes: bootstrapbandit loops for arms x horizon x simulations (times nr of agents).

sim <- simulation$run()

# plot the results

plot(sim, type = "cumulative", regret = FALSE, rate = TRUE,

legend_position = "topleft", ylim=c(0.48,0.87))

サンプル3：文脈付きバンディット問題でMovieLensのTop50の作品における評価の最適化

こちらのコードは、MovieLensのデータセットにおいて、特徴量として過去にユーザーが評価した映画のカテゴリーの割合を19カテゴリ分用意して、ユーザーの見た映画の評価を最も高めるという、文脈付きバンディット問題です。こちらは実行して、30分程度で処理が終わりました。先程のサンプルと同じで、LinUCBが累積の報酬が大きいようです。

library(contextual)
library(data.table)
library(splitstackshape)

# Movielens 100k ---------------------------------------------------------------------------------------------

# Info: https://d1ie9wlkzugsxr.cloudfront.net/data_movielens/ml-100k/ml-100k-README.txt

movies_dat      <- "http://d1ie9wlkzugsxr.cloudfront.net/data_movielens/ml-100k/u.item"
ratings_dat     <- "http://d1ie9wlkzugsxr.cloudfront.net/data_movielens/ml-100k/u.data"

# Import and merge files

movies_dat      <- fread(movies_dat, sep = "|", quote="")
setnames(movies_dat, c("V1", "V2"), c("MovieID", "Name"))
movies_dat[, (3:5)  := NULL ]
ratings_dat     <- fread(ratings_dat, quote="")
setnames(ratings_dat, c("V1", "V2", "V3", "V4"), c("UserID", "MovieID", "Rating", "Timestamp"))
all_movies      <- ratings_dat[movies_dat, on=c(MovieID = "MovieID")]

rm(movies_dat,ratings_dat)

# Data wrangling ---------------------------------------------------------------------------------------------

count_movies    <- all_movies[,.(MovieCount = .N), by = MovieID]
top_50          <- as.vector(count_movies[order(-MovieCount)][1:50]$MovieID)
not_50          <- as.vector(count_movies[order(-MovieCount)][51:nrow(count_movies)]$MovieID)
top_50_movies   <- all_movies[MovieID %in% top_50]

# User features: tags they've watched for non-top-50 movies normalized per user

user_features   <- all_movies[MovieID %in% not_50]
rm(all_movies)
user_features[, c("MovieID", "Rating", "Timestamp", "Name"):=NULL]
user_features   <- user_features[, lapply(.SD, sum, na.rm=TRUE), by=UserID ]
user_features[, total := rowSums(.SD, na.rm = TRUE), .SDcols = 2:20]
user_features[, 2:20 := lapply(.SD, function(x) x/user_features$total), .SDcols = 2:20]
user_features$total <- NULL

# Add user features to top50
top_50_movies      <- top_50_movies[user_features, on=c(UserID = "UserID")]
top_50_movies      <- na.omit(top_50_movies)

rm(user_features, not_50, top_50, count_movies)

top_50_movies[, choice := as.numeric(as.factor(MovieID))]
top_50_movies[, reward := ifelse(Rating <= 4, 0, 1)]

# Run simulation ---------------------------------------------------------------------------------------------

simulations <- 1
horizon     <- nrow(top_50_movies)

formula     <- formula("reward ~ choice | i.V6 + i.V7 + i.V8 +i.V9 + i.V10 + i.V11 + i.V12 + i.V13 + i.V14 +
                                          i.V15 + i.V16 + i.V17 + i.V18 + i.V19 + i.V20 + i.V21 + i.V22 +
                                          i.V23 + i.V24")

bandit      <- OfflineBootstrappedReplayBandit$new(formula = formula, data = top_50_movies)

agents      <-
  list(Agent$new(ThompsonSamplingPolicy$new(), bandit, "Thompson"),
       Agent$new(RandomPolicy$new(), bandit, "Random"),
       Agent$new(LinUCBDisjointOptimizedPolicy$new(2.05), bandit, "LinUCB Dis"))

simulation  <-
  Simulator$new(
    agents           = agents,
    simulations      = simulations,
    horizon          = horizon
  )

sim  <- simulation$run()

plot(sim,
     type = "cumulative",
     regret = FALSE,
     rate = TRUE,
     legend_position = "bottomright")

library(contextual)

library(data.table)

library(splitstackshape)

# Movielens 100k ---------------------------------------------------------------------------------------------

# Info: https://d1ie9wlkzugsxr.cloudfront.net/data_movielens/ml-100k/ml-100k-README.txt

movies_dat <- "http://d1ie9wlkzugsxr.cloudfront.net/data_movielens/ml-100k/u.item"

ratings_dat <- "http://d1ie9wlkzugsxr.cloudfront.net/data_movielens/ml-100k/u.data"

# Import and merge files

movies_dat <- fread(movies_dat, sep = "|", quote="")

setnames(movies_dat, c("V1", "V2"), c("MovieID", "Name"))

movies_dat[, (3:5) := NULL ]

ratings_dat <- fread(ratings_dat, quote="")

setnames(ratings_dat, c("V1", "V2", "V3", "V4"), c("UserID", "MovieID", "Rating", "Timestamp"))

all_movies <- ratings_dat[movies_dat, on=c(MovieID = "MovieID")]

rm(movies_dat,ratings_dat)

# Data wrangling ---------------------------------------------------------------------------------------------

count_movies <- all_movies[,.(MovieCount = .N), by = MovieID]

top_50 <- as.vector(count_movies[order(-MovieCount)][1:50]$MovieID)

not_50 <- as.vector(count_movies[order(-MovieCount)][51:nrow(count_movies)]$MovieID)

top_50_movies <- all_movies[MovieID %in% top_50]

# User features: tags they've watched for non-top-50 movies normalized per user

user_features <- all_movies[MovieID %in% not_50]

rm(all_movies)

user_features[, c("MovieID", "Rating", "Timestamp", "Name"):=NULL]

user_features <- user_features[, lapply(.SD, sum, na.rm=TRUE), by=UserID ]

user_features[, total := rowSums(.SD, na.rm = TRUE), .SDcols = 2:20]

user_features[, 2:20 := lapply(.SD, function(x) x/user_features$total), .SDcols = 2:20]

user_features$total <- NULL

# Add user features to top50

top_50_movies <- top_50_movies[user_features, on=c(UserID = "UserID")]

top_50_movies <- na.omit(top_50_movies)

rm(user_features, not_50, top_50, count_movies)

top_50_movies[, choice := as.numeric(as.factor(MovieID))]

top_50_movies[, reward := ifelse(Rating <= 4, 0, 1)]

# Run simulation ---------------------------------------------------------------------------------------------

simulations <- 1

horizon <- nrow(top_50_movies)

formula <- formula("reward ~ choice | i.V6 + i.V7 + i.V8 +i.V9 + i.V10 + i.V11 + i.V12 + i.V13 + i.V14 +

i.V15 + i.V16 + i.V17 + i.V18 + i.V19 + i.V20 + i.V21 + i.V22 +

i.V23 + i.V24")

bandit <- OfflineBootstrappedReplayBandit$new(formula = formula, data = top_50_movies)

agents <-

list(Agent$new(ThompsonSamplingPolicy$new(), bandit, "Thompson"),

Agent$new(RandomPolicy$new(), bandit, "Random"),

Agent$new(LinUCBDisjointOptimizedPolicy$new(2.05), bandit, "LinUCB Dis"))

simulation <-

Simulator$new(

agents = agents,

simulations = simulations,

horizon = horizon

)

sim <- simulation$run()

plot(sim,

type = "cumulative",

regret = FALSE,

rate = TRUE,

legend_position = "bottomright")

おわりに

2~3年前に、Tokyo Web Miningの懇親会でContextual Banditの論文いいぞとテラモナギさんが紹介していて、へー、そんなのあるんだと、「へー」の域を出なかったんですが、一歩前進した気がします。先人が切り開いた道を2~3年後に舗装されてから通るというのも遅いなと感じられるので、残業もっと減らして勉強時間増やしたいと思います。

参考情報

バンディット問題の理論とアルゴリズム (機械学習プロフェッショナルシリーズ)
Bandit Algorithms for Website Optimization: Developing, Deploying, and Debugging
Contextual package ~ Japan.R Shota Yasui
Package ‘contextual’
バンディットアルゴリズムの復習３：UCB(Upper Confidence Bound)

Rでオペレーションズ・リサーチ(OR)に関する情報をあさる/ コード付き

はじめに

私は基本的にデータ分析を生業としていますが、どうしても分析の案件が足りない時期は分析以外のものに手を染めることもあります。主に、RPAやクローリング、APIを用いたソーシャルリスニングなどです。今後も分析以外のことをやる時があるとしたら、レパートリーを増やしたいですよね。なので、ORについて調べてみることにしました。

ORとは

公益社団法人日本オペレーションズ・リサーチ学会による定義によると、

「現象を抽象化した数理モデルを構築し, モデル分析に基づいて種々の問題, とりわけ意思決定問題の解決を支援する方法論や技法の総称. 情報化社会の進展に伴って, 線形計画法に代表される最適化モデルや待ち行列理論に代表される確率的なモデル等, 多様なモデルに基づく分析が, 経営計画や生産・販売・財務等の企業意思決定や都市・公共システム等広く社会一般の問題解決に大きな役割を果たしている.」

とされています。うむ、お硬い感じの定義ですね。
他の記述にわかりやすい表現がありました。

「問題を科学的，つまり「筋のとおった方法」を用いて解決するための「問題解決学」であります」
これならわかりやすいです。

問題解決につながる、あらゆる科学的な手法を扱っているのがORだと考えてよいのだろうと思います。

ORの強みとしては、

大規模プロジェクトなどの遂行に役立つ
常識や過去の経験では判断が難しい問題に対する解の提供をしてくれる
経営、工学、医学、公共政策など幅広い分野での適用可能性がある

などがあげられています。

ORの手法

問題解決につながるなら何でもありということで、手法も幅広いようです。私の持っている参考書やOR学会のサイトの情報から判断すると少なくとも以下の手法が扱われているようです。

・数理最適化
・組合せ最適化
・シミュレーション
・待ち行列
・AHP(階層的意思決定法)
・DEA(包絡分析法)
・スケジューリング
・ゲーム理論
・ネットワーク理論
・データマイニング

データマイニングはもはや機械学習ブームなので特筆したものではないですが、最近データ分析を始めた人などは知らないことも多いのではないでしょうか。

仕事やプライベートでの使い所

社員のシフトを決めるタスク
ミーティングに参加する社員の移動距離が一番少ない、空いている会議室を見つけるタスク
コンビニのレジを何台おけば客の待ち時間が想定内になるのか
工程A、B、Cの全てを経る必要のある作業で、生産ラインを安定させるには今日はAとBとCのどの工程をどれだけすすめるべきか決めるタスク
旅行をする際に、予算を所与のもとで、どのスポットに立ち寄ることが効用が高いかを見つけるタスク

Rでの実践例

ようやく、この記事の本題です。Rでオペレーションズリサーチなどという書籍に出会えていないので、Rを用いてオペレーションズリサーチを行っている事例を集めてみようと思います。ブログを漁ればいろいろとありますね。

なお、答え合わせを兼ねて、登場する問題は『Ｅｘｃｅｌで学ぶＯＲ』の例題で表現を変えています。網羅性はないものの、出来るだけ取り上げてみようと思います。

数理最適化

シュークリーム専門店、pseudoカモノハシはシュークリームとパンケーキの生産をしているが、厨房があまりに狭すぎてシュークリームとパンケーキの同時生産ができない。また、労働基準法の観点から厨房の利用時間は40時間以内となる。

シュークリームの生産に関しての情報は以下の通りとします。

pseudoカモノハシはシュークリームとパンケーキをどれだけ生産すれば利益を最大にすることができるだろうかという問題です。

定式化すると以下のようになります。

RのlpSolveパッケージを用いてこの線形計画問題を解いてみます。

library(lpSolve)

# 目的関数の係数
f.obj <- c(18, 11) 

# 制約式の左辺の係数
f.con <- matrix (c(1, 1, 2, 0, 0, 1), ncol=2, byrow=TRUE)

# 制約式の等号・不等号
f.dir <- c("<=", "<=", "<=")

# 制約式の右辺
f.rhs <- c(40, 50, 25)

# 決定変数は非負と仮定されているので，非負条件の記述は不要
# 目的関数の最大値を返す
lp ("max", f.obj, f.con, f.dir, f.rhs)

#解を求める
lp("max", f.obj, f.con, f.dir, f.rhs)$solution

library(lpSolve)

# 目的関数の係数

f.obj <- c(18, 11)

# 制約式の左辺の係数

f.con <- matrix (c(1, 1, 2, 0, 0, 1), ncol=2, byrow=TRUE)

# 制約式の等号・不等号

f.dir <- c("<=", "<=", "<=")

# 制約式の右辺

f.rhs <- c(40, 50, 25)

# 決定変数は非負と仮定されているので，非負条件の記述は不要

# 目的関数の最大値を返す

lp ("max", f.obj, f.con, f.dir, f.rhs)

#解を求める

lp("max", f.obj, f.con, f.dir, f.rhs)$solution

実行するとこんな結果です。

> # 目的関数の最大値を返す
> lp ("max", f.obj, f.con, f.dir, f.rhs)
Success: the objective function is 615 
> # 解を求める
> lp("max", f.obj, f.con, f.dir, f.rhs)$solution
[1] 25 15

> # 目的関数の最大値を返す

> lp ("max", f.obj, f.con, f.dir, f.rhs)

Success: the objective function is 615

> # 解を求める

> lp("max", f.obj, f.con, f.dir, f.rhs)$solution

[1] 25 15

pseudoカモノハシは25個のシュークリーム、15個のパンケーキを生産することで615万円の売上をあげることができるということになります。

係数がどれくらい変わっても最適解が変化しないのかを知るための感度分析も簡単にできるようです。
係数がちょっと変わっただけで崩れる最適化とかだと実務で使う際に怖いので、大事な工程ですね。

# 感度分析
# 最適解が変化しない目的関数の係数の下限値
lp ("max", f.obj, f.con, f.dir, f.rhs, compute.sens=T)$sens.coef.from 
# 最適解が変化しない目的関数の係数の上限値
lp ("max", f.obj, f.con, f.dir, f.rhs, compute.sens=T)$sens.coef.to
# 最適解が変化しない目的関数の係数の上限値
lp ("max", f.obj, f.con, f.dir, f.rhs, compute.sens=T)

# 感度分析

# 最適解が変化しない目的関数の係数の下限値

lp ("max", f.obj, f.con, f.dir, f.rhs, compute.sens=T)$sens.coef.from

# 最適解が変化しない目的関数の係数の上限値

lp ("max", f.obj, f.con, f.dir, f.rhs, compute.sens=T)$sens.coef.to

# 最適解が変化しない目的関数の係数の上限値

lp ("max", f.obj, f.con, f.dir, f.rhs, compute.sens=T)

輸送問題（ネットワーク型の線形計画法）

シュークリーム専門店、pseudoカモノハシは事業拡大につき、3つの食料庫と4つの工房を持つに至った。シュークリームやパンケーキを生産するためには食料庫から工房までトラックで輸送をする必要がある。各々の食料庫から工房までの輸送コストは以下の表の1行1列目〜3行4列目までで表される。食料庫には置ける在庫が決まっており、表の5列目で与えられている。工房には客の注文ベースの生産ノルマが課されており、表の4行目で与えられている。

以上、pseudoカモノハシは輸送コストを最も下げて生産するにはどの食料庫からどの工房に材料を輸送すればよいか、と言う問題となります。

定式化すると以下のようになります。

library(lpSolve)

# 目的関数の係数
f.obj <- c(40, 48, 21, 15, 52, 35, 45, 60, 25, 43, 70, 85) 

# 制約式の左辺の係数
f.con <- matrix (c(1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0 ,0,
                   0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0,
                   0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1,
                   1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0,
                   0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0,
                   0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0,
                   0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1), ncol=12, byrow=TRUE)

# 制約式の等号・不等号
f.dir <- c("<=",
           "<=",
           "<=",
           "==",
           "==",
           "==",
           "==")

# 制約式の右辺
f.rhs <- c(25,
           35,
           40,
           15,
           20,
           35,
           30)

# 決定変数は非負と仮定されているので，非負条件の記述は不要
# 目的関数の最大値を返す
lp ("min", f.obj, f.con, f.dir, f.rhs)

# 解を求める
lp("min", f.obj, f.con, f.dir, f.rhs)$solution

library(lpSolve)

# 目的関数の係数

f.obj <- c(40, 48, 21, 15, 52, 35, 45, 60, 25, 43, 70, 85)

# 制約式の左辺の係数

f.con <- matrix (c(1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0 ,0,

0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0,

0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1,

1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0,

0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0,

0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0,

0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1), ncol=12, byrow=TRUE)

# 制約式の等号・不等号

f.dir <- c("<=",

"<=",

"==",

"==")

# 制約式の右辺

f.rhs <- c(25,

35,

40,

15,

20,

35,

30)

# 決定変数は非負と仮定されているので，非負条件の記述は不要

# 目的関数の最大値を返す

lp ("min", f.obj, f.con, f.dir, f.rhs)

# 解を求める

lp("min", f.obj, f.con, f.dir, f.rhs)$solution

このコードを実行すると以下のようになり、輸送コストの最小値は3610であることが示されている。

> # 目的関数の最大値を返す
> lp ("min", f.obj, f.con, f.dir, f.rhs)
Success: the objective function is 3610 
> # 解を求める
> lp("min", f.obj, f.con, f.dir, f.rhs)$solution
 [1]  0  0  0 25  0  0 30  5 15 20  5  0

> # 目的関数の最大値を返す

> lp ("min", f.obj, f.con, f.dir, f.rhs)

Success: the objective function is 3610

> # 解を求める

> lp("min", f.obj, f.con, f.dir, f.rhs)$solution

[1] 0 0 0 25 0 0 30 5 15 20 5 0

ナップサック問題

シュークリーム専門店、pseudoカモノハシは創業3周年記念としてご当地グルメとのコラボレーションを計画している。グルメ系コンサルが提案したプランA〜Jの10件のご当地グルメとのコラボには、それぞれ費用と便益がある。この施策への予算が2000万円だとして、pseudoカモノハシは総便益が最大になるようにどのコラボ企画を採用するべきか。

定式化すると以下のようになります。

こちらにRのコードがあったので拝借しております。

library(lpSolve)

#便益
benefit <- c(4, 5, 3, 6, 13, 23, 11, 7, 15, 9)
#費用
cost <- matrix(c(3, 4, 2, 5, 10, 15, 6, 4, 13, 7)    , nrow=1)
#予算
budget <- 20

ans <- lp(direction = "max",
          objective.in = benefit,
          const.mat = cost,
          const.dir = "<=",
          const.rhs = budget,
          all.bin = TRUE)

library(lpSolve)

#便益

benefit <- c(4, 5, 3, 6, 13, 23, 11, 7, 15, 9)

#費用

cost <- matrix(c(3, 4, 2, 5, 10, 15, 6, 4, 13, 7) , nrow=1)

#予算

budget <- 20

ans <- lp(direction = "max",

objective.in = benefit,

const.mat = cost,

const.dir = "<=",

const.rhs = budget,

all.bin = TRUE)

これを実行すると、以下の結果が得られます。最大便益が31であること、それを実現するプランとしてE、G、Hが選ばれることがそれぞれ示しています。ただし、31を実現する解は他にもあります。

> print(ans)
Success: the objective function is 31 
> print(ans$solution)
 [1] 0 0 0 0 1 0 1 1 0 0

> print(ans)

Success: the objective function is 31

> print(ans$solution)

[1] 0 0 0 0 1 0 1 1 0 0

混合整数計画

シュークリーム専門店、pseudoカモノハシは傘下のパンケーキチェーン店、ヒッグス・シングスの3店舗から、「繁忙につき代替生産をお願いしたい」と言われた。ヒッグス・シングスのパンケーキは3店舗それぞれ味が異なり、それによってパンケーキの製造コストなども異なる。pseudoカモノハシの工房でパンケーキは生産可能ではあるが、そこそこに忙しいので生産できて15個が限度だと考えられる。そこで、以下の表が与えられたもとで、pseudoカモノハシはヒッグス・シングスの3店舗それぞれのパンケーキをどれだけ生産すれば利益が最大化されるか。

定式化すると以下のようになります。

今回は調査の結果、Rglpkというパッケージがあることがわかったので、そちらを用います。

library(Rglpk)

obj <- c(6.5, 18, 4, -10, -50, -6)
mat <- matrix(c(1, 1, 0, 0,
                1, 0, 1, 0,
                1, 0, 0, 1,
                0, -8, 0, 0,
                0, 0, -6, 0,
                0, 0, 0, -20), nrow = 4)
dir <- c("<=",
         "<=",
         "<=",
         "<=")

rhs <- c(15,
         0,
         0,
         0)

types <- c("I", "I", "I", "B","B","B")

max <- TRUE

Rglpk_solve_LP(obj, mat, dir, rhs, types = types, max = max)

library(Rglpk)

obj <- c(6.5, 18, 4, -10, -50, -6)

mat <- matrix(c(1, 1, 0, 0,

1, 0, 1, 0,

1, 0, 0, 1,

0, -8, 0, 0,

0, 0, -6, 0,

0, 0, 0, -20), nrow = 4)

dir <- c("<=",

"<=",

"<=")

rhs <- c(15,

types <- c("I", "I", "I", "B","B","B")

max <- TRUE

Rglpk_solve_LP(obj, mat, dir, rhs, types = types, max = max)

このコードを実行すると、以下の結果が得られます。工房1と工房2で8個と6個の生産にコミットすることで、10万円の最大利益を実現できることが示されています。

> Rglpk_solve_LP(obj, mat, dir, rhs, types = types, max = max)
$optimum
[1] 100

$solution
[1] 8 6 0 1 1 0

$status
[1] 0

$solution_dual
[1] NA

$auxiliary
$auxiliary$primal
[1] 14  0  0  0

$auxiliary$dual
[1] NA

> Rglpk_solve_LP(obj, mat, dir, rhs, types = types, max = max)

$optimum

[1] 100

$solution

[1] 8 6 0 1 1 0

$status

[1] 0

$solution_dual

[1] NA

$auxiliary

$auxiliary$primal

[1] 14 0 0 0

$auxiliary$dual

[1] NA

ウェーバー問題（非線形計画）

ウェーバー問題というのは、ORWikiによると以下の定義とされています。

施設・顧客間の距離に需要量を乗じたものの総和を最小化するような単一の施設の配置を、平面上の任意の地点の中から決定する問題。

具体例で取り組んでみましょう。

シュークリーム専門店、pseudoカモノハシは新しく食料庫を設置したいと考えている。店舗の位置が(x,y)座標で与えられているものとし、各店舗からの距離が最も小さくなるような位置に食料庫を設けるとしたらどこになるか、という問題。

定式化すると、以下のようになります。

orlocaというパッケージを使えばウェーバー問題のような非線形計画問題を簡単に解くことができます。

library(orloca)

# A new unweighted loca.p object
loca <- loca.p(x = c(2, 0, 1, 6, 5, 5.5, 4, 5, 4.5, 5.5),
               y = c(5, 0, 3, 0, 5, 4, 4.5, 3.5, 4, 1.5))

# Compute the minimum
sol <- distsummin(loca)

# Show the result
sol

# Evaluation of the objective function at solution point
distsum(loca, sol[1], sol[2])

library(orloca)

# A new unweighted loca.p object

loca <- loca.p(x = c(2, 0, 1, 6, 5, 5.5, 4, 5, 4.5, 5.5),

y = c(5, 0, 3, 0, 5, 4, 4.5, 3.5, 4, 1.5))

# Compute the minimum

sol <- distsummin(loca)

# Show the result

sol

# Evaluation of the objective function at solution point

distsum(loca, sol[1], sol[2])

このコードを実行すると以下の結果が得られます。x座標が4.468528、y座標が3.843755の時に、距離の最小値が22.80259となることが示されています。

> # Show the result
> sol
[1] 4.468528 3.843755
> # Evaluation of the objective function at solution point
> distsum(loca, sol[1], sol[2])
[1] 22.80259

> # Show the result

> sol

[1] 4.468528 3.843755

> # Evaluation of the objective function at solution point

> distsum(loca, sol[1], sol[2])

[1] 22.80259

せっかくなので、いらすとやの画像を背景に結果などをプロットしてみます。赤い点が解となった点です。

#図にプロット
library(ggplot2)
library(png)

image <- readPNG("map_open.png")
dataset <- data.frame(x=loca@x, y=loca@y)

g <- ggplot(data = dataset, aes(x = x, y = y))
g <- g + annotation_raster(image, xmin = -2, xmax = 8, ymin = -2, ymax = 8)
g <- g + geom_point()
g <- g + geom_point(data = data.frame(x=sol[1],y=sol[2]),
                    aes(x = x,y = y),
                    colour = "red",
                    size = 3)
g

#図にプロット

library(ggplot2)

library(png)

image <- readPNG("map_open.png")

dataset <- data.frame(x=loca@x, y=loca@y)

g <- ggplot(data = dataset, aes(x = x, y = y))

g <- g + annotation_raster(image, xmin = -2, xmax = 8, ymin = -2, ymax = 8)

g <- g + geom_point()

g <- g + geom_point(data = data.frame(x=sol[1],y=sol[2]),

aes(x = x,y = y),

colour = "red",

size = 3)

待ち行列

シュークリーム専門店、pseudoカモノハシ本店における来店客の行列に関してシミュレーションするものとする。

pseudoカモノハシ本店は平均して1時間に50人来店し、それはポワソン分布に従うとされている。
pseudoカモノハシ本店の商品は商品ごとに提供するまでにかかる時間が異なり、客はどれか1品を選択するが、その選択確率がおおむね決まっている。（以下の表）

進め方としては、まず、一様乱数を生成させ、ポアソン分布の分布関数の逆関数を求め、その逆関数に乱数を入力することで顧客の来店時間の間隔を生成します。

次に、顧客の選択する商品をシミュレーションするために別で一様乱数を生成し、その乱数の取る値に応じて商品を割り当てます。
そして、商品の提供に時間がかかることから、開始時間に商品の提供時間を足して、終了時間を求めます。ただし、来店時間に前の顧客がまだ商品を受け取れていないと、待ち時間が発生するので、来店時間の時点で終了していない場合はその分だけ開始時間が遅れます。

以下のRコードを作成してみました。

library(tidyverse)

#乱数の生成
set.seed(101)
rand_uni <-  runif(500)
itm <- -1/(50/60)*log(1-rand_uni)
ggplot(data = data.frame(x=itm), aes(x = x)) + geom_histogram(bins = 10)

set.seed(5)
rand_uni_demand <-  runif(500)
simulation <- data.frame(interval = itm,
                         x=rand_uni_demand) %>%
                         mutate(menu=if_else(x <= 0.25, "シュークリーム",
                                      if_else(x <= 0.5,"パンケーキ",
                                        if_else(x <= 0.8,"パフェ","ずんだ餅"))),
                                time_required=if_else(x <= 0.25, 1.2,
                                             if_else(x <= 0.5,1.5,
                                                     if_else(x <= 0.8,2.1,0.7))))

simulation <- simulation %>% mutate(arrival = cumsum(interval),
                                    start = NA)

for (i in 1:nrow(simulation)) {
  if (i == 1) {
    simulation$start[i] <- simulation$arrival[i]
  }
  else {
    end_time <- simulation$time_required[i-1] + simulation$start[i-1]
    simulation$start[i] <- if_else( end_time > simulation$arrival[i], end_time, simulation$arrival[i])
  }
}

#終了時間
simulation <- simulation %>% mutate(end = time_required + start)

#待ち時間
simulation <- simulation %>% mutate(latency = start - arrival)

#開始待ち人数
simulation <- simulation %>% mutate(numberof_wait= NA)

for (i in 1:nrow(simulation)) {
  simulation$numberof_wait[i] <- sum(simulation$arrival < simulation$start[i])
}

#空き時間
simulation <- simulation %>% mutate(idle_time = start - lag(end))

library(tidyverse)

#乱数の生成

set.seed(101)

rand_uni <- runif(500)

itm <- -1/(50/60)*log(1-rand_uni)

ggplot(data = data.frame(x=itm), aes(x = x)) + geom_histogram(bins = 10)

set.seed(5)

rand_uni_demand <- runif(500)

simulation <- data.frame(interval = itm,

x=rand_uni_demand) %>%

mutate(menu=if_else(x <= 0.25, "シュークリーム",

if_else(x <= 0.5,"パンケーキ",

if_else(x <= 0.8,"パフェ","ずんだ餅"))),

time_required=if_else(x <= 0.25, 1.2,

if_else(x <= 0.5,1.5,

if_else(x <= 0.8,2.1,0.7))))

simulation <- simulation %>% mutate(arrival = cumsum(interval),

start = NA)

for (i in 1:nrow(simulation)) {

if (i == 1) {

simulation$start[i] <- simulation$arrival[i]

}

else {

end_time <- simulation$time_required[i-1] + simulation$start[i-1]

simulation$start[i] <- if_else( end_time > simulation$arrival[i], end_time, simulation$arrival[i])

}

#終了時間

simulation <- simulation %>% mutate(end = time_required + start)

#待ち時間

simulation <- simulation %>% mutate(latency = start - arrival)

#開始待ち人数

simulation <- simulation %>% mutate(numberof_wait= NA)

for (i in 1:nrow(simulation)) {

simulation$numberof_wait[i] <- sum(simulation$arrival < simulation$start[i])

}

#空き時間

simulation <- simulation %>% mutate(idle_time = start - lag(end))

まず、来店間隔のシミュレーションですが、以下のようになります。

続いて、来店時間に応じた、待ち時間の推移です。

来店者の到来とともに、ぐんぐんと伸びているのがわかります。私はシュークリームに50分も待てないですね。

続いて、任意の顧客が開始した時に、すでに待っている客の数の推移です。

最後に、空き時間です。

顧客が忍耐強く待つのであれば、ほとんどレジは休めていないという劣悪な労働環境になりそうですね。厨房の能力やレジの能力などを高める必要がありそうです。

最短路問題

シュークリーム専門店、pseudoカモノハシのオーナーSKUE氏は、いま店舗1にいるが、店舗6の店長との1on1があるため、向かおうと考えている。ついでに他の店舗の店長にも顔を出したいと考えており、他の店を経由して一番距離が短い経路を見つけたい。店舗と店舗の距離は以下のグラフのノード間のラベルで与えられているものとする。

Rのigraphパッケージに「Shortest (directed or undirected) paths between vertices」というグラフ間の最短経路を見つける関数があるので、そちらを使います。この関数は最短経路問題を解く際に効率的とされる、ダイクストラ（Dijkstra）法というアルゴリズムを採用しています。具体的に使った例が載っているブログとしてこちらを参考にしています。

このグラフから、以下のようなデータを作っておきます。

library(igraph)

shortpath_dataset <- read_csv(file = "shortpath_problem.csv")

net <- graph.data.frame(shortpath_dataset,directed=F)
E(net)$weight <- shortpath_dataset$distance
E(net)$label <- shortpath_dataset$distance
E(net)$width <- 13
V(net)$size <- 35
plot(net)

spv <- as.data.frame(shortest.paths(net))

#pに出発点、経由点、終着点を付値
p <- c("shop1","shop6")
#aにダミーのデータを入れてリスト形式のオブジェクトを作成
a <- list("test")
#ベクトル形式のオブジェクトbを作成
b <- 0

for(i in 1:(length(p)-1)){　　
  #spに最短距離のノードリストを付値
  sp <- get.shortest.paths(net, from=p[i], to=p[i+1])
  #最短距離のノードリストをaに格納
  a[[i]] <- V(net)$name[unlist(sp)]
  #距離行列の行と列を指定して距離を取り出すしてbに格納
  b[i] <- spv[p[i], p[i+1]]
}
print(a)
print(b)

library(igraph)

shortpath_dataset <- read_csv(file = "shortpath_problem.csv")

net <- graph.data.frame(shortpath_dataset,directed=F)

E(net)$weight <- shortpath_dataset$distance

E(net)$label <- shortpath_dataset$distance

E(net)$width <- 13

V(net)$size <- 35

plot(net)

spv <- as.data.frame(shortest.paths(net))

#pに出発点、経由点、終着点を付値

p <- c("shop1","shop6")

#aにダミーのデータを入れてリスト形式のオブジェクトを作成

a <- list("test")

#ベクトル形式のオブジェクトbを作成

b <- 0

for(i in 1:(length(p)-1)){　　

#spに最短距離のノードリストを付値

sp <- get.shortest.paths(net, from=p[i], to=p[i+1])

#最短距離のノードリストをaに格納

a[[i]] <- V(net)$name[unlist(sp)]

#距離行列の行と列を指定して距離を取り出すしてbに格納

b[i] <- spv[p[i], p[i+1]]

}

print(a)

print(b)

このコードを実行すると、以下のように、店舗1→店舗3→店舗5→店舗6の順番で店舗を巡ると最短経路である13が達成されることが示されます。

> print(a)
[[1]]
[1] "shop1" "shop3" "shop5" "shop6"

> print(b)
[1] 13

> print(a)

[[1]]

[1] "shop1" "shop3" "shop5" "shop6"

> print(b)

[1] 13

今回は店舗が少ないので、人間の目でも見つけることができますね。

巡回セールスマン問題

カワウソ急便の配達員が今、バレンタインデーの集配のためにpseudoカモノハシの店舗1（本店）にいる。配達員は他の全ての店舗の集配もする必要があり、集配後は確認のために本店に立ち寄る必要がある。どのようにして店舗を一度ずつ回れば移動距離が最も小さくなるかに関心がある。店舗ごとの位置は以下の図の通り。

各地点の座標は以下の表で与えられている。

このような問題を巡回セールスマン問題と呼ぶが、解く際は座標の点から各々店舗の距離を計算し、距離行列を作成し、その距離行列をもとに、1度しか通れないという制約条件を課しながら距離が最も小さくなる組み合わせを見つける。
RではTSP(Traveling Salesperson Problem)パッケージがあるので、proxyパッケージを使って距離行列を作ってしまえば、簡単に最適な順路を示してくれる。

#図にプロット
library(ggplot2)
library(png)
library(ggrepel)

image <- readPNG("map_open.png")
dataset <- data.frame(x = c(2, 4, 5, 3, 2),
                      y = c(2, 1, 4, 6, 5),
                      label = c(1, 2, 3, 4, 5))

g <- ggplot(data = dataset, aes(x = x, y = y, label=label))
g <- g + annotation_raster(image, xmin = 0, xmax = 7, ymin = -1, ymax = 8)
g <- g + geom_point(size = 3)
g <- g + geom_point(data = data.frame(x=2,y=2,label=1),
                    aes(x = x,y = y),
                    colour = "red",
                    size = 3) + geom_text_repel()
g


library(TSP)
library(proxy)

#距離行列の作成
data <- as.matrix(dist(dataset[1:2], method="Euclidean"))
tsp <- TSP(data)
tsp

## use some methods
n_of_cities(tsp)
labels(tsp)

## calculate a tour
#start=1を指定することで1から始まるTSPを解いてくれる。methodでは様々な計算アルゴリズムが選択できる。
tour <- solve_TSP(tsp, method = "nn",start=1)
tour[1:5]
tour_length(tour)

#図にプロット

library(ggplot2)

library(png)

library(ggrepel)

image <- readPNG("map_open.png")

dataset <- data.frame(x = c(2, 4, 5, 3, 2),

y = c(2, 1, 4, 6, 5),

label = c(1, 2, 3, 4, 5))

g <- ggplot(data = dataset, aes(x = x, y = y, label=label))

g <- g + annotation_raster(image, xmin = 0, xmax = 7, ymin = -1, ymax = 8)

g <- g + geom_point(size = 3)

g <- g + geom_point(data = data.frame(x=2,y=2,label=1),

aes(x = x,y = y),

colour = "red",

size = 3) + geom_text_repel()

library(TSP)

library(proxy)

#距離行列の作成

data <- as.matrix(dist(dataset[1:2], method="Euclidean"))

tsp <- TSP(data)

tsp

## use some methods

n_of_cities(tsp)

labels(tsp)

## calculate a tour

#start=1を指定することで1から始まるTSPを解いてくれる。methodでは様々な計算アルゴリズムが選択できる。

tour <- solve_TSP(tsp, method = "nn",start=1)

tour[1:5]

tour_length(tour)

このコードを実行すると、以下のようになり、1→2→3→4→5→1と巡回することで総距離が12.641で済むことがわかる。

> tour[1:5]
1 2 3 4 5 
1 2 3 4 5 
> tour_length(tour)
[1] 12.64099

> tour[1:5]

1 2 3 4 5

> tour_length(tour)

[1] 12.64099

これ以降は疲れ果てたので、実践というより紹介にとどまりますが、あしからず。

スケジューリングと集合被覆問題

Rでの事例が見つかりませんでした。Qiitaで以下のようなPythonによる実践があったので、それをもとにRで書き換えてみるのも良いと思います。（後日、追記したいと思います。）
組合せ最適化 – 典型問題 – 集合被覆問題
 組合せ最適化 – 典型問題 – 勤務スケジューリング問題

NPV、IRR

ファイナンスで基本のNPV（Net Present Value）やIRR(Internal Rate of Return)の計算ですが、Package ‘FinCal’パッケージで簡単に計算できます。これくらい自分で書いてもいい気もしますが。このパッケージに割引率とキャッシュフローのベクトルを入力したらNPVを返してくれるnpv関数やirr関数があります。

> FinCal::npv(r=0.12, cf=c(-5, 1.6, 2.4, 2.8))
[1] 0.3348214
> FinCal::irr(cf=c(-5, 1.6, 2.4, 2.8))
[1] 0.1551911

> FinCal::npv(r=0.12, cf=c(-5, 1.6, 2.4, 2.8))

[1] 0.3348214

> FinCal::irr(cf=c(-5, 1.6, 2.4, 2.8))

[1] 0.1551911

データ包絡分析法

先日、ブレインパッドさんがこちらの記事で公開していたデータ包絡分析法(Data Envelopment Analysis:DEA)ですが、Rのコードが付いてなかったので漁ったところ、早稲田大学の逆瀬川先生がこちらでコードを公開しているようです。

DEAは同質な複数の事業体の相対的な効率性評価のための方法と定義されています。ブレインパッドさんの例ではコストを入力として、出力としての売上が効率的かどうかを見るために使われています。

パッケージに関しては、Data Envelopment Analysisでググったら、Benchmarkingというパッケージを見つけました。このパッケージは、少なくともブレインパッドさんのブログで紹介されている、DRS(Decreasing Returns to Scale)とFDH(Free Disposal Hull)は引数で選択可能のようです。

library(Benchmarking)

x <- matrix(c(100,200,300,500,100,200,600),ncol=1)
y <- matrix(c(75,100,300,400,25,50,400),ncol=1)

dea.plot.frontier(x,y,RTS="fdh+",txt=TRUE)

library(Benchmarking)

x <- matrix(c(100,200,300,500,100,200,600),ncol=1)

y <- matrix(c(75,100,300,400,25,50,400),ncol=1)

dea.plot.frontier(x,y,RTS="fdh+",txt=TRUE)

このコードを実行すると、ブログと似たような図が描けるようです。

詳しくはドキュメントについている先行研究とかを見たいところです。

ポートフォリオ選択問題

これもファイナンスで基本となる、ポートフォリオ選択問題なのですが、投資家のリスク選好から無リスク資産（国債とか）とリスク資産（株式とか）をどれくらいの配分で持つのが効率的かを解くという問題になります。

どうやらtidyquantというパッケージを使うことで、効率的ポートフォリオを探索できるようです。こちらのドキュメント（Tidy_Portfoliomanagement_in_R）をまだ実践できていないですが、最終的に以下のような効率的ポートフォリオを描けるようです。

編集を終えて

まだまだ原理を深く理解できていないですが、いざ仕事の依頼が来た際の取っ掛かりとしては良いものが手に入った気がします。加えて、人間だと計算が厳しそうな問題を解けるというのは非常に面白いです。
「問題設定→定式化→コードに書き落とす」という一連の訓練を続けるとかなり力が付きそうな気がします。Rって統計学・機械学習以外にも本当に幅広く取り揃っていて飽きがこなくてよいですね。

参考情報

Ｅｘｃｅｌで学ぶＯＲ
 サルでもわかる待ち行列
 RでLinear Programming
問題解決の数理（’１７）
Sensitivity Analysis 感度分析もし○○○だったどうする？
Rでデータ解析を始めよう020 Rでナップサック問題を解いてみよう
 CRAN Task View: Optimization and Mathematical Programming
Rで数理計画
 ＲでＯＲ：待ち行列モデル
 ［R］ggplot2によるグラフィックスで、図にPNG形式の画像を貼る
 Rでクラスター分析〜距離行列の生成からクラスタリングまで
 RでTSPの練習

2018年に参加したデータ分析系の勉強会で得た知識の詰め合わせ

社内に分析チームがないことから、私は月に3~4件は刺激を求めて勉強会に足を運んでいます。新しい知見を得れることは然ることながら、社内だともらえないフィードバックをいただけたり、課題の共有などをできるのが良いと思います。

目の肥えた皆さんにとって新規性のある情報はあまりないかもしれませんが、詰め合わせた情報をお楽しみください。

統計学まわり

勉強会名
KDD論文読み会
- 会社名
  LINE
- 知見
  Winner’s Curse Bias（勝者の呪いバイアス）という問題を回避するための効果検証のアプローチで、ABテストに勝とうが負けようがテストごとに負債を負っていくモデルとされています。100回ABテストやって99回負けたら99回分の負債を計算した上で、累積の価値貢献額を見積もろうという仕組みのようです。
- 発表資料
  論文読んだ「Winner’s Curse: Bias Estimation for Total Effects of Features in Online Controlled Experiments 」

論文読んだ「Winner’s Curse: Bias Estimation for Total Effects of Features in Online Controlled Experiments 」

勉強会名
KDD論文読み会
- 会社名
  LINE
- 知見
  因果推論に関するチュートリアルの文献を知れた
  Tutorial on Causal Inference and Counterfactual Reasoning
  DoWhyライブラリについても語られていたので、関連リンクも。
  統計的因果推論のためのPythonライブラリDoWhyについて解説：なにができて、なにに注意すべきか

機械学習まわり

勉強会名
merpay×M3 機械学習 NIGHT
- 会社名
  M3
- 知見
  コンテンツをレコメンドする際のテクニックとして、MFとCNNの合わせ技について紹介されていました。訓練時には、アクセスログデータをもとにMFで潜在的な表現を抽出しそれのアイテム間の類似度を計算し、推薦時には、テキストのタイトルとキーワードなどをCNNで学習し訓練時と同じ次元になるようにアイテムのベクトルを出力する。そして、訓練時のものと近いアイテムを推薦することでCold-Start問題を克服するとのことでした。
- 発表資料
  Matrix Factorization と Text CNN による Cold Start Problem への取り組み

Matrix Factorization と Text CNN による Cold Start Problem への取り組み from masahiro nishiba

勉強会名
MLCT
- 会社名
  LINE
- 知見
  GBDTの木の結果として得られるノード自体を特徴量として、予測を行う。特徴量圧縮に繋がり、予測の精度に関しても遜色ないらしい。
- 発表資料
  The Road to Machine Learning Engineer from Data Scientistの32ページ目

The Road to Machine Learning Engineer from Data Scientist

勉強会名
NetaDashi Meetup
- 会社名
  NRI
- 知見
  Elmoを用いた文書分類。Word2Vecなどではできなかった、文脈を考慮して類似度などを算出できる。
  Elmoの多言語対応に関しては、このGitHubを参照すると良いらしい。
  https://github.com/HIT-SCIR/ELMoForManyLangs

勉強会名
NetaDashi Meetup
- 会社名
  ？？？
- 知見
  異常検知の評価指標についてのお話でした。訓練データは正常なデータだけで、テストが異常かどうかを知りたいというユースケースで、Lee-Liu metricと呼ばれる評価指標があるそうです。ベイズの定理を使っているらしい。
  それを計算するためのkenchiというPythonモジュールもある。
  https://kenchi.readthedocs.io/en/latest/_modules/kenchi/metrics.html
- 発表資料
  https://speakerdeck.com/yohrn/metrics-for-one-class-classification

異常検知の評価指標って何を使えばいいの？ / Metrics for one-class classification

勉強会名
グリー開発本部 Meetup #1 DataEngConf NYC報告会
- 会社名
  GREE
- 知見
  Contextual Banditについての紹介
  「Artwork Personalization at Netflix」という記事で2017年ごろに取り上げられていたようです。

エンジニアリングまわり

勉強会名
グリー開発本部 Meetup #1 DataEngConf NYC報告会
- 会社名
  GREE
- 知見
  LUIJI
  ・メリットとしては、少なくとも以下のものがあるそうな。
  　・Pythonで書ける
  　・エラーの途中で処理を止めて、それを解消したら、止めたポイントから開始できる
  　・様々なツール群と連携できる柔軟性
  　・10行程度でスクリプト書ける。
  　・複雑な依存関係も描ける。
- 発表資料
  https://www.slideshare.net/greetech/dataengconf-nyc18-1

DataEngConf NYC’18 セッションサマリー #1 from gree_tech

勉強会名
bq_sushi tokyo #9 2018総集編
- 会社名
  オープンハウス
- 知見
  BigQueryGIS
  BigQueryからGISの情報を扱うことが可能になったらしい。顧客の希望する物件の情報をレコメンドするために地理情報を扱うらしいです。
  ただ、基準とする測地系が国によって異なり、それらを考慮しないで推薦すると1~2kmはズレてしまうとのこと。家買う際にそんだけズレるとキツイですね。こちら（BigQueryGIS: Google und PostGIS　）はBigQueryGISに関連した情報を漁って見つけた記事ですが、BQで抽出した情報をそのままGoogleMapに表示できるのは面白いですね。

データ分析のツラミ系

勉強会名
merpay×M3 機械学習 NIGHT
- 会社名
  M3
- 知見
  メタデータの検索システムについて
  データセット名、テーブル名、カラム名、カラムのディスクリプションをキーワードで検索できる。
  日次でディスクリプションを取ってくるようにしている。どのドキュメントが一番見られているのかもモニタリングできるとのこと。似たような取り組みとして、リクルートがMetaLookingとかいう内製ツールを作っていたりしますね。私は各サービスごとのDBのテーブルの注意点などを適宜スプレッドシートに残す程度しかしていませんが、分析者がすぐにキャッチアップできる環境は重要ですね。

勉強会名
MLCT
- 会社名
  ？？？
- 知見
  事業計画書を作るリーンキャンバスの機械学習版とも言える、機械学習キャンバス0.1というものが質疑応答の際に紹介されていました。

機械学習キャンバス0.1 from nishio

勉強会名
グリー開発本部 Meetup #1 DataEngConf NYC報告会
- 会社名
  GREE
- 知見
  データリーク問題はどこも苦しんでいる？
  SalesForce社が顧客企業15万社の情報を活用して、機械学習モデルを構築しようとしたが、
  蓄積されたデータにおいては、ビジネスプロセスをやたらと予測できてしまうようなデータリーク問題が起きまくっていた。
  原因としては、データサイエンティスト不足（分析を前提としたデータ蓄積ができていない。）、手入力によるラベリングミスなどがあるらしい。
  どこの企業も苦しんでいると思うと、分析を前提にスナップショットを残し続けるという取り組みは競争優位性につながるのだろうか。
  SalesForce社は、訓練と検証の精度の差が大きいと注意したり時系列データを確認するなどして、データの信憑性に気をつけてモデルを作ったそうです。
  15万社にうまくフィットするモデルなので、精度は70~75%で満足できるものらしい。
- 発表資料
  https://www.slideshare.net/greetech/dataengconf-nyc18-1

DataEngConf NYC’18 セッションサマリー #1 from gree_tech

R Advent Calendar 2018 一発屋芸人の検索トレンドの分析

はじめに

昨年のR Advent Calendarはポケモンのデータをrvestでスクレイピングして、レアポケモンがどのような種族値における特徴があるのかを探ったり、経験値が必要な割に種族値が低い「コスパの悪いポケモン」などを見つけました。
今年のR Advent Calendarでは、年末年始といえば一発屋芸人のテレビなどでの露出が多くなることから、一発屋芸人の検索トレンドのデータを手に入れて分析してみたいと思います。

分析工程

・データの収集
・データの整形
・可視化
・分析

データの収集

こちらのサイト（流行した一発屋芸人一覧／年代流行）に一発屋の芸人さんが列挙されていました。私は普段テレビを見ないので大体の芸人さんがわからないです。

Googleトレンドから、芸人名に関するGoogle検索の時系列データを収集します。

非常に残酷なデータだなと思いました。

ただ、一つ弁護すると、Googleトレンドはレベルではなくピークを1として標準化した数値をデータとして提供してくれていますので、
ピークが著しく高ければ、今の水準が低くてもそこそこ検索されている可能性はあるとだけ言っておきます。

本当の検索回数が必要な場合は、Google Adwords（検索連動型広告）のアカウントの開設とともに検索ボリューム取得APIなどの申請が必要なので、正確なデータが必要な場合は会社として取り組んだほうが良いと思います。個人では厳しいです。

データの整形

各芸人さん（総勢21名）の検索トレンドデータのピークの6ヶ月前までのデータと6ヶ月後のデータまでの合計1年間の検索トレンドを各々抽出してみようと思います。
GoogleトレンドのデータはCSVでダウンロードできますので、そのCSVを読み込み、トレンドのデータを文字列から数値にし、ピークの前後12ヶ月ずつのデータを抽出します。
そうすることで、一発屋芸人のピークの前後に関するデータを作ります。（ただし、今朝、Google Trendのデータを取得できるgtrendsRというパッケージがR bloggerで紹介されていました。APIないはずなんですが、URLの工夫か裏でSelenium動かしていたりするんですかね。）

library(tidyverse)
library(directlabels)
library(TSclust)
library(gghighlight)

# データの準備 ------------------------------------------------------------------
trend_dataset <- data.frame()

for (i in 1:21) {
  #データの読込
  trend_data <- read_csv(file = paste0("multiTimeline_",i,".csv"),skip = 2)
  trend_data <- trend_data %>% mutate( gsub(": (日本)","",colnames(trend_data)[2]) )
  colnames(trend_data) <- c("month", "trend", "keyword")
  #1未満のデータをゼロにする
  trend_data <- trend_data %>% mutate(trend = as.numeric(replace(trend, trend=="1 未満", 0)))
  #ピークの月の前後12ヶ月を抽出
  trend_data <- trend_data[(which.max(trend_data$trend)-12):(which.max(trend_data$trend)+12),]
  trend_dataset <- trend_dataset %>% rbind(trend_data)
}

#キーワードごとにインデックスをふる

trend_dataset <- trend_dataset %>%
                        group_by(keyword) %>%
                        mutate(period = 1:n())

library(tidyverse)

library(directlabels)

library(TSclust)

library(gghighlight)

# データの準備 ------------------------------------------------------------------

trend_dataset <- data.frame()

for (i in 1:21) {

#データの読込

trend_data <- read_csv(file = paste0("multiTimeline_",i,".csv"),skip = 2)

trend_data <- trend_data %>% mutate( gsub(": (日本)","",colnames(trend_data)[2]) )

colnames(trend_data) <- c("month", "trend", "keyword")

#1未満のデータをゼロにする

trend_data <- trend_data %>% mutate(trend = as.numeric(replace(trend, trend=="1 未満", 0)))

#ピークの月の前後12ヶ月を抽出

trend_data <- trend_data[(which.max(trend_data$trend)-12):(which.max(trend_data$trend)+12),]

trend_dataset <- trend_dataset %>% rbind(trend_data)

}

#キーワードごとにインデックスをふる

trend_dataset <- trend_dataset %>%

group_by(keyword) %>%

mutate(period = 1:n())

可視化

作成したデータを実際にプロットしてみます。

old = theme_set(theme_gray(base_family="HiraKakuProN-W3"))

ggplot(data = trend_dataset,
            aes(x = period, y = trend, color=keyword))  + geom_line()

old = theme_set(theme_gray(base_family="HiraKakuProN-W3"))

ggplot(data = trend_dataset,

aes(x = period, y = trend, color=keyword)) + geom_line()

一発屋にも盛り上がり方に違いがあるようですね。

時系列クラスタリングの適用

多様な盛り上がり方があることから、TSclustというライブラリを使って時系列クラスタリングを行い、トレンドに関しての分類的なものを得たいと思います。
今回初めて使うのですが、参考文献によると様々な類似性指標を指定して、時系列ごとの類似性を計算するようです。ピアソン相関係数のようなシンプルなものもあれば、ユークリッド距離のものやFrechet距離とかいう聞いたことないものまで幅広く用意されています。今回はシンプルにピアソン相関係数にしてみます。そして、類似性指標を出してから、そのまま階層クラスタリングを行います。

trend_dataset_spread <- trend_dataset %>% select(-month) %>% tidyr::spread(key = keyword, value = trend)

# COR距離で距離行列を作成
d <- diss(trend_dataset_spread %>% select(-period), "COR")

#デフォルトの設定で階層クラスタリング
h <- hclust(d)

#階層クラスタリングの結果の可視化
par(cex=0.6)
par(family = "HiraKakuProN-W3")
plot(h, hang = -1)

trend_dataset_spread <- trend_dataset %>% select(-month) %>% tidyr::spread(key = keyword, value = trend)

# COR距離で距離行列を作成

d <- diss(trend_dataset_spread %>% select(-period), "COR")

#デフォルトの設定で階層クラスタリング

h <- hclust(d)

#階層クラスタリングの結果の可視化

par(cex=0.6)

par(family = "HiraKakuProN-W3")

plot(h, hang = -1)

こちらが、TSclustのdiss関数を用いて計算した時系列データごとの距離を、階層クラスタリングにより描いたデンドログラムです。

この分類だけ見ても、芸人さんを知らない私からすると何も共感がありませんので、先程のクラスタリング結果をもとに可視化をしてみます。
そこで、Tokyo.Rで知らない人はいないであろう、yutaniさんの作られたgghighlightを使ってみようと思います。

ただ、日本語のラベルの表示がうまくいかなかったので、芸人さんの名前をGoogleSpreadSheetのGoogle翻訳関数（GOOGLETRANSLATE）で英訳しておきます。

（Anyway bright YasumuraやThick slice Jasonは結構キャッチーなのでは？）

# クラスタ数は3とする
data.frame(cutree(h, 3))
clusters <- data.frame(cluster_number=cutree(h, 3))
clusters$keyword <- rownames(clusters)
rownames(clusters ) <- NULL
trend_dataset_withcluster <- trend_dataset %>% left_join(clusters, by = "keyword")

#英訳したデータの読み込みと結合
rename_keywordlist <- read_csv("rename_keywordlist.csv")
trend_dataset_withcluster <- trend_dataset_withcluster %>% left_join(rename_keywordlist, by = "keyword")


gghighlight_line(trend_dataset_withcluster, aes(period, trend, colour = keyword_en),
                 predicate = max(cluster_number) == 1)

gghighlight_line(trend_dataset_withcluster, aes(period, trend, colour = keyword_en),
                 predicate = max(cluster_number) == 2)

gghighlight_line(trend_dataset_withcluster, aes(period, trend, colour = keyword_en),
                 predicate = max(cluster_number) == 3)

# クラスタ数は3とする

data.frame(cutree(h, 3))

clusters <- data.frame(cluster_number=cutree(h, 3))

clusters$keyword <- rownames(clusters)

rownames(clusters ) <- NULL

trend_dataset_withcluster <- trend_dataset %>% left_join(clusters, by = "keyword")

#英訳したデータの読み込みと結合

rename_keywordlist <- read_csv("rename_keywordlist.csv")

trend_dataset_withcluster <- trend_dataset_withcluster %>% left_join(rename_keywordlist, by = "keyword")

gghighlight_line(trend_dataset_withcluster, aes(period, trend, colour = keyword_en),

predicate = max(cluster_number) == 1)

gghighlight_line(trend_dataset_withcluster, aes(period, trend, colour = keyword_en),

predicate = max(cluster_number) == 2)

gghighlight_line(trend_dataset_withcluster, aes(period, trend, colour = keyword_en),

predicate = max(cluster_number) == 3)

まずはクラスター1

比較的短期でピークに達し、すぐに検索されなくなる、一発屋の名に相違ない傾向を持ったクラスターのように思われます。「日本エレキテル連合」とか「楽しんご」とか「8.6秒バズーカ」とかです。

続いてクラスター2

急激にピークに達するものの、ややしぶとく残り続けるような一発屋のクラスターなのかなと思います。「レイザーラモンHG」とか「厚切りジェイソン」とか「ピコ太郎」とか「世界のナベアツ」です。

そしてクラスター3

3人の芸人さんしか属していないですね。クラスターの数は2個でもよかったかもしれない。段階的にピークに達し、一気に落とされるという一発屋のクラスターのようです。「とにかく明るい安村」とか「藤崎マーケット」とか「すぎちゃん」とかです。

様々な傾向の一発屋さんがいるのがわかりました。

トレンドの推定

今回扱っているデータは芸人さんの数×時点のデータの多変量時系列となります。都合の良いものはないかと考えていましたが、古典的なVARではサンプルサイズ的にかなり苦しいと思い、Stanによるダイナミックパネルデータ分析などの事例はないか漁っていましたが、なかなかありませんでした。

松浦さんの『StanとRでベイズ統計モデリング (Wonderful R)』の241pに書かれている、モデル式12-8や12-9が今回のものに適しているなと思いましたが、コードを上げている方は見当たらなかったです。よしそれならば作ろうかと思った矢先、logics-of-blueさんのStan Advent Calendarの投稿、「Stanで推定する多変量時系列モデル」がかなりどんぴしゃな内容でしたので、コードを拝借してこの一発屋データの推定をしてみようと思います。

まずは、stanのコード

data {
  int T;                       // データ取得期間の長さ
  int performer_num;           // 芸人さんの数
  matrix[T, performer_num] y;  // 観測値
}

parameters {
  vector[T] x;        // 状態の推定値
  vector[performer_num] r; // 芸人さん毎のランダム効果
  real<lower=0> s_w;  // 過程誤差の標準偏差
  real<lower=0> s_v;  // 観測誤差の標準偏差
  real<lower=0> s_r;  // ランダム効果の標準偏差
  //vector[performer_num] sigma; //芸人さん毎の標準偏差
}

model {
  // 状態方程式に従い、状態が遷移する
  for(i in 2:T) {
    x[i] ~ normal(x[i-1], s_w);
  }
  
  // ランダム効果
  r ~ normal(0, s_r);
  
  // 観測方程式に従い、観測値が得られる
  for(i in 1:T) {
    for(j in 1:performer_num) {
      y[i, j] ~ normal(x[i] + r[j], s_v);
    }
  }
}

data {

int T; // データ取得期間の長さ

int performer_num; // 芸人さんの数

matrix[T, performer_num] y; // 観測値

}

parameters {

vector[T] x; // 状態の推定値

vector[performer_num] r; // 芸人さん毎のランダム効果

real<lower=0> s_w; // 過程誤差の標準偏差

real<lower=0> s_v; // 観測誤差の標準偏差

real<lower=0> s_r; // ランダム効果の標準偏差

//vector[performer_num] sigma; //芸人さん毎の標準偏差

}

model {

// 状態方程式に従い、状態が遷移する

for(i in 2:T) {

x[i] ~ normal(x[i-1], s_w);

}

// ランダム効果

r ~ normal(0, s_r);

// 観測方程式に従い、観測値が得られる

for(i in 1:T) {

for(j in 1:performer_num) {

y[i, j] ~ normal(x[i] + r[j], s_v);

}

そしてキックして結果を可視化するためのRコード

library(rstan)
library(bayesplot)

T <- trend_dataset_spread %>% select(-period) %>% nrow()
performer_num <- trend_dataset_spread %>% select(-period) %>% ncol()

data <- list(T = T,
             performer_num = performer_num ,
             y = trend_dataset_spread %>% select(-period))

fit <- stan(file = 'multivariate_time_series.stan',
            data = data,
            seed = 1, 
            iter = 30000,
            warmup = 10000,
            thin = 10
            )

mcmc_rhat(rhat(fit))

# データの整形
stan_df_1 <- fit %>% 
             rstan::extract() %$% x %>% 
             apply(2, quantile, probs = c(0.025, 0.5, 0.975)) %>% 
             t() %>% 
             cbind(1:nrow(trend_dataset_spread)) %>% 
             data.frame
# 列名の変更
colnames(stan_df_1) <- c("lwr", "fit", "upr", "time")
# 結果
head(stan_df_1, n = 3)

ggplot(data = trend_dataset_withcluster) + 
  ggtitle("推定結果（ピコ太郎）") +
  geom_line(aes(x = period, y = trend, color = keyword_en)) + 
  gghighlight(keyword_en == "Pico Taro", use_group_by = FALSE) + 
  geom_line(data = stan_df_1, 
            aes(x = time, y = fit), size = 1.2) +
  geom_ribbon(data = stan_df_1, 
              aes(x = time, ymin = lwr, ymax = upr), alpha = 0.3)

library(rstan)

library(bayesplot)

T <- trend_dataset_spread %>% select(-period) %>% nrow()

performer_num <- trend_dataset_spread %>% select(-period) %>% ncol()

data <- list(T = T,

performer_num = performer_num ,

y = trend_dataset_spread %>% select(-period))

fit <- stan(file = 'multivariate_time_series.stan',

data = data,

seed = 1,

iter = 30000,

warmup = 10000,

thin = 10

)

mcmc_rhat(rhat(fit))

# データの整形

stan_df_1 <- fit %>%

rstan::extract() %$% x %>%

apply(2, quantile, probs = c(0.025, 0.5, 0.975)) %>%

t() %>%

cbind(1:nrow(trend_dataset_spread)) %>%

data.frame

# 列名の変更

colnames(stan_df_1) <- c("lwr", "fit", "upr", "time")

# 結果

head(stan_df_1, n = 3)

ggplot(data = trend_dataset_withcluster) +

ggtitle("推定結果（ピコ太郎）") +

geom_line(aes(x = period, y = trend, color = keyword_en)) +

gghighlight(keyword_en == "Pico Taro", use_group_by = FALSE) +

geom_line(data = stan_df_1,

aes(x = time, y = fit), size = 1.2) +

geom_ribbon(data = stan_df_1,

aes(x = time, ymin = lwr, ymax = upr), alpha = 0.3)

そのまんま実行して、一発屋の時系列の中央値を可視化したらこんな感じになりました。一発屋のトレンドをうまく抽出できているのかなと思います。

今後の改良としては、階層性を持たせ、芸人さんごとのハイパーパラメータを持たせるとかなのですが、正月にでも取り組みたいと思います。（芸人さん以外のデータでやりたい。）

一方で、他にも多変量時系列で何かないか漁っていたのですが、Applied Time Series Analysis for Fisheries and Environmental Sciences : Dynamic factor analysisで紹介されている、Dynamic Factor Analysisというものが面白そうだなと思いました。
bayesdfaというパッケージを用いて、多変量時系列データに存在するであろうトレンドをStanを用いて推定することができるようです。元となった論文には各エリアごとのノルウェーロブスターの個体数のトレンドを推定し、3つのトレンドが発見されたとしています。ただ、同時点間のデータではないという点から今回のデータへの適用は不適切です。

同時点間に観測されていないデータであるという問題を認識した上で、このパッケージを使ってどんなトレンドを抽出できるのか試してみようと思います。

mod_3 = bayesdfa::fit_dfa(y = trend_dataset_spread %>% select(-period) %>% t(), num_trends = 3)

rot = bayesdfa::rotate_trends(mod_3)
names(rot)

matplot(t(rot$trends_mean), type = "l", lwd = 2, ylab = "mean trend")

mod_3 = bayesdfa::fit_dfa(y = trend_dataset_spread %>% select(-period) %>% t(), num_trends = 3)

rot = bayesdfa::rotate_trends(mod_3)

names(rot)

matplot(t(rot$trends_mean), type = "l", lwd = 2, ylab = "mean trend")

徐々に増えてから一気に落ちるトレンドや、一気に増えてから徐々に落ちるトレンドなどがうまく捉えれている気がします。
さらなる試行として、AICのような情報量基準である、Leave One Out Information Criterion (LOOIC)が最も低くなるトレンドの数を探索してみます。

mod_1 = fit_dfa(y = trend_dataset_spread %>% select(-period) %>% t(), num_trends = 1)
mod_2 = fit_dfa(y = trend_dataset_spread %>% select(-period) %>% t(), num_trends = 2)
mod_3 = fit_dfa(y = trend_dataset_spread %>% select(-period) %>% t(), num_trends = 3)
mod_4 = fit_dfa(y = trend_dataset_spread %>% select(-period) %>% t(), num_trends = 4)
mod_5 = fit_dfa(y = trend_dataset_spread %>% select(-period) %>% t(), num_trends = 5)

mod_1 = fit_dfa(y = trend_dataset_spread %>% select(-period) %>% t(), num_trends = 1)

mod_2 = fit_dfa(y = trend_dataset_spread %>% select(-period) %>% t(), num_trends = 2)

mod_3 = fit_dfa(y = trend_dataset_spread %>% select(-period) %>% t(), num_trends = 3)

mod_4 = fit_dfa(y = trend_dataset_spread %>% select(-period) %>% t(), num_trends = 4)

mod_5 = fit_dfa(y = trend_dataset_spread %>% select(-period) %>% t(), num_trends = 5)

トレンド数を1から5まで指定して実行した結果、5の時が一番LOOICが低くなりました。

rot = bayesdfa::rotate_trends(mod_5)
names(rot)
matplot(t(rot$trends_mean), type = "l", lwd = 2, ylab = "mean trend")

rot = bayesdfa::rotate_trends(mod_5)

names(rot)

matplot(t(rot$trends_mean), type = "l", lwd = 2, ylab = "mean trend")

まぁ、適切な使い方ではないのですが、徐々に増えてから一気に落ちるトレンドや、一気に増えてから徐々に落ちるトレンドなどが引き続き捉えれているようです。

今後の課題

・Stanによる多変量時系列のモデリングをしてみる。（Dynamic Panel分析とかもできると良い。少なくともStanのドキュメントにはない。）
・Dynamic Factor Analysisの適切な事例での適用をしてみる。

それでは、どうか良い年末をお過ごし下さい！
メリークリスマス！

参考情報

Introduction to gghighlight: Highlight ggplot’s Lines and Points with Predicates
{TSclust} ではじめる時系列クラスタリング
 Applied Time Series Analysis for Fisheries and Environmental Sciences 9.7 Dynamic factor analysis
読了：Montero & Vilar (2014) RのTSclustパッケージで時系列クラスタリング

参加できなかった第74回TokyoRのキャッチアップと結婚式について

自身の結婚式があったため、参加できなかった第74回目のTokyoRについてキャッチアップするために、公開資料を読んだ際に感じたメモをここに記します。
そして、世のデータサイエンティストが結婚式の際に苦しまないように参考になりそうな情報を少し書きました。

今回はTogetterも初めて作ってみました。（第74回R勉強会@東京（#TokyoR）のタグが付いたものたち）

初心者セッション

初心者セッション1 – Data Import & Export –

資料なし

初心者セッション2 – Data Handling –

https://ymattu.github.io/TokyoR74/slide.html#/

いつもながら良い資料です。まだ社内でアクティブにRを広めてはいないですが、広めるならこの資料が良いですね。
上から下に直感的に書けるのは初学者には大事なので。
ただ、SQLの知識がないメンバーとかには補助教材が必要ですね。

登場するパッケージはdplyrとlubridateとstringrとforcatsとpurrr。

初心者セッション3 – Plot & Visualization –

資料なし

応用セッション

How LINE Corp Use R to Compete in a Data-Driven World

資料なし

LINEでのRの活用最前線の話なのでしょうが、資料がないのが寂しいです。
きっと、タイムラインで共有されていたこの記事をベースに話されているのではないか。

LINE の全社員が必要に応じて担当サービスのデータを分析できる環境の構築

総勢50名の機械学習エンジニア・データサイエンティスト・データプランナー・データエンジニアが活躍できる分析基盤を作っていますよと。
Hadoop クラスタのデータの全社公開&活用のために、エンドユーザ向けのWebインターフェース「OASIS」をゼロから新規に開発。
Spark, Spark SQL, PySpark, SparkRおよびPrestoクエリを叩くことができる。
約20のサービス・部署で利用され、月間利用者数は約200人（データ関連の人以外で150人は触っていることになる。すごく層が厚そう。）

あと、ブリスベンの写真がTwitterで写り込んでいたので、UseR!2018のお話をされているのだろうと思われます。
linerパッケージについても語られている模様。

UseR!2018に参加し、社内Rパッケージ「liner」の活用事例を紹介しました

これのことでしょう。写真も合致している。
「いらすとや」の画像がワールドワイドに使われているのがシュールでいいですね。

「データの取得、分析、レポーティング、そして結果の共有にいたるまで、様々な便利機能を提供」と書かれています。
図を見る限りは、

RStudioでの解析結果をDBにカジュアルに保存したり、通知したり、ドキュメント化したりできる
PrestoやらHiveQLなどもRStudio上で実行できる（Tab補完とかもしてくれるんでしょうか？そこは聞いてみないとわからない。）
コーポレートカラーを適用したggplot2を利用できる
A/Bテストの結果をShinyのアプリで確認できる

などの機能があるようです。50人もいる高単価な人達の時間を節約できるという点でも、すごくインパクトのある取り組みですね。

tidyeval入門以前

Yet Another Introduction to tidyeval from yutannihilation

(speakerdeckのembedに若干苦戦しましたｗ)

湯谷さんの考えるtidyevalについて英語で書かれています。
環境に応じた値をRは自動で引っ張ってくるけど、たまに干渉してしまうことがあり、実際dplyrとstatsは共にfilter()関数を持っていると。
どの優先順位をもたせるかのコントロールって難しそうですね。
干渉を避けるためにquosureやunquoteというのが説明されています。私の理解が追いついていないので、正直めちゃわかった感はないのですが、
様々な環境下で動くパッケージを作ろうと考える際は不可避な領域なのではないでしょうか。

LT

「うまい飯が作りたい」

introduction_of_recipes from Yutaka Kuroki

recipesパッケージの紹介です。面白い方なんですね。
tidymodelsパッケージの中に内包されているそうな。

recipesパッケージの使い方を丁寧に説明してくださっています。

目的変数と説明変数を最初に明示的に指定し、前処理の手法をパイプ演算子でつないでいくだけ。
前処理のステップは使い回せるとのことで、似たようなデータを扱う場合は使いまわして楽をできるとのこと。
これは試してみる価値がありそうですね。k-nn法を用いた欠損値補完なども関数として用意されているようです。

マジレスすると、モテるかどうかは容姿や性格によるところが大きいと思うので、
身なりを清潔に保つとか、連れて行くと喜ばれそうな場所や体験を提供するとかが近道な気がしますよね。

不連続回帰とrdrobustパッケージの紹介

https://www.slideshare.net/YusukeKaneko6/tokyor74rdd-122646880

計量経済学系のバックグラウンドをお持ちのkagglerの方のLTです。
イスラム教の社会が、女性の高校修了率に与える影響について、統計的因果推論をされています。
rdrobustというパッケージを用いて、不連続回帰（RDD）という手法を実践され、イスラム教政治がランダム割当になりやすい状況を作っています。
全データで推定すると負の影響が推定されていましたが、RDDを使うことで正の効果があるという結果となりました。これは他の学部領域での先行研究とも整合的とのこと。

PCAや対応分析で補完要素を使う

資料なし

CiNii API その2

資料なし

Rで健康体

資料は後ほど公開される模様
https://twitter.com/weda_654/status/1061193131335475201

データのみは公開されている
https://github.com/weda-654/my_health_log

Soccer × Attribution Analysis

u++さんのLTですね！アクセスログや第三者配信データを用いた分析でよく扱われるアトリビューション分析をサッカーの貢献度に利用したという話です。
ChannelAttributionパッケージを使われています。

マーケティングの業務で使う際は、これまで見逃していた意外なページや意外な参照元などがこの分析で見えてきたりします。
ただ、価値はあってもコストがかかっては元も子もないので、コストに関する記述があったのも実務で使われている方の視点だなと思いました。

サッカーはルールくらいしか知らないので誰が意外なのかはわからないですｗ

地理空間データの交差検証、正しくできていますか？

地理空間データにおける交差検証する際の手法として、Spacial Cross-ValidationとTarget-oriented cross-validationが挙げられています。
空間データ向けのパッケージとしては、sfやCASTが、学習周りではmlrやcaretが扱われていました。
Referenceがあるのが嬉しいですね。

reticulateパッケージとデータサイエンスフロー

資料なし

結婚式について

なぜ開くのか

大事な思い出づくりのため。

工程

Willing to Payの決定（全てはここ）
- 持ち込みし放題の式場
- 大学OB割引などがある式場
- 料理が美味しい（国賓を迎えたりしているか？）
- 荘厳な雰囲気
- アクセスが良い
教会式か人前式か神前式かの選定（教会式だとめっちゃ高かったりするところもある）
誘う対象の選定
- 共に勉強を頑張った仲
- 共に仕事を頑張った仲
- 親族
住所の聞き出し
- Googleフォームを活用
大量の切手の購入（送付用と返信用も）
- 郵便局は21時まで空いていたりするので助かった。（どこもそうなのかな？）
上司や友人へのスピーチの依頼、乾杯の依頼
- どういうオーディエンスなのかを事前に伝えておく。スピーチ作成者の負担を軽減する。
招待状の作成、発送
料理の試食
- いろんな種類を食べるのでお腹ぱんぱん。
- 料理は1000~2000円程度の予算アップはした方が面白い。
ウェディングケーキのデザインのすり合わせ
テーブルクロス、花、引き出物、ネームプレートなどのすり合わせ
司会者との打ち合わせ
自分のスピーチの作成（新郎新婦ともに）
- 意外性と感動と笑いを織り交ぜるのが良い。
メッセージカードの作成
- スプレッドシートに書きなぐり、それを手書きでひたすら書ききるのみ。
式場音楽の選定
イベントの進行のすり合わせ
DVDの作成（オープニング・プロフィール・エンディング）
大量のピン札の調達（交通費は全額支給しろという親の教え）
ウェルカムボードの作成（ダイソーで4~500円くらいで材料は手に入る）

工夫した点

DVDを3枚自作した
- お願いすると結構お金がかかるし、理想形に近づけるための試行錯誤の回数も限られる。
- ハイスペックPCを持っているデータサイエンティストなら動画の編集に耐えられるはず。
- マックのiMovieを使えばGUIで簡単に動画作成ができる。
十分な大きさのメッセージカード
- わざわざ来てくれた友人との思い出をとにかく書きまくる。

結果としてよかったこと

DVDの自作
- 3枚とも笑いを提供でき、上映後は拍手している人もいたようです。（裏手にいたので会場の様子は直接見れなかったですが）
定型文を避けたスピーチ
- 書く前によくある定型文を見たんですが、面白くないと思い、独自のスピーチにしました。結果として大爆笑を提供できました。
普段から写真や動画を撮りまくっていたことで、DVD制作の素材が潤沢にあった。
- Googleフォトに写真をアップしまくれば容量に悩まずに済むので、とにかく日常的に撮り続けましょう。

やはり、結婚式はある意味でエンターテインメントなので、自分が工夫できるところ（スピーチとムービー）は少し頑張ってみるといいのかな、と思いました。

反省点

結婚式前日は有給休暇を取るべき
- 仕事を19時であがったけど、その後の準備で疲れ果てた。目にクマが若干できた。

結婚を頑張るエンジニアやデータサイエンティストの皆さんへ

結婚式の準備は色々と時間がかかります。
世の優秀なデータサイエンティストが、そのようなことに時間を割きすぎるのは社会的な損失なので、この参考情報を元に少しでも楽に準備をしていただけると幸いですね。

学習済み分散表現を用いた文書分類に挑戦（一部再学習も）

はじめに

2018年9月のテキストアナリティクスシンポジウムに行った際に、学習済みの分散表現で事前学習したモデルを使って分類してうまくいく事例が紹介されていました。
全てのタスクにおいてうまくいくとは思えませんが、試すコストはあまりかからないので試してみます。

2017年のテキストアナリティクスシンポジウムにおいても、メルカリやGunosyでは分散表現を用いた手法が一番精度が高いと言われていましたし、今年の会ではNLP系の学会でも分散表現はデファクトスタンダードになっているという話も伺いました。
2013~14年はLDAを使った研究が多かった気がしますが、徐々にシフトしていっているんですね。

これまで（Word2Vecを用いて蒙古タンメン中本の口コミ評価を予測してみる）は4000件程度の蒙古タンメン中本の口コミの情報を元に分散表現を手に入れていましたが、学習済みの分散表現を用いたアプローチも有効かもしれないと思い、試してみようと思います。

分類タスク

某グルメ口コミサイトの蒙古タンメン中本の口コミのテキストから、3.5点以上の評価かどうかを予測するタスクを扱います。
本当は、ポケモン図鑑の説明文から水やら炎やらのタイプを予測するとかをしたいのですが、あいにく手元にデータがないので、以前集めた蒙古タンメン中本の口コミを使います。（実は後日、ポケモン図鑑のデータを集めたのですが、平仮名にまみれたデータな上に、データ数も800件しかなかったので、どのみち厳しかったです。）

学習済み分散表現

Word2Vecなどで大量の文書をもとに学習させた分散表現のことを指します。
大規模コーパスで分散表現を手に入れる際は、数十GBにも相当するテキストデータを数時間かけて推定するので、学習済みのモデルは非常にありがたいです。（4年前に会社のPCで計算した際は、12時間くらいかかったこともありました。）

無料で提供してくださっている分散表現については、すでにこちらのブログで紹介されています。そこで紹介されているものに少し付け足すと、日本語の分散表現に関しては以下のようなものがあります。

白ヤギコーポレーションのモデル：Gensim
東北大学乾・岡崎研究室のモデル：Gensim
Facebookの学習済みFastTextモデル：Gensim
NWJC から取得した単語の分散表現データ (nwjc2vec)：Gensim
NNLM embedding trained on Google News：TensorFlow

そこで、今回は各種学習済み分散表現と蒙古タンメン中本コーパスで求めた分散表現の文書分類の性能バトルをしてみたいと思います。
ただ、分散表現ではなく、単語の頻度をもとに特徴量を作ったものが一番精度が高いのですが、分散表現同士の比較でもってどの学習済み分散表現が中本の口コミ分類に役に立ちそうなのかを明らかにしようと思います。（本来は分析という観点から即でボツですが、見苦しくも比較していきます。）

前処理

前処理は以下の通りで、テキストデータを分かち書きして、数値や低頻度・高頻度語を除外しています。

import pandas as pd

# データの読み込み
corpus_data = pd.read_pickle("nakamoto_corpus.pickle").reset_index(drop=True)

import collections
import MeCab
import mojimoji
from string import digits

remove_digits = str.maketrans('', '', digits)

tagger = MeCab.Tagger("-d /usr/local/lib/mecab/dic/mecab-ipadic-neologd/")

def nouns_extract(line):
    keyword=[]
    node = tagger.parseToNode(line).next
    while node:
        if node.feature.split(",")[0] == "名詞":
            keyword.append(node.surface)
        node = node.next
    keyword = str(keyword).replace("', '"," ")
    keyword = keyword.replace("\'","")
    keyword = keyword.replace("[","")
    keyword = keyword.replace("]","")
    return keyword


#欠損データを除外する関数
def FilterNANData(dataset, column_name):
    result_drop_remove = copy.deepcopy(dataset)
    result_drop_remove = result_drop_remove[~result_drop_remove[column_name].isnull()].reset_index(drop=True)
    return result_drop_remove[column_name]

#形態素解析して名詞のみを抽出し単語の頻度を集計して降順で返す関数
def NounceSum(wordlist):
    m = MeCab.Tagger("-d /usr/local/lib/mecab/dic/mecab-ipadic-neologd/")
    noun_list = [] # 重複を含めた名詞のリスト
    
    for i in wordlist:
        for l in m.parse (i).splitlines():
            try:
                if l != 'EOS' and l.split('\t')[1].split(',')[0] == '名詞': # EOSを除き名詞のみ抽出
                    noun_list.append(l.split('\t')[0]) # 見出し追加
            except:
                pass

    noun_cnt = collections.Counter(noun_list) # 各名詞の数え上げ
    nouns_data = pd.DataFrame.from_dict(noun_cnt, orient='index').reset_index()
    nouns_data.columns = ['nouns', 'count']
    nouns_data = nouns_data.sort_values(by=["count"], ascending=None)
    return nouns_data

#ストップワードを取り除く
def stop_word(documents):
    texts = [word for word in documents.lower().split() if word not in stoplist]
    texts = " ".join(texts)
    return texts

#リストでもらったテキストからストップワードを取り除く
def stop_word_tolist(documents):
    document_total = ' '.join(documents.tolist())
    texts = [[word for word in document_total.lower().split() if word not in stoplist]]
    return texts


#全角を半角にする
corpus_data["text_fixed"] =  list(map(lambda text: mojimoji.zen_to_han(text, kana=False) , corpus_data.text))
#数字を除外する
corpus_data["text_fixed"] =  list(map(lambda text: text.translate(remove_digits) , corpus_data.text_fixed))
#形態素解析する
corpus_data["text_wakati"] = list(map(lambda text:nouns_extract(text) , corpus_data.text_fixed))

#単語の頻度の計算
term_freq = pd.DataFrame(NounceSum(FilterNANData(corpus_data, "text_fixed")).reset_index(drop=True))
term_freq["ratio"] = term_freq["count"]/term_freq["count"].sum()

#不要語の除去（指定したものや、頻出のもの、頻度の低すぎるものを除外）
stoplist = ["ーー", "HP", "http://", "https://"]
stoplist2 = term_freq.query(' 10 > count | count > 1000  ').nouns.tolist()
stoplist.extend(stoplist2)
stoplist = set(stoplist)

corpus_data["text_wakati_fixed"] = list(map(lambda text:stop_word(text) , corpus_data.text_wakati))

#3.5点以上であれば1そうでなければ0
corpus_data["label"] = np.where(corpus_data.rating >= 3.5, 1, 0)

import pandas as pd

# データの読み込み

corpus_data = pd.read_pickle("nakamoto_corpus.pickle").reset_index(drop=True)

import collections

import MeCab

import mojimoji

from string import digits

remove_digits = str.maketrans('', '', digits)

tagger = MeCab.Tagger("-d /usr/local/lib/mecab/dic/mecab-ipadic-neologd/")

def nouns_extract(line):

keyword=[]

node = tagger.parseToNode(line).next

while node:

if node.feature.split(",")[0] == "名詞":

keyword.append(node.surface)

node = node.next

keyword = str(keyword).replace("', '"," ")

keyword = keyword.replace("\'","")

keyword = keyword.replace("[","")

keyword = keyword.replace("]","")

return keyword

#欠損データを除外する関数

def FilterNANData(dataset, column_name):

result_drop_remove = copy.deepcopy(dataset)

result_drop_remove = result_drop_remove[~result_drop_remove[column_name].isnull()].reset_index(drop=True)

return result_drop_remove[column_name]

#形態素解析して名詞のみを抽出し単語の頻度を集計して降順で返す関数

def NounceSum(wordlist):

m = MeCab.Tagger("-d /usr/local/lib/mecab/dic/mecab-ipadic-neologd/")

noun_list = [] # 重複を含めた名詞のリスト

for i in wordlist:

for l in m.parse (i).splitlines():

try:

if l != 'EOS' and l.split('\t')[1].split(',')[0] == '名詞': # EOSを除き名詞のみ抽出

noun_list.append(l.split('\t')[0]) # 見出し追加

except:

pass

noun_cnt = collections.Counter(noun_list) # 各名詞の数え上げ

nouns_data = pd.DataFrame.from_dict(noun_cnt, orient='index').reset_index()

nouns_data.columns = ['nouns', 'count']

nouns_data = nouns_data.sort_values(by=["count"], ascending=None)

return nouns_data

#ストップワードを取り除く

def stop_word(documents):

texts = [word for word in documents.lower().split() if word not in stoplist]

texts = " ".join(texts)

return texts

#リストでもらったテキストからストップワードを取り除く

def stop_word_tolist(documents):

document_total = ' '.join(documents.tolist())

texts = [[word for word in document_total.lower().split() if word not in stoplist]]

return texts

#全角を半角にする

corpus_data["text_fixed"] = list(map(lambda text: mojimoji.zen_to_han(text, kana=False) , corpus_data.text))

#数字を除外する

corpus_data["text_fixed"] = list(map(lambda text: text.translate(remove_digits) , corpus_data.text_fixed))

#形態素解析する

corpus_data["text_wakati"] = list(map(lambda text:nouns_extract(text) , corpus_data.text_fixed))

#単語の頻度の計算

term_freq = pd.DataFrame(NounceSum(FilterNANData(corpus_data, "text_fixed")).reset_index(drop=True))

term_freq["ratio"] = term_freq["count"]/term_freq["count"].sum()

#不要語の除去（指定したものや、頻出のもの、頻度の低すぎるものを除外）

stoplist = ["ーー", "HP", "http://", "https://"]

stoplist2 = term_freq.query(' 10 > count | count > 1000 ').nouns.tolist()

stoplist.extend(stoplist2)

stoplist = set(stoplist)

corpus_data["text_wakati_fixed"] = list(map(lambda text:stop_word(text) , corpus_data.text_wakati))

#3.5点以上であれば1そうでなければ0

corpus_data["label"] = np.where(corpus_data.rating >= 3.5, 1, 0)

処理を施すとこのようなデータになります。

特徴量は、scikit-learnのCountVectorizerやTfidfVectorizer、分散表現の合計・平均・TF-IDFを求めたものを用意します。

from gensim.models.word2vec import Word2Vec
from gensim.models import word2vec
import matplotlib.pyplot as plt
import seaborn as sns

from tabulate import tabulate
from collections import Counter, defaultdict
from sklearn.feature_extraction.text import CountVectorizer
from sklearn.feature_extraction.text import TfidfVectorizer
from sklearn.ensemble import ExtraTreesClassifier
from sklearn.pipeline import Pipeline
from sklearn.metrics import accuracy_score
from sklearn.cross_validation import cross_val_score

#単語の分散表現の合計値を求めるクラスの定義
class SumEmbeddingVectorizer(object):
    def __init__(self, word2vec):
        self.word2vec = word2vec
        self.dim = word2vec.values()
        self.dim = next(iter(self.dim))
        self.dim = self.dim.size
    
    def fit(self, X, y):
        return self 
 
    def transform(self, X):
        return np.array([
            np.sum([self.word2vec[w] for w in words if w in self.word2vec] 
                    or [np.zeros(self.dim)], axis=0)
            for words in X
        ])


#単語の分散表現の平均値を求めるクラスの定義
class MeanEmbeddingVectorizer(object):
    def __init__(self, word2vec):
        self.word2vec = word2vec
        self.dim = word2vec.values()
        self.dim = next(iter(self.dim))
        self.dim = self.dim.size
    
    def fit(self, X, y):
        return self 
 
    def transform(self, X):
        return np.array([
            np.mean([self.word2vec[w] for w in words if w in self.word2vec] 
                    or [np.zeros(self.dim)], axis=0)
            for words in X
        ])

#TF-IDFで重み付けした分散表現を求めるクラスの定義
class TfidfEmbeddingVectorizer(object):
    def __init__(self, word2vec):
        self.word2vec = word2vec
        self.word2weight = None
        self.dim = word2vec.values()
        self.dim = next(iter(self.dim))
        self.dim = self.dim.size
        
    def fit(self, X, y):
        tfidf = TfidfVectorizer(analyzer=lambda x: x)
        tfidf.fit(X)
        max_idf = max(tfidf.idf_)
        self.word2weight = defaultdict(
            lambda: max_idf, 
            [(w, tfidf.idf_[i]) for w, i in tfidf.vocabulary_.items()])
    
        return self
    
    def transform(self, X):
        return np.array([
                np.mean([self.word2vec[w] * self.word2weight[w]
                         for w in words if w in self.word2vec] or
                        [np.zeros(self.dim)], axis=0)
                for words in X
            ])

#入力変数と出力変数の指定
X, y = np.array(corpus_data.text_wakati_fixed), np.array(corpus_data.label)

#gensimで読み込むための形式にする
sentences = [token.split(" ") for token in corpus_data.text_wakati]

from gensim.models.word2vec import Word2Vec

from gensim.models import word2vec

import matplotlib.pyplot as plt

import seaborn as sns

from tabulate import tabulate

from collections import Counter, defaultdict

from sklearn.feature_extraction.text import CountVectorizer

from sklearn.feature_extraction.text import TfidfVectorizer

from sklearn.ensemble import ExtraTreesClassifier

from sklearn.pipeline import Pipeline

from sklearn.metrics import accuracy_score

from sklearn.cross_validation import cross_val_score

#単語の分散表現の合計値を求めるクラスの定義

class SumEmbeddingVectorizer(object):

def __init__(self, word2vec):

self.word2vec = word2vec

self.dim = word2vec.values()

self.dim = next(iter(self.dim))

self.dim = self.dim.size

def fit(self, X, y):

return self

def transform(self, X):

return np.array([

np.sum([self.word2vec[w] for w in words if w in self.word2vec]

or [np.zeros(self.dim)], axis=0)

for words in X

])

#単語の分散表現の平均値を求めるクラスの定義

class MeanEmbeddingVectorizer(object):

def __init__(self, word2vec):

self.word2vec = word2vec

self.dim = word2vec.values()

self.dim = next(iter(self.dim))

self.dim = self.dim.size

def fit(self, X, y):

return self

def transform(self, X):

return np.array([

np.mean([self.word2vec[w] for w in words if w in self.word2vec]

or [np.zeros(self.dim)], axis=0)

for words in X

])

#TF-IDFで重み付けした分散表現を求めるクラスの定義

class TfidfEmbeddingVectorizer(object):

def __init__(self, word2vec):

self.word2vec = word2vec

self.word2weight = None

self.dim = word2vec.values()

self.dim = next(iter(self.dim))

self.dim = self.dim.size

def fit(self, X, y):

tfidf = TfidfVectorizer(analyzer=lambda x: x)

tfidf.fit(X)

max_idf = max(tfidf.idf_)

self.word2weight = defaultdict(

lambda: max_idf,

[(w, tfidf.idf_[i]) for w, i in tfidf.vocabulary_.items()])

return self

def transform(self, X):

return np.array([

np.mean([self.word2vec[w] * self.word2weight[w]

for w in words if w in self.word2vec] or

[np.zeros(self.dim)], axis=0)

for words in X

])

#入力変数と出力変数の指定

X, y = np.array(corpus_data.text_wakati_fixed), np.array(corpus_data.label)

#gensimで読み込むための形式にする

sentences = [token.split(" ") for token in corpus_data.text_wakati]

蒙古タンメン中本の口コミ4000件から作成した分散表現：Gensim

まず、以前のブログで紹介した蒙古タンメン中本の分散表現ですが、以下のように推定しています。

#Word2Vecを実行する。 
model = Word2Vec(sentences,
                 sg=1,
                 size=50,
                 window=5, 
                 min_count=5,
                 workers=2,
                 seed=123)

#単語ごとの分散表現を手に入れる。
w2v = {w: vec for w, vec in zip(model.wv.index2word, model.wv.syn0)}

#Word2Vecを実行する。

model = Word2Vec(sentences,

sg=1,

size=50,

window=5,

min_count=5,

workers=2,

seed=123)

#単語ごとの分散表現を手に入れる。

w2v = {w: vec for w, vec in zip(model.wv.index2word, model.wv.syn0)}

Pipelineを用いてExtraTreesClassifierによる学習をします。特徴量は先程あげた、テキストベースのCountVectorizerやTfidfVectorizer、分散表現の合計・平均・TF-IDFで、評価指標はAUCのクロスバリデーションスコアとします。

#ベースラインとなる既存手法のモデルの準備
#etree
etree = Pipeline([("count_vectorizer",CountVectorizer(analyzer=lambda x: x)), 
                    ("extra trees", ExtraTreesClassifier(n_estimators=200))])
#etreeのTF-IDF版
etree_tfidf = Pipeline([("tfidf_vectorizer", TfidfVectorizer(analyzer=lambda x: x)),
                          ("extra trees", ExtraTreesClassifier(n_estimators=200))])

#Word2Vecを特徴量としてExtraTreesによる分類器を準備する。
etree_w2v_sum = Pipeline([("word2vec vectorizer", SumEmbeddingVectorizer(w2v)), 
                        ("extra trees", ExtraTreesClassifier(n_estimators=200))])

etree_w2v = Pipeline([("word2vec vectorizer", MeanEmbeddingVectorizer(w2v)), 
                        ("extra trees", ExtraTreesClassifier(n_estimators=200))])
 
etree_w2v_tfidf = Pipeline([("word2vec vectorizer", TfidfEmbeddingVectorizer(w2v)), 
                        ("extra trees", ExtraTreesClassifier(n_estimators=200))])
 
#各モデルを実行し、クロスバリデーションスコアを計算し、出力させる。
all_models = [
    ("etree",etree),
    ("etree_tfidf",etree_tfidf),
    ("w2v_sum", etree_w2v_sum),    
    ("w2v", etree_w2v),
    ("w2v_tfidf", etree_w2v_tfidf)
]

scores = sorted([(name, cross_val_score(model, X, y, cv=5, scoring = "roc_auc").mean()) 
                 for name, model in all_models], 
                key = lambda x:x[0])

print(tabulate(scores, floatfmt=".4f", headers=("model", 'score')))

#ベースラインとなる既存手法のモデルの準備

#etree

etree = Pipeline([("count_vectorizer",CountVectorizer(analyzer=lambda x: x)),