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

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

先行研究

海外のブログ(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

階層ベイズモデルの直帰率分析への適用 with rstan

松浦さんの『StanとRでベイズ統計モデリング』の8章の階層ベイズがすごくわかりやすいなぁと思いつつも、自分の持っているデータで試していなかったので、これを機に実践してみようと思います。
やや変数を追加しているくらいで大した変更点はありませんが、題材としては当ブログのアクセスログにおける直帰率に関するデータで、どのような要素が直帰率に影響を与えるのかを分析します。

目次
・モデル概要
・前処理
・推定
・結果（非階層モデルとの比較）
・参考文献

モデル概要

モデルは8章のロジスティック回帰の階層モデルに一部変数を追加していますが、ほぼそのままです。記事ごとのパラメータやリファラーごとのパラメータを想定しています。
Nは記事数でnはそのインデックス、Cはリファラーの数でcはそのインデックス、Iはログとして残っているセッションの数でiはそのインデックスとなっています。hatebuは記事のはてぶ数、stringlineは記事の行数、holidayは休日or祝日ダミー変数、daytimeは12:00~18:00なら1をとるダミー変数、revisitedは再訪問ユーザーなら1を取るダミー変数となっています。記事ごと・リファラーごとに直帰のしやすさが違う（パラメータが従う正規分布のパラメータがそれぞれ異なる）という仮定のもとに立ったモデルとなります。

$$x[i] = b_{1} + x_{記事}[記事ID[i]] \\ + x_{リファラー}[リファラーID[i]] + x_{セッション}[i] $$

$$q[i] = inverselogit(x[i]) $$

$$Y[i] \sim Bernoulli(q[i]) $$

$$x_{記事}[n] = b_{2}hatebu + b_{3}stringline[n] \\ + b_{記事間の差}[n] $$

$$b_{記事間の差}[n] \sim Normal(0, \sigma_{記事番号}) $$

$$x_{リファラー}[c] = b_{リファラー間の差}[c] $$

$$b_{リファラー間の差}[c] \sim Normal(0, \sigma_{リファラー番号}) $$

$$x_{セッション}[i] = b_{4}holiday[i] 　+ b_{5}divice[i] \\ + b_{6}daytime[i] + b_{7}revisited[i] $$

前処理

GAのAPIからデータを取得して1セッション1記事になるようにデータを作成しています。数ヶ月で25000件ほどデータがあったのですが、計算に時間がかかるので、データ数を2400件くらいにサンプリングしています。

library(RGA)
library(tidyverse)
library(Nippon)

authorize()
prof <-list_profiles()

start_date <- "2017-04-01"
end_date   <- "2017-10-20"

accesslogdata <- get_ga(profileId = prof$id[2],
                        start.date = start_date,
                        end.date = end_date,
                        dimensions = "ga:pagePath,
                        ga:dateHourMinute,
                        ga:deviceCategory,
                        ga:userType,
                        ga:referralPath,
                        ga:fullReferrer",
                        sort = "-ga:sessions",
                        metrics = "ga:sessions,ga:bounces",
                        fetch.by = "day")

#reshape_url
accesslogdata$pagePath <- vapply(strsplit(accesslogdata$pagePath,"\\?"), `[`, 1, FUN.VALUE=character(1))
accesslogdata <- accesslogdata %>% filter(grepl(x = pagePath,"/archives/[0-9]+$"))
accesslogdata$pagePath <- gsub(accesslogdata$pagePath,pattern = "/archives/",replacement = "article_")

accesslogdata_filtered <- accesslogdata %>% filter(sessions == 1)
accesslogdata_filtered <- accesslogdata_filtered %>% 
                          mutate(fullReferrer = ifelse(grepl(x = fullReferrer,'t.co/'),'twitter', fullReferrer)) %>% 
                          mutate(fullReferrer = ifelse(grepl(x = fullReferrer,'b.hatena.ne.jp/'),'b.hatena.ne.jp', fullReferrer)) %>% 
                          mutate(fullReferrer = ifelse(grepl(x = fullReferrer,'d.hatena.ne.jp/'),'d.hatena.ne.jp', fullReferrer)) %>% 
                          mutate(fullReferrer = ifelse(grepl(x = fullReferrer,'facebook.com/'),'facebook.com', fullReferrer))

#refferer list
refferer_cat <- c("google", "(direct)", "twitter", 
                  "yahoo", "b.hatena.ne.jp", "facebook.com",
                  "bing")

daytime_cat <- c("12","13","14","15","16","17","18")

#making weekday data and daytime data
accesslogdata_filtered <- accesslogdata_filtered %>% 
                            filter(fullReferrer %in% refferer_cat) %>% 
                            mutate(date=as.Date(format(substr(accesslogdata_filtered$dateHourMinute,start = 1,stop = 12), format="%Y%m%d%"),format = "%Y%m%d")) %>% 
                            mutate(hourminutes = substr(accesslogdata_filtered$dateHourMinute,start = 9,stop = 10 )) %>% 
                            mutate(holiday = ifelse(is.jholiday(date), 1, 0)) %>% 
                            mutate(daytime = ifelse(hourminutes %in% daytime_cat, 1, 0)) %>% 
                            mutate(device = ifelse(deviceCategory == "desktop", 1, 0)) %>% 
                            mutate(revisited = ifelse(userType == "Returning Visitor", 1, 0))

selected_dataset <- accesslogdata_filtered %>% 
                      select(pagePath,fullReferrer,
                             device,revisited,holiday,daytime,bounces) %>% 
                      filter(!(pagePath %in% c("article_10", "article_53") )) %>% 
                      bind_cols(rand =runif(nrow(selected_dataset), min = 0, max = 1)) %>% 
                      filter(rand <= 2400/nrow(selected_dataset))

referer_cat <- selected_dataset %>% 
                  select(fullReferrer) %>% 
                  distinct() %>% 
                  mutate(referer_id = 1:n())

article_cat <- selected_dataset %>% 
                  select(pagePath) %>% 
                  distinct() %>% 
                  mutate(article_id = 1:n())

selected_dataset <- selected_dataset %>% 
                      left_join(article_cat, by="pagePath") %>% 
                      left_join(referer_cat, by="fullReferrer") %>% 
                      select(-pagePath, -fullReferrer)

article_data <- read_csv(file = "kamonohashiperry_text.csv")
article_data$url <- gsub(article_data$url,
                         pattern = "http://kamonohashiperry.com/archives/",replacement = "article_")

article_data <- article_data %>% left_join(article_cat, by=c("url"="pagePath"))
article_data <- article_data[!(is.na(article_data$article_id)),]

library(RGA)

library(tidyverse)

library(Nippon)

authorize()

prof <-list_profiles()

start_date <- "2017-04-01"

end_date <- "2017-10-20"

accesslogdata <- get_ga(profileId = prof$id[2],

start.date = start_date,

end.date = end_date,

dimensions = "ga:pagePath,

ga:dateHourMinute,

ga:deviceCategory,

ga:userType,

ga:referralPath,

ga:fullReferrer",

sort = "-ga:sessions",

metrics = "ga:sessions,ga:bounces",

fetch.by = "day")

#reshape_url

accesslogdata$pagePath <- vapply(strsplit(accesslogdata$pagePath,"\\?"), `[`, 1, FUN.VALUE=character(1))

accesslogdata <- accesslogdata %>% filter(grepl(x = pagePath,"/archives/[0-9]+$"))

accesslogdata$pagePath <- gsub(accesslogdata$pagePath,pattern = "/archives/",replacement = "article_")

accesslogdata_filtered <- accesslogdata %>% filter(sessions == 1)

accesslogdata_filtered <- accesslogdata_filtered %>%

mutate(fullReferrer = ifelse(grepl(x = fullReferrer,'t.co/'),'twitter', fullReferrer)) %>%

mutate(fullReferrer = ifelse(grepl(x = fullReferrer,'b.hatena.ne.jp/'),'b.hatena.ne.jp', fullReferrer)) %>%

mutate(fullReferrer = ifelse(grepl(x = fullReferrer,'d.hatena.ne.jp/'),'d.hatena.ne.jp', fullReferrer)) %>%

mutate(fullReferrer = ifelse(grepl(x = fullReferrer,'facebook.com/'),'facebook.com', fullReferrer))

#refferer list

refferer_cat <- c("google", "(direct)", "twitter",

"yahoo", "b.hatena.ne.jp", "facebook.com",

"bing")

daytime_cat <- c("12","13","14","15","16","17","18")

#making weekday data and daytime data

accesslogdata_filtered <- accesslogdata_filtered %>%

filter(fullReferrer %in% refferer_cat) %>%

mutate(date=as.Date(format(substr(accesslogdata_filtered$dateHourMinute,start = 1,stop = 12), format="%Y%m%d%"),format = "%Y%m%d")) %>%

mutate(hourminutes = substr(accesslogdata_filtered$dateHourMinute,start = 9,stop = 10 )) %>%

mutate(holiday = ifelse(is.jholiday(date), 1, 0)) %>%

mutate(daytime = ifelse(hourminutes %in% daytime_cat, 1, 0)) %>%

mutate(device = ifelse(deviceCategory == "desktop", 1, 0)) %>%

mutate(revisited = ifelse(userType == "Returning Visitor", 1, 0))

selected_dataset <- accesslogdata_filtered %>%

select(pagePath,fullReferrer,

device,revisited,holiday,daytime,bounces) %>%

filter(!(pagePath %in% c("article_10", "article_53") )) %>%

bind_cols(rand =runif(nrow(selected_dataset), min = 0, max = 1)) %>%

filter(rand <= 2400/nrow(selected_dataset))

referer_cat <- selected_dataset %>%

select(fullReferrer) %>%

distinct() %>%

mutate(referer_id = 1:n())

article_cat <- selected_dataset %>%

select(pagePath) %>%

distinct() %>%

mutate(article_id = 1:n())

selected_dataset <- selected_dataset %>%

left_join(article_cat, by="pagePath") %>%

left_join(referer_cat, by="fullReferrer") %>%

select(-pagePath, -fullReferrer)

article_data <- read_csv(file = "kamonohashiperry_text.csv")

article_data$url <- gsub(article_data$url,

pattern = "http://kamonohashiperry.com/archives/",replacement = "article_")

article_data <- article_data %>% left_join(article_cat, by=c("url"="pagePath"))

article_data <- article_data[!(is.na(article_data$article_id)),]

推定

stanコードはこちらになります。

data {
  int N; #the number of article
  int C; #the number of referer
  int I; #the number of log
  int<lower=0> hatebu[N]; 
  real<lower=0> stringline[N];
  int<lower=1, upper=N> article_id[I];
  int<lower=1, upper=C> referer_id[I];
  real<lower=0, upper=1> holiday[I];
  real<lower=0, upper=1> device[I];
  real<lower=0, upper=1> daytime[I];
  real<lower=0, upper=1> revisited[I];
  int<lower=0, upper=1> Y[I];
}

parameters {
  real b[7];
  real b_P[N];
  real b_C[C];
  real<lower=0> s_P;
  real<lower=0> s_C;
}

transformed parameters {
  real x_P[N];
  real x_C[C];
  real x_J[I];
  real x[I];
  real q[I];
  for (n in 1:N)
    x_P[n] = b[2]*stringline[n] + b[3]*hatebu[n] + b_P[n];
  for (c in 1:C)
    x_C[c] = b_C[c];
  for (i in 1:I) {
    x_J[i] = b[4]*holiday[i] + b[5]*device[i] + b[6]*daytime[i] + b[7]*revisited[i];
    x[i] = b[1] + x_P[article_id[i]] + x_C[referer_id[i]] + x_J[i];
    q[i] = inv_logit(x[i]);
  }
}

model {
  for (n in 1:N)
    b_P[n] ~ normal(0, s_P);
  for (c in 1:C)
    b_C[c] ~ normal(0, s_C);
  for (i in 1:I)
    Y[i] ~ bernoulli(q[i]);
}

data {

int N; #the number of article

int C; #the number of referer

int I; #the number of log

int<lower=0> hatebu[N];

real<lower=0> stringline[N];

int<lower=1, upper=N> article_id[I];

int<lower=1, upper=C> referer_id[I];

real<lower=0, upper=1> holiday[I];

real<lower=0, upper=1> device[I];

real<lower=0, upper=1> daytime[I];

real<lower=0, upper=1> revisited[I];

int<lower=0, upper=1> Y[I];

}

parameters {

real b[7];

real b_P[N];

real b_C[C];

real<lower=0> s_P;

real<lower=0> s_C;

}

transformed parameters {

real x_P[N];

real x_C[C];

real x_J[I];

real x[I];

real q[I];

for (n in 1:N)

x_P[n] = b[2]*stringline[n] + b[3]*hatebu[n] + b_P[n];

for (c in 1:C)

x_C[c] = b_C[c];

for (i in 1:I) {

x_J[i] = b[4]*holiday[i] + b[5]*device[i] + b[6]*daytime[i] + b[7]*revisited[i];

x[i] = b[1] + x_P[article_id[i]] + x_C[referer_id[i]] + x_J[i];

q[i] = inv_logit(x[i]);

}

model {

for (n in 1:N)

b_P[n] ~ normal(0, s_P);

for (c in 1:C)

b_C[c] ~ normal(0, s_C);

for (i in 1:I)

Y[i] ~ bernoulli(q[i]);

}

rstanを用いたstan実行用のRコードです。ヴァイオリンプロットで主要な係数の分布を見る処理も書かれています。

library(rstan)

N <- nrow(article_cat)
C <- nrow(referer_cat)
I <- nrow(selected_dataset)

data <- list(N = N,
             C = C,
             hatebu = article_data$hatebu,
             stringline = article_data$str_lines,
             article_id = selected_dataset$article_id,
             referer_id = selected_dataset$referer_id,
             holiday = selected_dataset$holiday,
             device = selected_dataset$device,
             daytime = selected_dataset$daytime,
             revisited = selected_dataset$revisited,
             Y = selected_dataset$bounces)

fit <- stan(file = "model/model8-8_access_analysis.stan",
            data = data,
            pars = c("b", "b_P", "b_C", "s_P", "s_C", "q"),
            seed = 1234)

source('../common.R')

ms <- rstan::extract(fit)
N_mcmc <- length(ms$lp__)

param_names <- c('mcmc', paste0('b', 1:7), 's_P', 's_C')
d_est <- data.frame(1:N_mcmc, ms$b, ms$s_P, ms$s_C)
colnames(d_est) <- param_names
d_qua <- data.frame.quantile.mcmc(x=param_names[-1], y_mcmc=d_est[,-1])
d_melt <- reshape2::melt(d_est, id=c('mcmc'), variable.name='X')
d_melt$X <- factor(d_melt$X, levels=rev(levels(d_melt$X)))

p <- ggplot()
p <- p + theme_bw(base_size=18)
p <- p + coord_flip()
p <- p + geom_violin(data=d_melt, aes(x=X, y=value), fill='white', color='grey80', size=2, alpha=0.3, scale='width')
p <- p + geom_pointrange(data=d_qua, aes(x=X, y=p50, ymin=p2.5, ymax=p97.5), size=1)
p <- p + labs(x='parameter', y='value')
p <- p + scale_y_continuous(breaks=seq(from=-2, to=6, by=2))
p

library(rstan)

N <- nrow(article_cat)

C <- nrow(referer_cat)

I <- nrow(selected_dataset)

data <- list(N = N,

C = C,

hatebu = article_data$hatebu,

stringline = article_data$str_lines,

article_id = selected_dataset$article_id,

referer_id = selected_dataset$referer_id,

holiday = selected_dataset$holiday,

device = selected_dataset$device,

daytime = selected_dataset$daytime,

revisited = selected_dataset$revisited,

Y = selected_dataset$bounces)

fit <- stan(file = "model/model8-8_access_analysis.stan",

data = data,

pars = c("b", "b_P", "b_C", "s_P", "s_C", "q"),

seed = 1234)

source('../common.R')

ms <- rstan::extract(fit)

N_mcmc <- length(ms$lp__)

param_names <- c('mcmc', paste0('b', 1:7), 's_P', 's_C')

d_est <- data.frame(1:N_mcmc, ms$b, ms$s_P, ms$s_C)

colnames(d_est) <- param_names

d_qua <- data.frame.quantile.mcmc(x=param_names[-1], y_mcmc=d_est[,-1])

d_melt <- reshape2::melt(d_est, id=c('mcmc'), variable.name='X')

d_melt$X <- factor(d_melt$X, levels=rev(levels(d_melt$X)))

p <- ggplot()

p <- p + theme_bw(base_size=18)

p <- p + coord_flip()

p <- p + geom_violin(data=d_melt, aes(x=X, y=value), fill='white', color='grey80', size=2, alpha=0.3, scale='width')

p <- p + geom_pointrange(data=d_qua, aes(x=X, y=p50, ymin=p2.5, ymax=p97.5), size=1)

p <- p + labs(x='parameter', y='value')

p <- p + scale_y_continuous(breaks=seq(from=-2, to=6, by=2))

結果

係数を見る限りは、符号の向きが確かなのはb5(PCダミー)とb7(再訪問ユーザーダミー)なので、PCの方が直帰しにくく、再訪問ユーザーの方が直帰しにくいという傾向があると考えることができます。

教科書ではAUCを非階層モデルと比較していましたので、比較してみようと思います。
AUCの計算を行うためのコードもGithubに載っていましたのでそちらを使います。

#calculate auc
library(pROC)

ms <- rstan::extract(fit)
N_mcmc <- length(ms$lp__)
spec <- seq(from=0, to=1, len=201)
probs <- c(0.1, 0.5, 0.9)

auces <- numeric(N_mcmc)
m_roc <- matrix(nrow=length(spec), ncol=N_mcmc)
for (i in 1:N_mcmc) {
  roc_res <- roc(selected_dataset$bounces, ms$q[i,])
  auces[i] <- as.numeric(roc_res$auc)
  m_roc[,i] <- coords(roc_res, x=spec, input='specificity', ret='sensitivity')
}
quantile(auces, prob=probs)

#calculate auc

library(pROC)

ms <- rstan::extract(fit)

N_mcmc <- length(ms$lp__)

spec <- seq(from=0, to=1, len=201)

probs <- c(0.1, 0.5, 0.9)

auces <- numeric(N_mcmc)

m_roc <- matrix(nrow=length(spec), ncol=N_mcmc)

for (i in 1:N_mcmc) {

roc_res <- roc(selected_dataset$bounces, ms$q[i,])

auces[i] <- as.numeric(roc_res$auc)

m_roc[,i] <- coords(roc_res, x=spec, input='specificity', ret='sensitivity')

}

quantile(auces, prob=probs)

＃ロジスティック回帰の階層ベイズ推定でのAUC
> quantile(auces, prob=probs)
      10%       50%       90% 
0.6683084 0.6770603 0.6846594


#ロジスティック回帰のベイズ推定でのAUC
> quantile(auces_non_hiral, prob=probs_non_hiral)
      10%       50%       90% 
0.5463287 0.5623697 0.5776898

＃ロジスティック回帰の階層ベイズ推定でのAUC

> quantile(auces, prob=probs)

10% 50% 90%

0.6683084 0.6770603 0.6846594

#ロジスティック回帰のベイズ推定でのAUC

> quantile(auces_non_hiral, prob=probs_non_hiral)

10% 50% 90%

0.5463287 0.5623697 0.5776898

80%が良いとされているAUCには程遠いですが、記事やリファラーごとの差を考慮しない非階層のものよりもAUCが高いと言えます。
ちなみに、教科書の例のAUCは80%ほどでした。
Webマーケのデータ分析においてロジットは汎用性が高いで、今回のコードを土台に色々と業務で試していこうと思います。

参考文献

StanとRでベイズ統計モデリング (Wonderful R)
ベイズ統計モデリング: R,JAGS, Stanによるチュートリアル原著第2版

確率的プログラミングライブラリ「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入門

RのFactoRizationMachinesパッケージを用いたFMのページビューデータへの適用

Googleタグマネージャーで集めたアクセスログデータを用いて、前回と同様に記事のレコメンドにチャレンジしてみようと思います。FactoRizationMachinesパッケージという便利そうなパッケージの存在も知れたことから、今回は以前から気になっていたFactorization Machineを扱います。

【目次】
・Factorization Machine(FM)の概要
・パッケージ紹介とインストール
・サンプルデータの構造把握と前処理
・FMの実行
・結果
・参考文献

Factorization Machine(FM)の概要

組み合わせ特徴量を扱う教師あり学習モデル。行列分解とSVMを合体させた手法。
スパースになりやすいデータの予測問題で扱う。
1ユーザーのある商品に対しての評価を、1評価1行として表して、ユーザーとアイテムの交互作用の特徴ベクトルを扱う。
相互作用項に関して、時間や文脈などを自由に入れられる。
相互作用項を次元圧縮する際の要素数を事前に決める必要がある。
Matrix Factorizationよりも精度が良いとされている。特徴量エンジニアリングなどで使われているようです。（Click-Through Rate Prediction）

パッケージ紹介とインストール

FactoRizationMachinesパッケージは線形SVMと2次のFMと高次のFMを実行することができ、引数で正則化項も加えることができます。現段階においては回帰のみで分類問題への適用は今後の開発となるようです。CRANから普通に
install.packages(‘FactoRizationMachines’)
でインストールします。libFMexeパッケージの場合は、libFMをインストールしてパスを指定しておく必要がありますが、このパッケージに関しては不要となります。

サンプルデータの構造把握と前処理

FactoRizationMachinesパッケージのサンプルコードにおいては、MovieLensのデータがサンプルデータとして載せられていました。ユーザーのID(整数)、映画のID(整数)、評価(整数、5段階)、日時(整数)からなるデータに対して、sparseMatrixに変換していました。

今回は、前回の投稿非負値行列因子分解(NMF)でブログ記事のレコメンドをしてみると同じデータを使って、アクセスログデータに適用しようと思います。FactoRizationMachinesの形式に合わせるために、このブログのアクセスログも、クッキーのIDを整数に、記事のIDを整数に、閲覧回数を5段階(5以上を5に変換)に、日時を整数に変更しています。

library(RGA)
library(tidyverse)
library(reshape2)

authorize()
prof <-list_profiles()

start_date <- "2017-07-04"
end_date   <- "2017-10-07"

#demention until 7
accesslogdata <- get_ga(profileId = prof$id[2],
                        start.date = start_date,
                        end.date = end_date,
                        dimensions = "ga:pagePath,
                        ga:dimension1,
                        ga:dateHourMinute,
                        ga:deviceCategory,
                        ga:userType,
                        ga:fullReferrer",
                        sort = "-ga:sessions",
                        metrics = "ga:sessions,ga:goal1Completions,ga:pageviews",
                        fetch.by = "day")

accesslogdata$pagePath <- vapply(strsplit(accesslogdata$pagePath,"\\?"), `[`, 1, FUN.VALUE=character(1))
accesslogdata <- accesslogdata %>% filter(grepl(x = pagePath,"/archives/[0-9]+$"))
accesslogdata$pagePath <- gsub(accesslogdata$pagePath,pattern = "/archives/",replacement = "article_")

#ユーザーIDの整数化
user <- unique(accesslogdata$dimension1)
user <- data.frame(dimension1=user) %>% mutate(user_number=1:n())

#記事IDの整数化
articles <- unique(accesslogdata$pagePath)
articles <- data.frame(pagePath=articles) %>% mutate(article_number=1:n())

#データの結合
accesslogdata <- accesslogdata %>% mutate(date=as.numeric(as.POSIXlt(as.Date(format(substr(accesslogdata$dateHourMinute,start = 1,stop = 12), format="%Y%m%d%"),format = "%Y%m%d"))))
accesslogdata <- accesslogdata %>% left_join(articles,by="pagePath")
accesslogdata <- accesslogdata %>% left_join(user,by="dimension1")

#最終接触日の抽出
last_date <- accesslogdata %>% select(user_number,article_number,date) %>% 
                arrange(user_number,desc(date))
last_date <- last_date[!duplicated(last_date[c("user_number","article_number")]),]

#各ユーザーの各記事に対するページビュー数の集計
accesslogdata_sum <- accesslogdata %>% group_by(article_number,user_number) %>% summarise(pageview=n()) %>% 
                        select(user_number,article_number,pageview)
accesslogdata_sum <- accesslogdata_sum %>% mutate(pageview=ifelse(pageview < 5,pageview,5))

accesslogdata_sum <- accesslogdata_sum %>% left_join(last_date, by = c("user_number" = "user_number", "article_number" = "article_number"))

library(RGA)

library(tidyverse)

library(reshape2)

authorize()

prof <-list_profiles()

start_date <- "2017-07-04"

end_date <- "2017-10-07"

#demention until 7

accesslogdata <- get_ga(profileId = prof$id[2],

start.date = start_date,

end.date = end_date,

dimensions = "ga:pagePath,

ga:dimension1,

ga:dateHourMinute,

ga:deviceCategory,

ga:userType,

ga:fullReferrer",

sort = "-ga:sessions",

metrics = "ga:sessions,ga:goal1Completions,ga:pageviews",

fetch.by = "day")

accesslogdata$pagePath <- vapply(strsplit(accesslogdata$pagePath,"\\?"), `[`, 1, FUN.VALUE=character(1))

accesslogdata <- accesslogdata %>% filter(grepl(x = pagePath,"/archives/[0-9]+$"))

accesslogdata$pagePath <- gsub(accesslogdata$pagePath,pattern = "/archives/",replacement = "article_")

#ユーザーIDの整数化

user <- unique(accesslogdata$dimension1)

user <- data.frame(dimension1=user) %>% mutate(user_number=1:n())

#記事IDの整数化

articles <- unique(accesslogdata$pagePath)

articles <- data.frame(pagePath=articles) %>% mutate(article_number=1:n())

#データの結合

accesslogdata <- accesslogdata %>% mutate(date=as.numeric(as.POSIXlt(as.Date(format(substr(accesslogdata$dateHourMinute,start = 1,stop = 12), format="%Y%m%d%"),format = "%Y%m%d"))))

accesslogdata <- accesslogdata %>% left_join(articles,by="pagePath")

accesslogdata <- accesslogdata %>% left_join(user,by="dimension1")

#最終接触日の抽出

last_date <- accesslogdata %>% select(user_number,article_number,date) %>%

arrange(user_number,desc(date))

last_date <- last_date[!duplicated(last_date[c("user_number","article_number")]),]

#各ユーザーの各記事に対するページビュー数の集計

accesslogdata_sum <- accesslogdata %>% group_by(article_number,user_number) %>% summarise(pageview=n()) %>%

select(user_number,article_number,pageview)

accesslogdata_sum <- accesslogdata_sum %>% mutate(pageview=ifelse(pageview < 5,pageview,5))

accesslogdata_sum <- accesslogdata_sum %>% left_join(last_date, by = c("user_number" = "user_number", "article_number" = "article_number"))

FMの実行

デフォルトの設定c(1, 10)では線形のウェイトが有効で、2次の項の要素数が10で正則化項なしのFMを実行することになります。引数に関する詳しい情報はPackage ‘FactoRizationMachines’に書かれています。今回はサンプルを参考に正則化項ありでモデルを実行します。まず、アクセスログデータに対して、ユーザーのIDからなる整数ベクトル、記事のIDからなる整数ベクトル、セッションのあった日時のデータからなる整数ベクトルを作成し、sparseMatrix関数を用いて元データを変形し、80%のデータをトレーニングに、20%のデータをテストに割り当てます。さらに、テストデータに関して、予測値との平均二乗誤差を計算します。

library(FactoRizationMachines)
library(Matrix)

user = accesslogdata_sum$user_number
items = accesslogdata_sum$article_number + max(user)
wdays = (as.POSIXlt(accesslogdata_sum$date,origin="2017-07-04")$wday+1)+max(items)

# Transform access log to feature form
data = sparseMatrix(i=rep(1:nrow(accesslogdata_sum),3),j=c(user,items,wdays),giveCsparse=F)
target = accesslogdata_sum$pageview

# Subset data to training and test data
set.seed(123)
subset = sample.int(nrow(data),nrow(data)*.8)
data.train = data[subset,]
data.test = data[-subset,]
target.train = target[subset]
target.test = target[-subset]

# Predict ratings with Support Vector Machine with linear kernel
model = SVM.train(data.train,target.train)

# RMSE resulting from test data prediction
sqrt(mean((predict(model,data.test)-target.test)^2))

# Predict ratings with second-order Factorization Machine
# with second-order 10 factors (default) and regularization
model = FM.train(data.train,target.train,regular=0.1)
# RMSE resulting from test data prediction
sqrt(mean((predict(model,data.test)-target.test)^2))

library(FactoRizationMachines)

library(Matrix)

user = accesslogdata_sum$user_number

items = accesslogdata_sum$article_number + max(user)

wdays = (as.POSIXlt(accesslogdata_sum$date,origin="2017-07-04")$wday+1)+max(items)

# Transform access log to feature form

data = sparseMatrix(i=rep(1:nrow(accesslogdata_sum),3),j=c(user,items,wdays),giveCsparse=F)

target = accesslogdata_sum$pageview

# Subset data to training and test data

set.seed(123)

subset = sample.int(nrow(data),nrow(data)*.8)

data.train = data[subset,]

data.test = data[-subset,]

target.train = target[subset]

target.test = target[-subset]

# Predict ratings with Support Vector Machine with linear kernel

model = SVM.train(data.train,target.train)

# RMSE resulting from test data prediction

sqrt(mean((predict(model,data.test)-target.test)^2))

# Predict ratings with second-order Factorization Machine

# with second-order 10 factors (default) and regularization

model = FM.train(data.train,target.train,regular=0.1)

# RMSE resulting from test data prediction

sqrt(mean((predict(model,data.test)-target.test)^2))

結果

各モデルについての平均二乗誤差を計算しています。
線形モデルや高次元モデルよりも、2次の項を持つFMが精度が高いようです。

> model = SVM.train(data.train,target.train)
> sqrt(mean((predict(model,data.test)-target.test)^2))
[1] 1.514266
> model = FM.train(data.train,target.train,regular=0.1)
> sqrt(mean((predict(model,data.test)-target.test)^2))
[1] 1.369983
> model　=　HoFM.train(data.train,target.train,c(1,3,1),regular=0.1)
> sqrt(mean((predict(model,data.test)-target.test)^2))
[1] 1.388247

> model = SVM.train(data.train,target.train)

> sqrt(mean((predict(model,data.test)-target.test)^2))

[1] 1.514266

> model = FM.train(data.train,target.train,regular=0.1)

> sqrt(mean((predict(model,data.test)-target.test)^2))

[1] 1.369983

> model　=　HoFM.train(data.train,target.train,c(1,3,1),regular=0.1)

> sqrt(mean((predict(model,data.test)-target.test)^2))

[1] 1.388247

こちらは、この中で性能の良かった2次の項を持つFMの予測結果とテストデータの結果をプロットしたものです。4点を超える値をあまり予測できていないようです。今回はサンプルを回しただけなので、本来であれば次元の数kや正則化のセッティングをいろいろいじったり、相互作用項を新しく追加するなどして精度を高めることが必要です。

結果の比較だけでは仕事で使えないので、実際に予測した結果を取り出したいと思います。
実際に運用するとなると、ページIDを所与として、ページビュー数を0とおいて（型をそろえるため。NULLだとエラーになった）、任意のタイミング(date)を想定して、モデルにデータを適用し、評価の高いものをサジェストするスタイルになるのではないでしょうか。

この結果だと、ユーザー98に記事1を見せることに対して4.02点が与えられています。

参考文献

Factorization Machinesを今更読みました
 Factorization Machines
High-order factorization machines with R #tokyor 61
Factorization Machinesのおはなし。
libFMexeを動かすまで (R Wrapper for the libFM Executable参照記事)
一歩Matrix Factorization、二歩Factorization Machines、三歩Field-aware Factorization Machines…『分解、三段突き！！』
[論文] Factorization Machines (ICDM 2010) 読んだ 22:41
Factorization machines with r
Factorization Machinesについて調べてみた

非負値行列因子分解(NMF)でブログ記事のレコメンドをしてみる

概要

2017年7月より当ブログにおいて、Google AnalyticsのカスタムディメンジョンにUser IDを取得して残せるようにしました。これにより、ユーザーの一人一人が閲覧した記事の分析を行うことが可能となっています。そこで今回は、記事の数も50近くなってきたことから、レコメンド手法の適用を行ってみたいと思います。Exampleデータでレコメンドをしても、ドメイン知識がないと、どれくらい良さそうな推薦なのかの判断が付きにくいと考えていたので、自分で運営しているサイトのデータであれば評価をつけやすいのではないかと思われます。
手法としてはNMF（non-Negative Matrix Factorization）を適用します。

【目次】
・前処理
・データ確認
・行列因子分解（Matrix Factorization）
・レコメンド結果の確認
・感想
・参考文献

前処理

ユーザー単位の閲覧した記事データを作るためにGAのAPIで取得したデータを整形し、集計します。このようなデータを手に入れるには自身のブログにおいて、GTM(Google Tag Maneger)を設置して、GTMのコンテナにカスタムJS変数などを格納する必要があります。詳しくはこちらのブログを参照ください。（ 2016年の新定番！ユーザーエクスプローラーをもっと活用するための簡単な方法）

#パッケージの読み込み
library(RGA)
library(tidyverse)

#GAの認証
authorize()
prof <-list_profiles()

#取得期日の指定
start_date <- "2017-07-04"
end_date   <- "2017-09-16"


#GAのAPIからデータ取得（ディメンションの数は7つまで）
accesslogdata <- get_ga(profileId = prof$id[2],
              start.date = start_date,
              end.date = end_date,
              dimensions = "ga:pagePath,
                            ga:dimension1,
                            ga:dateHourMinute,
                            ga:deviceCategory,
                            ga:userType,
                            ga:fullReferrer",
              sort = "-ga:sessions",
              metrics = "ga:sessions,ga:goal1Completions,ga:pageviews",
              fetch.by = "day")

#記事のURLだけに絞り込む
accesslogdata$pagePath <- vapply(strsplit(accesslogdata$pagePath,"\\?"), `[`, 1, FUN.VALUE=character(1))
accesslogdata <- accesslogdata %>% filter(grepl(x = pagePath,"/archives/[0-9]+$"))

#ユーザーごとに重複した記事の情報を削除する
df_dups <- accesslogdata[c("pagePath","dimension1")]
accesslogdata <- accesslogdata[!duplicated(df_dups),]
accesslogdata$pagePath <- gsub(accesslogdata$pagePath,pattern = "/archives/",replacement = "article_")

articles <- unique(accesslogdata$pagePath)

#記事ごとの閲覧有無に関するカラムを生成する
for(i in 1:length(articles)){
  accesslogdata <- accesslogdata %>% mutate(variable=if_else(pagePath==articles[i],1,0))
  accesslogdata <- eval(parse(text=paste0("accesslogdata %>% rename( ",articles[i],"=variable)")))
}

#クッキー単位で接触した記事ごとに集計する
accesslogdata <- accesslogdata %>% select(-pagePath,-dateHourMinute,-deviceCategory,-userType,-fullReferrer,-sessions,-goal1Completions,-dateHour,-pageviews) %>% 
                    group_by(dimension1) %>% summarise_all(funs(sum))

#パッケージの読み込み

library(RGA)

library(tidyverse)

#GAの認証

authorize()

prof <-list_profiles()

#取得期日の指定

start_date <- "2017-07-04"

end_date <- "2017-09-16"

#GAのAPIからデータ取得（ディメンションの数は7つまで）

accesslogdata <- get_ga(profileId = prof$id[2],

start.date = start_date,

end.date = end_date,

dimensions = "ga:pagePath,

ga:dimension1,

ga:dateHourMinute,

ga:deviceCategory,

ga:userType,

ga:fullReferrer",

sort = "-ga:sessions",

metrics = "ga:sessions,ga:goal1Completions,ga:pageviews",

fetch.by = "day")

#記事のURLだけに絞り込む

accesslogdata$pagePath <- vapply(strsplit(accesslogdata$pagePath,"\\?"), `[`, 1, FUN.VALUE=character(1))

accesslogdata <- accesslogdata %>% filter(grepl(x = pagePath,"/archives/[0-9]+$"))

#ユーザーごとに重複した記事の情報を削除する

df_dups <- accesslogdata[c("pagePath","dimension1")]

accesslogdata <- accesslogdata[!duplicated(df_dups),]

accesslogdata$pagePath <- gsub(accesslogdata$pagePath,pattern = "/archives/",replacement = "article_")

articles <- unique(accesslogdata$pagePath)

#記事ごとの閲覧有無に関するカラムを生成する

for(i in 1:length(articles)){

accesslogdata <- accesslogdata %>% mutate(variable=if_else(pagePath==articles[i],1,0))

accesslogdata <- eval(parse(text=paste0("accesslogdata %>% rename( ",articles[i],"=variable)")))

}

#クッキー単位で接触した記事ごとに集計する

accesslogdata <- accesslogdata %>% select(-pagePath,-dateHourMinute,-deviceCategory,-userType,-fullReferrer,-sessions,-goal1Completions,-dateHour,-pageviews) %>%

group_by(dimension1) %>% summarise_all(funs(sum))

このような形で行にユーザーが、列に記事が、各要素にアクセスの有無が入るようになっています。本来は評価データを用いますが、評価データがないので、アクセスしたかどうかの1-0データとして扱います。

データ確認

#ユーザーごとの接触ページ数を集計し、ヒストグラムで描く
session <- data.frame(PageView=rowSums(accesslogdata[,-1]))
ggplot(data = session, aes(x = PageView)) + geom_histogram(binwidth = 1) +
    ylab("User Count") + ggtitle("Article Page View Each Users") + 
    theme(plot.title = element_text(hjust = 0.5))

#ユーザーごとの接触ページ数を集計し、ヒストグラムで描く

session <- data.frame(PageView=rowSums(accesslogdata[,-1]))

ggplot(data = session, aes(x = PageView)) + geom_histogram(binwidth = 1) +

ylab("User Count") + ggtitle("Article Page View Each Users") +

theme(plot.title = element_text(hjust = 0.5))

以下より、行列の要素のうち、値が入っている要素は全体の2.9%であることがわかります。参考文献は2.6%をスパースと言っていたことから、今回のデータもスパースであると考えることができます。

> X <- as(accesslogdata[,-1],"matrix")
> table(X)
X
     0      1 
122741   3777 
> table(X)[2]/(table(X)[1]+table(X)[2])
         1 
0.02985346

> X <- as(accesslogdata[,-1],"matrix")

> table(X)

0 1

122741 3777

> table(X)[2]/(table(X)[1]+table(X)[2])

0.02985346

行列因子分解（Matrix Factorization）

行列因子分解はユーザーごとのアイテム購買履歴などのデータからなる行列の因子分解を行うもので、レコメンドにおいては因子分解の結果として、ゼロ以上の要素を持つ場合に推薦対象とするなどの活用がされています。非負値行列因子分解は要素が非負となるように制約を付けて行列を分解するもので、正の値を取ることから推薦に使いやすいとされています。

今回はNMFパッケージを用いて、ほぼデフォルトな設定で結果を眺めてみますが。本来であれば以下のような考慮が必要となります。最適化などをがっつり進めようとするとクロスバリデーションなどを行ったり、計算リソースが結構かかったりすると思われます。

アルゴリズムの選択
- brunet,KL,lee,Frobenius,offset,nsNMF,ls-nmf,pe-nmf,siNMF,snmf/r,snmf/lがある。
- デフォルトではbrunetとなっている。
シーディングメソッドの選択（初期値の決め方）
- 固定値(none),ランダム(random),ica,nndsvdがある。
- デフォルトではランダムとなっている。
実行回数の設定
- １００〜200回を選ぶことが多い。
並列処理
- 処理量が多いので最適化などを考えると並列処理も重要。
因子分解ランクの推定
- ランクの数が多いと残差は減るが、過学習になりやすい。クロスバリデーションを行う必要がある。

library(NMF)

#NMFの実行（6次元以上の次元を削除）
res <- nmf(accesslogdata[,-1], 5, seed = 1234, .options = "t")

#試行回数と損失関数としてのカルバックライブラー情報量の可視化
plot(res)

#基底行列の取得
w <- basis(res)

#係数行列の取得
h <- coef(res)

#レコメンデーション評価値を得る
accesslogdata.nmf <- w %*% h

#確認するユーザーの行番号
userid <- c(377,987,2584)

#RAWデータとレコメンデーション結果の比較
recommend_result <- data.frame(round(t(accesslogdata.nmf[userid,]), 2),
           t(accesslogdata[userid,-1]))

#記事タイトルの読み込み
article_title <- read_csv(file = "article_titles.csv")

recommend_result <- left_join(rownames_to_column(recommend_result), article_title, by = c("rowname" = "article_id")) %>% 
                        mutate(rowname=title) %>% select(-title)

library(NMF)

#NMFの実行（6次元以上の次元を削除）

res <- nmf(accesslogdata[,-1], 5, seed = 1234, .options = "t")

#試行回数と損失関数としてのカルバックライブラー情報量の可視化

plot(res)

#基底行列の取得

w <- basis(res)

#係数行列の取得

h <- coef(res)

#レコメンデーション評価値を得る

accesslogdata.nmf <- w %*% h

#確認するユーザーの行番号

userid <- c(377,987,2584)

#RAWデータとレコメンデーション結果の比較

recommend_result <- data.frame(round(t(accesslogdata.nmf[userid,]), 2),

t(accesslogdata[userid,-1]))

#記事タイトルの読み込み

article_title <- read_csv(file = "article_titles.csv")

recommend_result <- left_join(rownames_to_column(recommend_result), article_title, by = c("rowname" = "article_id")) %>%

mutate(rowname=title) %>% select(-title)

レコメンド結果の確認

各ユーザーに対して推薦されたアイテムのうち、比較的値の大きなものに関してコメントをしたいと思います。

ユーザー1
「XGBoostやパラメータチューニングの仕方に関する調査」を閲覧したユーザー1に対して、「XGBoostのパラメータチューニング実践 with Python」を推薦できている一方で、「OpenCV&Pythonで画像の類似度を計算させる〜イケメンの顔比較」など関係性の薄そうな記事も推薦してしまっています。

ユーザー2
「人工知能学会全国大会2017のWebマーケティングで参考になりそうな研究9選」を閲覧したユーザー2に対して、「Kaggleで使われた特徴量エンジニアリングとアルゴリズムまとめ」を推薦できている。

ユーザー3
「洋楽の歌詞データでDoc2vecを実行してみる」を閲覧したユーザー3に対して、「LDA（潜在的ディリクレ配分法）まとめ手法の概要と試行まで」を推薦できている一方で、「ベイジアンネットワークをRのbnlearnパッケージで推定して予測してみる」と関係性の薄そうな記事も推薦している。本来、近いと思われる、「Word2Vecでクラシックの楽曲情報をコーパスとして類似度を出してみる」を推薦できていない。

感想

結果に関して、納得のいくものといかないもので半々な印象でデフォルトの設定では満足できない。
NMFはクロスバリデーションなどにより過学習を避けた形での最適なランク数やアルゴリズムの選定、初期値の設定などが必要と思われるので、そこらへんの取り組みも今後行いたい。
FM(Factorization Machine)は回帰項を追加できるので、今後挑戦してみたい。
1-0のデータではなく評価のデータが欲しい。

参考文献

岩波データサイエンス Vol.5
非負値行列因子分解（NMF）によるレコメンドのちょっとした例
 An introduction to NMF package
非負ッ値行列・ファクトリゼーション～半群じゃないから難しくないもん！

人工知能学会全国大会2017のWebマーケティングで参考になりそうな研究9選

はじめに

今更ではありますが、2017年5月に開かれた人工知能学会全国大会2017の公開資料の中から、私の本業であるWebマーケティングで役に立ちそうな研究を独断と偏見で9本ほど集めてみました。思っていたよりもWebマーケティングぽい領域が多かったので、社内での分析業務におけるアイデアに直結しそうな気がしています。

ちなみに、全ての資料はこちらにあります。

基本的に各研究ごとに
・目的と結果
・対象となるデータ
・手法の概要
・PDFのリンク
について記していきます。

紹介する研究

今回は以下の研究を紹介したいと思います。
「オンライン広告におけるスパース性と遅れコンバージョンを考慮した予測モデル」「テキストマイニングを用いた転職サイトの会員離脱予測」「インターネット広告におけるスパースなユーザー行動ベクトルからのユーザー特徴抽出モデル」「ユーザーのページビュー系列からのコンバージョン予測」「SNSにおけるユーザの行動のモチベーション要因の分析」「状況に応じた楽曲推薦に向けたソーシャルデータ分析」「テキストマイニングを用いた口コミ分析による点数評価の信頼性確認手法」「不動産仲介マーケティングのためのユーザ行動予測」「SUUMOでの不動産データ活用の取り組みと未来」

オンライン広告におけるスパース性と遅れコンバージョンを考慮した予測モデル

目的と結果

オンライン広告のCVR予測の精度向上を目的としている。
これまでの課題としては、
「質的変数による特徴量のスパース性」
「広告クリックからのCVまでの期間の長さ（期間が空いてからCVしたりするケースがあるので正例だったはずが負例と扱われる可能性がある）」
などがあった。提案手法により従来手法に比べ高い精度でCVR予測が可能であることが示された。

対象となるデータ

過去3週間分の広告配信データでテスト期間は7日間

手法の概要

次元圧縮
- Factorization Machine(FM)を用いて、従来のロジスティック回帰では考慮できていなかった交互作用をモデリングし、低次元のベクトル内積に分解して表現している。
遅れCV問題の対応
- 遅れてCVが発生する事象に関しては、指数分布を用いて遅れてCVする事象をモデリングするDelayed Feedback Modelを構築することで、本当は事後的にCVしそうなユーザーの推定を行っている。
予測
- これらの手法を合わせた、Delayed Feedback Factorization Machineを構築し、SGD(確率的勾配降下法)やAdaGradを用いてパラメータを推定している。

PDFのリンク

https://kaigi.org/jsai/webprogram/2017/pdf/110.pdf

テキストマイニングを用いた転職サイトの会員離脱予測

目的と結果

登録情報や履歴書のテキスト情報を用いて特徴量とし、転職サイトの登録会員の離脱予測を行っている。ブラックボックスでは無い形での結果の解釈ができ、予測精度も高いことが示された。

対象となるデータ

2009年10月1日〜2016年9月28日までの履歴書、職務経歴書、自己PR書、希望条件
離脱会員の定義は「登録して1ヶ月以内に一度も応募しない会員」としている。

手法の概要

STEP1：Pythonのmojimojiというライブラリを利用してカナを全角に、数字とアルファベットは半角に変換。
STEP2：MeCabを用いて名詞のみ抽出
STEP3：「職務」「経歴」「業務」「内容」などを不用語（stop word）とした。
STEP4：Grid Searchでパラメータチューニングしたランダムフォレストを用いて変数重要度（Feature Importance）を出し、重要度の上位90%を占める単語以外を削除
STEP5：登録情報のデータ（年齢、転職回数、住所）を数値データに変換
STEP6：選定した特徴量を用いて決定木アルゴリズムを適用

PDFのリンク

https://kaigi.org/jsai/webprogram/2017/pdf/949.pdf

インターネット広告におけるスパースなユーザー行動ベクトルからのユーザー特徴抽出モデル

目的と結果

広告のCTR予測を行う際の特徴量が非常にスパースであることから、ユーザーをクラスタリングして広告配信に効果的なユーザーを抽出することを目的としている。研究ではCVRの高いユーザー群の抽出を可能にしている。

対象となるデータ

2016年6月1日〜2016年6月21日の間のクリックやコンバージョンのデータ（株式会社アイモバイルが提供）
クリック実績のある広告数、コンバージョン実績のある広告数、クリックされた広告番号、コンバージョンが発生した広告番号、ユーザー番号など

手法の概要

高次元スパースデータを低次元に落とし込むために、オートエンコーダを用いている。
オートエンコーダにより得られた特徴量をもとにユーザーのクラスタリングを行っている。(オートエンコーダは入力された情報をなるべく欠損しないように、効率的に圧縮する構造を持たせたニューラルネットワークのこと。)
データにおけるゼロの成分を確率的に欠損とみなして計算処理を早めている。
学習モデルはAdaGrad(学習率を学習に合わせて最適化する手法)
訓練データ（特徴量の数がクリックで23339個、コンバージョンで5619個）にスパースデータ許容型オートエンコーダで500次元まで圧縮し、さらに通常のオートエンコーダで146次元まで圧縮している。
次元圧縮して獲得したデータに対してk-means法でクラスタリングを行い、CVRなどを比較している。

PDFのリンク

https://kaigi.org/jsai/webprogram/2017/pdf/12.pdf

ユーザーのページビュー系列からのコンバージョン予測

目的と結果

Web広告において、ユーザーの閲覧履歴からコンバージョンを予測することを目的に、Recurrent Neural Network(RNN)を用いた結果、非時系列のSVMよりもわずかばかりかの精度向上が観察された。データ数が少ない場合はSVMに軍配があがる結果となっている。

対象となるデータ

2016年7月1日〜2016年10月31日の広告の接触履歴データなど（広告の業種別データやユーザーのページビュー数）を特徴量に用いて、2016年11月1日〜11月30日までの期間を予測対象とした。データは株式会社マイクロアドが提供している模様。

手法の概要

Long Short-Term Memory RNN(LSTM)
- ソフトアテンションモデル
  - 可変長の系列に対して内積で重みを算出し、足し合わせることで一部から情報を取り出すモデル。途中の層まで複数業種のデータを用いることで、複数の業種などの情報を共有することができる。（最後の隠れ層には業種別の特徴量を用いている。）
  - 勾配の計算に関してはAdamを用いて、実装はKerasを用いている。活性化関数はTanh、ドロップアウトはGaussianDropout
    としている。
- RNNで用いる特徴量としては、特定のWebページを640種類に分類して、1日のうちに各トピックにページビューがあったかなかったかを1-0の変数にするなどの前処理を行っている。

PDFのリンク

https://kaigi.org/jsai/webprogram/2017/pdf/765.pdf

SNSにおけるユーザの行動のモチベーション要因の分析

目的と結果

SNSユーザーのモチベーションの要因（一ヶ月先に投稿するかどうか）をユーザー同士の関係性（staticな関係）と、そのユーザーの周りのインタラクション（dynamicな関係）という観点から分析し、ユーザー同士の関係性とユーザー周りのインタラクションを同時に用いたモデルにおいて、それらがモチベーションの要因につながりうる傾向が観察された。

対象となるデータ

RoomClipという部屋の家具やインテリアの様子の写真を投稿するSNSのデータ
- 気に入った写真に「いいね」を付けたり、コメントを送れる。
2015年3月の一ヶ月間を学習期間とし、2015年4月の投稿の有無を予測対象とする。

手法の概要

ネットワーク特徴量
- 各ユーザーに対して、
  - 中心性（次数中心性、媒介中心性、近接中心性）
  - クラスタ係数、PageRank、フォロー数とフォロワー数の比率、フォローしているユーザのアクティブ度
フィードバック特徴量
- 特定期間中にユーザが周りのユーザから自身の投稿に対してもらった「いいね」「コメント」に対して
  - 反応の数
  - 反応を行ったユーザの数
  - 反応を行ったユーザのうち、相互フォロー関係にあるユーザからのものの割合
予測器としてSVMを利用

PDFのリンク

https://kaigi.org/jsai/webprogram/2017/pdf/534.pdf

状況に応じた楽曲推薦に向けたソーシャルデータ分析

目的と結果

ユーザへの楽曲推薦システムの構築を目的に、楽曲とその再生時の状況の関係を収集したソーシャルデータより分析している。取得したデータのうちの7%ほどではあるが、ツィート情報から楽曲と再生時の状況との関係を抽出可能であるということが示された。

対象となるデータ

2016年8月10日〜2016年9月3日の間に投稿されたツィート（66879件）

手法の概要

ハッシュタグ#nowplaying付きのツィートを収集して、アーティスト名や楽曲名や楽曲再生時の状況（同一ユーザの5分前後の投稿内容）を収集
取得したツィートからアーティスト名、楽曲名を除去して再生時の状況を取得
取得したデータからデータベースを作成

PDFのリンク

https://kaigi.org/jsai/webprogram/2017/pdf/102.pdf

テキストマイニングを用いた口コミ分析による点数評価の信頼性確認手法

目的と結果

製品などの口コミの評価レビューに関する信頼性を確かめることを目的としている。レビューデータに対して感性評価を行い、さらにその結果に対して総合評価点数別に感性表現（Goodなど）の伴った表現の頻出ランキングを作成することで、点数が高い理由を把握しやすくなった。

対象となるデータ

ソニーの製品サイトにおける、ユーザの製品に対する評価とコメントのデータ（1406件）

手法の概要

テキストマイニングスタジオによる感性評価（ポジティブ、中立、ネガティブ）
総合評価点数別（5段階）の「ポジティブ、中立、ネガティブ」の割合の調査
総合評価4の口コミかつ中立な表現にフィルターをかけて、感性情報を有する単語と係り受け関係となった単語の頻度ランキングのTOP10を抽出

PDFのリンク

https://kaigi.org/jsai/webprogram/2017/pdf/1115.pdf

不動産仲介マーケティングのためのユーザ行動予測　Prediction of user behaviors for real estate brokerage marketing

目的と結果

不動産紹介業において、顧客の属性情報と来店初期の行動から、来店者が最終的に商品を購買するかどうかの予測モデル(不均衡データにおける予測タスク)を構築することでマーケティングの広告効果測定に適用することを目的としている。結果として、顧客の来店後30分間という短い時間における行動だけでも識別力のある予測が可能であることが示されている。

対象となるデータ

オンライン不動産賃貸仲介店舗iettyの顧客データベース
- 希望物件の情報
  - 引っ越し希望日
  - 住みたいエリアの駅
  - 間取り
  - 広さ
  - 家賃
- 顧客属性
  - 在住の地域
  - 性別
  - 年齢
- 行動データ
  - 登録時
    - フォームの記入にかかった時間
    - フリーコメントの長さ
    - 利用デバイス
  - 登録後
    - 初回チャット発言までの経過時間
    - 初回物件評価までの経過時間

手法の概要

STEP1：顧客の属性および登録初期(30分以内)の行動から特徴量を作成し、2ヶ月以内に契約すれば正例、それ以外を負例とする。（正例は619件、負例は33212件）
STEP2：検証用データを無作為に20%取り置く
STEP3：残り80%のデータに対して負例に対する、ランダムアンダーサンプリングと、正例に対するSMOTE(Synthetic Minority Over-sampling Technique)を用いたオーバーサンプリングによりサンプル数を調整する。
STEP4：ランダムフォレストを用いて、無効な特徴量の取捨選択や離散化のレンジ調整する。
STEP5：全特徴量（83次元）を二値化して、RBFカーネル(Gaussian カーネル)のSVMを適用し、様々なCやγの組み合わせを試行して、二値分類を行う。
STEP6：評価を予測確率の平均とAUCの2種類で行う。

PDFのリンク

https://kaigi.org/jsai/webprogram/2017/pdf/466.pdf

SUUMOでの不動産データ活用の取り組みと未来

目的と結果

物件検索ユーザーが比較検討を行いやすくすることを目的に様々な取り組みをしている。

対象となるデータ

不動産ポータルサイトSUUMOのデータ
- 行動履歴(Webサイト、アプリ)
- 物件情報(不動産情報)
- 街情報(独自で収集)
- アンケート

手法の概要

実活用
- 相場推定
- 単純な加算平均ではなく、外れ値や時系列変動、データ量などを考慮して推定している。
- リアルタイム物件レコメンド
- ランキング学習により推定したモデルを利用し、リアルタイムに計算を行っている。個人ごとにパーソナライズしたランキングも用いている。
- バンディットアルゴリズムによる配信施策の最適化
- メールやプッシュ通知などで件名やデザイン、配信時間など様々な選択肢があり、その最適なパターンを見つけるために、バンディットアルゴリズムを適用している。
実験段階
- 間取り画像の分析
- ユーザーの求めている間取りが何なのかを画像から抽出しようとしている。
- ユーザーの検索行動、嗜好の推定
- ユーザーの検索行動から、嗜好性を掴み取り、レコメンドする際の理由付けができるようにモデルを作ろうとしている。

PDFのリンク

https://kaigi.org/jsai/webprogram/2017/pdf/743.pdf

感想

特徴量の削減のためにランダムフォレストを適用する流れ
高次元スパースデータに対するニューラルネットワークの適用
CVの時系列予測のためにRNNの適用
SNS分析におけるネットワーク特徴量とフィードバック特徴量という概念
口コミデータ分析の作法(点数別に感性表現の係り受けを見る)
など勉強になりました。Kaggleと違って日本語なのが読みやすかったですね。

参考文献

初めてのディープラーニング –オープンソース”Caffe”による演習付き
 SMOTE で不均衡データの分類

Kaggleで使われた特徴量エンジニアリングとアルゴリズムまとめ

Kaggleにおいて、人によってはソリューションを書いてくれているものがあります。特徴量エンジニアリングやアルゴリズムなど業務での参考になるだろうと思われるので、仕事で関わりそうなデータという観点から4つほどですが事例を集めてみました。このような事例のまとまった本とかあったらいいのですが。

基本的に各事例ごとに
・目的
・どんなデータセット
・どんな特徴量エンジニアリング
・どんなアルゴリズム
・リンク
を記していきます。

Walmart Recruiting – Store Sales Forecasting

目的

ウォルマートの部門ごとの売上予測

どんなデータセット

外部のデータ
・気温
・物価指数
・スーパーボウルの時期
など

どんな特徴量エンジニアリング

を行った程度のようです。

どんなアルゴリズム

機械学習の手法
3. Random Forest
4. Linear Regression
5. K nearest regression
6. Principle Component Regression

以上の6つのモデルから平均値をとって予測し、部門ごとにモデルを作ったようです。

リンク

6 bad models make 1 good model: Power of Ensemble Learning

Algorithmic Trading Challenge

目的

大規模なトレーディングにおけるマーケットの反応を予測

どんなデータセット

どんな特徴量エンジニアリング

どんなアルゴリズム

リンク

Winning the Kaggle Algorithmic Trading Challenge with the Composition of Many Models and Feature Engineering

Predict Closed Questions on Stack Overflow

目的

質問が閉じられるかどうかを予測

どんなデータセット

どんな特徴量エンジニアリング

どんなアルゴリズム

Vowpal WabbitというMicrosoftやYahoo!が出資している機械学習ツールを使い、多クラス分類を行ったようです。カテゴリカルな変数をone-hotに変換する必要がなく、楽に分析できるようです。

リンク

Predicting closed questions on Stack Overflow

vowpal_wabbit

Basic Multiclass Classification

Outbrain Click Prediction

目的

ユーザーにレコメンドするコンテンツのクリック予測

どんなデータセット

どんな特徴量エンジニアリング

どんなアルゴリズム

リンク

2nd place solution | team brain-afk

感想

調べてみて、複数の時系列モデルの予測結果の平均値で予測する手法、特徴量の選択をstepAICのようにランダムフォレストで行う手法、Vowpal Wabbit（今回の2つのソリューションで扱われていた）などを新たに知れました。Kaggleには他にもKernelという手法をシェアする場があるので、その情報も今後キャッチアップしていきたいと思います。

参考情報

機械学習コンペティションの進展と今後の展開

ベイジアンネットワークをRのbnlearnパッケージで推定して予測してみる

ベイジアンネットワークの知見が無かったので、調べた情報をまとめています。一応、載せているスクリプトでRを用いて予測するということができます。

【目次】
・ベイジアンネットワークとは
・ベイジアンネットワークの用途
・ベイジアンネットワークの推定のステップ
・Rでの実行例
・おまけ：Webサービスへの応用例
・参考文献

ベイジアンネットワークとは

・複数の確率変数の間の定性的な依存関係をグラフ構造によって表し、個々の変数の間の定量的な関係を条件付確率で表した確率モデル。前提として、有向非循環（グルグルと回らないグラフ）となっているグラフの構造を持つものに限定している。
・入力となる変数と出力となる変数はモデルの中では区別されない。
・時間という明確な因果関係などをモデルに組み込みやすいので、系列データなどを扱うケースが多い。モデル設計者がデータが生成されるプロセスを考慮しながらモデルを組んでいける。
・非循環性とd分離の仮定のみによって導かれる現在考えられる最も自然な離散モデルであり、現在の様々なモデルの中でも、最も表現力と予測力を持つモデルとされる。

ベイジアンネットワークの用途

・物体追跡
・ジェスチャ認識
・Webサイトなどのレコメンデーションサービス
・広告配信
・メルマガなどの配信最適化
など画像処理問題から消費者行動問題までいろいろな分野で活用されているようです。

ベイジアンネットワークの推定のステップ

モデルの選択や推定は学習 (learning)と呼ばれ、以下の2段階のステップを踏みます。

Step1：構造学習(structure learning)…データからネットワーク構造を学習する。
Step2：パラメータ学習(parameter learning)…Step1で学習した構造によって意味付けした分布のパラメータを学習する。

ベイジアンネットワークの推定は条件付き確率やMAP(Maximum a posteriori)推定量などを用います。構造学習の際に、事前知識を導入することができます。

Rでの実行例

・データ
何千もの主要な免疫系細胞から取り出した11種類のリン酸タンパク質とリン脂質の測定値からなるデータです。以下のサイトよりダウンロードできます。
Supporting Online Material Causal Protein-Signaling Networks Derived from Multiparameter Single-Cell Data

・パッケージ
ベイジアンネットワークを実行できるパッケージはいろいろありますが、参考文献に従いbnlearnを用います。加えて、ネットワーク図の可視化のためにRgraphvizを用います。Rgraphvizがcranでは対応したものがなかったので、こちらを参考にして、インストールしました。
Provides plotting capabilities for R graph objects

・ネットワークの図示

#ネットワークの図示
library(Rgraphviz)
library(bnlearn)

spec = paste("[PKC][PKA|PKC][praf|PKC:PKA][pmek|PKC:PKA:praf]",
             "[p44.42|pmek:PKA][pakts473|p44.42:PKA][P38|PKC:PKA]",
             "[pjnk|PKC:PKA][plcg][PIP3|plcg][PIP2|plcg:PIP3]")

net = model2network(spec)
class(net)

graphviz.plot(net, shape = "ellipse")

#ネットワークの図示

library(Rgraphviz)

library(bnlearn)

spec = paste("[PKC][PKA|PKC][praf|PKC:PKA][pmek|PKC:PKA:praf]",

"[p44.42|pmek:PKA][pakts473|p44.42:PKA][P38|PKC:PKA]",

"[pjnk|PKC:PKA][plcg][PIP3|plcg][PIP2|plcg:PIP3]")

net = model2network(spec)

class(net)

graphviz.plot(net, shape = "ellipse")

・データ確認

#データ確認
> protein_dataset <-  read.csv(file = "protain_dataset.csv",as.is = TRUE)
> head(protein_dataset,5)
   praf pmek  plcg PIP2  PIP3 p44.42 pakts473  PKA   PKC  P38 pjnk
1 39.20  8.9  4.26 9.22  1.88   25.9     46.6  223  3.02 31.6 65.5
2 42.60 20.2 12.10 1.45 40.00   30.2     51.9  264  5.94 14.1 16.1
3  5.88 17.9  3.43 2.19 18.10  105.0    300.0 1928  6.98 20.5 12.9
4 34.30 18.8  3.65 1.70 16.50  352.0    519.0 1263 18.30 59.4 35.9
5 24.60 16.0  6.92 9.39  8.06   15.7     53.8  346  9.31 18.1 12.5

#データ確認

> protein_dataset <- read.csv(file = "protain_dataset.csv",as.is = TRUE)

> head(protein_dataset,5)

praf pmek plcg PIP2 PIP3 p44.42 pakts473 PKA PKC P38 pjnk

1 39.20 8.9 4.26 9.22 1.88 25.9 46.6 223 3.02 31.6 65.5

2 42.60 20.2 12.10 1.45 40.00 30.2 51.9 264 5.94 14.1 16.1

3 5.88 17.9 3.43 2.19 18.10 105.0 300.0 1928 6.98 20.5 12.9

4 34.30 18.8 3.65 1.70 16.50 352.0 519.0 1263 18.30 59.4 35.9

5 24.60 16.0 6.92 9.39 8.06 15.7 53.8 346 9.31 18.1 12.5

・構造学習
データセットから、ネットワーク構造を推定します。推定するための手法は様々あるようですが、今回はヒルクライムアルゴリズムを用います。

#構造学習の実行
> print(hc(protein_dataset))

  Bayesian network learned via Score-based methods

  model:
   [PIP2][p44.42][PKC][PIP3|PIP2][pakts473|p44.42][pjnk|PKC][plcg|PIP3]
   [PKA|p44.42:pakts473][P38|PKC:pjnk][pmek|plcg:P38][praf|pmek]
  nodes:                                 11 
  arcs:                                  11 
    undirected arcs:                     0 
    directed arcs:                       11 
  average markov blanket size:           2.18 
  average neighbourhood size:            2.00 
  average branching factor:              1.00 

  learning algorithm:                    Hill-Climbing 
  score:                                 BIC (Gauss.) 
  penalization coefficient:              3.348517 
  tests used in the learning procedure:  185 
  optimized:                             TRUE

#構造学習の実行

> print(hc(protein_dataset))

Bayesian network learned via Score-based methods

model:

[PIP2][p44.42][PKC][PIP3|PIP2][pakts473|p44.42][pjnk|PKC][plcg|PIP3]

[PKA|p44.42:pakts473][P38|PKC:pjnk][pmek|plcg:P38][praf|pmek]

nodes: 11

arcs: 11

undirected arcs: 0

directed arcs: 11

average markov blanket size: 2.18

average neighbourhood size: 2.00

average branching factor: 1.00

learning algorithm: Hill-Climbing

score: BIC (Gauss.)

penalization coefficient: 3.348517

tests used in the learning procedure: 185

optimized: TRUE

構造学習より得られたネットワーク構造を可視化します。因果関係がダメな感じのネットワークが出来ていないかチェックする必要がありますね。機械的に作るだけでなく、背景的知識も考慮するという進め方が推奨されています。

#推定したネットワークを可視化
spec.estimated = paste("[PIP2][p44.42][PKC][PIP3|PIP2][pakts473|p44.42][pjnk|PKC][plcg|PIP3]",
                        "[PKA|p44.42:pakts473][P38|PKC:pjnk][pmek|plcg:P38][praf|pmek]")
net.estimated = model2network(spec.estimated)
class(net.estimated)
graphviz.plot(net.estimated, shape = "ellipse")

#推定したネットワークを可視化

spec.estimated = paste("[PIP2][p44.42][PKC][PIP3|PIP2][pakts473|p44.42][pjnk|PKC][plcg|PIP3]",

"[PKA|p44.42:pakts473][P38|PKC:pjnk][pmek|plcg:P38][praf|pmek]")

net.estimated = model2network(spec.estimated)

class(net.estimated)

graphviz.plot(net.estimated, shape = "ellipse")

・パラメータ学習
先ほど、データセットから求めた構造について、以下のコードでパラメータ推定します。パラメータの推定や構造推定は試行錯誤するところなので、こんな簡単に済む訳では無いようです。

#パラメータ学習
bn.fit(学習した構造, データセット,method = "mle")

1 2	#パラメータ学習 bn.fit(学習した構造, データセット,method = "mle")

・予測
テストデータを用意して、データセットから求めた構造・パラメータを用いて、任意の変数の予測を行います。予測したい変数の値は欠損していたらエラーで予測してくれないですが、全部0にしておけば回ります。運用上は予測したい変数の値はわからないはずなので。

#ベイジアンネットワークによる予測
library(forecast)
training.set = protein_dataset[1:405, ] #パラメータ推定用の訓練用データ
test.set = protein_dataset[406:810, ]  #テスト用データ
baysian_structure = hc(training.set)   #訓練用データでベイジアンネットワークの構造を学習
fitted = bn.fit(baysian_structure, training.set,method = "mle")     #パラメータの学習
test.set2 <- test.set
test.set2$PKC <- 0 #本番想定でこの値を知らないということにしておく。
pred = predict(fitted, "PKC", test.set2)  #テストデータが与えられたもとでのPKCの予測
head(cbind(pred, test.set[, "PKC"]))      #予測値と実績値の比較
accuracy(f = pred, x = test.set[, "PKC"]) #正確度の算出


> head(cbind(pred, test.set[, "PKC"]))      #予測値と実績値の比較
         pred     
[1,] 19.76000 22.3
[2,] 21.66300 26.4
[3,] 11.82080 27.9
[4,] 23.85982 53.3
[5,] 21.82902 34.3
[6,] 17.29319 16.8
> accuracy(f = pred, x = test.set[, "PKC"]) #正確度の算出
              ME     RMSE      MAE       MPE     MAPE
Test set 0.47376 10.02839 7.479115 -90.60837 114.6982

#ベイジアンネットワークによる予測

library(forecast)

training.set = protein_dataset[1:405, ] #パラメータ推定用の訓練用データ

test.set = protein_dataset[406:810, ] #テスト用データ

baysian_structure = hc(training.set) #訓練用データでベイジアンネットワークの構造を学習

fitted = bn.fit(baysian_structure, training.set,method = "mle") #パラメータの学習

test.set2 <- test.set

test.set2$PKC <- 0 #本番想定でこの値を知らないということにしておく。

pred = predict(fitted, "PKC", test.set2) #テストデータが与えられたもとでのPKCの予測

head(cbind(pred, test.set[, "PKC"])) #予測値と実績値の比較

accuracy(f = pred, x = test.set[, "PKC"]) #正確度の算出

> head(cbind(pred, test.set[, "PKC"])) #予測値と実績値の比較

pred

[1,] 19.76000 22.3

[2,] 21.66300 26.4

[3,] 11.82080 27.9

[4,] 23.85982 53.3

[5,] 21.82902 34.3

[6,] 17.29319 16.8

> accuracy(f = pred, x = test.set[, "PKC"]) #正確度の算出

ME RMSE MAE MPE MAPE

Test set 0.47376 10.02839 7.479115 -90.60837 114.6982

#予測結果の可視化
library(ggplot2)
df <- data.frame(prediction=pred,actual=test.set[, "PKC"])
p <- ggplot(df,aes(x = prediction,y = actual)) + geom_point() 
p <- p + geom_line(data = data.frame(x = c(0,40), y = c(0,40)),aes(x = x, y = y), colour = "red")
p

#予測結果の可視化

library(ggplot2)

df <- data.frame(prediction=pred,actual=test.set[, "PKC"])

p <- ggplot(df,aes(x = prediction,y = actual)) + geom_point()

p <- p + geom_line(data = data.frame(x = c(0,40), y = c(0,40)),aes(x = x, y = y), colour = "red")

全然45度線上に乗っていないので、あまり精度は高くないようです。

cpquery関数と言って、条件付き確率を予測できる関数もあるようです。
Perform conditional probability queries

おまけ：Webサービスへの応用例

参考文献には、Webサイトの閲覧データをもとに、ユーザーの行動を予測するというモデルの紹介がなされていました。ただ、Webサイトの場合、ページ数や商品の数・ユーザーの数も膨大なことから、ネットワーク構造の推定が難しいようです。そこで、ネットワークを作るに際して、変数を確率的潜在意味解析（pLSA）で扱いやすい数に絞るなどの工夫をされていました。自社のサイトに関しても適用する際に、変数の多さは待った無しだと思うので、いざ適用する際は情報圧縮技術を駆使したいですね。

参考文献

ベイジアンネットワーク技術ユーザ・顧客のモデル化と不確実性推論

確率的グラフィカルモデル

Learning Bayesian Networks in R an Example in Systems Biology

グラフィカルモデル入門

ベイジアンネットワークを用いたWeb レコメンデーションシステムの開発

PRML8章

洋楽の歌詞データで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

先行研究

コード

参考文献

モデル概要

前処理

推定

結果

参考文献

Edwardとは

Edwardでできること

参考スライド

参考文献

Factorization Machine(FM)の概要

パッケージ紹介とインストール

サンプルデータの構造把握と前処理

FMの実行

結果

参考文献

概要

前処理

データ確認

行列因子分解（Matrix Factorization）

レコメンド結果の確認

感想

参考文献

はじめに

紹介する研究

オンライン広告におけるスパース性と遅れコンバージョンを考慮した予測モデル

目的と結果

対象となるデータ

手法の概要

PDFのリンク

テキストマイニングを用いた転職サイトの会員離脱予測

目的と結果

対象となるデータ

手法の概要

PDFのリンク

インターネット広告におけるスパースなユーザー行動ベクトルからのユーザー特徴抽出モデル

目的と結果

対象となるデータ

手法の概要

PDFのリンク

ユーザーのページビュー系列からのコンバージョン予測

目的と結果

対象となるデータ

手法の概要

PDFのリンク

SNSにおけるユーザの行動のモチベーション要因の分析

目的と結果

対象となるデータ

手法の概要

PDFのリンク

状況に応じた楽曲推薦に向けたソーシャルデータ分析

目的と結果

対象となるデータ

手法の概要

PDFのリンク

テキストマイニングを用いた口コミ分析による点数評価の信頼性確認手法

目的と結果

対象となるデータ

手法の概要

PDFのリンク

不動産仲介マーケティングのためのユーザ行動予測 Prediction of user behaviors for real estate brokerage marketing

目的と結果

対象となるデータ

手法の概要

PDFのリンク

SUUMOでの不動産データ活用の取り組みと未来

目的と結果

対象となるデータ

手法の概要

PDFのリンク

感想

参考文献

不動産仲介マーケティングのためのユーザ行動予測　Prediction of user behaviors for real estate brokerage marketing