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

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

はじめに

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)),

("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')))

汗に関してコンテキストの似ている単語を抽出しています。

結果は、以下の通りで、分散表現を使わない方がAUCが高いです。ただ、w2v_tfidf（分散表現のTF-IDFを特徴量にしたもの）が分散表現の中でAUCが高いようです。今回はこの60.5%をベースラインに比較していこうと思います。

#中本の分散表現
model          score
-----------  -------
etree         0.6340
etree_tfidf   0.6499
w2v           0.5955
w2v_sum       0.5765
w2v_tfidf     0.6051

#中本の分散表現

model score

----------- -------

etree 0.6340

etree_tfidf 0.6499

w2v 0.5955

w2v_sum 0.5765

w2v_tfidf 0.6051

白ヤギコーポレーションのモデル：Gensim

こちらのリンク、「word2vecの学習済み日本語モデルを公開します」から、ダウンロードしてそのまま以下のコードでモデルを扱えます。

from gensim.models.word2vec import Word2Vec

model_path = "latest-ja-word2vec-gensim-model/word2vec.gensim.model"
model_shiroyagi = Word2Vec.load(model_path)
w2v_shiroyagi = {w: vec for w, vec in zip(model_shiroyagi.wv.index2word, model_shiroyagi.wv.syn0)}

from gensim.models.word2vec import Word2Vec

model_path = "latest-ja-word2vec-gensim-model/word2vec.gensim.model"

model_shiroyagi = Word2Vec.load(model_path)

w2v_shiroyagi = {w: vec for w, vec in zip(model_shiroyagi.wv.index2word, model_shiroyagi.wv.syn0)}

汗の関連語を抽出していますが、中国の歴史の何かですか？可汗とかいう単語は聞いたことあるかも。

まずは白ヤギさんの分散表現をそのまま使って予測してみます。（コードは先程のものとほぼ重複するので省略しています。）
残念ながら、ベースラインの60.5%には至りませんでした。

#学習済み分散表現をそのまま用いたもの
model          score
-----------  -------
w2v           0.6007
w2v_sum       0.5648
w2v_tfidf     0.5744

#学習済み分散表現をそのまま用いたもの

model score

----------- -------

w2v 0.6007

w2v_sum 0.5648

w2v_tfidf 0.5744

hogehoge.modelというフルモデル形式の場合は、再学習が可能です。詳しくはこちら（models.word2vec – Word2vec embeddings model）に書かれています。

train(sentences=None,
          corpus_file=None,
          total_examples=None, 
          total_words=None, 
          epochs=None, 
          start_alpha=None, 
          end_alpha=None, 
          word_count=0,
          queue_factor=2,
          report_delay=1.0,
          compute_loss=False,
          callbacks=())

train(sentences=None,

corpus_file=None,

total_examples=None,

total_words=None,

epochs=None,

start_alpha=None,

end_alpha=None,

word_count=0,

queue_factor=2,

report_delay=1.0,

compute_loss=False,

callbacks=())

今回は、白ヤギさんの分散表現に対して、追加で蒙古タンメン中本のテキストを食わせて再学習させます。

from gensim.models.word2vec import Word2Vec

model_path = "latest-ja-word2vec-gensim-model/word2vec.gensim.model"
model_shiroyagi = Word2Vec.load(model_path)

#リストの重複を除く
def remove_duplicates(l):
    return list(set(l))

#ユニーク単語のリストを作る
flat_list = [item for sublist in sentences for item in sublist]

model_shiroyagi.train(sentences=sentences,
                    total_examples=len(sentences),
                    total_words=len(remove_duplicates(flat_list)) ,
                    word_count= len(model_shiroyagi.wv.index2word),
                    epochs=4)

w2v_shiroyagi = {w: vec for w, vec in zip(model_shiroyagi.wv.index2word, model_shiroyagi.wv.syn0)}

from gensim.models.word2vec import Word2Vec

model_path = "latest-ja-word2vec-gensim-model/word2vec.gensim.model"

model_shiroyagi = Word2Vec.load(model_path)

#リストの重複を除く

def remove_duplicates(l):

return list(set(l))

#ユニーク単語のリストを作る

flat_list = [item for sublist in sentences for item in sublist]

model_shiroyagi.train(sentences=sentences,

total_examples=len(sentences),

total_words=len(remove_duplicates(flat_list)) ,

word_count= len(model_shiroyagi.wv.index2word),

epochs=4)

w2v_shiroyagi = {w: vec for w, vec in zip(model_shiroyagi.wv.index2word, model_shiroyagi.wv.syn0)}

ベースラインの60.5%よりも下回り、さきほどの白ヤギさんのもともとの分散表現よりも下回りました。

model          score
-----------  -------
w2v           0.5731
w2v_sum       0.5515
w2v_tfidf     0.5774

model score

----------- -------

w2v 0.5731

w2v_sum 0.5515

w2v_tfidf 0.5774

再学習してもかえって精度が下がったりすることから、簡単に精度が出るわけではなさそうです。まぁ、理想はその適用領域での大量のテキストデータがあることで、Wikipediaを元に作成した分散表現に強く依存しても駄目なのだろうと思われます。

東北大学乾・岡崎研究室のモデル：Gensim

日本語 Wikipedia エンティティベクトルからダウンロードした学習済み分散表現を用います。ダウンロード後は普通に「gzip -d file.txt.gz」みたいにターミナル上で解凍します。以下のコードを実行すればすぐに使うことができます。
ただし、KeyedVectors形式のものは白ヤギさんのように再学習ができません。（Why use KeyedVectors instead of a full model?）

import gensim.models.keyedvectors as word2vec_for_txt

#学習済み分散表現の読み込み
file_name = "jawiki.all_vectors.100d.txt"
model_tonpei = word2vec_for_txt.Word2VecKeyedVectors.load_word2vec_format(file_name)

＃単語の分散表現の取得
w2v_tonpei = {w: vec for w, vec in zip(model_tonpei.wv.index2word, model_tonpei.wv.syn0)}

import gensim.models.keyedvectors as word2vec_for_txt

#学習済み分散表現の読み込み

file_name = "jawiki.all_vectors.100d.txt"

model_tonpei = word2vec_for_txt.Word2VecKeyedVectors.load_word2vec_format(file_name)

＃単語の分散表現の取得

w2v_tonpei = {w: vec for w, vec in zip(model_tonpei.wv.index2word, model_tonpei.wv.syn0)}

汗の類似語に関しては、難しい単語が高めに出ているようです。

results = model_tonpei.most_similar(positive="汗", topn=10)

for result in results:
    print(result[0], '\t', result[1])

曷多   0.6981213092803955
拭い   0.6970207691192627
びっしょり    0.6847031116485596
垢    0.6723113059997559
喝    0.6681551933288574
裴羅   0.662636935710907
拭く   0.6584706902503967
洗う   0.6570826768875122
汚    0.6568220257759094
小便   0.6558053493499756

results = model_tonpei.most_similar(positive="汗", topn=10)

for result in results:

print(result[0], '\t', result[1])

曷多 0.6981213092803955

拭い 0.6970207691192627

びっしょり 0.6847031116485596

垢 0.6723113059997559

喝 0.6681551933288574

裴羅 0.662636935710907

拭く 0.6584706902503967

洗う 0.6570826768875122

汚 0.6568220257759094

小便 0.6558053493499756

残念ながら、ベースラインの60.5%には至りませんでした。

model          score
-----------  -------

w2v           0.5901
w2v_sum       0.5478
w2v_tfidf     0.5786

model score

----------- -------

w2v 0.5901

w2v_sum 0.5478

w2v_tfidf 0.5786

Facebookの学習済みFastTextモデル：Gensim

FastTextはGoogleにいたTomas Mikolov氏がFacebookに転職されて作られた分散表現を求めるためのモデルです。Gensimでも呼び出せます。学習済みのものはこちらのGitHub（Pre-trained word vectors）にあるのですが、NEologdで形態素解析したものをベースに学習し公開されている方がいるとのことで、こちら（fastTextの学習済みモデルを公開しました）からダウンロードしたものを使わせていただきました。

file_name = "model.vec"
model_fasttext = word2vec_for_txt.KeyedVectors.load_word2vec_format(file_name, binary=False)

w2v_fasttext = {w: vec for w, vec in zip(model_fasttext.wv.index2word, model_fasttext.wv.syn0)}

file_name = "model.vec"

model_fasttext = word2vec_for_txt.KeyedVectors.load_word2vec_format(file_name, binary=False)

w2v_fasttext = {w: vec for w, vec in zip(model_fasttext.wv.index2word, model_fasttext.wv.syn0)}

何だこれはレベルの結果が返ってきました。中国の歴史上の人物か何かなんでしょうか。

results = model_fasttext.most_similar(positive="汗", topn=10)

for result in results:
    print(result[0], '\t', result[1])

処羅可汗     0.6946621537208557
突利   0.6881076097488403
乙息記可汗    0.6672964096069336
都藍可汗     0.6537224650382996
頡利   0.6425886154174805
葉護可汗     0.6401641368865967
忠貞可汗     0.6395537257194519
撅    0.6240820288658142
默啜   0.6238714456558228
テュルギシュ   0.61381596326828

results = model_fasttext.most_similar(positive="汗", topn=10)

for result in results:

print(result[0], '\t', result[1])

処羅可汗 0.6946621537208557

突利 0.6881076097488403

乙息記可汗 0.6672964096069336

都藍可汗 0.6537224650382996

頡利 0.6425886154174805

葉護可汗 0.6401641368865967

忠貞可汗 0.6395537257194519

撅 0.6240820288658142

默啜 0.6238714456558228

テュルギシュ 0.61381596326828

若干ですがベースラインの60.5%よりも良い結果が得られましたが、誤差の範囲な気がします。

model          score
-----------  -------
w2v           0.6064
w2v_sum       0.5701
w2v_tfidf     0.5977

model score

----------- -------

w2v 0.6064

w2v_sum 0.5701

w2v_tfidf 0.5977

NWJC から取得した単語の分散表現データ (nwjc2vec)：Gensim

国立国語研究所の収集されたテキストデータを元に学習した分散表現が提供されています。ただし、利用するためには申請する必要があります。申請が受理されたらこちら（NWJC から取得した単語の分散表現データ (nwjc2vec) を頒布）からダウンロードして使えます。

file_name = "nwjc_word_skip_300_8_25_0_1e4_6_1_0_15.txt.vec"
model_bonten = word2vec_for_txt.KeyedVectors.load_word2vec_format(file_name, binary=False)

w2v_bonten = {w: vec for w, vec in zip(model_bonten.wv.index2word, model_bonten.wv.syn0)}

file_name = "nwjc_word_skip_300_8_25_0_1e4_6_1_0_15.txt.vec"

model_bonten = word2vec_for_txt.KeyedVectors.load_word2vec_format(file_name, binary=False)

w2v_bonten = {w: vec for w, vec in zip(model_bonten.wv.index2word, model_bonten.wv.syn0)}

汗の関連語ですが、うまく関連付けれているように思われます。少なくとも中国史ぽくはありません。しかしながら、顔文字まで学習していたとは。

results = model_bonten.most_similar(positive="汗", topn=10)

for result in results:
    print(result[0], '\t', result[1])

大汗   0.8739632368087769
滝汗   0.7782824039459229
冷汗   0.7470424175262451
苦笑   0.7295455932617188
￣∀￣;)    0.6987063884735107
笑    0.6849342584609985
冷や汗      0.6749577522277832
脂汗   0.6701192259788513
泣    0.6603907942771912
ダクダク     0.6596677303314209

results = model_bonten.most_similar(positive="汗", topn=10)

for result in results:

print(result[0], '\t', result[1])

大汗 0.8739632368087769

滝汗 0.7782824039459229

冷汗 0.7470424175262451

苦笑 0.7295455932617188

￣∀￣;) 0.6987063884735107

笑 0.6849342584609985

冷や汗 0.6749577522277832

脂汗 0.6701192259788513

泣 0.6603907942771912

ダクダク 0.6596677303314209

ベースラインの60.5%よりも1%ポイントほど高い結果となりました。

model          score
-----------  -------
w2v           0.6163
w2v_sum       0.5839
w2v_tfidf     0.5916

model score

----------- -------

w2v 0.6163

w2v_sum 0.5839

w2v_tfidf 0.5916

NNLM embedding trained on Google News：TensorFlow

こちら（tensorflow-hubで超簡単にテキスト分類モデルが作成できる）で紹介されているように、GoogleがTensorFlowでGoogleニュースのテキストをもとに学習した分散表現が提供されています。

こちらのGitHub（NNLM embedding trained on Google News）から、Japaneseのnnlm-ja-dim50、nnlm-ja-dim50-with-normalizationなどが使えます。分散表現の説明についてはこちらのドキュメント（Token based text embedding trained on Japanese Google News 6B corpus.）にあります。

import tensorflow_hub as hub
import tensorflow as tf
from sklearn.utils import shuffle

embed = hub.Module("https://tfhub.dev/google/nnlm-ja-dim50/1")

df = copy.deepcopy(corpus_data)
df['category_id'] = df.label.factorize()[0]
df = shuffle(df)


train_input_fn = tf.estimator.inputs.pandas_input_fn(
    df[:3500], df[:3500]["category_id"], num_epochs=None, shuffle=True)

predict_test_input_fn = tf.estimator.inputs.pandas_input_fn(
    df[3500:], df[3500:]["category_id"], shuffle=False)

embedded_text_feature_column = hub.text_embedding_column(
    key="text_wakati", 
    module_spec="https://tfhub.dev/google/nnlm-ja-dim50/1")

estimator = tf.estimator.DNNClassifier(
    hidden_units=[500, 100],
    feature_columns=[embedded_text_feature_column],
    n_classes=2,
    optimizer=tf.train.AdagradOptimizer(learning_rate=0.003))

estimator.train(input_fn=train_input_fn, steps=1000);
test_eval_result = estimator.evaluate(input_fn=predict_test_input_fn)

print("Test set accuracy: {accuracy}".format(**test_eval_result))

import tensorflow_hub as hub

import tensorflow as tf

from sklearn.utils import shuffle

embed = hub.Module("https://tfhub.dev/google/nnlm-ja-dim50/1")

df = copy.deepcopy(corpus_data)

df['category_id'] = df.label.factorize()[0]

df = shuffle(df)

train_input_fn = tf.estimator.inputs.pandas_input_fn(

df[:3500], df[:3500]["category_id"], num_epochs=None, shuffle=True)

predict_test_input_fn = tf.estimator.inputs.pandas_input_fn(

df[3500:], df[3500:]["category_id"], shuffle=False)

embedded_text_feature_column = hub.text_embedding_column(

key="text_wakati",

module_spec="https://tfhub.dev/google/nnlm-ja-dim50/1")

estimator = tf.estimator.DNNClassifier(

hidden_units=[500, 100],

feature_columns=[embedded_text_feature_column],

n_classes=2,

optimizer=tf.train.AdagradOptimizer(learning_rate=0.003))

estimator.train(input_fn=train_input_fn, steps=1000);

test_eval_result = estimator.evaluate(input_fn=predict_test_input_fn)

print("Test set accuracy: {accuracy}".format(**test_eval_result))

AUCが65%となっているものの、先程のsklearnでのクロスバリデーションのものとの比較ではないので、なんとも言えないですが、Googleニュースのデータだし結構精度が出そうな可能性を感じますね。
今後、TensorFlowでクロスバリデーションによるAUCスコアの出し方を調べてみて、順当に比較できるようにしたいです。（Kerasを使って計算している事例は見つけた。）

{'accuracy': 0.62835246,
 'accuracy_baseline': 0.578544,
 'auc': 0.6501204,
 'auc_precision_recall': 0.577037,
 'average_loss': 0.64989257,
 'global_step': 1000,
 'label/mean': 0.42145595,
 'loss': 67.848785,
 'precision': 0.5855263,
 'prediction/mean': 0.4055166,
 'recall': 0.40454546}

{'accuracy': 0.62835246,

'accuracy_baseline': 0.578544,

'auc': 0.6501204,

'auc_precision_recall': 0.577037,

'average_loss': 0.64989257,

'global_step': 1000,

'label/mean': 0.42145595,

'loss': 67.848785,

'precision': 0.5855263,

'prediction/mean': 0.4055166,

'recall': 0.40454546}

比較

今回の分類タスクはそもそも分散表現では精度が出なかったのですが、学習済み分散表現の中で序列を作るとすると、梵天が一番良く、FastTextが少しだけ良かったです。
TensorFlowをほぼ業務で使わないので、Googleニュースの分散表現を今回の比較対象にできなかったのですが、後日比較できるようにしたいと思います。

あと、今回の口コミの点数を当てるタスクよりも、分散表現にとって相性がいいタスクがあるかもしれないので、今回の結果で諦めることなく色々と試して行きたいです。

おわりに

様々なシンポジウムなどでスタンダードとなってきた分散表現ですが、学習済み分散表現をそのまま使って分類問題で役に立つのかを見てきました。残念ながら、口コミの評価予測タスクにおいては全然効果がなさそうでした。ただ、分散表現の中でもタスクによって相性の良い学習済み分散表現がありそうです。
先程も述べたように、理想は大量のテキストデータで学習した分散表現を求め、それを予測に使うことなので大量のテキストデータを集めて再チャレンジしたいです。どれくらいのテキストデータがあれば十分なのかの規模感もわからないので、実践あるのみなんですかね。

参考情報

Word Embeddingだけで文書分類する
 tensorflow-hubで超簡単にテキスト分類モデルが作成できる
 Error: ”Word2vec’ object has no attribute index2word
Word2vec Tutorial Online training / Resuming training
Word Embeddingモデル再訪
 Googleの事前学習済みモデルを手軽に利用出来るTensorFlow Hub
ゼロから作るDeep Learning ❷ ―自然言語処理編

pyreaperで音声データのピッチを掴むためのメモ書き（F0の抽出）

はじめに

最近、仕事で音声データを扱うことがあるのですが、アプローチの仕方を色々と調べています。

今のところ検討しているのが、

クラウドAPIなどによる文字起こしによるテキストアナリティクス
音声の波形に関する分析
音声のピッチに関する分析

の3つですが、1つ目と2つ目については比較的すぐにググってチャレンジできましたし、ビジネス上で役に立ちそうな発見もありました。

そこで今回は3つ目の「音声のピッチ」（音程）について扱おうと思います。

ピッチがわかると何が嬉しいかというと、顧客対応などの音声でピッチがどのようなものであれば顧客満足度が高いのかを分析できるようになり、担当者へのフィードバックに使える可能性があります。
今までなんとなく蓄積していた音声データが宝の山になるかどうか、それをジャッジするための一つの手段と考えることができます。

後輩からの聞き込み

大学院で動画の圧縮について研究していた元後輩社員くんがいたので、きっと知っているはずだろうと思い、その方からピッチに関する研究を進める上でのキーワードを教えてもらいました。
彼曰く、F1やF0というキーワードで探していけばピッチの抽出ができるそうなので、早速調べて実践できないかチャレンジします。

F0とF1

Wikipediaで調べて結果を以下に記します。

F0：基本周波数
信号を正弦波の合成で表したときの最も低い周波数成分の周波数を基本周波数と呼ぶ。
F1：フォルマント周波数
フォルマント：言葉を発している人の音声スペクトルを観察すると分かる、時間的に移動している複数のピークのこと。
周波数の低い順に、第1フォルマント、第2フォルマントというように数字を当てて呼び、それぞれF1、F2とも表記する。
観察の仕方としては、縦軸に周波数、横軸に時間を置くものとする。

ピッチはどうすれば求めることができるのか？

断創録さんのブログによると、ピッチとは声の高さのことで、基本周波数やF0とも呼ばれているそうです。
そのため、今回はF0を求めれば目的を果たすことができそうです。

今回扱うモジュール／今回のツール

今回は、ピッチ検出を行うためのツール、REAPERのPythonラッパーであるpyreaperを使います。

pip install pyreaper

1	pip install pyreaper

一発で入るはずです。

続いて、音声ファイルをWAVファイルに変換するためのツールとして、FFmpegをインストールします。
Macだと、

brew install ffmpeg

1	brew install ffmpeg

で一発で入るはずです。
使い方としては、

ffmpeg -i ファイル名.m4a -ar 16000 ファイル名.wav

1	ffmpeg -i ファイル名.m4a -ar 16000 ファイル名.wav

のようにターミナル上で使います。

データ

今回は、無料の音声素材でなにか無いかと探していたのですが、NHKの百人一首の読み上げデータにしました。
句に関してはこだわりはないので、適当に第5句にしています。

NHKの百人一首第5句
奥山に紅葉踏み分け鳴く鹿の　声聞く時ぞ秋は悲しき／猿丸大夫

もう秋ですし、風情があって良いですね。

プログラム

実行プログラムはご丁寧にpyreaperのPyPIのページに良いものがありましたので、こちらをまるごと使います。

まず、百人一首の音声がm4a形式だったので、wav形式になるように、
先程インストールしたFFmpegを用いて、サンプルレート16000、チャネル数1のWAVファイルに変換します。

ffmpeg -i speech_example.m4a -ar 16000 -map_channel 0.0.0 speech_example.wav

1	ffmpeg -i speech_example.m4a -ar 16000 -map_channel 0.0.0 speech_example.wav

なお、電話などの話者が2名いるステレオのデータの場合、FFmpegを用いて左側成分、右側成分だけを抽出することができますので、分けて解析すると良いと思います。左右の分割に関してはFFmpegの参考文献で紹介されています。そうすれば、片方の話者の声だけを抽出することができますので、解析が捗ります。

%matplotlib inline
import pysptk
import pyreaper
from scipy.io import wavfile
from IPython.display import Audio

file = "D0002110105_00000_A_001.wav"
fs, x = wavfile.read(file)

plot(x)
Audio(x, rate=fs)

#pyreaperモジュールを使ってピッチマークやF0や相関係数を計算する。
pm_times, pm, f0_times, f0, corr = pyreaper.reaper(x, fs)

plot(pm_times, pm, linewidth=3, color="red", label="Pitch mark")
legend(fontsize=15)

plot(f0_times, f0, linewidth=3, color="green", label="F0 contour")
legend(fontsize=15)

plot(f0_times, corr, linewidth=3, color="blue", label="Correlations")
legend(fontsize=15)

%matplotlib inline

import pysptk

import pyreaper

from scipy.io import wavfile

from IPython.display import Audio

file = "D0002110105_00000_A_001.wav"

fs, x = wavfile.read(file)

plot(x)

Audio(x, rate=fs)

#pyreaperモジュールを使ってピッチマークやF0や相関係数を計算する。

pm_times, pm, f0_times, f0, corr = pyreaper.reaper(x, fs)

plot(pm_times, pm, linewidth=3, color="red", label="Pitch mark")

legend(fontsize=15)

plot(f0_times, f0, linewidth=3, color="green", label="F0 contour")

legend(fontsize=15)

plot(f0_times, corr, linewidth=3, color="blue", label="Correlations")

legend(fontsize=15)

実行結果

百人一首の音声の波形です。

ピッチマークです。ピッチマークとは、「単位波形を抽出する基準点」とされているとの説明があるのですが、あまり良くわかりません。無音のときに生成されていないので、ピッチ抽出の確信度という意味なんでしょうか。

F0です。音程とみなすことができるので、こちらを会話の分析の際に使ってみるのが良いのだと思います。

相関係数です。こちらもピッチ検出に使う指標のようです。

今後の使いみち

F0を集計してある水準のものが多い、時系列で見ると前半と後半でF0の水準が違っているというような情報と、その会話におけるラベル（会話の結果としての商談の成否）について分析することで、商談などがうまくいく際のヒントが得られるかもしれないですね。

追記

数十分の音声に関してF0を抽出しようとしたら、処理に非常に時間がかかり、結局最初の数分とかのF0の抽出にとどまりました。予感はしていましたが、結構計算リソースが必要な領域なのですね。

参考情報

IPythonデータサイエンスクックブック ―対話型コンピューティングと可視化のためのレシピ集
 音声生成の基礎と音声学
 Is there a way to convert audio files in Mac OS X or the command line without using iTunes?
やる夫で学ぶディジタル信号処理東北大学大学院情報科学研究科鏡慎吾
 [vDSP][信号処理]オーディオ・音声分析への道7 FFT 自己相関関数ピッチ検出
 REAPER: Robust Epoch And Pitch EstimatoR
FFmpeg wiki:AudioChannelManipulation
基本周波数抽出器 REAPERによる抽出実験
 SPTKの使い方 (3) ピッチ抽出

蒙古タンメン中本コーパスに対してのLDAの適用とトピック数の探索

モチベーション

前回の記事では、Webスクレイピングにより入手した、蒙古タンメン中本の口コミデータに関して、Word2Vecを適用した特徴量エンジニアリングの事例を紹介しました。
今回はせっかく興味深いデータがあるので、どのようなトピックがあるのかをLDAを適用したいと思います。加えて、これまで記事で扱ってきたLDAの事例では評価指標であるPerplexityやCoherenceを扱ってこなかったことから、トピック数がどれくらいであるべきなのか、考察も含めて行いたいと思います。以前扱った階層ディリクレ過程であれば、トピック数を事前に決める必要が無いのですが、今回は扱わないものとします。

環境

・MacBook Pro
・Python3.5
・R version 3.4.4

Gensimで行うLDA

今回もPythonのGensimライブラリを用いて行います。

パープレキシティ
- テストデータに対して計算
- 負の対数尤度で、低いほどよい。
  - パープレキシティが低いと、高い精度で予測できるよい確率モデルと見なされる。汎化能力を表す指標。
  - トピックの数をいくらでも増やせばパープレキシティは下がる傾向が出ている。
  - 教科書でのパープレキシティの事例に関しては、トピック数を増やせば低くなるという傾向が出ている。

以下のコードでパープレキシティを計算します。

import pandas as pd
from gensim import corpora, models
import  gensim as  gensim

import matplotlib.pyplot as plt
import numpy as np
from ipywidgets import FloatProgress
from IPython.display import display, clear_output

nakamoto_corpus = pd.read_pickle("nakamoto_corpus.pickle")

output_list = pd.DataFrame(columns=["y","x"])
plt.ion()
fig = plt.figure()
axe = fig.add_subplot(111)

for iterations in range(2, 100):
                        clear_output(wait = True)

                        nakamoto_test = nakamoto_corpus.sample(frac=0.1, replace=True)
                        nakamoto_train= nakamoto_corpus[~nakamoto_corpus.index.isin(nakamoto_test.index)]

                        texts = [ ]
                        for line in nakamoto_train.text_wakati:
                            texts.append(line.split())

                        texts_test = [ ]
                        for line in nakamoto_test.text_wakati:
                            texts_test.append(line.split())

                        # 辞書作成
                        dictionary = corpora.Dictionary(texts)
                        dictionary.filter_extremes(no_below=20, no_above=0.3)

                        # コーパスを作成
                        corpus = [dictionary.doc2bow(text) for text in texts]
                        corpus_test = [dictionary.doc2bow(text) for text in texts_test]

                        # LDA の計算
                        topic_N = 1 + iterations
                        parameters = topic_N
                        lda = gensim.models.ldamodel.LdaModel(
                            corpus=corpus,
                            alpha='auto',
                            num_topics=topic_N,
                            id2word=dictionary
                        )

                        d = {'y':lda.log_perplexity(chunk=corpus_test),'x':parameters}
                        df = pd.DataFrame(data=d,index=[0])
                        output_list = output_list.append(df)
                        
                        axe.plot(output_list.x,output_list.y)
                        fig.set_size_inches(8, 8)
                        display(fig)
                        axe.cla()
                        
                        for i in range(topic_N):
                            print('TOPIC:', i, '__', lda.print_topic(i))

import pandas as pd

from gensim import corpora, models

import gensim as gensim

import matplotlib.pyplot as plt

import numpy as np

from ipywidgets import FloatProgress

from IPython.display import display, clear_output

nakamoto_corpus = pd.read_pickle("nakamoto_corpus.pickle")

output_list = pd.DataFrame(columns=["y","x"])

plt.ion()

fig = plt.figure()

axe = fig.add_subplot(111)

for iterations in range(2, 100):

clear_output(wait = True)

nakamoto_test = nakamoto_corpus.sample(frac=0.1, replace=True)

nakamoto_train= nakamoto_corpus[~nakamoto_corpus.index.isin(nakamoto_test.index)]

texts = [ ]

for line in nakamoto_train.text_wakati:

texts.append(line.split())

texts_test = [ ]

for line in nakamoto_test.text_wakati:

texts_test.append(line.split())

# 辞書作成

dictionary = corpora.Dictionary(texts)

dictionary.filter_extremes(no_below=20, no_above=0.3)

# コーパスを作成

corpus = [dictionary.doc2bow(text) for text in texts]

corpus_test = [dictionary.doc2bow(text) for text in texts_test]

# LDA の計算

topic_N = 1 + iterations

parameters = topic_N

lda = gensim.models.ldamodel.LdaModel(

corpus=corpus,

alpha='auto',

num_topics=topic_N,

id2word=dictionary

)

d = {'y':lda.log_perplexity(chunk=corpus_test),'x':parameters}

df = pd.DataFrame(data=d,index=[0])

output_list = output_list.append(df)

axe.plot(output_list.x,output_list.y)

fig.set_size_inches(8, 8)

display(fig)

axe.cla()

for i in range(topic_N):

print('TOPIC:', i, '__', lda.print_topic(i))

実際に、中本コーパスで計算したトピック数に対してのパープレキシティは以下のように推移しました。

Ldaのモデル選択におけるperplexityの評価によると、
”複数のトピック数で比べて、Perplexityが最も低いものを選択する。」という手法は人間にとって有益なモデルを選択するのに全く役に立たない可能性がある。”と記されています。

『トピックモデルによる統計的潜在意味解析』には、”識別問題の特徴量として使う場合は識別問題の評価方法で決定すればよい”とあるので、目的によってはパープレキシティにこだわらなくても良いと思われます。

今回のケースだと、パープレキシティだけだと、決めかねてしまいますね。

コヒーレンス
- トピックごとの単語間類似度の平均
- トピック全体のコヒーレンスが高ければ、良い学習アルゴリズムとみなす。

以下のコードでコヒーレンスを計算します。

output_list = pd.DataFrame(columns=["y","x"])
plt.ion()
fig = plt.figure()
axe = fig.add_subplot(111)

for iterations in range(2, 100):
                        clear_output(wait = True)

                        nakamoto_test = nakamoto_corpus.sample(frac=0.1, replace=True)
                        nakamoto_train= nakamoto_corpus[~nakamoto_corpus.index.isin(nakamoto_test.index)]

                        texts = [ ]
                        for line in nakamoto_train.text_wakati:
                            texts.append(line.split())

                        texts_test = [ ]
                        for line in nakamoto_test.text_wakati:
                            texts_test.append(line.split())

                        # 辞書作成
                        dictionary = corpora.Dictionary(texts)
                        dictionary.filter_extremes(no_below=20, no_above=0.3)

                        # コーパスを作成
                        corpus = [dictionary.doc2bow(text) for text in texts]
                        corpus_test = [dictionary.doc2bow(text) for text in texts_test]

                        # LDA の計算
                        topic_N = 1 + iterations
                        parameters = topic_N
                        lda = gensim.models.ldamodel.LdaModel(
                            corpus=corpus,
                            alpha='auto',
                            num_topics=topic_N,
                            id2word=dictionary
                        )

                        
                        cm = models.coherencemodel.CoherenceModel(model=lda, corpus=corpus, coherence='u_mass')  # tm is the trained topic model                        
                        d = {'y':cm.get_coherence(),'x':parameters}
                        df = pd.DataFrame(data=d,index=[0])
                        output_list = output_list.append(df)
                        
                        axe.plot(output_list.x,output_list.y)
                        fig.set_size_inches(8, 8)
                        display(fig)
                        axe.cla()
                        
                        for i in range(topic_N):
                            print('TOPIC:', i, '__', lda.print_topic(i))

output_list = pd.DataFrame(columns=["y","x"])

plt.ion()

fig = plt.figure()

axe = fig.add_subplot(111)

for iterations in range(2, 100):

clear_output(wait = True)

nakamoto_test = nakamoto_corpus.sample(frac=0.1, replace=True)

nakamoto_train= nakamoto_corpus[~nakamoto_corpus.index.isin(nakamoto_test.index)]

texts = [ ]

for line in nakamoto_train.text_wakati:

texts.append(line.split())

texts_test = [ ]

for line in nakamoto_test.text_wakati:

texts_test.append(line.split())

# 辞書作成

dictionary = corpora.Dictionary(texts)

dictionary.filter_extremes(no_below=20, no_above=0.3)

# コーパスを作成

corpus = [dictionary.doc2bow(text) for text in texts]

corpus_test = [dictionary.doc2bow(text) for text in texts_test]

# LDA の計算

topic_N = 1 + iterations

parameters = topic_N

lda = gensim.models.ldamodel.LdaModel(

corpus=corpus,

alpha='auto',

num_topics=topic_N,

id2word=dictionary

)

cm = models.coherencemodel.CoherenceModel(model=lda, corpus=corpus, coherence='u_mass') # tm is the trained topic model

d = {'y':cm.get_coherence(),'x':parameters}

df = pd.DataFrame(data=d,index=[0])

output_list = output_list.append(df)

axe.plot(output_list.x,output_list.y)

fig.set_size_inches(8, 8)

display(fig)

axe.cla()

for i in range(topic_N):

print('TOPIC:', i, '__', lda.print_topic(i))

実際に推定してみたところ、トピック数が20を超えたあたりからコヒーレンスが下がる傾向があるので、
それ以上のトピック数は追い求めない方が良いのかもしれません。

Rでもやってみる

Rでトピック数を決める良い方法がないか調べてみたところ、ldatuningとかいうパッケージがあることがわかりました。複数の論文（Griffiths2004, CaoJuan2009, Arun2010,Deveaud2014）で扱われている手法を元に、適切なトピック数を探れるようです。このパッケージを紹介しているブログの事例では、90から140の範囲で最適なトピック数となることが示されています。詳しくはこちらを見てください。
Select number of topics for LDA model

以下のコードで実行しました。一部、驚異のアニヲタさんのコードを拝借しております。なお、ldaパッケージのlexicalize関数を用いることで、ldatuningに入力するデータを作成することができます。

library(ldatuning)
library(topicmodels)
library(tidyverse)

nakamoto_data <- read_csv(file = "nakamoto_dataset.csv")

TFIDF <- function(corpus, progress=FALSE){ 
  res <- matrix(0, nr=length(corpus$vocab), nc=4)
  dimnames(res) <- list(corpus$vocab, c("documents", "count", "freq", "score"))
  res[, "documents"] <- length(corpus$documents)
  wordset <- mapply(function(x) x[1,], corpus$documents) 
  allfreq <- matrix(unlist(corpus$documents), nr=2)
  wordfreq <- tapply(allfreq[2,], allfreq[1,], sum) 
  for(v in seq(corpus$vocab)){ 
    count_docs <- sum(sapply(lapply(wordset, "==", v-1), any)) #
    res[v, "freq"] <- count_docs
    if(progress){ 
      pb <- txtProgressBar(min=1, max=length(corpus$vocab), style=3)
      setTxtProgressBar(pb, v)
    }
  }
  res[, "count"] <- wordfreq
  res[, "score"] <- log(res[, "count"]) * log(res[, "documents"]/res[, "freq"])
  return(as.data.frame(res))
}

lex1 <- lda::lexicalize(nakamoto_data$text_wakati_remove_freq)

s0 <- TFIDF(lex1, TRUE)
term1 <- rownames(s0)[s0$score > 0]
lex2 <- list(documents=lda::lexicalize(nakamoto_data$text_wakati_remove_freq, vocab=term1), vocab=term1)
dtm2 <- ldaformat2dtm(lex2$documents, lex2$vocab)

result <- FindTopicsNumber(
  dtm2,
  topics = seq(from = 2, to = 100, by = 1),
  metrics = c("Griffiths2004", "CaoJuan2009", "Arun2010", "Deveaud2014"),
  method = "Gibbs",
  control = list(seed = 77),
  mc.cores = 2L,
  verbose = TRUE
)

FindTopicsNumber_plot(values = result)

library(ldatuning)

library(topicmodels)

library(tidyverse)

nakamoto_data <- read_csv(file = "nakamoto_dataset.csv")

TFIDF <- function(corpus, progress=FALSE){

res <- matrix(0, nr=length(corpus$vocab), nc=4)

dimnames(res) <- list(corpus$vocab, c("documents", "count", "freq", "score"))

res[, "documents"] <- length(corpus$documents)

wordset <- mapply(function(x) x[1,], corpus$documents)

allfreq <- matrix(unlist(corpus$documents), nr=2)

wordfreq <- tapply(allfreq[2,], allfreq[1,], sum)

for(v in seq(corpus$vocab)){

count_docs <- sum(sapply(lapply(wordset, "==", v-1), any)) #

res[v, "freq"] <- count_docs

if(progress){

pb <- txtProgressBar(min=1, max=length(corpus$vocab), style=3)

setTxtProgressBar(pb, v)

}

res[, "count"] <- wordfreq

res[, "score"] <- log(res[, "count"]) * log(res[, "documents"]/res[, "freq"])

return(as.data.frame(res))

}

lex1 <- lda::lexicalize(nakamoto_data$text_wakati_remove_freq)

s0 <- TFIDF(lex1, TRUE)

term1 <- rownames(s0)[s0$score > 0]

lex2 <- list(documents=lda::lexicalize(nakamoto_data$text_wakati_remove_freq, vocab=term1), vocab=term1)

dtm2 <- ldaformat2dtm(lex2$documents, lex2$vocab)

result <- FindTopicsNumber(

dtm2,

topics = seq(from = 2, to = 100, by = 1),

metrics = c("Griffiths2004", "CaoJuan2009", "Arun2010", "Deveaud2014"),

method = "Gibbs",

control = list(seed = 77),

mc.cores = 2L,

verbose = TRUE

)

FindTopicsNumber_plot(values = result)

これを見る限りは、60〜70個の辺りに落ち着くのでしょうか。

トピックの吐き出し

Rでの結果から、60個程度のトピックで推定し、各記事に割り当てが最大のトピックを付与して、トピック別の口コミ評価をみてみようと思います。

以下のコードではトピック別の口コミ評価のしやすさからtopicmodelsパッケージを用いた推定となっています。

nbo_topics <- 60
lda_estimate <- topicmodels::LDA(dtm2,control=list(verbose=1, alpha = 0.1),
                                 k = nbo_topics,
                                 method = "Gibbs")

#トピックの上位5単語を確認する
terms_each_topics <- data.frame(terms(lda_estimate,5))
topic_keywords <- data.frame(topic_keywords_5 = apply(t(terms_each_topics),1,paste,collapse=","))
topic_keywords <- topic_keywords %>% mutate(topic_id=1:n())

#割り振られた最大の確率のトピックを抽出し、口コミデータと統合する
topics_each_document <- data.frame(topic_id=topics(lda_estimate,1))
topics_each_document <- nakamoto_data %>% bind_cols(topics_each_document)

#トピックごとの口コミ評価を計算する
topic_rating_summary <- topics_each_document %>% 
                            group_by(topic_id) %>% 
                            summarise(average_rating = mean(rating),
                                      count=n()) %>% 
                            left_join(topic_keywords, by = "topic_id") %>% 
                            arrange(desc(average_rating))

nbo_topics <- 60

lda_estimate <- topicmodels::LDA(dtm2,control=list(verbose=1, alpha = 0.1),

k = nbo_topics,

method = "Gibbs")

#トピックの上位5単語を確認する

terms_each_topics <- data.frame(terms(lda_estimate,5))

topic_keywords <- data.frame(topic_keywords_5 = apply(t(terms_each_topics),1,paste,collapse=","))

topic_keywords <- topic_keywords %>% mutate(topic_id=1:n())

#割り振られた最大の確率のトピックを抽出し、口コミデータと統合する

topics_each_document <- data.frame(topic_id=topics(lda_estimate,1))

topics_each_document <- nakamoto_data %>% bind_cols(topics_each_document)

#トピックごとの口コミ評価を計算する

topic_rating_summary <- topics_each_document %>%

group_by(topic_id) %>%

summarise(average_rating = mean(rating),

count=n()) %>%

left_join(topic_keywords, by = "topic_id") %>%

arrange(desc(average_rating))

口コミ評価の点数が上位のトピックはこんな感じです。

口コミ評価の点数が下位のトピックはこんな感じです。

中本は社会人2〜3年目で新規メディアの立ち上げのストレス解消で数回行きましたが、北極の赤さは異常だと思います。北極を食べたり、トッピングする余裕のある人、ましてや辛さを倍にするという時点で口コミ評価も高くなると考えるのは自然なのかもしれません。

参考情報

トピックモデル (機械学習プロフェッショナルシリーズ)
トピックモデルによる統計的潜在意味解析 (自然言語処理シリーズ)
models.ldamodel – Latent Dirichlet Allocation
Ldaのモデル選択におけるperplexityの評価
 pythonでgensimを使ってトピックモデル(LDA)を行う
 gensim0.8.6のチュートリアルをやってみた【コーパスとベクトル空間】
LDA 実装の比較
 Jupyter notebookにMatplotlibでリアルタイムにチャートを書く
 Inferring the number of topics for gensim’s LDA – perplexity, CM, AIC, and BIC
Select number of topics for LDA model
47の心得シリーズをトピックモデルで分類する。 – 驚異のアニヲタ社会復帰への道

Word2Vecを用いて蒙古タンメン中本の口コミ評価を予測してみる

はじめに

word2vecを用いた分類は以前からやってみたいと思っていたのですが、関心を持てるテキストデータがなかったのでなかなか手を出していませんでした。
ある時、ふとしたことから某グルメ系口コミサイトから蒙古タンメン中本の口コミと評価点を抽出して、その評価をword2vecでやってみるのは面白いだろうと思いついたので、さっそくやってみます。
こういう時にはじめて、データ分析だけでなくクローリング屋としても業務をやっていて良かったなと思うところですね。
コードは以前見つけて紹介した「分散表現を特徴量として文書分類するための方法について調べてみた」のものを再利用します。

目的

某グルメ系口コミサイトの口コミを収集し、個々人の口コミの内容から個々人の店に対する評価が高いか低いかを予測する。

データ収集

BeautifulSoupで収集しており、各店舗あわせて数千件ほど集めました。（実行コードはこちらでは紹介しません。）

このようなデータが手に入っている前提で以下の分析を進めていきます。

形態素解析

文書を形態素解析して、名詞のみを抽出するためのコードを用意します。

import MeCab
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

import MeCab

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

先ほどのデータフレームに対して以下のように実行すれば、名詞のみの分かち書きを行ったカラムが手に入ります。

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

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

集計

点数のヒストグラム

3.5点から4点の間が最も評価が多いようです。1点台をつける人はほとんどいないことがわかります。

単語数のヒストグラム

大体の口コミで100単語未満のようです。

単語数と点数の散布図

どうやら口コミにおいて500語を超える記述をしている人は評価が3点を下回ることはないようですが、文字数と点数でキレイに傾向が出ているわけではないですね。

形態素解析結果の集計、単語ランキング

名詞の抽出に関して非常に便利なMeCab Neologdを用いています。蒙古タンメンもきちんと捉えることができています。

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

noun_list = [] # 重複を含めた名詞のリスト
for i in corpus_data.text:
    for l in tagger.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) # 各名詞の数え上げ

word_count_df = pd.DataFrame.from_dict(noun_cnt, orient='index').reset_index()
word_count_df.columns = ["word", "count"]
word_count_df = word_count_df[~word_count_df.word.str.contains("^[あ-ん]$")]
word_count_df.sort_values("count",ascending=False).head(10)

import collections

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

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

for i in corpus_data.text:

for l in tagger.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) # 各名詞の数え上げ

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

word_count_df.columns = ["word", "count"]

word_count_df = word_count_df[~word_count_df.word.str.contains("^[あ-ん]$")]

word_count_df.sort_values("count",ascending=False).head(10)

味噌よりも北極の方が出現しているようです。北極は言わずもがな、極端に辛い罰ゲームレベルの一品。味噌タンメンは辛さが抑えめのラーメンで、知人の間では最もおいしいのがこのレベルだという合意があったりしますね。

分散表現とは

単語の意味を低次元の密な実数値ベクトルで表現したもの。
入力層から中間層への重み自体が各単語の分散表現となっている。
2017年9月のテキストアナリティクスシンポジウムにてメルカリとGunosyが特徴量として分散表現を活用しており性能が出ているとの発言があった。

word2vecについて

単語の分散表現を作ることを目的としている。

CBOW(Continuous Bag-of-Words)
注目している単語の前後N単語を文脈と呼び、その文脈をBag-of-Words表現として入力し、注目している単語を出力するというニューラルネットワークを学習する。入力層から隠れ層への結合は単語の位置を問わず同じとし、隠れ層の活性化関数をただの恒等関数としている。
Skip-gram
文脈のBOWを突っ込むCBOWとは異なり、入力層に1単語だけを入れる。1単語を入力し、正解データとして他の単語を入れることを繰り返して学習し、ある単語の入力に対して、どの単語の出現確率が高いかどうかを計算する。正解確率が上がるようにニューラルネットワークの重みを調整する。深層学習で使われる自己符号化器と似たような構造とされている。

gensimのword2vecの引数

gensimのword2vecには数多くの引数が存在します。gensimのドキュメントに英語で書かれていますが、せっかくなのでこちらで紹介します。

class gensim.models.word2vec.Word2Vec(sentences=None,
　size=100,
　alpha=0.025,
 window=5,
 min_count=5,
 max_vocab_size=None,
 sample=0.001, 
 seed=1,
 workers=3,
 min_alpha=0.0001,
 sg=0,
 hs=0,
 negative=5,
 cbow_mean=1,
 hashfxn=<built-in function hash>,
 iter=5,
 null_word=0,
 trim_rule=None,
 sorted_vocab=1,
 batch_words=10000,
 compute_loss=False,
 callbacks=())

class gensim.models.word2vec.Word2Vec(sentences=None,

　size=100,

　alpha=0.025,

window=5,

min_count=5,

max_vocab_size=None,

sample=0.001,

seed=1,

workers=3,

min_alpha=0.0001,

sg=0,

hs=0,

negative=5,

cbow_mean=1,

hashfxn=<built-in function hash>,

iter=5,

null_word=0,

trim_rule=None,

sorted_vocab=1,

batch_words=10000,

compute_loss=False,

callbacks=())

sentences
解析に使う1行1センテンスで書かれた文書。日本語の場合はLineSentenceフォーマットを使えばうまくいった。単語が空白文字で区切られていて、文章は改行で区切られていれば問題ない。
sg
{1,0}の整数で訓練アルゴリズムを設定できる。 1を選べばskip-gramで、0ならばCBOWを使う。
size
特徴ベクトルの次元を設定する。
window
文書内における現在の単語と予測した単語の間の距離の最大値を設定する。言い換えると、文脈の最大単語数を設定する。
alpha
学習率の初期値を設定する。
min_alpha
訓練の過程で徐々に落ちていく学習率の最小値を設定する。
seed
乱数を生成する際のシード番号を設定する。
min_count
一定の頻度以下の単語を除外する際の値を設定する。
max_vocab_size
語彙ベクトルを構築している際のメモリ制限を設定する。
sample
(0, 1e-5)の範囲で、頻度語がランダムに削除される閾値を設定する。高速化と精度向上を狙っており、自然言語処理においても高頻度語はストップワードとして除去するなどの対応が取られている。
workers
モデルを訓練するために多くのワーカースレッドを利用するかどうか設定する。（並列化関連）
hs
{1,0}の整数で、1であれば階層的ソフトマックスがモデルの訓練で用いられ、0であり引数negativeがnon-zeroであればネガティヴサンプリングが設定できる。全部計算することが大変なので、階層的なグループに分けて各グループごとに学習するというのがモチベーション。
negative
0よりも大きければネガティブサンプリングが用いられる。5〜20などを選び、どれだけノイズワードが描かれているかを識別する。0であればネガティブサンプリングが適用されない。ネガティブサンプリングは計算高速化を目的に出力層で正解ニューロン以外のニューロンを更新しないように学習する手法。
cbow_mean
{1,0}の整数で、0であれば単語ベクトルの合計を用い、1であればCBOWが用いられた際の平均が用いられる。
hashfxn
訓練の再現性のためにランダムに初期値のウエイト付けできる。
iter
コーパスにおける繰り返し回数（エポック数）を設定できる。
trim_rule
ある単語を語彙に含めるべきかどうかを識別する、語彙のトリミングルールを設定する。
sorted_vocab
{1,0}の整数で、1であれば頻度の降順で語彙を並べ替える。
batch_words
ワーカースレッドにわたすバッチの大きさを指定する。（並列化関連）
compute_loss
Trueであれば損失関数の計算と蓄積を行う。
callbacks
訓練時の特定の段階で実行する際に必要なコールバックのリストを指定できる。

word2vecによる文書分類の適用

口コミの点数が4点以上であれば1、そうでなければ0を取る変数を作成し、それをラベルとして文書分類を行います。
以前、紹介したブログ同様に、scikit-learnのExtraTreesClassifierを用いてCountVectorizerとTfidfVectorizerを特徴量としたものをベースラインとして、同様の手法に対してword2vecで作成した分散表現を特徴量として用いたものとを比較します。評価指標はクロスバリデーションスコア(5-folds)とします。

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
from sklearn.cross_validation import StratifiedShuffleSplit

#単語の分散表現の平均値を求めるクラスの定義
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
            ])

#ベースラインとなる既存手法のモデルの準備
#etree
etree = Pipeline([("count_vectorizer",CountVectorizer(analyzer=lambda x: x)), 
                    ("linear svc", SVC(kernel="linear"))])
#etreeのTF-IDF版
etree_tfidf = Pipeline([("tfidf_vectorizer", TfidfVectorizer(analyzer=lambda x: x)),
                          ("extra trees", ExtraTreesClassifier(n_estimators=200))])

#Word2Vecを特徴量としてExtraTreesによる分類器を準備する。
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))])

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

from sklearn.cross_validation import StratifiedShuffleSplit

#単語の分散表現の平均値を求めるクラスの定義

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

])

#ベースラインとなる既存手法のモデルの準備

#etree

etree = Pipeline([("count_vectorizer",CountVectorizer(analyzer=lambda x: x)),

("linear svc", SVC(kernel="linear"))])

#etreeのTF-IDF版

etree_tfidf = Pipeline([("tfidf_vectorizer", TfidfVectorizer(analyzer=lambda x: x)),

("extra trees", ExtraTreesClassifier(n_estimators=200))])

#Word2Vecを特徴量としてExtraTreesによる分類器を準備する。

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))])

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

#scikit-learnで使う入力とラベルの設定
X, y = np.array(corpus_data.text_wakati), np.array(corpus_data.label)

sentences2 = [token.split(" ") for token in corpus_data.text_wakati]


#Word2Vecを実行する。 
model = Word2Vec(sentences2, 
                 size=100,
                 window=5, 
                 min_count=5,
                 workers=2)

w2v = {w: vec for w, vec in zip(model.index2word, model.syn0)}

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

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

#scikit-learnで使う入力とラベルの設定

X, y = np.array(corpus_data.text_wakati), np.array(corpus_data.label)

sentences2 = [token.split(" ") for token in corpus_data.text_wakati]

#Word2Vecを実行する。

model = Word2Vec(sentences2,

size=100,

window=5,

min_count=5,

workers=2)

w2v = {w: vec for w, vec in zip(model.index2word, model.syn0)}

分類の前に、せっかくword2vecを使ったので、任意の単語に類似した単語を見てみます。

まずは初心者向けの味噌ラーメン

続いて、中級者向けの蒙古タンメン

そして、上級者向けの北極ラーメン

最後に、誰もが経験する翌日という単語。

どれも関連性の高いと思われる単語が抽出できているように思われます。

それでは分類モデルの学習を以下のコードで行います。
scikit-learnを使えば、データさえあれば非常に短いコードで書けてしまいます。

#各モデルを実行し、クロスバリデーションスコアを計算し、出力させる。
all_models = [
    ("etree",etree),
    ("etree_tfidf",etree_tfidf),
    ("w2v", etree_w2v),
    ("w2v_tfidf", etree_w2v_tfidf)
]
scores = sorted([(name, cross_val_score(model, X, y, cv=5).mean()) 
                 for name, model in all_models], 
                key = lambda x:x[0])
print(tabulate(scores, floatfmt=".4f", headers=("model", 'score')))
 
#棒グラフでクロスバリデーションスコアを比較する。
plt.figure(figsize=(15, 6))
sns.barplot(x=[name for name, _ in scores], y=[score for _, score in scores])

#各モデルを実行し、クロスバリデーションスコアを計算し、出力させる。

all_models = [

("etree",etree),

("etree_tfidf",etree_tfidf),

("w2v", etree_w2v),

("w2v_tfidf", etree_w2v_tfidf)

]

scores = sorted([(name, cross_val_score(model, X, y, cv=5).mean())

for name, model in all_models],

key = lambda x:x[0])

print(tabulate(scores, floatfmt=".4f", headers=("model", 'score')))

#棒グラフでクロスバリデーションスコアを比較する。

plt.figure(figsize=(15, 6))

sns.barplot(x=[name for name, _ in scores], y=[score for _, score in scores])

model          score
-----------  -------
etree         0.6959
etree_tfidf   0.8031
w2v           0.8046
w2v_tfidf     0.8038

model score

----------- -------

etree 0.6959

etree_tfidf 0.8031

w2v 0.8046

w2v_tfidf 0.8038

一応、ベースラインよりもword2vecを特徴量としたものの方がスコアが高いのですが、わずかです。TF-IDFベースで特徴量を作成したモデルは十分に性能が出ているようです。
word2vecを用いることによる旨味はそれほどなさそうですが、パラメータを試行錯誤していけばよくなるかもしれません。

終わりに

蒙古タンメン中本のテキストをWebスクレイピングし、その口コミ情報をコーパスとして口コミ評価の二値分類に挑戦しましたが、TF-IDFよりもわずかに優秀な特徴量になりうるという結果になりました。もっと劇的な向上を夢見ていたのですが、パラメータの試行錯誤を今後の宿題としようと思います。

参考情報

Chainer v2による実践深層学習
 word2vecによる自然言語処理
 models.word2vec – Deep learning with word2vec
Python で「老人と海」を word2vec する
 Python3 – MeCabで日本語文字列の名詞出現数の出し方について
 Transform a Counter object into a Pandas DataFrame

LIMEで赤ワインのデータをいじってみる with Python

はじめに

2018年1月のTokyoR（ TokyoR67に行ってきました）で機械学習結果の解釈可能性について多く語られていたので、Hello World的な試行を赤ワインのデータで行ってみたいと思います。コードは論文を書いた人のGitHubのものを拝借しました。

営業やマーケティングのメンバーに機械学習手法を提案する際に、ひとつひとつの予測結果において、なんでその予測結果になったのか理由が知りたいという要望を受けることが多いです。
ある講演で、実装用のモデルにランダムフォレストやSVMを使い、マネジャーに説明する用に決定木の結果を見せたとかいう話もありました。わざわざモデルを二つ作って説明のための工数を取らなくてよくなると思うと非常にありがたい技術です。事業側の解釈可能性を重視するあまり、シンプルな手法を取らざるを得ない現場には朗報ですね。

目次
・LIMEとは
・データ
・コード
・結果の解釈
・参考情報

LIMEとは

Local Interpretable Model-agnostic Explanationsの頭文字をとったもので、機械学習によって構築したモデルに関して、その予測結果を人間が解釈しやすくする技術です。
流れとしては、

まずランダムフォレストなりXGBoostなりで分類器を作る。
任意のデータxを取り出し、解釈可能バージョンのx’（x’∈{0,1}で、xの非ゼロ要素を1としている。）を用意する。
x’の周辺のデータをサンプリングする。
サンプリングしたデータを使って、元のモデルを近似するために、モデルの距離を目的関数とした最適化問題を解く。
　（K-LASSOという線形モデルの手法を使うことで、最終的に近似するモデルの変数の数を決めている。）
学習したモデルの偏回帰係数を確認して、予測結果への影響度を見る。

という流れのようです。間違っているかもなので、参考情報をご覧になってください。
紹介動画も作られているようです。そういえば、Stanも紹介動画ありましたね。

データ

ワインデータから赤ワインのデータのみを利用します。データの詳細はこちらにあります。
http://archive.ics.uci.edu/ml/machine-learning-databases/wine-quality/winequality.names

このようなデータセットでデータサイズは1599です。

ワインの質に関するヒストグラムはこんな感じです。

今回の分析の目的

赤ワインの評価値が7以上であれば良いワイン（Y=1）、そうでなければ悪いワイン（Y=0）として、評価値が7を超えるような赤ワインの分類器を学習させ、任意のテストデータを取り出し、そのテストデータが良いワインである要因を探るものとします。

進め方

・LIMEのインストール（pip install lime で一発）
・scikit-learnによる機械学習（ランダムフォレスト）
・LIMEを用いた予測結果の解釈の提示

コード

今回のコードはこちらにもあります。
（ kamonohashiperry.com/lime_study/Lime_With_Wine Data.ipynb ）

import sklearn
import sklearn.ensemble
import numpy as np
import lime
import lime.lime_tabular
import pandas as pd
from __future__ import print_function
np.random.seed(1)

# データのインポート
wine_dataset = pd.read_csv("https://archive.ics.uci.edu/ml/machine-learning-databases/wine-quality/winequality-red.csv", sep=";")


#ワインの質に関する数値
Y = wine_dataset.quality.values
#質に関するデータを落としている。
wine_dataset = wine_dataset.drop('quality', axis =1)
#7よりも小さかったら0、それ以外は1とする。
Y = np.asarray([1 if  i >= 7 else 0 for i in Y])
Y_label = np.asarray(['good' if  i>=1 else 'bad' for i in Y])

# 訓練データ、テストデータの作成
train, test, labels_train, labels_test = sklearn.model_selection.train_test_split(wine_dataset.values,
                                                                                  Y,
                                                                                  train_size=0.80)

# ランダムフォレストの学習
rf = sklearn.ensemble.RandomForestClassifier(n_estimators=500,
                                             class_weight="balanced" )
rf.fit(train, labels_train)


sklearn.metrics.accuracy_score(labels_test, rf.predict(test))

# LIMEの計算
explainer = lime.lime_tabular.LimeTabularExplainer(train,
                                                   feature_names = list(wine_dataset.columns.values),
                                                   class_names = np.unique(Y_label.astype('<U10')),
                                                   discretize_continuous=True)

import sklearn

import sklearn.ensemble

import numpy as np

import lime

import lime.lime_tabular

import pandas as pd

from __future__ import print_function

np.random.seed(1)

# データのインポート

wine_dataset = pd.read_csv("https://archive.ics.uci.edu/ml/machine-learning-databases/wine-quality/winequality-red.csv", sep=";")

#ワインの質に関する数値

Y = wine_dataset.quality.values

#質に関するデータを落としている。

wine_dataset = wine_dataset.drop('quality', axis =1)

#7よりも小さかったら0、それ以外は1とする。

Y = np.asarray([1 if i >= 7 else 0 for i in Y])

Y_label = np.asarray(['good' if i>=1 else 'bad' for i in Y])

# 訓練データ、テストデータの作成

train, test, labels_train, labels_test = sklearn.model_selection.train_test_split(wine_dataset.values,

train_size=0.80)

# ランダムフォレストの学習

rf = sklearn.ensemble.RandomForestClassifier(n_estimators=500,

class_weight="balanced" )

rf.fit(train, labels_train)

sklearn.metrics.accuracy_score(labels_test, rf.predict(test))

# LIMEの計算

explainer = lime.lime_tabular.LimeTabularExplainer(train,

feature_names = list(wine_dataset.columns.values),

class_names = np.unique(Y_label.astype('<U10')),

discretize_continuous=True)

結果を解釈する以前に、二値分類モデルとして精度が低かったら元も子もないので、精度やAUCを確認してみます。

# accuracyの計算
sklearn.metrics.accuracy_score(labels_test, rf.predict(test))
0.93125

# confusion_matrixの確認
sklearn.metrics.confusion_matrix(labels_test, rf.predict(test))
array([[284,   4],
       [ 18,  14]])

# AUCの計算
fpr, tpr, thresholds = sklearn.metrics.roc_curve(labels_test, rf.predict(test), pos_label=1)
sklearn.metrics.auc(fpr, tpr)
0.7118055555555556

# accuracyの計算

sklearn.metrics.accuracy_score(labels_test, rf.predict(test))

0.93125

# confusion_matrixの確認

sklearn.metrics.confusion_matrix(labels_test, rf.predict(test))

array([[284, 4],

[ 18, 14]])

# AUCの計算

fpr, tpr, thresholds = sklearn.metrics.roc_curve(labels_test, rf.predict(test), pos_label=1)

sklearn.metrics.auc(fpr, tpr)

0.7118055555555556

scikit-learnの引数でclass_weight=”balanced”にしているので、少しは不均衡データに対応できているようです。AUCは8割を超えたかったですが、いったんこれで進めます。

結果の解釈

こちらのコードで、ランダムにインスタンスを選んで、その予測結果とその予測結果に影響を与えている変数を見てみます。

i = np.random.randint(0, test.shape[0])
exp = explainer.explain_instance(test[i], rf.predict_proba, num_features=2, top_labels=1)

exp.show_in_notebook(show_table=True, show_all=False)

i = np.random.randint(0, test.shape[0])

exp = explainer.explain_instance(test[i], rf.predict_proba, num_features=2, top_labels=1)

exp.show_in_notebook(show_table=True, show_all=False)

どうやら、このインスタンスを良いワインと予測しており、sulphates（硫酸）が0.74を超えていた、volatile acidity（酢酸とその派生物質）が0.39よりも小さいというのが理由のようです。画像では結果が消えていますが、exp.as_list()で結果を抽出できます。

著者のコードには続きがあって、このインスタンスに対して、値を足したりすることで分類確率がどのように変わるのかが記されていました。今回はアルコールを3%ポイント増やして、総亜硫酸濃度を30増やしてみるものとします。

feature_index = lambda x: list(wine_dataset.columns.values).index(x)

print('Increasing alcohol')
temp = test[i].copy()
print('P(Bad) before:', rf.predict_proba(temp.reshape(1,-1))[0,0])
temp[feature_index('alcohol')] = test[i][feature_index('alcohol')] + 3
print('P(Bad) after:', rf.predict_proba(temp.reshape(1,-1))[0,0])
print ()
print('Increasing total sulfur dioxide')
temp = test[i].copy()
print('P(Bad) before:', rf.predict_proba(temp.reshape(1,-1))[0,0])
temp[feature_index('total sulfur dioxide')] = test[i][feature_index('total sulfur dioxide')] + 30   
print('P(Bad) after:', rf.predict_proba(temp.reshape(1,-1))[0,0])
print()
print('Increasing both')
temp = test[i].copy()
print('P(Bad) before:', rf.predict_proba(temp.reshape(1,-1))[0,0])
temp[feature_index('alcohol')] = test[i][feature_index('alcohol')] + 3
temp[feature_index('total sulfur dioxide')] = test[i][feature_index('total sulfur dioxide')] + 30
print('P(Bad) after:', rf.predict_proba(temp.reshape(1,-1))[0,0])

feature_index = lambda x: list(wine_dataset.columns.values).index(x)

print('Increasing alcohol')

temp = test[i].copy()

print('P(Bad) before:', rf.predict_proba(temp.reshape(1,-1))[0,0])

temp[feature_index('alcohol')] = test[i][feature_index('alcohol')] + 3

print('P(Bad) after:', rf.predict_proba(temp.reshape(1,-1))[0,0])

print ()

print('Increasing total sulfur dioxide')

temp = test[i].copy()

print('P(Bad) before:', rf.predict_proba(temp.reshape(1,-1))[0,0])

temp[feature_index('total sulfur dioxide')] = test[i][feature_index('total sulfur dioxide')] + 30

print('P(Bad) after:', rf.predict_proba(temp.reshape(1,-1))[0,0])

print()

print('Increasing both')

temp = test[i].copy()

print('P(Bad) before:', rf.predict_proba(temp.reshape(1,-1))[0,0])

temp[feature_index('alcohol')] = test[i][feature_index('alcohol')] + 3

temp[feature_index('total sulfur dioxide')] = test[i][feature_index('total sulfur dioxide')] + 30

print('P(Bad) after:', rf.predict_proba(temp.reshape(1,-1))[0,0])

これを見る限り、アルコールを3%ポイント増やすと、悪いワインの確率が6.2%ポイントあがり、総亜硫酸濃度を30増やすと、悪いワインの確率が49.8%ポイントあがることが示されています。
ワインあまり飲まないので結果の解釈以前に変数の解釈ができていないのは今回のオチになるんでしょうか。このコードをもとに、仕事現場で使ってみようと思います。きっと解釈できるはず。

参考情報

“Why Should I Trust You?” Explaining the Predictions of Any Classifier
lime/doc/notebooks/Tutorial – continuous and categorical features.ipynb
3.2.4.3.1. sklearn.ensemble.RandomForestClassifier
機械学習と解釈可能性 by Sinhrks
機械学習モデルの予測結果を説明するための力が欲しいか…？
LIMEで機械学習の予測結果を解釈してみる
 ワインの味（美味しさのグレード）は予測できるか？（１）

分散表現を特徴量として文書分類するための方法について調べてみた

以前、テキストアナリティクスシンポジウムに参加した際に登壇者が機械学習のタスクにおいて分散表現を特徴量に使ったと言っていて、実務で使えるようにしたいと思ったので、調べた手法について記します。

先行研究

海外のブログ(Text Classification With Word2Vec)では、以下の手順で文書分類のための特徴量としてWord2Vecが用いられていました。

STEP1:ロイターのコーパスをもとにWord2Vecを求める。（先行研究では100次元に圧縮しています。）
STEP2:テキストごとに属する単語に対して、STEP1で求めた100次元の分散表現の平均値をとる。あるいは、TF-IDFのスコアで重み付けしたものを作る。これで100次元の分散表現が各テキストごとに用意できる。
STEP3：STEP2で作った100次元の分散表現を用いて、テキストごとのラベルについてExtra Treesによる分類器を学習させる。

ざっとこんな感じのアプローチになっていました。”ベクトルの足し算が意味の足し算に対応する「加法構成性」”という考え方からのアプローチと言えるのでしょうか。

そのあと、ベースラインとの比較をしていましたが、結論としては、ラベル付きのトレーニングデータが非常に少ない場合において有利である傾向があったものの、SVMなどの既存手法より精度が出ているとは言えないような結果でした。

今回の事例におけるタスクではうまくいかなかったようですが、テキストアナリティクスシンポジウムにおいてはアツいとされていたので、自社のデータで適用するなどしておきたいです。

今回取り上げたブログやGithubにはPython2.X系のコードしかなかったため、3.X系で動くように一部変更するなどをしました。
以下で3系で動き、同様の結果が出ているものを載せています。

コード

コードの概要としては
・訓練データの作成
・モデルの定義
・モデルのベンチマーク
・結果の描写
です。

事前にテキストデータをダウンロードして手に入れておく必要がありますが、それに関してはターミナルで実行できるものなので、ここでは載せていません。参考文献のblog_stuff/classification_w2v/benchmarking.ipynbを参照ください。

#ライブラリの読み込み
from tabulate import tabulate
%matplotlib inline
import matplotlib.pyplot as plt
import seaborn as sns
import pandas as pd
import numpy as np
from gensim.models.word2vec import Word2Vec
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.naive_bayes import BernoulliNB, MultinomialNB
from sklearn.pipeline import Pipeline
from sklearn.svm import SVC
from sklearn.metrics import accuracy_score
from sklearn.cross_validation import cross_val_score
from sklearn.cross_validation import StratifiedShuffleSplit

TRAIN_SET_PATH = "r8-no-stop.txt"

#行ごとにテキストをTabで分割し、ラベルデータとテキストデータに分ける。
X, y = [], []
with open(TRAIN_SET_PATH, "rb") as infile:
    for line in infile:
        label, text = line.decode('utf-8').split("\t")
        X.append(text.split())
        y.append(label)
X, y = np.array(X), np.array(y)
print("total examples %s" % len(y))

#Word2Vecを実行する。
model = Word2Vec(X, size=100, window=5, min_count=5, workers=2)
w2v = {w: vec for w, vec in zip(model.index2word, model.syn0)}

#ベースラインとなる既存手法のモデルの準備
mult_nb = Pipeline([("count_vectorizer", CountVectorizer(analyzer=lambda x: x)), ("multinomial nb", MultinomialNB())])
bern_nb = Pipeline([("count_vectorizer", CountVectorizer(analyzer=lambda x: x)), ("bernoulli nb", BernoulliNB())])
mult_nb_tfidf = Pipeline([("tfidf_vectorizer", TfidfVectorizer(analyzer=lambda x: x)), ("multinomial nb", MultinomialNB())])
bern_nb_tfidf = Pipeline([("tfidf_vectorizer", TfidfVectorizer(analyzer=lambda x: x)), ("bernoulli nb", BernoulliNB())])
svc = Pipeline([("count_vectorizer", CountVectorizer(analyzer=lambda x: x)), ("linear svc", SVC(kernel="linear"))])
svc_tfidf = Pipeline([("tfidf_vectorizer", TfidfVectorizer(analyzer=lambda x: x)), ("linear svc", SVC(kernel="linear"))])

#単語の分散表現の平均値を求めるクラスの定義
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
            ])

#Word2Vecを特徴量としてExtraTreesによる分類器を準備する。
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 = [
    ("mult_nb", mult_nb),
    ("mult_nb_tfidf", mult_nb_tfidf),
    ("bern_nb", bern_nb),
    ("bern_nb_tfidf", bern_nb_tfidf),
    ("svc", svc),
    ("svc_tfidf", svc_tfidf),
    ("w2v", etree_w2v),
    ("w2v_tfidf", etree_w2v_tfidf)
]
scores = sorted([(name, cross_val_score(model, X, y, cv=5).mean()) 
                 for name, model in all_models], 
                key = lambda x:x[0])
print(tabulate(scores, floatfmt=".4f", headers=("model", 'score')))


#棒グラフでクロスバリデーションスコアを比較する。
plt.figure(figsize=(15, 6))
sns.barplot(x=[name for name, _ in scores], y=[score for _, score in scores])

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

#ライブラリの読み込み

from tabulate import tabulate

%matplotlib inline

import matplotlib.pyplot as plt

import seaborn as sns

import pandas as pd

import numpy as np

from gensim.models.word2vec import Word2Vec

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.naive_bayes import BernoulliNB, MultinomialNB

from sklearn.pipeline import Pipeline

from sklearn.svm import SVC

from sklearn.metrics import accuracy_score

from sklearn.cross_validation import cross_val_score

from sklearn.cross_validation import StratifiedShuffleSplit

TRAIN_SET_PATH = "r8-no-stop.txt"

#行ごとにテキストをTabで分割し、ラベルデータとテキストデータに分ける。

X, y = [], []

with open(TRAIN_SET_PATH, "rb") as infile:

for line in infile:

label, text = line.decode('utf-8').split("\t")

X.append(text.split())

y.append(label)

X, y = np.array(X), np.array(y)

print("total examples %s" % len(y))

#Word2Vecを実行する。

model = Word2Vec(X, size=100, window=5, min_count=5, workers=2)

w2v = {w: vec for w, vec in zip(model.index2word, model.syn0)}

#ベースラインとなる既存手法のモデルの準備

mult_nb = Pipeline([("count_vectorizer", CountVectorizer(analyzer=lambda x: x)), ("multinomial nb", MultinomialNB())])

bern_nb = Pipeline([("count_vectorizer", CountVectorizer(analyzer=lambda x: x)), ("bernoulli nb", BernoulliNB())])

mult_nb_tfidf = Pipeline([("tfidf_vectorizer", TfidfVectorizer(analyzer=lambda x: x)), ("multinomial nb", MultinomialNB())])

bern_nb_tfidf = Pipeline([("tfidf_vectorizer", TfidfVectorizer(analyzer=lambda x: x)), ("bernoulli nb", BernoulliNB())])

svc = Pipeline([("count_vectorizer", CountVectorizer(analyzer=lambda x: x)), ("linear svc", SVC(kernel="linear"))])

svc_tfidf = Pipeline([("tfidf_vectorizer", TfidfVectorizer(analyzer=lambda x: x)), ("linear svc", SVC(kernel="linear"))])

#単語の分散表現の平均値を求めるクラスの定義

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

])

#Word2Vecを特徴量としてExtraTreesによる分類器を準備する。

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 = [

("mult_nb", mult_nb),

("mult_nb_tfidf", mult_nb_tfidf),

("bern_nb", bern_nb),

("bern_nb_tfidf", bern_nb_tfidf),

("svc", svc),

("svc_tfidf", svc_tfidf),

("w2v", etree_w2v),

("w2v_tfidf", etree_w2v_tfidf)

]

scores = sorted([(name, cross_val_score(model, X, y, cv=5).mean())

for name, model in all_models],

key = lambda x:x[0])

print(tabulate(scores, floatfmt=".4f", headers=("model", 'score')))

#棒グラフでクロスバリデーションスコアを比較する。

plt.figure(figsize=(15, 6))

sns.barplot(x=[name for name, _ in scores], y=[score for _, score in scores])

このようにSVCに負けている結果が出ます。

model            score
-------------  -------
bern_nb         0.7954
bern_nb_tfidf   0.7954
mult_nb         0.9467
mult_nb_tfidf   0.8615
svc             0.9562
svc_tfidf       0.9656
w2v             0.9526
w2v_tfidf       0.9549

model score

------------- -------

bern_nb 0.7954

bern_nb_tfidf 0.7954

mult_nb 0.9467

mult_nb_tfidf 0.8615

svc 0.9562

svc_tfidf 0.9656

w2v 0.9526

w2v_tfidf 0.9549

比較の棒グラフはこんな感じ。

自社のテキストデータで是非とも試したい。

参考文献

岩波データサイエンス Vol.2
blog_stuff/classification_w2v/benchmarking.ipynb
分散表現(単語埋め込み)
Text Classification With Word2Vec

確率的プログラミングライブラリ「Edward」まとめ

Edwardで何ができるのか知らなかったので、忘備録として残しておきます。

目次
・Edwardとは
・Edwardでできること
・参考スライド
・参考文献

Edwardとは

・LDAで有名なコロンビア大のBlei先生の研究室で、2016年より開発されている確率的プログラミング((　プログラミング言語の変数をモデルの構成要素として使うプログラミング))のPythonライブラリ。
・積み木のように明快な形で確率的モデリングを行うことができる。（モデル→推論→評価　を一括でできる。）
・ベイズ統計と機械学習、深層学習、確率的プログラミングを融合させている。
・計算の際にTensorFlowを用いている。TensorBoardを可視化の際に用いることもできる。
・計算速度がStanやPyMC3よりも速い。GPUを用いた高速化も可能。((　pip install tensorflow-gpuでGPU版のTensolFlowを入れておく必要がある。))
・統計学者のGeorge Edward Pelham Boxから名前を取っている。

Edwardでできること

一般的なベイズ推定は当然ながら、深層学習向けのベイズ適用系の事例が豊富なようです。

・ベイズ線形回帰( Supervised Learning (Regression) )
・バッチトレーニング（巨大なデータセットにおける学習で用いる）( Batch Training )
・Tensorboardを用いた可視化（ Tensorboard ）
・Automated Transformations( Automated Transformations )
・線形混合効果モデル( Linear Mixed Effects Models )
・教師あり学習による分類( Supervised Learning (Classification) )
・教師なし学習( Unsupervised Learning )
・ニューラルネットワークの潜在空間モデル( Latent Space Models for Neural Data )
・混合密度ネットワーク( Mixture Density Networks )
・GAN( Generative Adversarial Networks )
・確率的デコーダー( Probabilistic Decoder )
・ネットワークの推論( Inference Networks )
・ベイジアンニューラルネットワーク( Bayesian Neural Network )
・確率的PCA(主成分分析)( Probabilistic PCA )

jupyterのコードたちはblei-lab/edwardのnotebookに載っています。

2層のニューラルネットワークへのベイズ推定の適用(
Bayesian Deep Learning with Edward (and a trick using Dropout) – Andrew Rowan – PyData London 2017)
Gounosyの方のブログによると、CTR予測などで扱うことができるようです。

参考スライド

確率的プログラミングライブラリEdward from Yuta Kashino

ベイズ統計の数理と深層学習 @db analytics showcase Sapporo 2017 from Takeshi Yoneda

（DL hacks輪読）Bayesian Neural Network from Masahiro Suzuki

参考文献

Edward: A library for probabilistic modeling, inference, and criticism

DEEP PROBABILISTIC PROGRAMMING

EdwardでBayesian DNN＋Variational Inferenceをやってみた話

DEEP PROBABILISTIC PROGRAMMING —”深層学習＋ベイズ”のライブラリ— Edwardの紹介

【Edward】MCMCの数学的基礎からStochastic Gradient Langevin Dynamicsの実装まで

Hello, world! Stan, PyMC3, and Edward
stanの開発者の方がstanとPyMC3とEdwardを比較しています。

Pythonで体験するベイズ推論 PyMCによるMCMC入門

洋楽の歌詞データでDoc2vecを実行してみる

はじめに

仕事で記事間の類似度計算などがあったりするんですが、TF-IDFにしてCOS類似度を計算するなどの方法で行っているのが現状です。そろそろ他の手法にも移行したいので、それに変わる類似度計算の手段としてDoc2vecを試してみたいと思います。

データ

以前より収集している洋楽の歌詞データを用います。Billboardのランキングに登場した楽曲の歌詞データを961曲分集めたものとなります。英語なので、日本語のように形態素解析は不要ですが、ストップワードを除去するなどの処理を施したコーパスを用います。Rのtmパッケージによるストップワードの除去についてはBillboard100位以内の楽曲の歌詞情報にLDAを適用してみたをご覧ください。日本語でのDoc2Vecの適用は参考文献において紹介しています。

類似度計算

TaggedLineDocumentを用いて、doc2vecで扱えるオブジェクトを作成します。TaggedLineDocumentに指定するファイルは主にtxtファイルで、その満たすべき条件は「1行につき1文書」「単語がスペースで区切られている」などです。あとは、doc2vecを実行するだけです。パラメータなどの細かい指定については今後の宿題としたいです。

#ライブラリの読み込み
from gensim.models.doc2vec import LabeledSentence
from collections import namedtuple
import pandas as pd

#データの読み込みとtxtファイルの生成
billboard_list = pd.read_csv("lyrics_data_preprocessed.csv", delimiter=",", encoding='utf-8')
billboard_list.columns = ['title','artist', 'lyrics']
billboard_songs = billboard_list['lyrics']
billboard_songs.to_csv("billboard_data.txt",header=None,index=None)

#1行1ドキュメントとしてdoc2vecで扱えるオブジェクトに変換
sentences = doc2vec.TaggedLineDocument("billboard_data.txt")

#doc2vecの実行
model = models.Doc2Vec(sentences, dm=0, size=300, window=15, alpha=.025,
        min_alpha=.025, min_count=1, sample=1e-6)

#トレーニングの開始
print('\nStart Training')
for epoch in range(20):
    print('Epoch: {}'.format(epoch + 1))
    model.train(sentences)
    model.alpha -= (0.025 - 0.0001) / 19
    model.min_alpha = model.alpha

#ライブラリの読み込み

from gensim.models.doc2vec import LabeledSentence

from collections import namedtuple

import pandas as pd

#データの読み込みとtxtファイルの生成

billboard_list = pd.read_csv("lyrics_data_preprocessed.csv", delimiter=",", encoding='utf-8')

billboard_list.columns = ['title','artist', 'lyrics']

billboard_songs = billboard_list['lyrics']

billboard_songs.to_csv("billboard_data.txt",header=None,index=None)

#1行1ドキュメントとしてdoc2vecで扱えるオブジェクトに変換

sentences = doc2vec.TaggedLineDocument("billboard_data.txt")

#doc2vecの実行

model = models.Doc2Vec(sentences, dm=0, size=300, window=15, alpha=.025,

min_alpha=.025, min_count=1, sample=1e-6)

#トレーニングの開始

print('\nStart Training')

for epoch in range(20):

print('Epoch: {}'.format(epoch + 1))

model.train(sentences)

model.alpha -= (0.025 - 0.0001) / 19

model.min_alpha = model.alpha

類似度の算出

早速、気になる楽曲に関して、類似度の高い楽曲を抽出してみたいと思います。

#推定したモデルの保存
model.save('doc2vec.model')
model = models.Doc2Vec.load('doc2vec.model')

#調べたい楽曲のインデックスの確認
billboard_list[billboard_list['title'].str.contains("Radioactive")]

#推定したモデルの保存

model.save('doc2vec.model')

model = models.Doc2Vec.load('doc2vec.model')

#調べたい楽曲のインデックスの確認

billboard_list[billboard_list['title'].str.contains("Radioactive")]

どうやら、Radioactiveという曲はkings of leonというグループも歌っているようですが、私はimagine dragonsの方の楽曲に関心がありますので、インデックスを409にして歌詞情報の近い楽曲を抽出します。

#最も近い歌詞の楽曲情報を抽出
mostsimilarlyrics = model.docvecs.most_similar(409)
billboard_list['title'][mostsimilarlyrics[0][0]]

'Made In America'

#最も近い歌詞の楽曲情報を抽出

mostsimilarlyrics = model.docvecs.most_similar(409)

billboard_list['title'][mostsimilarlyrics[0][0]]

'Made In America'

どうやら、toby keithのMade In Americaという楽曲が最も近いようです。類似度は35%程度ですが、全然単語が被っていないので本当に近いのか納得がいかないです。

次に、lady gagaのBorn This Wayに近い楽曲を出してみます。Dancing Queenという非常に懐かしい曲が選ばれていますが、類似度は49%と先ほどよりも高いです。queenやgirlやcanやrightなど共通の単語が含まれているので、先ほどの結果よりは近いのかなぁと思います。

正しく推定できているのか不安だったので、類似度が90%と非常に高かった、Just The Way You Areという楽曲の最も近い楽曲を見てみます。

調べたところ、同じ楽曲のカバー版のようです。近いものは、ちゃんと近いと見なせるようです。近いかどうかの基準をどの水準に置くのかは難しい判断ですね。

参考情報

models.doc2vec – Deep learning with paragraph2vec
Doc2Vecの仕組みとgensimを使った文書類似度算出チュートリアル
 Pythonによるデータ分析入門 ―NumPy、pandasを使ったデータ処理

Billboard100位以内の楽曲の歌詞情報にLDAを適用してみた

目次

・はじめに
・データ収集
・Rによる分析
・LDAの結果
・参考文献

はじめに

前回の投稿でBillboardの週次洋楽ランキングデータをWebスクレイピングで取得し、楽曲の消費サイクルのような順位の挙動を確かめることができました。（某洋楽ヒットチャートの週次ランキングデータをBeautiful Soupで集めてみた）今回は、歌詞の情報を用いて順位データとつなぐことにより、どのような単語の入っている洋楽がBillboardにおいてTop10に入る傾向があるのかをLDAを行うことで確かめたいと思います。

データ収集

残念なことに、Billboardのサイトに歌詞の情報は載っていません。そこで、洋楽の歌詞が取り上げられている某サイトをPython(3系)でWebスクレイピングし、名寄せを頑張って順位データと歌詞データを繋ぎます。

幸いなことに某サイトのURLに規則性があったので、アーティスト名からなるURLを生成し、そのURLをWebスクレイピングして楽曲のリストを集め、今回のBillboardのランキングに入った楽曲のみに絞ります。

#アーティストの楽曲一覧の取得

import urllib
from bs4 import BeautifulSoup
from urllib.request import urlopen
import requests
from urllib.error import HTTPError, URLError
import csv, re, time
from requests.exceptions import ConnectionError

f = open('artist_url_list.csv', 'r')
dataReader = csv.reader(f)

#結果の出力用のデータフレームを作る。
data01 =[] #URL
data02 =[] #曲名
data03 =[] #link

for row in dataReader:
       for url in row:
            time.sleep(10.0) #sleep(秒指定)
            try:
                    r = requests.get(url)
                    soup =  BeautifulSoup(r.content, 'html.parser')
                    
                    for body in soup.findAll("td",{'class':'colfirst'}):
                        for link in body.findAll("a"):
                                data01.append(url)
                                data02.append(''.join(link.findAll(text=True)))
                                data03.append(link.get("href"))
                                data = zip(data01,data02,data03)
                                #CSV出力
                                with open('artistpage_result.csv','wt',errors='backslashreplace') as fout:
                                    writecsv = csv.writer(fout,lineterminator='\n')
                                    writecsv.writerows(data)                                   
                                
                                    
            except HTTPError as e:
                print(e.code)
                
            except URLError:
                print("URLError")

#アーティストの楽曲一覧の取得

import urllib

from bs4 import BeautifulSoup

from urllib.request import urlopen

import requests

from urllib.error import HTTPError, URLError

import csv, re, time

from requests.exceptions import ConnectionError

f = open('artist_url_list.csv', 'r')

dataReader = csv.reader(f)

#結果の出力用のデータフレームを作る。

data01 =[] #URL

data02 =[] #曲名

data03 =[] #link

for row in dataReader:

for url in row:

time.sleep(10.0) #sleep(秒指定)

try:

r = requests.get(url)

soup = BeautifulSoup(r.content, 'html.parser')

for body in soup.findAll("td",{'class':'colfirst'}):

for link in body.findAll("a"):

data01.append(url)

data02.append(''.join(link.findAll(text=True)))

data03.append(link.get("href"))

data = zip(data01,data02,data03)

#CSV出力

with open('artistpage_result.csv','wt',errors='backslashreplace') as fout:

writecsv = csv.writer(fout,lineterminator='\n')

writecsv.writerows(data)

except HTTPError as e:

print(e.code)

except URLError:

print("URLError")

楽曲をランキングに含まれるもののみに絞ったら、歌詞詳細ページを取得します。

#歌詞詳細の歌詞該当部分のみ抽出

import urllib
from bs4 import BeautifulSoup
from urllib.request import urlopen
import requests
from urllib.error import HTTPError, URLError
import csv, re, time
from requests.exceptions import ConnectionError

f = open('song_detail_url.csv', 'r')
dataReader = csv.reader(f)

#結果の出力用のデータフレームを作る。
data01 =[] #URL
data02 =[] #歌詞

for row in dataReader:
       for url in row:
            time.sleep(10.0) #sleep(秒指定)
            try:
                    r = requests.get(url)
                    soup =  BeautifulSoup(r.content, 'html.parser')
                    
                    for body in soup.findAll("div",{'id':'content_h'}):
                                data01.append(url)
                                data02.append(''.join(body.findAll(text=True)))
                                data = zip(data01,data02)
                                #CSV出力
                                with open('lyrics_result.csv','wt',errors='backslashreplace') as fout:
                                    writecsv = csv.writer(fout,lineterminator='\n')
                                    writecsv.writerows(data)                                   
                                
                                    
            except HTTPError as e:
                print(e.code)
                           
            except URLError:
                print("URLError")

#歌詞詳細の歌詞該当部分のみ抽出

import urllib

from bs4 import BeautifulSoup

from urllib.request import urlopen

import requests

from urllib.error import HTTPError, URLError

import csv, re, time

from requests.exceptions import ConnectionError

f = open('song_detail_url.csv', 'r')

dataReader = csv.reader(f)

#結果の出力用のデータフレームを作る。

data01 =[] #URL

data02 =[] #歌詞

for row in dataReader:

for url in row:

time.sleep(10.0) #sleep(秒指定)

try:

r = requests.get(url)

soup = BeautifulSoup(r.content, 'html.parser')

for body in soup.findAll("div",{'id':'content_h'}):

data01.append(url)

data02.append(''.join(body.findAll(text=True)))

data = zip(data01,data02)

#CSV出力

with open('lyrics_result.csv','wt',errors='backslashreplace') as fout:

writecsv = csv.writer(fout,lineterminator='\n')

writecsv.writerows(data)

except HTTPError as e:

print(e.code)

except URLError:

print("URLError")

うまいこと歌詞情報を手に入れることができました。ざっと947曲です。

	url	lyrics
0	http://www.lyricsfreak.com/e/eminem/love+the+w...	Just gonna stand there and watch me burnBut th...
1	http://www.lyricsfreak.com/t/taio+cruz/dynamit...	I came to dance-dance-dance-dance (Yeah)I hate...
2	http://www.lyricsfreak.com/t/taylor+swift/mine...	Oh, oh, ohOh, oh, ohYou were in college, worki...
3	http://www.lyricsfreak.com/e/enrique+iglesias/...	One life, one loveEnrique Iglesias, PitbullY'a...
4	http://www.lyricsfreak.com/b/bob/airplanes_208...	Can we pretend that airplanesIn the night sky ...
5	http://www.lyricsfreak.com/m/mike+posner/coole...	If I could write you a song,And make you fall ...
6	http://www.lyricsfreak.com/j/jason+derulo/ridi...	Yea yeah, yeah, yeah, yeah,I'm feeling like a ...
7	http://www.lyricsfreak.com/t/travie+mccoy/bill...	I wanna be a billionaire so freakin' badBuy al...
8	http://www.lyricsfreak.com/d/drake/find+your+l...	I'm more than just an option (hey, hey, hey) R...
9	http://www.lyricsfreak.com/u/usher/omg_2087748...	Oh my goshBaby let meDid it again, so Imma let...
10	http://www.lyricsfreak.com/b/bob/magic_2087969...	I got the magic in meEvery time I touch that t...
11	http://www.lyricsfreak.com/n/nicki+minaj/your+...	[Chorus]Shawty I'm a only tell you this once, ...
12	http://www.lyricsfreak.com/m/maroon+5/misery_2...	Oh yeahOh yeahSo scared of breaking itThat you...
13	http://www.lyricsfreak.com/t/train/hey+soul+si...	Hey, hey, heyYour lipstick stains on the front...
14	http://www.lyricsfreak.com/b/bruno+mars/just+t...	Oh her eyes, her eyesMake the stars look like ...
15	http://www.lyricsfreak.com/l/lady+gaga/alejand...	I know that we are young,And I know that you m...
16	http://www.lyricsfreak.com/l/la+roux/bulletpro...	Been there, done that, messed aroundI'm having...
17	http://www.lyricsfreak.com/f/flo+rida/club+can...	You know I know howTo make 'em stop and stare ...
18	http://www.lyricsfreak.com/s/shontelle/impossi...	I remember years agoSomeone told me I should t...
19	http://www.lyricsfreak.com/p/paramore/the+only...	When I was youngerI saw my daddy cryAnd curse ...
20	http://www.lyricsfreak.com/u/usher/there+goes+...	Yeah, Right,Usher baby, OKYeah man, rightThere...

url lyrics

0 http://www.lyricsfreak.com/e/eminem/love+the+w... Just gonna stand there and watch me burnBut th...

1 http://www.lyricsfreak.com/t/taio+cruz/dynamit... I came to dance-dance-dance-dance (Yeah)I hate...

2 http://www.lyricsfreak.com/t/taylor+swift/mine... Oh, oh, ohOh, oh, ohYou were in college, worki...

3 http://www.lyricsfreak.com/e/enrique+iglesias/... One life, one loveEnrique Iglesias, PitbullY'a...

4 http://www.lyricsfreak.com/b/bob/airplanes_208... Can we pretend that airplanesIn the night sky ...

5 http://www.lyricsfreak.com/m/mike+posner/coole... If I could write you a song,And make you fall ...

6 http://www.lyricsfreak.com/j/jason+derulo/ridi... Yea yeah, yeah, yeah, yeah,I'm feeling like a ...

7 http://www.lyricsfreak.com/t/travie+mccoy/bill... I wanna be a billionaire so freakin' badBuy al...

8 http://www.lyricsfreak.com/d/drake/find+your+l... I'm more than just an option (hey, hey, hey) R...

9 http://www.lyricsfreak.com/u/usher/omg_2087748... Oh my goshBaby let meDid it again, so Imma let...

10 http://www.lyricsfreak.com/b/bob/magic_2087969... I got the magic in meEvery time I touch that t...

11 http://www.lyricsfreak.com/n/nicki+minaj/your+... [Chorus]Shawty I'm a only tell you this once, ...

12 http://www.lyricsfreak.com/m/maroon+5/misery_2... Oh yeahOh yeahSo scared of breaking itThat you...

13 http://www.lyricsfreak.com/t/train/hey+soul+si... Hey, hey, heyYour lipstick stains on the front...

14 http://www.lyricsfreak.com/b/bruno+mars/just+t... Oh her eyes, her eyesMake the stars look like ...

15 http://www.lyricsfreak.com/l/lady+gaga/alejand... I know that we are young,And I know that you m...

16 http://www.lyricsfreak.com/l/la+roux/bulletpro... Been there, done that, messed aroundI'm having...

17 http://www.lyricsfreak.com/f/flo+rida/club+can... You know I know howTo make 'em stop and stare ...

18 http://www.lyricsfreak.com/s/shontelle/impossi... I remember years agoSomeone told me I should t...

19 http://www.lyricsfreak.com/p/paramore/the+only... When I was youngerI saw my daddy cryAnd curse ...

20 http://www.lyricsfreak.com/u/usher/there+goes+... Yeah, Right,Usher baby, OKYeah man, rightThere...

Rによる分析

ここから、Rにてテキストマイニングを行いたいと思います。まず、tmパッケージを用いて、不要語（stop word）を除去します。具体的にはtheとかyouとかを除外しています。

#歌詞データの読み込み
lyrics_dataset <- read.csv(file = "lyrics_result.csv",as.is = TRUE,header = FALSE)
colnames(lyrics_dataset) <- c("link","lyrics")

library(tm)

#歌詞を小文字にする
document_dataset$lyrics <- tolower(document_dataset$lyrics)

#stop wordを除去する
stopwords_regex = paste(stopwords('en'), collapse = '\\b|\\b')
stopwords_regex = paste0('\\b', stopwords_regex, '\\b')
document_dataset$lyrics <- stringr::str_replace_all(document_dataset$lyrics, stopwords_regex, '')

#歌詞データの読み込み

lyrics_dataset <- read.csv(file = "lyrics_result.csv",as.is = TRUE,header = FALSE)

colnames(lyrics_dataset) <- c("link","lyrics")

library(tm)

#歌詞を小文字にする

document_dataset$lyrics <- tolower(document_dataset$lyrics)

#stop wordを除去する

stopwords_regex = paste(stopwords('en'), collapse = '\\b|\\b')

stopwords_regex = paste0('\\b', stopwords_regex, '\\b')

document_dataset$lyrics <- stringr::str_replace_all(document_dataset$lyrics, stopwords_regex, '')

続いて、LDAを実行できるtopicmodelsパッケージで扱えるようにするために、テキストデータに以下の処理を施します。

#文書単語行列の作成

corpus <- Corpus(VectorSource(document_dataset$lyrics))
inspect(corpus)
dtm <- DocumentTermMatrix(corpus)
findFreqTerms(dtm)

#文書単語行列の作成

corpus <- Corpus(VectorSource(document_dataset$lyrics))

inspect(corpus)

dtm <- DocumentTermMatrix(corpus)

findFreqTerms(dtm)

あとは以下のコードでLDAを実行するだけです。トピック数はアドホックに20としています。研究者の方、いい加減ですみません。

#LDAの実行

library(topicmodels)

nbo_topics <- 20
lda <- LDA(dtm,control=list(verbose=1), k = nbo_topics,method = "Gibbs")

#LDAの実行

library(topicmodels)

nbo_topics <- 20

lda <- LDA(dtm,control=list(verbose=1), k = nbo_topics,method = "Gibbs")

LDAの結果

まずは推定されたトピックごとの上位10単語をみてみます。トピック1はラブソングとかでしょうか。トピック17にパリピっぽい単語が、トピック18にスラングが含まれていますね。

#トピックの上位10単語を確認する
terms_each_topics <- data.frame(terms(lda,10))

1 2	#トピックの上位10単語を確認する terms_each_topics <- data.frame(terms(lda,10))

> terms_each_topics
   Topic.1 Topic.2 Topic.3 Topic.4 Topic.5 Topic.6 Topic.7 Topic.8 Topic.9 Topic.10 Topic.11   Topic.12 Topic.13 Topic.14 Topic.15 Topic.16 Topic.17 Topic.18 Topic.19　Topic.20
1     love    know    feel  better   wanna     hey     one  chorus    back     yeah    gonna      never     like      got      ’m     keep     stop      ain      let      get
2     like     now   heart   world    want    said     ooh     way    time     baby      low       will      new     like   don’t     even     just     like      can     good
3     make    just    life    whoa    take     old    call   verse    know     girl     hear      still     city     back     ’re  getting    hands     shit      say    night
4    touch    need    away     run    rock   every   cause     can    come   little     mean       ever     high    right    young      one      put     fuck      believe      ain
5     know   think    just   light     see     woo  gettin     pre    long     like    sound       eyes     bout     wish      ’s     give    party      got      fly     kind
6   nobody   cause   break  things    kiss    left    born    tell    like     just    shaky     always     ride     hold     ’ll     lose    crazy     hook      fall   really
7     baby    much    real    find    come    told     day     got    best    right     just everything     yeah     boom      que   please     live      gon      made    sleep
8    cause    give  enough    show    make nothing   makes    like    home      get     solo      leave      fun     know   ain’t   youand   lights    nigga      first champion
9     name    mind   every     see    body   daddy    came    used     til      can  tonight       hope     know      one     para     just     play   niggas      words    catch
10  loving  really    find waiting tonight   sweet   stand    made alright     look   wicked       stay      get     come   can’t  without      see    money      lonely      til

> terms_each_topics

Topic.1 Topic.2 Topic.3 Topic.4 Topic.5 Topic.6 Topic.7 Topic.8 Topic.9 Topic.10 Topic.11 Topic.12 Topic.13 Topic.14 Topic.15 Topic.16 Topic.17 Topic.18 Topic.19　Topic.20

1 love know feel better wanna hey one chorus back yeah gonna never like got ’m keep stop ain let get

2 like now heart world want said ooh way time baby low will new like don’t even just like can good

3 make just life whoa take old call verse know girl hear still city back ’re getting hands shit say night

4 touch need away run rock every cause can come little mean ever high right young one put fuck believe ain

5 know think just light see woo gettin pre long like sound eyes bout wish ’s give party got fly kind

6 nobody cause break things kiss left born tell like just shaky always ride hold ’ll lose crazy hook fall really

7 baby much real find come told day got best right just everything yeah boom que please live gon made sleep

8 cause give enough show make nothing makes like home get solo leave fun know ain’t youand lights nigga first champion

9 name mind every see body daddy came used til can tonight hope know one para just play niggas words catch

10 loving really find waiting tonight sweet stand made alright look wicked stay get come can’t without see money lonely til

見ずらいので、行を一つにまとめて、トピックにidを割り振ります。

#トピックの上位10単語をまとめて、トピックにidをふる
topic_keywords <- data.frame(topic_keywords_10=apply(t(terms_each_topics),1,paste,collapse=","))
topic_keywords <- topic_keywords %>% mutate(topic_id=1:n())

#トピックの上位10単語をまとめて、トピックにidをふる

topic_keywords <- data.frame(topic_keywords_10=apply(t(terms_each_topics),1,paste,collapse=","))

topic_keywords <- topic_keywords %>% mutate(topic_id=1:n())

> topic_keywords
                                              topic_keywords_10 topic_id
1       love,like,make,touch,know,nobody,baby,cause,name,loving        1
2          know,now,just,need,think,cause,much,give,mind,really        2
3        feel,heart,life,away,just,break,real,enough,every,find        3
4      better,world,whoa,run,light,things,find,show,see,waiting        4
5          wanna,want,take,rock,see,kiss,come,make,body,tonight        5
6          hey,said,old,every,woo,left,told,nothing,daddy,sweet        6
7           one,ooh,call,cause,gettin,born,day,makes,came,stand        7
8              chorus,way,verse,can,pre,tell,got,like,used,made        8
9           back,time,know,come,long,like,best,home,til,alright        9
10           yeah,baby,girl,little,like,just,right,get,can,look       10
11     gonna,low,hear,mean,sound,shaky,just,solo,tonight,wicked       11
12 never,will,still,ever,eyes,always,everything,leave,hope,stay       12
13               like,new,city,high,bout,ride,yeah,fun,know,get       13
14             got,like,back,right,wish,hold,boom,know,one,come       14
15        ’m,don’t,’re,young,’s,’ll,que,ain’t,para,can’t       15
16   keep,even,getting,one,give,lose,please,youand,just,without       16
17         stop,just,hands,put,party,crazy,live,lights,play,see       17
18           ain,like,shit,fuck,got,hook,gon,nigga,niggas,money       18
19         let,can,say,believe,fly,fall,made,first,words,lonely       19
20      get,good,night,ain,kind,really,sleep,champion,catch,til       20

> topic_keywords

topic_keywords_10 topic_id

1 love,like,make,touch,know,nobody,baby,cause,name,loving 1

2 know,now,just,need,think,cause,much,give,mind,really 2

3 feel,heart,life,away,just,break,real,enough,every,find 3

4 better,world,whoa,run,light,things,find,show,see,waiting 4

5 wanna,want,take,rock,see,kiss,come,make,body,tonight 5

6 hey,said,old,every,woo,left,told,nothing,daddy,sweet 6

7 one,ooh,call,cause,gettin,born,day,makes,came,stand 7

8 chorus,way,verse,can,pre,tell,got,like,used,made 8

9 back,time,know,come,long,like,best,home,til,alright 9

10 yeah,baby,girl,little,like,just,right,get,can,look 10

11 gonna,low,hear,mean,sound,shaky,just,solo,tonight,wicked 11

12 never,will,still,ever,eyes,always,everything,leave,hope,stay 12

13 like,new,city,high,bout,ride,yeah,fun,know,get 13

14 got,like,back,right,wish,hold,boom,know,one,come 14

15 ’m,don’t,’re,young,’s,’ll,que,ain’t,para,can’t 15

16 keep,even,getting,one,give,lose,please,youand,just,without 16

17 stop,just,hands,put,party,crazy,live,lights,play,see 17

18 ain,like,shit,fuck,got,hook,gon,nigga,niggas,money 18

19 let,can,say,believe,fly,fall,made,first,words,lonely 19

20 get,good,night,ain,kind,really,sleep,champion,catch,til 20

最後に、BillboardでTop10に入ったかどうかのデータを作っておき、そのデータと各歌詞を繋ぎ、各歌詞ごとに割りふられた確率が最大のトピックで集計をします。

#上位10位に入ったことがあるかどうかのダミーを作成
top_10_songs <- merge_dataset %>% group_by(link.y) %>% summarise(top_10=sum(top_10))
top_10_songs <- top_10_songs %>% mutate(top_10_dummy=ifelse(top_10>0,1,0))

#割り振られた最大の確率のトピックを抽出し、歌詞データと統合する
topics_each_document <- data.frame(topic_id=topics(lda,1))
topics_each_document <- cbind(link.y=document_dataset$link.y,topics_each_document)

#上位10位ダミーを繋ぎ、トピックの上位10単語の表現も繋ぐ
topics_each_document <- topics_each_document %>% left_join(top_10_songs,by="link.y")
topics_each_document <- topics_each_document %>% left_join(topic_keywords,by="topic_id")

#トピックごとのBillboardのTop10ランクイン割合をもとめる
topics_each_document %>% group_by(topic_keywords_10) %>% summarise(top_10_dummy=mean(top_10_dummy),count=n())

# A tibble: 20 × 3
                                              topic_keywords_10 `mean(top_10_dummy)` count
                                                         <fctr>                <dbl> <int>
1         ’m,don’t,’re,young,’s,’ll,que,ain’t,para,can’t           0.12195122    41
2            ain,like,shit,fuck,got,hook,gon,nigga,niggas,money           0.14705882    68
3           back,time,know,come,long,like,best,home,til,alright           0.21666667    60
4      better,world,whoa,run,light,things,find,show,see,waiting           0.25000000    40
5              chorus,way,verse,can,pre,tell,got,like,used,made           0.14814815    54
6        feel,heart,life,away,just,break,real,enough,every,find           0.10416667    48
7       get,good,night,ain,kind,really,sleep,champion,catch,til           0.22857143    35
8      gonna,low,hear,mean,sound,shaky,just,solo,tonight,wicked           0.22857143    35
9              got,like,back,right,wish,hold,boom,know,one,come           0.13888889    36
10         hey,said,old,every,woo,left,told,nothing,daddy,sweet           0.09756098    41
11   keep,even,getting,one,give,lose,please,youand,just,without           0.18918919    37
12         know,now,just,need,think,cause,much,give,mind,really           0.21333333    75
13         let,can,say,believe,fly,fall,made,first,words,lonely           0.15217391    46
14               like,new,city,high,bout,ride,yeah,fun,know,get           0.14705882    34
15      love,like,make,touch,know,nobody,baby,cause,name,loving           0.17021277    47
16 never,will,still,ever,eyes,always,everything,leave,hope,stay           0.22222222    72
17          one,ooh,call,cause,gettin,born,day,makes,came,stand           0.28205128    39
18         stop,just,hands,put,party,crazy,live,lights,play,see           0.25000000    48
19         wanna,want,take,rock,see,kiss,come,make,body,tonight           0.10810811    37
20           yeah,baby,girl,little,like,just,right,get,can,look           0.14814815    54

#上位10位に入ったことがあるかどうかのダミーを作成

top_10_songs <- merge_dataset %>% group_by(link.y) %>% summarise(top_10=sum(top_10))

top_10_songs <- top_10_songs %>% mutate(top_10_dummy=ifelse(top_10>0,1,0))

#割り振られた最大の確率のトピックを抽出し、歌詞データと統合する

topics_each_document <- data.frame(topic_id=topics(lda,1))

topics_each_document <- cbind(link.y=document_dataset$link.y,topics_each_document)

#上位10位ダミーを繋ぎ、トピックの上位10単語の表現も繋ぐ

topics_each_document <- topics_each_document %>% left_join(top_10_songs,by="link.y")

topics_each_document <- topics_each_document %>% left_join(topic_keywords,by="topic_id")

#トピックごとのBillboardのTop10ランクイン割合をもとめる

topics_each_document %>% group_by(topic_keywords_10) %>% summarise(top_10_dummy=mean(top_10_dummy),count=n())

# A tibble: 20 × 3

topic_keywords_10 `mean(top_10_dummy)` count

1 ’m,don’t,’re,young,’s,’ll,que,ain’t,para,can’t 0.12195122 41

2 ain,like,shit,fuck,got,hook,gon,nigga,niggas,money 0.14705882 68

3 back,time,know,come,long,like,best,home,til,alright 0.21666667 60

4 better,world,whoa,run,light,things,find,show,see,waiting 0.25000000 40

5 chorus,way,verse,can,pre,tell,got,like,used,made 0.14814815 54

6 feel,heart,life,away,just,break,real,enough,every,find 0.10416667 48

7 get,good,night,ain,kind,really,sleep,champion,catch,til 0.22857143 35

8 gonna,low,hear,mean,sound,shaky,just,solo,tonight,wicked 0.22857143 35

9 got,like,back,right,wish,hold,boom,know,one,come 0.13888889 36

10 hey,said,old,every,woo,left,told,nothing,daddy,sweet 0.09756098 41

11 keep,even,getting,one,give,lose,please,youand,just,without 0.18918919 37

12 know,now,just,need,think,cause,much,give,mind,really 0.21333333 75

13 let,can,say,believe,fly,fall,made,first,words,lonely 0.15217391 46

14 like,new,city,high,bout,ride,yeah,fun,know,get 0.14705882 34

15 love,like,make,touch,know,nobody,baby,cause,name,loving 0.17021277 47

16 never,will,still,ever,eyes,always,everything,leave,hope,stay 0.22222222 72

17 one,ooh,call,cause,gettin,born,day,makes,came,stand 0.28205128 39

18 stop,just,hands,put,party,crazy,live,lights,play,see 0.25000000 48

19 wanna,want,take,rock,see,kiss,come,make,body,tonight 0.10810811 37

20 yeah,baby,girl,little,like,just,right,get,can,look 0.14814815 54

#トピックごとのランクイン割合を高い順に棒グラフで描写

library(ggplot2)
ggplot(topics_and_top10, aes(x=reorder(topic_keywords_10,top_10_dummy), y=top_10_dummy)) +
  geom_bar(stat='identity') + coord_flip() + xlab("topics")

#トピックごとのランクイン割合を高い順に棒グラフで描写

library(ggplot2)

ggplot(topics_and_top10, aes(x=reorder(topic_keywords_10,top_10_dummy), y=top_10_dummy)) +

geom_bar(stat='identity') + coord_flip() + xlab("topics")

BillboardのTop10ランクイン割合の高いトピックTop3
「one,ooh,call,cause,gettin,born,day,makes,came,stand」
「better,world,whoa,run,light,things,find,show,see,waiting」・・・明るい感じ？
「stop,just,hands,put,party,crazy,live,lights,play,see」・・・パリピぽい

BillboardのTop10ランクイン割合の低いトピックTop3
「wanna,want,take,rock,see,kiss,come,make,body,tonight」・・・欲求系？
「feel,heart,life,away,just,break,real,enough,every,find」・・・癒し系？
「hey,said,old,every,woo,left,told,nothing,daddy,sweet」

あまり洋楽を聴かないので、得られたトピックの解釈が中々できないのがもどかしいです。ただ、スラングの歌詞を含む歌詞はそんなにランクイン割合が悪いわけではなさそうですね。洋楽をもっと聴いて、前処理などもう少し工夫してリベンジしたいですね。

参考文献

トピックモデルによる統計的潜在意味解析 (自然言語処理シリーズ)

Pythonクローリング&スクレイピング -データ収集・解析のための実践開発ガイド-

モダンなRによるテキスト解析 topicmodels: An R Package for Fitting Topic Models

はじめに

Uplift Modelingとはなにか

Uplift Modelingの卑近な例え話

Uplift Modelingのサンプルデータ

tools4upliftについて

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

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

おわりに

参考文献

はじめに

分類タスク

学習済み分散表現

前処理

蒙古タンメン中本の口コミ4000件から作成した分散表現：Gensim

白ヤギコーポレーションのモデル：Gensim

東北大学 乾・岡崎研究室のモデル：Gensim

Facebookの学習済みFastTextモデル：Gensim

NWJC から取得した単語の分散表現データ (nwjc2vec)：Gensim

NNLM embedding trained on Google News：TensorFlow

比較

おわりに

参考情報

はじめに

後輩からの聞き込み

F0とF1

ピッチはどうすれば求めることができるのか？

今回扱うモジュール／今回のツール

データ

プログラム

実行結果

今後の使いみち

追記

参考情報

モチベーション

環境

Gensimで行うLDA

Rでもやってみる

トピックの吐き出し

参考情報

はじめに

目次

目的

データ収集

形態素解析

集計

点数のヒストグラム

単語数のヒストグラム

単語数と点数の散布図

形態素解析結果の集計、単語ランキング

分散表現とは

word2vecについて

gensimのword2vecの引数

word2vecによる文書分類の適用

終わりに

参考情報

はじめに

LIMEとは

データ

今回の分析の目的

進め方

コード

結果の解釈

参考情報

先行研究

コード

参考文献

Edwardとは

Edwardでできること

参考スライド

参考文献

東北大学乾・岡崎研究室のモデル：Gensim