Stanで順序プロビット(Ordered Probit)の推定のためのメモ書き

最近はBayesian Statistics and Marketingという本に関心があって、そこで取り上げられているモデルをStanに落とし込めないか模索しています。
そこで順序プロビット（Ordered Probit）の推定が必要であることがわかったため、Stanでの適用事例を漁っていました。まだマーケティング事例への適用はうまくいってないですが、いったん順序プロビットを簡単にまとめて今後の備忘録としておきます。

順序プロビットとは

被説明変数yが連続潜在変数y∗に対応していると考えるとする。
潜在変数は観察できないが、被説明変数yは観察でき、これらの2つ変数の関係は次のように表される。
（今回扱うデータは3から8までの順序データのため、以下のような表記になる。）

この対応関係は閾値メカニズムと呼ばれている。
各被説明変数をとる確率は以下のように記され、プロビットでは正規分布を扱うため以下のようになる。

これらの選択確率からなる尤度関数を最大にしたものが順序プロビットの推定となる。（c0=−∞でc6=∞とする。σは1とする。）

このように、潜在的な順序関係を想定し、それを満たすように閾値とパラメータを推定する点において、潜在変数を用いたモデルの柔軟性の高さが感じられる。

なぜ順序プロビットを使うのか

マーケティング業務において扱うデータにおいて、NPSやアンケートなど順序尺度の質的変数が多いので、それらのデータを二値データに落とし込んだり、そのまま基数データとして扱うのではなく、適切に扱いたいというモチベーションがあります。加えて、順序尺度の質的変数をもとに予測する際は普通のOLSだと、今回のケースで3を下回ったり、8を超えたりする可能性があり、予測結果として使いにくいです。
アンケートの点数をそのまま被説明変数として回帰しているケースは、データアナリティクスにこだわりの無いメンバーとかであればままあることなので、順序プロビットの民主化というか、布教していきたいと思います。

今回扱うデータ

勝手ながら大人のirisだと思っているワインデータです。今回は赤ワインに絞って、品質に関する順序変数を被説明変数として、各変数との相関を見ていきます。
まずはGGallyパッケージのggpairs()関数を適用して傾向を掴みます。見にくいので是非コードを回して確かめてください。

データに関する説明はワインの味（美味しさのグレード）は予測できるか？（１）で丁寧になされていますので、ご確認ください。

モデル

データセットに含まれる全部を含めて順序プロビットで回帰してみようと思います。
つまり、「酒石酸濃度、酢酸濃度、クエン酸濃度
残留糖分濃度、塩化ナトリウム濃度、遊離亜硫酸濃度
総亜硫酸濃度、密度、ｐＨ、硫酸カリウム濃度、アルコール度数」
の全てを使って赤ワインの質への影響を見ていきます。

Stanコード

最初に、Stanのユーザーガイド2.17の138ページにあるOrdered Probitのサンプルコードを使ってみたのですが、
収束しなかったので、初期値を設定するか弱情報事前分布を導入するかの判断が必要となります。
そこで、jabranhamさんが係数が平均0で分散10の正規分布に従うとするサンプルコードを書かれていたので、そちらを使って推定します。
書き換えているところはデータの制約くらいです。

data{
  int<lower=1> N; // number of obs
  int<lower=3> J; // number of categories
  int<lower=2> K; // num of predictors
  int<lower=0,upper=10> y[N]; // outcome var 
  matrix[N, K] x; // predictor vars 
}

parameters{
  ordered[J-1] tau; // thresholds
  vector[K] beta; // beta coefficients 
}

model{
  vector[J] theta;
  vector[N] xB;
  beta ~ normal(0, 10);
  xB <- x*beta;
  for(n in 1:N){
    theta[1] <- 1 - Phi(xB[n]-tau[1]);
    for(j in 2:J-1)
      theta[j] <- Phi(xB[n]-tau[j-1]) - Phi(xB[n]-tau[j]);
    theta[J] <- Phi(xB[n] - tau[J-1]);
    y[n] ~ categorical(theta);
  }
}

data{

int<lower=1> N; // number of obs

int<lower=3> J; // number of categories

int<lower=2> K; // num of predictors

int<lower=0,upper=10> y[N]; // outcome var

matrix[N, K] x; // predictor vars

}

parameters{

ordered[J-1] tau; // thresholds

vector[K] beta; // beta coefficients

}

model{

vector[J] theta;

vector[N] xB;

beta ~ normal(0, 10);

xB <- x*beta;

for(n in 1:N){

theta[1] <- 1 - Phi(xB[n]-tau[1]);

for(j in 2:J-1)

theta[j] <- Phi(xB[n]-tau[j-1]) - Phi(xB[n]-tau[j]);

theta[J] <- Phi(xB[n] - tau[J-1]);

y[n] ~ categorical(theta);

}

StanをキックするためのRコード

推定結果の可視化を行うためのcommon.Rは松浦さんのGitHubにあるものになります。

library(tidyverse)
library(rstan)
library(GGally)
library(shinystan)

rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())

wine_dataset <- read.csv("http://ieor.berkeley.edu/~ieor265/homeworks/winequality-red.csv", sep=";" )

y <- wine_dataset$quality - 2
x <- as.matrix(wine_dataset %>% select(-quality))
x <- scale(x)

# Visualization -----------------------------------------------------------
ggpairs(wine_dataset)


# Estimation --------------------------------------------------------------
stanmodel <- stan_model(file = "orderedprobit.stan")


N <- nrow(x)
J <- 6L
K <- ncol(x)

data_customer_list_test <- list(N=N,
                                J=J,
                                K=K,
                                y=y,
                                x=x)

fit <- stan(file = "orderedprobit.stan",
            data = data_customer_list_test,
            iter = 1000,
            chains = 4)

summary(fit)
traceplot(fit)

# Convergence Check -------------------------------------------------------
launch_shinystan(fit)

# Result Plot -------------------------------------------------------------
source('common.R')

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


param_names <- c('mcmc', colnames(wine_dataset %>% select(-quality)))

d_est <- data.frame(1:N_mcmc, ms$b)
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(tidyverse)

library(rstan)

library(GGally)

library(shinystan)

rstan_options(auto_write = TRUE)

options(mc.cores = parallel::detectCores())

wine_dataset <- read.csv("http://ieor.berkeley.edu/~ieor265/homeworks/winequality-red.csv", sep=";" )

y <- wine_dataset$quality - 2

x <- as.matrix(wine_dataset %>% select(-quality))

x <- scale(x)

# Visualization -----------------------------------------------------------

ggpairs(wine_dataset)

# Estimation --------------------------------------------------------------

stanmodel <- stan_model(file = "orderedprobit.stan")

N <- nrow(x)

J <- 6L

K <- ncol(x)

data_customer_list_test <- list(N=N,

J=J,

K=K,

y=y,

x=x)

fit <- stan(file = "orderedprobit.stan",

data = data_customer_list_test,

iter = 1000,

chains = 4)

summary(fit)

traceplot(fit)

# Convergence Check -------------------------------------------------------

launch_shinystan(fit)

# Result Plot -------------------------------------------------------------

source('common.R')

ms <- rstan::extract(fit)

N_mcmc <- length(ms$lp__)

param_names <- c('mcmc', colnames(wine_dataset %>% select(-quality)))

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

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

結果

まず、MCMCが収束したかどうかの判断ですが、ShinyStanに従うものとします。

ShinyStanによる収束診断をクリアできています。

続いて、推定したパラメータです。

係数の符号がはっきりと分かれている、赤ワインの品質に影響を与えそうな変数としては、volatile acidity（酢酸濃度）、chlorides（塩化ナトリウム濃度）、total.sulfur.dioxide（総亜硫酸濃度）、sulphates（硫酸カリウム濃度）、alcohol（アルコール度数）のようです。

最後に、推定した閾値です。

1と2の閾値が近く、2と3の開きが大きく、あとは比較的均等のようです。

比較

質的変数をそのまま重回帰した際の結果ですが、符号やその大小はあまり変わらないです。やはり予測の際にどの順序尺度の値に対応するかがわかるのが使う利点だと思います。

> dataset <- data.frame(cbind(y, x))
> multiple_regresion_model <- lm(formula = y ~ .,data = dataset)
> summary(multiple_regresion_model)

Call:
lm(formula = y ~ ., data = dataset)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.68911 -0.36652 -0.04699  0.45202  2.02498 

Coefficients:
                     Estimate Std. Error t value Pr(>|t|)    
(Intercept)           3.63602    0.01621 224.372  < 2e-16 ***
fixed.acidity         0.04351    0.04518   0.963   0.3357    
volatile.acidity     -0.19403    0.02168  -8.948  < 2e-16 ***
citric.acid          -0.03556    0.02867  -1.240   0.2150    
residual.sugar        0.02303    0.02115   1.089   0.2765    
chlorides            -0.08821    0.01973  -4.470 8.37e-06 ***
free.sulfur.dioxide   0.04562    0.02271   2.009   0.0447 *  
total.sulfur.dioxide -0.10739    0.02397  -4.480 8.00e-06 ***
density              -0.03375    0.04083  -0.827   0.4086    
pH                   -0.06386    0.02958  -2.159   0.0310 *  
sulphates             0.15533    0.01938   8.014 2.13e-15 ***
alcohol               0.29433    0.02822  10.429  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.648 on 1587 degrees of freedom
Multiple R-squared:  0.3606,    Adjusted R-squared:  0.3561 
F-statistic: 81.35 on 11 and 1587 DF,  p-value: < 2.2e-16

> dataset <- data.frame(cbind(y, x))

> multiple_regresion_model <- lm(formula = y ~ .,data = dataset)

> summary(multiple_regresion_model)

Call:

lm(formula = y ~ ., data = dataset)

Residuals:

Min 1Q Median 3Q Max

-2.68911 -0.36652 -0.04699 0.45202 2.02498

Coefficients:

Estimate Std. Error t value Pr(>|t|)

(Intercept) 3.63602 0.01621 224.372 < 2e-16 ***

fixed.acidity 0.04351 0.04518 0.963 0.3357

volatile.acidity -0.19403 0.02168 -8.948 < 2e-16 ***

citric.acid -0.03556 0.02867 -1.240 0.2150

residual.sugar 0.02303 0.02115 1.089 0.2765

chlorides -0.08821 0.01973 -4.470 8.37e-06 ***

free.sulfur.dioxide 0.04562 0.02271 2.009 0.0447 *

total.sulfur.dioxide -0.10739 0.02397 -4.480 8.00e-06 ***

density -0.03375 0.04083 -0.827 0.4086

pH -0.06386 0.02958 -2.159 0.0310 *

sulphates 0.15533 0.01938 8.014 2.13e-15 ***

alcohol 0.29433 0.02822 10.429 < 2e-16 ***

---

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.648 on 1587 degrees of freedom

Multiple R-squared: 0.3606, Adjusted R-squared: 0.3561

F-statistic: 81.35 on 11 and 1587 DF, p-value: < 2.2e-16

考察

マーケティングにおいて、順序尺度の質的変数を扱う際に順序プロビットを積極的に使っていきたいと思いますが、アンケート分析を行う際は、ユーザーごとの評点の癖が点数に影響を与えている可能性があります。
そのため、点数の付き方がユーザーごとに違うとする階層モデルへの拡張を今後行っていくのが面白いと思いますし、実際に研究されている論文があります。

参考文献

StanとRでベイズ統計モデリング (Wonderful R)
stan-examples/limited-dv/oprobit.stan
stan-dev/stan users-guide-2.18.0.pdf
Stanによる順序ロジット回帰
 第9章順序選択モデル：年金投資選択問題
 Wine Quality Data Set
ワインの味（美味しさのグレード）は予測できるか？（１）
世界一簡単な収束[シナイ]Stanコード
 RStanとShinyStanによるベイズ統計モデリング入門

おまけ

数式をブログに載せる際は、こちら
Online LaTeX Equation Editor – create, integrate and download
でインタラクティブに数式を作成し、その結果を
QuickLaTex Publish Math on the Web without compromising quality
に貼り付けて画像を出力しています。

Bayesian Statistics and Marketing – 混合ガウス×階層モデルのマーガリン購買データへの適用

前回の分析では、価格への反応係数の事前分布が正規分布を仮定したモデルを用いていましたが、事後分布から多峰性が観察されました。そこで今回は、各個人の価格への反応係数の事前分布が混合ガウス分布に従うとした場合の事例を扱いたいと思います。

データのおさらい

データ自体は前回のブログと同じですが、先日のTokyo.Rで松浦さんがオススメしていたGGallyパッケージのggpairs関数を用いて、今回扱うデータを可視化してみます。

まず、購買したマーガリンのブランド選択(6ブランド分)と、購入した価格（ドル）からなるデータセットの可視化をすると以下のようになります。

1つ目のブランドが最も多く選択されているようです。購入価格の分布は、ブランドによって多峰性がありそうです。

続いて、家計ごとの属性（家族構成、学歴、職位、退職の有無など）からなるデータセットの可視化をすると以下のようになります。

ほとんどダミー変数なので面白みには欠いていますね。ホワイトカラーの家計が多く、退職していない家計が多く、学歴が低い人が多いようです。年収に関しては対数正規分布に従ってそうです。

モデル

前回と同様に、「Bayesian Statistics and Marketing」の5章に載っている混合ガウスを想定した階層ベイズモデルを扱います。

yi ∼ MNL(Xi, βi) は購買レコードごとの意思決定が多項ロジスティック回帰モデルに従うということを意味し、βiは(説明変数の数×1)のベクトルとなります。
βi = Z∆[i,] + ui は購買レコードごとの価格にかかってくる係数で、その係数が家計の属性データに係数Δ(家計の数×属性データからなる説明変数の数)をかけ合わせたものと潜在的な項の和となります。なお、ここでの属性データには定数項を含めていません。係数Δは平均deltabar、分散A_delta^-1の多変量正規分布に従います。一方で、潜在的な項は平均µind、分散Σindの多変量正規分布に従います。この平均や分散に振られているindが、混合正規分布のパラメータの通し番号となり、多項分布に従います。この多項分布の割当確率がディリクレ分布に従います。最後に混合正規分布の各パラメータは、平均は正規分布に分散は逆ウィシャート分布に従います。

前回との大きな違いは、価格の係数の一部である、潜在的な項において、混合正規分布が仮定されているところになります。
ちなみに、モデルに関しての詳細はbayesmパッケージのマニュアルの61ページ目に記載されていました。

今回のモデルのDAG

Pythonのdaftというモジュールを使うことで、非常に簡単に今回のモデルのDAG(有向非巡回グラフ)を描くことができます。

今回はこちらのPythonコードで描けました。

import daft
from matplotlib import rc

pgm = daft.PGM(shape=[6,6])

# Nodes
pgm.add_node(daft.Node("pvec", r"pvec",1, 2)) # 名前　ラベル　座標(横、縦)
pgm.add_node(daft.Node("ind", r"ind",2, 2))
pgm.add_node(daft.Node("sigma", r"$\Sigma_k$",3, 1))
pgm.add_node(daft.Node("mu", r"$\mu_k$",3, 3))
pgm.add_node(daft.Node("theta", r"$\theta_i$",4, 2))
pgm.add_node(daft.Node("delta", r"$\Delta$",4, 0.5))
pgm.add_node(daft.Node("y", r"$y_i$",5, 2,observed=True))

# Edges
pgm.add_edge("pvec", "ind")
pgm.add_edge("ind", "sigma")
pgm.add_edge("ind", "mu")
pgm.add_edge("sigma", "mu")
pgm.add_edge("sigma", "theta")
pgm.add_edge("mu", "theta")
pgm.add_edge("delta", "theta")
pgm.add_edge("theta", "y")

pgm.render()
pgm.figure.savefig("Hierachical_Models_ForMixtureNormal.png")

import daft

from matplotlib import rc

pgm = daft.PGM(shape=[6,6])

# Nodes

pgm.add_node(daft.Node("pvec", r"pvec",1, 2)) # 名前　ラベル　座標(横、縦)

pgm.add_node(daft.Node("ind", r"ind",2, 2))

pgm.add_node(daft.Node("sigma", r"$\Sigma_k$",3, 1))

pgm.add_node(daft.Node("mu", r"$\mu_k$",3, 3))

pgm.add_node(daft.Node("theta", r"$\theta_i$",4, 2))

pgm.add_node(daft.Node("delta", r"$\Delta$",4, 0.5))

pgm.add_node(daft.Node("y", r"$y_i$",5, 2,observed=True))

# Edges

pgm.add_edge("pvec", "ind")

pgm.add_edge("ind", "sigma")

pgm.add_edge("ind", "mu")

pgm.add_edge("sigma", "mu")

pgm.add_edge("sigma", "theta")

pgm.add_edge("mu", "theta")

pgm.add_edge("delta", "theta")

pgm.add_edge("theta", "y")

pgm.render()

pgm.figure.savefig("Hierachical_Models_ForMixtureNormal.png")

stanコード

今回扱ったstanコードとなります。誤りがある場合はお知らせしていただけると幸いです。

data{
  int<lower=0> N_x; // 購買レコードの数
  int<lower=0> N_z; // 家計の数
  int<lower=0> p_x; // 購買レコードの項目数
  int<lower=0> p_z; // 家計の属性データの項目数
  
  int y[N_x]; // 選択肢
  matrix[N_x, p_x] X; // 説明変数
  matrix[N_z, p_z] Z; // 家計の属性データ
  int<lower=0> hhid[N_x];  // 家計ID

  int K; // 混合ガウス分布の要素数
}

transformed data{
  real nu;
  matrix[p_x, p_x] I; // 購買レコードの説明変数の数の正方行列
  matrix[p_x, p_x] J; // 属性データの説明変数の数の正方行列
  
  nu = p_x + 3; // 説明変数の項に3を足す
  I = diag_matrix(rep_vector(1, p_x)); // 1を繰り返しp_x個並べた対角行列を作成
  J = diag_matrix(rep_vector(1, p_x));
}

parameters{
  vector[p_x] theta_ast[N_z]; // 説明変数の数だけある、購買ごとのパラメータ
  matrix[p_z, p_x] Delta; // 属性データの説明変数の数×購買データの説明変数の数だけのパラメータ
  vector[p_x] u[N_z]; // 購買レコードごとのパラメータ
  vector[p_x] mu[K]; // 混合分布を構成する平均値
  cov_matrix[p_x] Sigma[K]; // 共分散行列
  simplex[K] pi[N_z]; // シンプレックス（各要素が[0,1]の範囲で合計が1という条件を満たす。）
}

transformed parameters{
  vector[p_x] theta[N_z]; #家計の数だけの係数ベクトル

  for(i in 1:N_z){ // 家計の数だけ繰り返す
    theta[i] = theta_ast[i] + Delta' * Z[i]'; // 係数は家計属性ごとの特徴に異質なDeltaとbeta_astの和で決まる 
  }
}

model{
  real ps[K]; // 混合正規分布の対数

  for(i in 1:N_x){ // 購買レコードの数だけ繰り返す
    y[i] ~ categorical(softmax(theta[hhid[i]] .* to_vector(X[i]))); //カテゴリカル分布にsoftmaxを組み合わせて多項ロジスティック回帰を行う

    pi[hhid[i]] ~ dirichlet(rep_vector(2, K));
    for(k in 1:K){ // 混合分布の構成要素の数だけ繰り返す
      ps[k] = log(pi[hhid[i]][k]) + multi_normal_lpdf(theta_ast[hhid[i]] | mu[k], Sigma[k]);
    }
    target += log_sum_exp(ps); // 離散パラメータを消去した形で対数尤度を表現する際に必要な計算（周辺化消去）。  
  }

  for(k in 1:K){ // 混合パラメータの数だけ繰り返す
    mu[k] ~ multi_normal(rep_vector(0, p_x), 100*Sigma[k]);
    Sigma[k] ~ inv_wishart(nu, nu*I); 
  }

  for(i in 1:p_z){ // 家計の属性データの数だけ繰り返す
    Delta[i] ~ multi_normal_cholesky(rep_vector(0, p_x), 100*J); 
  }
}

data{

int<lower=0> N_x; // 購買レコードの数

int<lower=0> N_z; // 家計の数

int<lower=0> p_x; // 購買レコードの項目数

int<lower=0> p_z; // 家計の属性データの項目数

int y[N_x]; // 選択肢

matrix[N_x, p_x] X; // 説明変数

matrix[N_z, p_z] Z; // 家計の属性データ

int<lower=0> hhid[N_x]; // 家計ID

int K; // 混合ガウス分布の要素数

}

transformed data{

real nu;

matrix[p_x, p_x] I; // 購買レコードの説明変数の数の正方行列

matrix[p_x, p_x] J; // 属性データの説明変数の数の正方行列

nu = p_x + 3; // 説明変数の項に3を足す

I = diag_matrix(rep_vector(1, p_x)); // 1を繰り返しp_x個並べた対角行列を作成

J = diag_matrix(rep_vector(1, p_x));

}

parameters{

vector[p_x] theta_ast[N_z]; // 説明変数の数だけある、購買ごとのパラメータ

matrix[p_z, p_x] Delta; // 属性データの説明変数の数×購買データの説明変数の数だけのパラメータ

vector[p_x] u[N_z]; // 購買レコードごとのパラメータ

vector[p_x] mu[K]; // 混合分布を構成する平均値

cov_matrix[p_x] Sigma[K]; // 共分散行列

simplex[K] pi[N_z]; // シンプレックス（各要素が[0,1]の範囲で合計が1という条件を満たす。）

}

transformed parameters{

vector[p_x] theta[N_z]; #家計の数だけの係数ベクトル

for(i in 1:N_z){ // 家計の数だけ繰り返す

theta[i] = theta_ast[i] + Delta' * Z[i]'; // 係数は家計属性ごとの特徴に異質なDeltaとbeta_astの和で決まる

}

model{

real ps[K]; // 混合正規分布の対数

for(i in 1:N_x){ // 購買レコードの数だけ繰り返す

y[i] ~ categorical(softmax(theta[hhid[i]] .* to_vector(X[i]))); //カテゴリカル分布にsoftmaxを組み合わせて多項ロジスティック回帰を行う

pi[hhid[i]] ~ dirichlet(rep_vector(2, K));

for(k in 1:K){ // 混合分布の構成要素の数だけ繰り返す

ps[k] = log(pi[hhid[i]][k]) + multi_normal_lpdf(theta_ast[hhid[i]] | mu[k], Sigma[k]);

}

target += log_sum_exp(ps); // 離散パラメータを消去した形で対数尤度を表現する際に必要な計算（周辺化消去）。

}

for(k in 1:K){ // 混合パラメータの数だけ繰り返す

mu[k] ~ multi_normal(rep_vector(0, p_x), 100*Sigma[k]);

Sigma[k] ~ inv_wishart(nu, nu*I);

}

for(i in 1:p_z){ // 家計の属性データの数だけ繰り返す

Delta[i] ~ multi_normal_cholesky(rep_vector(0, p_x), 100*J);

}

stanをキックするためのコードです。先人が書かれた混合ガウスのスクリプトをHMCで実行した際に26時間ほどかかったので、今回はより複雑なモデルであることから、変分ベイズ法による推定を行ってみることにしました。松浦さんの教科書にあるように、vb関数を用いて変分ベイズ推論を行っています。

library(bayesm)
library(tidyverse)
library(rstan)
library(GGally)

rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())

data("margarine")

#1,2,3,4,5,7の商品に関してデータを抽出し、家計IDごとにカウントし、5件以上のものに絞る。
hhid_selected <- margarine$choicePrice %>% 
                                  filter(choice %in% c(1,2,3,4,5,7)) %>% 
                                  group_by(hhid) %>% 
                                  summarise(purc_cnt = n()) %>% 
                                  filter(purc_cnt >= 5)

#今回扱う商品のカラムだけを抽出し、先ほど絞ったユーザーのリストに合致するデータでフィルターする。
choicePrice.selected <- margarine$choicePrice %>% 
                                             filter(choice %in% c(1,2,3,4,5,7) & hhid %in% hhid_selected$hhid)
#並べにくいので7を6に置き換える。
choicePrice.selected$choice[choicePrice.selected$choice == 7] <- 6

#家計ごとに関する属性データの抽出
demos.selected <- margarine$demos %>% filter(hhid %in% hhid_selected$hhid)

#データサイズ
N <- nrow(choicePrice.selected)

#選択肢の数（特に使っているデータではない。）
p <- n_distinct(choicePrice.selected$choice)

#被説明変数
y <- choicePrice.selected$choice

#説明変数
X <- choicePrice.selected %>% select(3,4,5,6,7,9)

#家計の属性データから家計IDを除く
Z <- demos.selected %>% 
                                   select(-hhid)

#可視化
ggpairs(choicePrice.selected %>% select(2,3,4,5,6,7,9))
ggpairs(Z)


#家計の属性データから家計IDを抽出し、1から行数までのインデックスを付与する。
hhid_index <- demos.selected %>%
                                                   select(hhid) %>% 
                                                   mutate(ind = seq(1,nrow(demos.selected)))

#購買データの家計IDを抽出し、先ほどのインデックスとjoinする
hhid_x <- choicePrice.selected %>% 
                                                    select(hhid) %>% 
                                                    left_join(hhid_index)

#stanで扱うデータリストの作成
d.dat <- list(N_x=nrow(X), N_z=nrow(Z), 
              p_x=ncol(X), p_z=ncol(Z),
              y=y, X=X, Z=Z,
              hhid = hhid_x$ind,
              K = 3)

#推定
stanmodel <- stan_model(file = "Hierarchical_Mixture.stan")
d.fit_vb <- vb(stanmodel, data=d.dat, seed=123)

library(bayesm)

library(tidyverse)

library(rstan)

library(GGally)

rstan_options(auto_write = TRUE)

options(mc.cores = parallel::detectCores())

data("margarine")

#1,2,3,4,5,7の商品に関してデータを抽出し、家計IDごとにカウントし、5件以上のものに絞る。

hhid_selected <- margarine$choicePrice %>%

filter(choice %in% c(1,2,3,4,5,7)) %>%

group_by(hhid) %>%

summarise(purc_cnt = n()) %>%

filter(purc_cnt >= 5)

#今回扱う商品のカラムだけを抽出し、先ほど絞ったユーザーのリストに合致するデータでフィルターする。

choicePrice.selected <- margarine$choicePrice %>%

filter(choice %in% c(1,2,3,4,5,7) & hhid %in% hhid_selected$hhid)

#並べにくいので7を6に置き換える。

choicePrice.selected$choice[choicePrice.selected$choice == 7] <- 6

#家計ごとに関する属性データの抽出

demos.selected <- margarine$demos %>% filter(hhid %in% hhid_selected$hhid)

#データサイズ

N <- nrow(choicePrice.selected)

#選択肢の数（特に使っているデータではない。）

p <- n_distinct(choicePrice.selected$choice)

#被説明変数

y <- choicePrice.selected$choice

#説明変数

X <- choicePrice.selected %>% select(3,4,5,6,7,9)

#家計の属性データから家計IDを除く

Z <- demos.selected %>%

select(-hhid)

#可視化

ggpairs(choicePrice.selected %>% select(2,3,4,5,6,7,9))

ggpairs(Z)

#家計の属性データから家計IDを抽出し、1から行数までのインデックスを付与する。

hhid_index <- demos.selected %>%

select(hhid) %>%

mutate(ind = seq(1,nrow(demos.selected)))

#購買データの家計IDを抽出し、先ほどのインデックスとjoinする

hhid_x <- choicePrice.selected %>%

select(hhid) %>%

left_join(hhid_index)

#stanで扱うデータリストの作成

d.dat <- list(N_x=nrow(X), N_z=nrow(Z),

p_x=ncol(X), p_z=ncol(Z),

y=y, X=X, Z=Z,

hhid = hhid_x$ind,

K = 3)

#推定

stanmodel <- stan_model(file = "Hierarchical_Mixture.stan")

d.fit_vb <- vb(stanmodel, data=d.dat, seed=123)

実行結果

今回はK={1,3,5,10}の混合要素数で推定を行いました。

K=1のケース

K=3のケース

K=5のケース

K=10のケース

Kが小さいと散らばりが比較的小さそうに見えます。逆に、Kが大きくなると散らばりが出てくるようです。

おわりに

　２回に渡って、マーケティングデータを用いたベイズ統計モデリングを学んできましたが、数式からstanコードに落とし込む作業の際に、stanの関数をある程度知らないとやりにくいという至極当然なことを実感しました。これまでに使ってきたstanの関数は限られたものしか扱っていなかったと言えます。特に今回は離散パラメータをstanで扱うパーツがあったので、先行研究や松浦さんの本を読みながらの手探りが多かったです。
　あと、複雑な階層ベイズモデルを扱う際に、頭の中を整理しないと手が止まってしまう感じがあったので、数式と対応するコードを横に並べながら進めました。
　マーケティングにおいては、顧客の属性ごとに多峰性のあるような事例を扱うことが多く、かつ各々のサンプルサイズも期待できないことが多いので、今回の学びを分析業務で試してみたいと思います。

参考情報

Bayesian Statistics and Marketing (Wiley Series in Probability and Statistics)
Package ‘bayesm’
Multivariate Gaussian Mixture Model done properly
StanとRでベイズ統計モデリング (Wonderful R)
機械学習スタートアップシリーズベイズ推論による機械学習入門 (KS情報科学専門書)
daftでグラフィカルモデル

Bayesian Statistics and Marketingの5章 – 家計の異質性を考慮した階層ベイズモデル

はじめに

ゴールデンウィークで実家に持ち込む本としてチョイスしたのが、2005年出版の「Bayesian Statistics and Marketing」です。大学院のときに購入して、ちょっとしか読んでませんでした。

この本は、字面の通りマーケティング関連の分析に関してベイズ統計を使ってアプローチするというもので、この書籍のために作られた、Rのbayesmというパッケージの紹介もあり、理論だけでなくRで実践することもできます。1章から7章までの全ての分析事例に対して実行可能な関数が用意されています。（CRANにあるdocumentも120p程度と割と大きめのパッケージです。）

和書で言うと、東北大学の照井先生の「ベイズモデリングによるマーケティング分析」などがありますが、その82pでもBayesian Statistics and Marketingとbayesmパッケージが紹介されています。

今回は、5章に載っている階層ベイズモデルを用いた、家計の異質性を考慮したブランド選択モデルの分析を紹介します。加えて、GitHubでstanによる再現を試みている方がいらっしゃったので、その方のコードの紹介も行います。

最近はこれまで以上にベイズ統計が流行ってきていますが、マーケティング×ベイズの書籍は限られている印象なので、少しでもリサーチのお役に立てれば幸いです。

目的

マーガリンの購買データから、ブランドごと、家計ごとのマーガリン価格に対しての反応の違いを明らかにしたい。

データ

bayesmパッケージにある、margarineデータ。data(margarine)で呼び出せ、詳細はcranのドキュメントに載っています。

Household Panel Data on Margarine Purchasesには、516家計の購買データと、家計ごとのデモグラフィック情報が収められています。1991年の論文のデータとなるので、かなり昔のデータです。

購買データは価格（USドル）と選択したブランドのID（10種類）
デモグラフィック情報はfamily size（家族構成）、学歴、職位、退職の有無などのダミー変数

今回の事例では、5回以上購買した家計に限定して分析しているため、
313家計・3405の購買レコードからなるデータセットとなります。

モデル

家計ごとに異なる、マーガリン価格に対する反応を想定。各マーガリンのブランドの価格に対するパラメータの数は家計の数だけある。
価格に対する反応は家計の属性によっても決まる。
という前提に立ち、以下のセッティングで推論していきます。
6つのブランド選択に関する多項ロジスティックモデル（カテゴリカル分布とsoftmax関数の適用）
1階層目はブランドごとの価格を説明変数とし、価格に対する反応係数をかけ合わせたものを多項ロジスティックモデルの入力とする。
2階層目はブランドの価格に対する反応係数が家計ごとの定数項と属性データに属性ごとの係数をかけ合わせたものからなる。
家計ごとの定数項は平均0、分散V_betaの正規分布に従う。
属性ごとの係数は平均vec(delta_bar)、分散V_betaクロネッカーのデルタA^(-1)の正規分布に従う。
分散V_betaは平均υ、分散Vの逆ウィシャート分布に従う。
A = 0.01、υ = 6 + 3 = 9、V = υI(Iは単位行列)

$\textit{ l } $ \beta_i | y_i , X_i $ [ Multinomial Logit ]$

$B = Z \varDelta + U , u_i \sim N(0,V_\beta)$

$vec(\varDelta | V_\beta ) \sim N(vec(\bar{\varDelta}), V_\beta \otimes A^{-1} )$

$V_\beta \sim IW(\upsilon ,V )$

コード

kefitsさんがいくつかの章に登場するbayesmでの実践例をstanに書き直されているようですので、そちらのコードで学ばせていただこうと思います。
https://github.com/kefits/Bayesian-Statistics-and-Marketing

以下が、stanのコードとなっています。ここでは、Hierarchical_MNL.stanとして保存します。

data{
  int<lower=0> N_x; // 購買レコードの数
  int<lower=0> N_z; // 家計の数
  int<lower=0> p_x; // 購買レコードの項目数
  int<lower=0> p_z; // 家計の属性データの項目数
  
  int y[N_x]; // 選択肢
  matrix[N_x, p_x] X; // 説明変数
  matrix[N_z, p_z] Z; // 家計の属性データ
  int<lower=0> hhid[N_x];  // 家計ID
}

transformed data{
  real nu;
  matrix[p_x, p_x] I; // 説明変数の数の正方行列
  
  nu = p_x + 3; // 説明変数の項に3を足す
  I = diag_matrix(rep_vector(1, p_x)); // 1を繰り返しp_x個並べた対角行列を作成
}

parameters{
  vector[p_x] beta_ast[N_z]; // 説明変数の数だけある、家計ごとのパラメータ
  matrix[p_z, p_x] Delta; // 属性データの説明変数の数×購買データの説明変数の数だけのパラメータ
  cov_matrix[p_x] V_b; // 共分散行列
}

transformed parameters{
  vector[p_x] beta[N_z]; #家計の数だけの係数ベクトル
  matrix[p_x, p_x] L_b; #共分散行列（beta(家計ごとの係数の共分散)）
  matrix[p_x, p_x] L_d; #共分散行列（delta(属性データの係数の共分散)）
  
  L_b = cholesky_decompose(V_b); // 共分散行列のコレスキー因子をもとめる
  L_d = cholesky_decompose(100*V_b); // 共分散行列に0.01で割ったもののコレスキー因子をもとめる
  for(i in 1:N_z){
    beta[i] = beta_ast[i] + Delta' * Z[i]'; // 係数は家計属性ごとの特徴に異質なDeltaとbeta_astの和で決まる 
  }
}

model{
  for(i in 1:N_x){
    y[i] ~ categorical(softmax(beta[hhid[i]] .* to_vector(X[i]))); //カテゴリカル分布にsoftmaxを組み合わせて多項ロジスティック回帰を行う
  }
  for(i in 1:p_z){
    Delta[i] ~ multi_normal_cholesky(rep_vector(0, p_x), L_d); // コレスキー因子（L_d）を引数にとる多変量正規分布(推定の高速化のために用いることがある。)
  }
  beta_ast ~ multi_normal_cholesky(rep_vector(0, p_x), L_b); // コレスキー因子（L_b）を引数にとる多変量正規分布
  V_b ~ inv_wishart(nu, nu*I); // 正規分布の共分散行列の共役事前分布として逆ウィシャート分布を利用
}

data{

int<lower=0> N_x; // 購買レコードの数

int<lower=0> N_z; // 家計の数

int<lower=0> p_x; // 購買レコードの項目数

int<lower=0> p_z; // 家計の属性データの項目数

int y[N_x]; // 選択肢

matrix[N_x, p_x] X; // 説明変数

matrix[N_z, p_z] Z; // 家計の属性データ

int<lower=0> hhid[N_x]; // 家計ID

}

transformed data{

real nu;

matrix[p_x, p_x] I; // 説明変数の数の正方行列

nu = p_x + 3; // 説明変数の項に3を足す

I = diag_matrix(rep_vector(1, p_x)); // 1を繰り返しp_x個並べた対角行列を作成

}

parameters{

vector[p_x] beta_ast[N_z]; // 説明変数の数だけある、家計ごとのパラメータ

matrix[p_z, p_x] Delta; // 属性データの説明変数の数×購買データの説明変数の数だけのパラメータ

cov_matrix[p_x] V_b; // 共分散行列

}

transformed parameters{

vector[p_x] beta[N_z]; #家計の数だけの係数ベクトル

matrix[p_x, p_x] L_b; #共分散行列（beta(家計ごとの係数の共分散)）

matrix[p_x, p_x] L_d; #共分散行列（delta(属性データの係数の共分散)）

L_b = cholesky_decompose(V_b); // 共分散行列のコレスキー因子をもとめる

L_d = cholesky_decompose(100*V_b); // 共分散行列に0.01で割ったもののコレスキー因子をもとめる

for(i in 1:N_z){

beta[i] = beta_ast[i] + Delta' * Z[i]'; // 係数は家計属性ごとの特徴に異質なDeltaとbeta_astの和で決まる

}

model{

for(i in 1:N_x){

y[i] ~ categorical(softmax(beta[hhid[i]] .* to_vector(X[i]))); //カテゴリカル分布にsoftmaxを組み合わせて多項ロジスティック回帰を行う

}

for(i in 1:p_z){

Delta[i] ~ multi_normal_cholesky(rep_vector(0, p_x), L_d); // コレスキー因子（L_d）を引数にとる多変量正規分布(推定の高速化のために用いることがある。)

}

beta_ast ~ multi_normal_cholesky(rep_vector(0, p_x), L_b); // コレスキー因子（L_b）を引数にとる多変量正規分布

V_b ~ inv_wishart(nu, nu*I); // 正規分布の共分散行列の共役事前分布として逆ウィシャート分布を利用

}

以下はstanをキックするためのRコードです。

library(bayesm)
library(dplyr)
library(rstan)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())

data("margarine")

#1,2,3,4,5,7の商品に関してデータを抽出し、家計IDごとにカウントし、5件以上のものに絞る。
hhid_selected <- margarine$choicePrice %>% 
                  filter(choice %in% c(1,2,3,4,5,7)) %>% 
                  group_by(hhid) %>% 
                  summarise(purc_cnt = n()) %>% 
                  filter(purc_cnt >= 5)

#今回扱う商品のカラムだけを抽出し、先ほど絞ったユーザーのリストに合致するデータでフィルターする。
choicePrice.selected <- margarine$choicePrice %>% 
                          filter(choice %in% c(1,2,3,4,5,7) & hhid %in% hhid_selected$hhid)
#並べにくいので7を6に置き換える。
choicePrice.selected$choice[choicePrice.selected$choice == 7] <- 6

#家計ごとに関する属性データの抽出
demos.selected <- margarine$demos %>% filter(hhid %in% hhid_selected$hhid)

#データサイズ
N <- nrow(choicePrice.selected)

#選択肢の数（特に使っているデータではない。）
p <- n_distinct(choicePrice.selected$choice)

#被説明変数
y <- choicePrice.selected$choice

#説明変数
X <- choicePrice.selected %>% select(3,4,5,6,7,9)

#家計の属性データから家計IDを除く
Z <- demos.selected %>% 
        select(-hhid)

#定数項を1列目に追加する
Z <- data.frame(intercept = rep(1, nrow(Z))) %>% 
        bind_cols(Z)

#家計の属性データから家計IDを抽出し、1から行数までのインデックスを付与する。
hhid_index <- demos.selected %>%
                select(hhid) %>% 
                mutate(ind = seq(1,nrow(demos.selected)))

#購買データの家計IDを抽出し、先ほどのインデックスとjoinする
hhid_x <- choicePrice.selected %>% 
            select(hhid) %>% 
            left_join(hhid_index)

#stanで扱うデータリストの作成
d.dat <- list(N_x=nrow(X), N_z=nrow(Z), 
              p_x=ncol(X), p_z=ncol(Z),
              y=y, X=X, Z=Z,
              hhid = hhid_x$ind)

#推定
d.fit <- stan("../Chapter5/Hierarchical_MNL.stan",
              data = d.dat,
              iter = 500,
              chains = 4)

library(bayesm)

library(dplyr)

library(rstan)

rstan_options(auto_write = TRUE)

options(mc.cores = parallel::detectCores())

data("margarine")

#1,2,3,4,5,7の商品に関してデータを抽出し、家計IDごとにカウントし、5件以上のものに絞る。

hhid_selected <- margarine$choicePrice %>%

filter(choice %in% c(1,2,3,4,5,7)) %>%

group_by(hhid) %>%

summarise(purc_cnt = n()) %>%

filter(purc_cnt >= 5)

#今回扱う商品のカラムだけを抽出し、先ほど絞ったユーザーのリストに合致するデータでフィルターする。

choicePrice.selected <- margarine$choicePrice %>%

filter(choice %in% c(1,2,3,4,5,7) & hhid %in% hhid_selected$hhid)

#並べにくいので7を6に置き換える。

choicePrice.selected$choice[choicePrice.selected$choice == 7] <- 6

#家計ごとに関する属性データの抽出

demos.selected <- margarine$demos %>% filter(hhid %in% hhid_selected$hhid)

#データサイズ

N <- nrow(choicePrice.selected)

#選択肢の数（特に使っているデータではない。）

p <- n_distinct(choicePrice.selected$choice)

#被説明変数

y <- choicePrice.selected$choice

#説明変数

X <- choicePrice.selected %>% select(3,4,5,6,7,9)

#家計の属性データから家計IDを除く

Z <- demos.selected %>%

select(-hhid)

#定数項を1列目に追加する

Z <- data.frame(intercept = rep(1, nrow(Z))) %>%

bind_cols(Z)

#家計の属性データから家計IDを抽出し、1から行数までのインデックスを付与する。

hhid_index <- demos.selected %>%

select(hhid) %>%

mutate(ind = seq(1,nrow(demos.selected)))

#購買データの家計IDを抽出し、先ほどのインデックスとjoinする

hhid_x <- choicePrice.selected %>%

select(hhid) %>%

left_join(hhid_index)

#stanで扱うデータリストの作成

d.dat <- list(N_x=nrow(X), N_z=nrow(Z),

p_x=ncol(X), p_z=ncol(Z),

y=y, X=X, Z=Z,

hhid = hhid_x$ind)

#推定

d.fit <- stan("../Chapter5/Hierarchical_MNL.stan",

data = d.dat,

iter = 500,

chains = 4)

実行結果

Core i5、8GBメモリのMacBook Proで40分ほどかかりました。

traceplot(d.fit)で以下のように4回の試行結果が描かれますが、収束しているようです。

summary関数を使えばわかりますが、3913行ものパラメータたちのサマリーが得られます。

313家計の家計ごとのブランドに対するパラメータ(1878個)
313家計の家計ごとのブランドに対する潜在パラメータ(1878個)
6ブランドの係数の共分散行列(36個)
6ブランドの係数の分散のハイパーパラメータの行列(36個)
6ブランドの属性データ(8つ)に対する係数(48個)
6ブランドの属性データに対する係数の共分散行列(36個)
lp(log posterior(確率密度の和でモデル比較で扱う。))(1個)

64番目の家計の各ブランドの価格に対する係数の分布を確認すると、4番目・5番目のブランドの係数が他のブランドに比べて小さいことがわかります。

続いて、家計ごとの係数に関して集計し、係数ごとの相関係数を見てみると、各ブランドごとに正の相関、負の相関がありそうです。

#トレースプロット
traceplot(d.fit)

#係数のサマリー
summary_table <- summary(d.fit)$summary

draws <- extract(d.fit)
beta <- as.data.frame(draws$beta)
Delta <- as.data.frame(draws$Delta)
V_b <- as.data.frame(draws$V_b)

hhid_info <- inner_join(hhid_index, hhid_selected)

# 1000行*313列のデータを313000行*1列のデータにしたい。
for (i in 1:6) {
  nam <- paste("beta", i, sep = "")
  assign(nam, beta[,(1+313*(i-1)):(313*(i))] %>% tidyr::gather(key, value))
}

beta_matrix <- beta1 %>% bind_cols(beta2,beta3,beta4,beta5,beta6)
beta_matrix <- beta_matrix %>% select(-starts_with("key"))

#相関係数
cor(beta_matrix)

            value     value1     value2     value3     value4    value5
value   1.0000000  0.5902734  0.4864998 -0.1798877 -0.4781558 0.3188025
value1  0.5902734  1.0000000  0.6134343 -0.2484286 -0.4441728 0.2002954
value2  0.4864998  0.6134343  1.0000000  0.1336322 -0.4512474 0.3663121
value3 -0.1798877 -0.2484286  0.1336322  1.0000000  0.6149186 0.2671819
value4 -0.4781558 -0.4441728 -0.4512474  0.6149186  1.0000000 0.1591574
value5  0.3188025  0.2002954  0.3663121  0.2671819  0.1591574 1.0000000

#トレースプロット

traceplot(d.fit)

#係数のサマリー

summary_table <- summary(d.fit)$summary

draws <- extract(d.fit)

beta <- as.data.frame(draws$beta)

Delta <- as.data.frame(draws$Delta)

V_b <- as.data.frame(draws$V_b)

hhid_info <- inner_join(hhid_index, hhid_selected)

# 1000行*313列のデータを313000行*1列のデータにしたい。

for (i in 1:6) {

nam <- paste("beta", i, sep = "")

assign(nam, beta[,(1+313*(i-1)):(313*(i))] %>% tidyr::gather(key, value))

}

beta_matrix <- beta1 %>% bind_cols(beta2,beta3,beta4,beta5,beta6)

beta_matrix <- beta_matrix %>% select(-starts_with("key"))

#相関係数

cor(beta_matrix)

value value1 value2 value3 value4 value5

value 1.0000000 0.5902734 0.4864998 -0.1798877 -0.4781558 0.3188025

value1 0.5902734 1.0000000 0.6134343 -0.2484286 -0.4441728 0.2002954

value2 0.4864998 0.6134343 1.0000000 0.1336322 -0.4512474 0.3663121

value3 -0.1798877 -0.2484286 0.1336322 1.0000000 0.6149186 0.2671819

value4 -0.4781558 -0.4441728 -0.4512474 0.6149186 1.0000000 0.1591574

value5 0.3188025 0.2002954 0.3663121 0.2671819 0.1591574 1.0000000

最後に、家計ごとに集計した、ブランドに対する価格反応係数の事後分布を描きます。

~~多峰性などはなく、正規分布に従っているようです。他のブランドと比較して、5番目の係数が小さいようです。~~

というのは誤りで、一週間後に気づいたのですが、家計ごとのブランドごとの係数の事後分布の平均値をプロットするべきでした。
正しくはこちらです。

事前情報として正規分布を仮定していましたが、係数に関して正規分布に従っていません。
そのため、事前情報として対称性のあるような正規分布を扱うのは適切ではなさそうです。

おわりに

2005年の本とは言え、十分に使いみちのある本だと思いました。まだまだ扱いきれていないですが、引き続き勉強していきます。
この本にはケーススタディが5つほどあるのですが、それのstanコード化などをしていけばかなり力がつくような気がします。

マーケティングの部署で働くデータアナリストにとって、マーケティング×ベイズの話は非常にモチベーションの上がるところなので、こういう文献を今後も見つけていきたい。

参考文献

Bayesian Statistics and Marketing (Wiley Series in Probability and Statistics)
Bayesian Statistics and Marketingのサポートサイト
 ベイズモデリングによるマーケティング分析
 StanとRでベイズ統計モデリング (Wonderful R)
RStanのおさらいをしながら読む岩波DS 1 Shinya Uryu
Stanのlp__とは何なのかうなどん
 ‘LP__’ IN STAN OUTPUT
Package ‘bayesm’

蒙古タンメン中本コーパスに対しての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

「NOSQLの基礎知識」を読んで基礎知識を養う

同じ職場のエンジニアに『NOSQLの基礎知識』はいいぞ〜と勧められて読もうとしているのですが、マーケターレベルだと、2章と4章を読めば良いとのことです。それ以外は実際に手を動かす人が必要な知識らしい。なので、今回は1、2、4章に書かれている内容を簡単にまとめてみようと思います。

NOSQLとは(1章)

Not Only SQLの略で、SQLではないと言う意味。
RDBは単体のハードウェア上で利用するには最適だが、複数台並べてデータを分散して管理するのが不得意だがNOSQLはそれが得意。
SQLは膨大かつ多種多様で激しく変化するデータの塊を制限なしに保存することが難しいが、NOSQLはそれを可能にする。

NOSQLの特徴(2章)

機能がRDBほど豊富ではない
Joinなどの演算子もサポートされていない。（今はMongoDBではlookupとかいう機能が付いている模様）
データの整合性がゆるい
結果整合性で良いという考え
一時的なデータの不整合は問題としない。
NOSQLは正規化、更新書き込み量の最小化、メモリキャッシュの利用などの専門性が不要。
SQLとNOSQLは補完関係にある。
SQLは豊富な機能とデータの整合性に強いが、ビッグデータに弱い。NOSQLはその逆。

NOSQLの種類(2章)

世界に100種類以上存在
ビッグデータに関連するカテゴリは以下の4つ
・キーバリュー
・カラム指向
・ドキュメント指向
・グラフ型
一般的なRDBなどのSQL
将来どのようなデータが発生するか想定の範囲内で設計が行われる。
キーバリュー
・キーとバリューの対が行として格納される。（データモデルがシンプル）
・キーを指定すればバリューを探し出せるのでサーバが何台増えても対応でき、スケールアウトに適している。
・複数のサーバに複数のバリューを複製した場合、複数のバージョンが存在する可能性がある。
・具体的なNOSQL：Dynamo,Voldemort,Riak,Hibari,Redis,Scalaris,Tokyo Cabinet/Tyrant
カラム指向
・キーの付された行が複数のカラムを持つことを許容している。
・あるカラムにおいて全てにデータが格納されていなくても良い。
・カラム数を後から増やしても良い。
・具体的なNOSQL：Bigtable,Cassandra,HBase,Hypertable
ドキュメント指向
・JSONやXMLといったデータ記述形式で記述されている。
・各ドキュメントは階層構造を持たず、相互の関係が横並びに管理されている。
・具体的なNOSQL：CouchDB,MongoDB
グラフ型
・データとデータ間のつながりを管理できる。
・個々のデータの属性に依存することなく、相互のつながり方だけを独立した形で保存できる。
・ノード（頂点）、リレーションシップ（関係の有無、関係の方向）、プロパティ（属性）の３つの要素から構成。
・属性は、期間や有効/無効のデータなどを持たせることができる。
・データ間の関係性を示すグラフを作成するという目的のために開発され、最適化されている。
・具体的なNOSQL：InfiniteGraph,Neo4j

HadoopはNOSQLではない(4章)

Hadoopは複数のソフトウェアを包含するフレームワーク
HDFSは大量のデータを複数のサーバに分散して格納する点でNOSQLのように見えるが、MapReduceによるバッチ処理前提となっており、リアルタイムでの応答の速さを求められる対話型のNOSQLとは分けて捉える必要がある。

MapReduce処理(4章)

MapReduceは大規模なデータを複数のノードで並列分散処理するためのプログラミングパターン
データ処理をMapタスク（膨大な元データを分解し必要な情報を抽出し有用な形に変換して出力）とReduceタスク（抽出された情報を集約し一塊のデータとして出力）の2段階のフェーズで分けて行う。
多くのNOSQLにもMapReduceは実装されている。

さいごに

これまでR、Python、SQLしか心得ていなかった自分の知識の狭さを、エンジニアと一緒に働くことで実感しました。一冊読むくらいでエンジニアの話に少しはついていけると思うと読書の幅は広げておきたいですね。今後、データを活用したマーケティング支援ツールの開発を先導する際に、技術選定などの意思決定に役立つかもしれません。データアナリストという職種である者からすると、マーケティング・エンジニアリング・データサイエンスを緩急つけて学んでくというスタンスが重要なのかもと思う今日この頃です。

参考文献

NOSQLの基礎知識 (ビッグデータを活かすデータベース技術)

データアナリストもLinuxについて学んでみる

はじめに

もともとデータアナリストとして働いているのですが、最近、エンジニアと一緒に仕事をすることが増えたので、良いとされている教科書を読むなどしています。その中で、『新しいLinuxの教科書』がわかりやすく書かれていたので、今回はそこで出てきたコマンドたちを箇条書きにして忘備録としておきたいと思います。なのでいつもよりショボいです。
普段、Rをメインとして機械学習・クローリング目的でPythonを扱うものとしては、エンジニア寄りの知識が薄いなぁと課題に感じていたので、これを良い機会と思ってCLIへの耐性を付けていきたいと思います。

感想

分析業務において、仮想環境でCentOSやUbuntuを扱うことはあるのですが、分析環境を作った上で、クローリング業務・RPAシステムなどのスケジュール実行をするくらいで基礎からちゃんと勉強したわけではありませんでした。そのため、教科書には知らないコマンドもありました。
この本はLinuxを使って開発する際の作法（移植性大事、GUI使うな、十字キー使うな）から、歴史的なコマンドの変遷、シェルスクリプトを用いた作業の効率化、Gitなどのチーム開発ツールの紹介まで幅広く扱われていて十分な内容だと思いました。図による解説もあるので、普段R使っているような非エンジニアでも読めるはずです。これからエンジニアと一緒に分析業務をしていくデータアナリストの方には強くおすすめできる本だと思います。Rでの当たり前の機能が、Linuxのシェルで普通にあるんだと色々と学びがありました。

コマンドたち

date

現在の日時の表示

echo

文字列の表示

pwd

現在のディレクトリを確認

cd

カレントディレクトリの変更

ls

ディレクトリ内のファイルの表示、オプションでパターン指定も可能

mkdir

ディレクトリを作成する、pオプションで深い階層も一気に作れる

touch

ファイルを作成する

rm

ファイル・ディレクトリを削除する

rmdir

空のディレクトリを削除する

cat

ファイルを表示する

less

ファイル内容をスクロール表示する

cp

ファイル・ディレクトリをコピーする

mv

ファイル・ディレクトリを移動する

ln

ファイルにリンクを張る

find

ディレクトリツリーからファイルを探す

locate

ファイル名データベースからファイルを探す（高速）、何日か前のファイルを探すの適している。

which

コマンドのフルパスを表示する

alias

エイリアスの設定（わかりやすい別名を与える）

set

bashのオプションを切り替える

shopt

bashのオプションを切り替える（種類豊富）

PATH

コマンド検索パス

printenv

環境変数の表示

export

環境変数の設定

sudo

コマンドをスーパーユーザとして実行

ps

プロセスの表示（x, ux, ax, aux, auxww）

jobs

現在のジョブ一覧を表示

fg

ジョブをフォアグラウンドにする（入力受付状態にする）

bg

ジョブをバックグラウンドにする

kill

ジョブ・プロセスを終了させる

|

パイプライン（処理と処理をつなぐ）
R使いからするとRのdplyrでパイプ処理を知ったので、Linuxに普通にあるのを知れて新鮮な気分になった。

head

先頭の部分を表示
こちらもRでよく使うものですね。

du

指定したファイルやディレクトリの使用容量を表示する

wc

入力ファイルの行数・単語数・バイト数を数える

sort

行単位でテキストを並び替える

uniq

同じ内容の行を省く

cut

入力の一部を切り出して出力する（元データから特定の部分を抽出したい際など）

tr

入力の文字を置き換える・削除もできる（-d）

tail

ファイルも末尾部分を表示する

diff

二つのファイルの差分を表示する

grep

文字列を検索する

sed

テキストの編集を元ファイルを変更せずに行える

awk

テキストの検索や抽出・加工などの編集操作を行う（sedより高機能で、スクリプトは「パターン」と「アクション」からなる）
・列選択もできる
　・CSVファイルからのスコア集計などもできる。
　・あんちべさんの本『データ解析の実務プロセス入門』でも紹介されていた。

xargs

標準入力として引数のリストを与える（findやgrepとの組み合わせで使うことが多いらしい）

tar

ファイルをアーカイブファイルにまとめたり、アーカイブファイルから元のファイルを取り出したりする

gzip

ファイルを圧縮・展開する（一つのファイルしか圧縮できないのでtarとの併用を行う）

bzip2

gzipよりも高い圧縮率で圧縮を行うが、gzipよりも圧縮・展開に時間がかかる

zip

アーカイブと圧縮を同時に行う（多くのディストリビューションで標準インストールされていない。unzipも必要）

git

ファイル変更履歴を保存し管理するツール
　・git init リポジトリを作成
　・git add コミット対象として登録
　・git commit コミットする
　・git status ワークツリーの状態を表示
　・git diff 差分を表示する
　・git log 履歴を表示する
　・git revert コミットを取り消す
　・git branch ブランチの一覧表示、新しくブランチを作成
　・git checkout ブランチを切り替える
　・git merge 指定したブランチを現在のブランチにマージする
　・git branch -d ブランチを削除する
　・git push 変更の履歴を送る
　・git clone リポジトリを複製する
　・git remote add リポジトリパスに別名を設定する
　・git fetch 他のリポジトリの変更を取り込む

yum

パッケージ管理（CentOS）
　・yum install パッケージのインストール
　・yum erase/remove パッケージのアンインストール、削除
　・yum search パッケージの検索
　・yum info パッケージの詳細情報を表示

apt

パッケージ管理（Ubuntu）
　・sudo apt-get install パッケージのインストール
　・sudo apt-get remove パッケージの削除
　・sudo apt-get purge 設定ファイルも含めた完全削除
　・apt-cache search パッケージを検索

ssh

リモートログインする

info

オンラインマニュアルを表示する（索引や階層構造を設定できるなどmanより充実）

参考文献

新しいLinuxの教科書

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で機械学習の予測結果を解釈してみる
 ワインの味（美味しさのグレード）は予測できるか？（１）

深層学習のマーケティング適用事例をarXivから漁る

はじめに

深層学習が流行りだして久しいですが、マーケティングにおいて適用することが全くないため、学習のモチベーションを維持する事ができておりません。そろそろアカデミアの領域でマーケティング事例での深層学習の盛り上がりがあるのではないかと期待して、調べてみようと思います。

基本的に各研究ごとに
・目的と結果
・対象となるデータ
・手法の概要
・PDFのリンク
について記していきます。リンクは下の方にある参考情報を参照ください。

調べ方

いったん、皆さんの大好きなarXiv.orgに掲載されている論文の中から、マーケティングに関わりそうなものを見つけます。見つけ方としては、Google検索において、「marketing deep learning site:https://arxiv.org」するものとします。site:コマンドを使うことでarXiv.org以外の情報を弾いて検索できます。

今回紹介する研究

Predicting online user behaviour using deep learning algorithms
Churn analysis using deep convolutional neural networks and autoencoders
Customer Lifetime Value Prediction Using Embeddings

Predicting online user behaviour using deep learning algorithms

目的と結果

大規模ECサイトのユーザー行動からの24時間以内での購買意図の予測を目的としており、ベースラインである既存の手法と比べて、深層学習の手法（Deep Belief Networks・Stacked Denoising auto-Encoders）が予測において勝り、不均衡データの分類においても秀でているという結果となった。

対象となるデータ

6ヶ月間のECサイトのユーザー行動データ
- ユーザーID
- タイムスタンプ
- イベントタイプ
  - 商品のPageView（価格などの詳細データが伴う。商品数は25,000個）
  - カゴのPageView（価格などの詳細データが伴う。）
  - 購入
  - 広告クリック（今回の分析では無視）
  - 広告ビュー（今回の分析では無視）

手法の概要

前処理
・購入までのセッション期間
・ユーザーごとの購入に対するクリックの比率
・購入までのセッション数の中央値
・商品の説明文
・商品の価格
・商品に費やされたセッション時間の合計
・セッションが発生した時刻
・セッション単位でのクリック数
・セッション単位での購入した商品の平均価格
・過去24時間でのページビュー数
・先週のページビュー数
- 購入アイテムの説明文に対して、word2vecにより50次元の分散表現を作成しその算術平均をとった。
- サイトにおいてクリック数が10回以下のユーザーのデータを除外。
- 特徴量の次元圧縮のためにNMF(Non-Negative Matrix Factorization)を適用。
アルゴリズム（評価方法はAUC、10分割交差検証。Kerasで実行。）
- Deep Belief Networks
  2006年にHintonが提案した教師なし学習アルゴリズム。制約付きボルツマンマシン（RBM）を多数重ねて、最上位層だけ無向でそれ以外は下層への有向エッジを伴うネットワーク。
  RBMよりも複雑なデータ生成の仕組みをモデル化することができるとされている。
- Stacked Denoising auto-Encoders
  AutoEncoderは出力データが入力データをそのまま再現するニューラルネットワークで、入力層から中間層への変換器をencoder、中間層から出力層への変換器をdecoderと呼ぶ。
  そのAutoEncoderを高次元に写像できるように改良したものがDenoising auto-Encodersと言われています。
  さらに、そのDenoising auto-Encodersを層ごとに学習し積み重ねて教師ありの学習器を繋いだものをStacked Denoising auto-Encodersと呼びます。
- 決定木
- ランダムフォレスト
- Drop outによる正則化
  過学習を避けるために、ランダムに隠れ層や入力層を除いている。

PDFのリンク

Predicting online user behaviour using deep learning algorithms

Churn analysis using deep convolutional neural networks and autoencoders

目的と結果

電話通信業におけるユーザーの解約予測を行うことを目的とし、深層学習を用いた予測の結果、既存手法よりも精度が高いことがわかった。
また、解約予測だけでなく、解約理由を探るに際しての可視化においてもAutoencoderを用いることで理解が深まることが示されている。

対象となるデータ

ユーザー情報（おそらく携帯電話会社であるタイのトゥルー・コーポレーションのデータ）
- データ使用量
  モバイルのプリペイド補充金額
  モバイルのプリペイド補充頻度
  音声通話
  通話時間
  SMSメッセージ
  など

手法の概要

ユーザーの情報を画像のように扱う。各々の行が日にちを表し、各々の列がユーザー行動の種類を表す。
二値変数の定義
解約の定義は30日間行動が観察されなかったかどうか。（解約率は3.57%）
最後の通話は過去14日間での通話として、最後の14日間に記録のなかったユーザーは分析対象から除外している。プロアクティブな解約予測をするに祭して、すでに解約してしまったと思われる人を除外するということとなる。
ユーザー行動を画像のように変換したデータに対してCNN(Convolutional Neural Network)を実行しています。
CNNの構成としては、
・2連続のConvolutional Layer
・2×1の最大Pooling Layer
・128個のFully Connected Layer
・バイナリーな結果を返すsoftmax output
となっています。
訓練データとテストデータは80:20で分割
Docker上で環境を構築
オープンソースライブラリとしてはTheano、TensorFlow、Kerasを利用。
ユーザー理解のためにAutoencoderを適用している。

PDFのリンク

Churn analysis using deep convolutional neural networks and autoencoders

Customer Lifetime Value Prediction Using Embeddings

目的と結果

手探りで特徴量を構築してきた、既存のCustomer Lifetime Value（CLTV）の算出手法に勝るモデルを作ることを目的とし、十分なデータがあれば、既存の特徴量よりも予測パフォーマンスにおいて秀でていることが示されている。しかしながら、ランダムフォレストを超えるレベルの精度を追求するとニューロンの数を大幅に増やす必要があり、計算コストがかかるため商業的なバリューにおいて懸念がある。

対象となるデータ

ユーザー情報（アパレルECサイトのASOS.comのデータ）
注文数
注文日の標準偏差
直近四半期のセッション数
売り上げ
など

手法の概要

特徴量の生成
word2vecを行うに際して、単語の代わりとしてユーザーの商品ビューデータを入力としてニューラルネットワークによる分散表現の生成を行なっている。出力層において、全てのユーザーへのノードがあり、ユーザー数が1億人近くいるため計算が膨大になってしまうが、Negative Sampling(NE)によるSkipGramを行うことで、計算の高速化（近似による計算コストの低減）が図られている。分散表現の次元は32-128次元までが望ましい傾向を示している。
Hybridモデルの適用
既存の特徴量に分散表現を加えた、ロジスティック回帰とのHybridなDNNによるモデルの適用を行なっている。
CLTVの研究ながらも、予測し検証する目的変数は解約（12ヶ月購買なし）の有無。
検証用のデータは5万人のユーザーデータ
ベースラインとの比較となる評価指標はAUCが使われている。ベースラインはロジスティック回帰とランダムフォレスト。
オープンソースライブラリとしてはTensorFlowを利用。GPUはTesla K80を利用。

PDFのリンク

Customer Lifetime Value Prediction Using Embeddings

まとめ

ユーザーの行動ログに対して、縦を時間、横にカテゴリとして画像データのように扱うアプローチを知れた。
データ数が十分に多そうな企業の適用例が多かったが、扱うデータの規模感が少しわかった。
画像や音声ではなく、ユーザーデータにおいてもDNNでランダムフォレストを超える精度を叩き出せる可能性はあるが、計算コストがかかってしまい商業的なバリューが下がってしまう可能性もある。
ログや閲覧したページのコンテンツベースに生成した分散表現を特徴量とするアプローチが3例中、2例あった。

データ量の観点からまだまだ社内での適用は視野に入っていないですが、私のDNN理解も浅いので、今回のような分析事例をモチベーションにして勉強を進め適用可能性を探っていきたいと思います。

参考情報

Chainer v2による実践深層学習
 深層学習 (機械学習プロフェッショナルシリーズ)
Stacked Denoising Autoencoderを用いた語義判別
 Dropoutの実装と重みの正則化
 RBMから考えるDeep Learning　～黒魔術を添えて～
Convolutional Neural Networkとは何なのか

R Advent Calendar 2017 rvestを用いてポケモンデータをスクレイピング&分析してみた

R Advent Calendar 2017の11日目を担当するMr_Sakaueです。
今回はrvestパッケージを用いて、友人がハマっているポケモンの情報を集めてみようと思います。
もっとも、業務でWebスクレイピングする際はPythonでBeautifulSoupやSeleniumを使うことがほとんどなのですが、たまにはRでやってみようと思います。

目次
・やりたいこと
・rvestについて
・データの取得と集計と可視化と分析
・まとめ
・参考情報

やりたいこと

今回はポケモンたちのデータを集めた上で、以下の内容を行いたいと思います。

ポケモンのサイトから種族値を取得
ポケモンの種族値を標準化して再度ランキング
ポケモンのレア度や経験値に関する情報を取得
レア度や経験値と相関しそうな種族値を探る

今回扱った全てのコードはこちらに載せております。
https://github.com/KamonohashiPerry/r_advent_calendar_2017/tree/master

※種族値はゲームにおける隠しパラメータとして設定されている、ポケモンの能力値とされている。

rvestについて

rvestはRでWebスクレイピングを簡単に行えるパッケージです。ここでの説明は不要に思われますが、今回はread_html()、html_nodes()、html_text()、html_attr()の4つ関数を用いました。

基本的に以下の3ステップでWebの情報を取得することができます。

STEP1
read_html()でHTMLからソースコードを取得する。(Pythonでいう、requestとBeautifulSoup)
STEP2
html_nodes()でソースコードから指定した要素を抽出する。(PythonでいうところのfindAll)
STEP3
html_text()やhtml_attr()で抽出した要素からテキストやリンクを抽出する。(Pythonでいうところのget(‘href’)など)

データの取得と集計と可視化

検索エンジンで検索してだいたい1位のサイトがあったので、そちらのWebサイトに載っているポケモンの種族値の一覧をスクレイピング対象とさせていただきます。

ポケモンのサイトから種族値を取得

library(rvest)
library(tidyverse)
library(magrittr)
library(reshape2)

# htmlソースコードを読み込む
pokemon_ranking <- read_html("https://yakkun.com/sm/status_list.htm")

# class属性がtdタグのノードを抽出
node_extracted <- html_nodes(pokemon_ranking, "td")

# ノードからテキストを抽出して行列にして、名前以外を数値に型変換して、変数名を変更する。
pokemon_data <- data.frame(matrix(html_text(node_extracted),
                       ncol = 9,byrow = TRUE),stringsAsFactors = FALSE) %>% 
                set_colnames(c('id', 'name', 'Hit_Points', 'Attack',
                               'Defense', 'Special_Attack',
                               'Special_Defense', 'Speed', 'Total')) %>%
                mutate_at(vars(-name), as.numeric)


# ポケモン別のページを取得するためのURLの取得
pokemon_link <- pokemon_ranking %>% html_nodes("td") %>% html_nodes('a') %>% html_attr('href')
pokemon_link <- gsub(x = pokemon_link,
                     pattern = './zukan',
                     replacement = "https://yakkun.com/sm/zukan")
# ポケモン別のURLを先ほどのデータに加える
pokemon_data <- pokemon_data %>% mutate(url = pokemon_link)

library(rvest)

library(tidyverse)

library(magrittr)

library(reshape2)

# htmlソースコードを読み込む

pokemon_ranking <- read_html("https://yakkun.com/sm/status_list.htm")

# class属性がtdタグのノードを抽出

node_extracted <- html_nodes(pokemon_ranking, "td")

# ノードからテキストを抽出して行列にして、名前以外を数値に型変換して、変数名を変更する。

pokemon_data <- data.frame(matrix(html_text(node_extracted),

ncol = 9,byrow = TRUE),stringsAsFactors = FALSE) %>%

set_colnames(c('id', 'name', 'Hit_Points', 'Attack',

'Defense', 'Special_Attack',

'Special_Defense', 'Speed', 'Total')) %>%

mutate_at(vars(-name), as.numeric)

# ポケモン別のページを取得するためのURLの取得

pokemon_link <- pokemon_ranking %>% html_nodes("td") %>% html_nodes('a') %>% html_attr('href')

pokemon_link <- gsub(x = pokemon_link,

pattern = './zukan',

replacement = "https://yakkun.com/sm/zukan")

# ポケモン別のURLを先ほどのデータに加える

pokemon_data <- pokemon_data %>% mutate(url = pokemon_link)

以上のコードを実行すれば、こんな感じでポケモンの種族値一覧を得る事ができます。

とりあえず、種族値合計（Total Tribal Value 以下、TTV）のランキングの上位を確認してみます。知らないんですが、メガミュウツーとかいうイカつそうなポケモンが上位にいるようです。昭和の世代には縁のなさそうなポケモンばかりですねぇ。

■TTVランキング

取得した種族値を項目別に集計したり、Boxプロットを描いてみます。どうやら、攻撃の平均が高く、ヒットポイントや素早さの平均は低いようです。

# 集計
pokemon_data_melt <- melt(pokemon_data %>% select(-url), id.vars = 'name')
pokemon_data_melt %>% 
    group_by(variable) %>% 
    summarise(mean = mean(value),
              median = median(value),
              sd = sd(value),
              max = max(value),
              min = min(value),
              cv = sd/mean)

# 集計

pokemon_data_melt <- melt(pokemon_data %>% select(-url), id.vars = 'name')

pokemon_data_melt %>%

group_by(variable) %>%

summarise(mean = mean(value),

median = median(value),

sd = sd(value),

max = max(value),

min = min(value),

cv = sd/mean)

■種族値のサマリー

# Box-Plotを描く
ggplot(data = pokemon_data_melt %>% filter(!(variable %in% c('id','Total'))),
       aes(x = variable, y = value)) +
   geom_boxplot() + ggtitle("Tribal Value") + 
   theme(plot.title = element_text(hjust = 0.5)) + coord_flip()

# Box-Plotを描く

ggplot(data = pokemon_data_melt %>% filter(!(variable %in% c('id','Total'))),

aes(x = variable, y = value)) +

geom_boxplot() + ggtitle("Tribal Value") +

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

■種族値のBoxプロット

ポケモンの種族値を標準化して再度ランキング

さて、攻撃の平均が高かったり、ヒットポイントと素早さの平均が低かったりしたので、各々の項目を標準化した上で、再度ランキングを作ってみたいと思います。

pokemon_data_standardized <- pokemon_data
pokemon_data_standardized <- pokemon_data_standardized %>% 
                              mutate_at(vars(Hit_Points,
                                             Attack,
                                             Defense,
                                             Special_Attack,
                                             Special_Defense,
                                             Speed),funs(scale(.) %>% as.vector))

pokemon_data_standardized <- pokemon_data_standardized %>% 
                              mutate(Total = rowSums(select(.,c(3:8))))

pokemon_data_standardized_ranking <- pokemon_data_standardized %>% arrange(desc(Total))
pokemon_data_standardized_ranking <- pokemon_data_standardized_ranking %>% mutate(standardized_ranking = 1:n())


pokemon_data_standardized_melt <- melt(pokemon_data_standardized %>% select(-url), id.vars = 'name')
pokemon_data_standardized_melt %>% 
  group_by(variable) %>% 
  summarise(mean = mean(value),
            median = median(value),
            sd = sd(value),
            max = max(value),
            min = min(value),
            cv = sd/mean)

pokemon_data_standardized <- pokemon_data

pokemon_data_standardized <- pokemon_data_standardized %>%

mutate_at(vars(Hit_Points,

Attack,

Defense,

Special_Attack,

Special_Defense,

Speed),funs(scale(.) %>% as.vector))

pokemon_data_standardized <- pokemon_data_standardized %>%

mutate(Total = rowSums(select(.,c(3:8))))

pokemon_data_standardized_ranking <- pokemon_data_standardized %>% arrange(desc(Total))

pokemon_data_standardized_ranking <- pokemon_data_standardized_ranking %>% mutate(standardized_ranking = 1:n())

pokemon_data_standardized_melt <- melt(pokemon_data_standardized %>% select(-url), id.vars = 'name')

pokemon_data_standardized_melt %>%

group_by(variable) %>%

summarise(mean = mean(value),

median = median(value),

sd = sd(value),

max = max(value),

min = min(value),

cv = sd/mean)

■標準化した種族値のサマリー

平均0、分散1にできているようです。

# Box-Plotを描く
ggplot(data = pokemon_data_standardized_melt %>% filter(!(variable %in% c('id','Total','Total_standardized'))),
       aes(x = variable, y = value)) +
  geom_boxplot() + ggtitle("Tribal Value") + 
  theme(plot.title = element_text(hjust = 0.5)) + coord_flip()

# Box-Plotを描く

ggplot(data = pokemon_data_standardized_melt %>% filter(!(variable %in% c('id','Total','Total_standardized'))),

aes(x = variable, y = value)) +

geom_boxplot() + ggtitle("Tribal Value") +

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

■標準化した種族値のBoxプロット

他よりも低かったヒットポイントと、高かった攻撃がならされていることが確認できます。

■標準化前後でのTTVランキングのギャップが大きかったものをピックアップ

ラッキーが144位ほど出世しています。攻撃が低く、ヒットポイントの高いラッキーが標準化により優遇されるようになったと考える事ができます。ポケモン大会の上位ランカーである後輩社員もラッキーは手強いですと言っていたのでまんざらでもないのでしょう。

ポケモンのレア度や経験値に関する情報を取得

今回のサイトには、個別にポケモン別のページが用意されており、そちらから、ゲットしやすさや経験値に関する情報を抽出します。

# ポケモンの個別ページの情報を格納するデータフレームの作成
pokemon_detail_database <- data.frame(url = as.character(),
                                      name = as.character(),
                                      rarity = as.integer(),
                                      experience = as.integer())

# ポケモン別のURLからゲットしやすさなどを抽出するための関数
Pokemon_Detail_Get <- function(pokemon_url){
  pokemon_detail <- read_html(pokemon_url)
  
  # XPathで名前とゲットしやすさと経験値タイプを取得
  node_extracted_pokemon_name <- pokemon_detail %>% html_nodes(xpath="//tr[1]") %>% html_text()
  node_extracted_pokemon_name <- node_extracted_pokemon_name[1]
  
  node_extracted_pokemon_get <- pokemon_detail %>% html_nodes(xpath="//tr[24]/td[2]") %>% html_text()
  node_extracted_pokemon_get <- as.integer(gsub(x = node_extracted_pokemon_get[1], pattern = "\u00A0", replacement = ""))
  
  node_extracted_pokemon_exp <- pokemon_detail %>% html_nodes(xpath="//tr[26]/td[2]") %>% html_text()
  node_extracted_pokemon_exp <- as.integer(gsub(x = node_extracted_pokemon_exp[1], pattern = "万", replacement = "0000"))
  
  pokemon_detail_data <- data.frame(url = pokemon_url,
                                    name = node_extracted_pokemon_name,
                                    rarity = node_extracted_pokemon_get,
                                    experience = node_extracted_pokemon_exp)
  return(pokemon_detail_data)
  
  Sys.sleep(30)
}

# ポケモン別のページをスクレイピングする
pokemon_detail_database <- map_dfr(pokemon_link ,
                               ~Pokemon_Detail_Get(.))

# 重複したURLを削除する
pokemon_detail_database <- pokemon_detail_database %>% distinct(url, .keep_all = TRUE)

# 種族値のデータとゲットしやすさなどのデータを繋ぎこむ
pokemon_data_standardized <- pokemon_data_standardized %>% left_join(pokemon_detail_database %>% select(-name), by ="url")

# ポケモンの個別ページの情報を格納するデータフレームの作成

pokemon_detail_database <- data.frame(url = as.character(),

name = as.character(),

rarity = as.integer(),

experience = as.integer())

# ポケモン別のURLからゲットしやすさなどを抽出するための関数

Pokemon_Detail_Get <- function(pokemon_url){

pokemon_detail <- read_html(pokemon_url)

# XPathで名前とゲットしやすさと経験値タイプを取得

node_extracted_pokemon_name <- pokemon_detail %>% html_nodes(xpath="//tr[1]") %>% html_text()

node_extracted_pokemon_name <- node_extracted_pokemon_name[1]

node_extracted_pokemon_get <- pokemon_detail %>% html_nodes(xpath="//tr[24]/td[2]") %>% html_text()

node_extracted_pokemon_get <- as.integer(gsub(x = node_extracted_pokemon_get[1], pattern = "\u00A0", replacement = ""))

node_extracted_pokemon_exp <- pokemon_detail %>% html_nodes(xpath="//tr[26]/td[2]") %>% html_text()

node_extracted_pokemon_exp <- as.integer(gsub(x = node_extracted_pokemon_exp[1], pattern = "万", replacement = "0000"))

pokemon_detail_data <- data.frame(url = pokemon_url,

name = node_extracted_pokemon_name,

rarity = node_extracted_pokemon_get,

experience = node_extracted_pokemon_exp)

return(pokemon_detail_data)

Sys.sleep(30)

}

# ポケモン別のページをスクレイピングする

pokemon_detail_database <- map_dfr(pokemon_link ,

~Pokemon_Detail_Get(.))

# 重複したURLを削除する

pokemon_detail_database <- pokemon_detail_database %>% distinct(url, .keep_all = TRUE)

# 種族値のデータとゲットしやすさなどのデータを繋ぎこむ

pokemon_data_standardized <- pokemon_data_standardized %>% left_join(pokemon_detail_database %>% select(-name), by ="url")

以上のコードを実行すれば、やや時間がかかりますが、全ポケモンのゲットしやすさや経験値のデータを抽出する事ができます。それらの情報がゲットできたら、まずは可視化します。

# ゲットしやすさのヒストグラム
ggplot(data = pokemon_data_standardized, aes(x = rarity)) + geom_histogram() 

# 経験値のヒストグラム
ggplot(data = pokemon_data_standardized, aes(x = experience)) + geom_histogram()

# ゲットしやすさのヒストグラム

ggplot(data = pokemon_data_standardized, aes(x = rarity)) + geom_histogram()

# 経験値のヒストグラム

ggplot(data = pokemon_data_standardized, aes(x = experience)) + geom_histogram()

■ゲットしやすさのヒストグラム

ゲットのしやすさは、小さいほど捕まえる難易度が高くなっています。難易度の高いポケモンである0が多過ぎるので、このデータは欠損値が0になっているのではないかと疑われます。

■経験値のヒストグラム

経験値は、レベル100になるまでに要する経験値をさしています。ほとんどが100万程度となっているようです。

■ゲットしやすさと標準化TTVの散布図

やはり、ゲットしやすさに関してはデータに不備があるようで、コラッタ（アローラの姿）のような雑魚ポケのゲットのしやすさが0だったり、伝説のポケモンであるネクロズマが255だったりします。ただ、上限と下限のデータを間引けば右下がりの傾向が見られそうです。

■経験値と標準化TTVの散布図

経験値が多く必要にも関わらず、TTVが低い集団があります。どうやらこの集団に属するのは、「キノガッサ」・「マクノシタ」・「イルミーゼ」・「ゴクリン」・「シザリガー」などで、一回しか進化しないポケモンのようです。これらのポケモンは育てにくく、TTVの低い、コスパの悪そうなポケモンと考えることができるのではないでしょうか。（技や特性によってはバリューあるかもしれませんが。）

レア度や経験値と相関しそうな種族値を探る

先ほどのレア度に関しては、データがおかしそうだったので、レア度0と255に関しては除外してみます。

# おかしそうなレア度0と255のデータを除外する。
pokemon_data_standardized_filtered <- pokemon_data_standardized %>% filter(rarity > 0, rarity < 255)

# ゲットのしやすさと標準化TTV
ggplot(data = pokemon_data_standardized_filtered, aes(x = rarity, y = Total)) + 
  geom_point() + ylab('Total Tribal Value')

# おかしそうなレア度0と255のデータを除外する。

pokemon_data_standardized_filtered <- pokemon_data_standardized %>% filter(rarity > 0, rarity < 255)

# ゲットのしやすさと標準化TTV

ggplot(data = pokemon_data_standardized_filtered, aes(x = rarity, y = Total)) +

geom_point() + ylab('Total Tribal Value')

■ゲットしやすさと標準化TTVの散布図

やはり除外する事で、理想的な右下がりの傾向を示す散布図が得られたと思います。
さて、各種族値がレア度にどれだけ相関しているのかを分析したいのですが、その前にレア度を表す二項変数を作成します。

■ゲットしやすさが50以下であれば1、それ以外を0にする変数を作成

pokemon_data_standardized_filtered <- pokemon_data_standardized_filtered %>% mutate(y = ifelse(rarity <= 50, 1, 0))

1	pokemon_data_standardized_filtered <- pokemon_data_standardized_filtered %>% mutate(y = ifelse(rarity <= 50, 1, 0))

続いて、各種族値を説明変数として、レア度を目的変数としたロジスティック回帰モデルの推定をrstanで実行させます。

■stanコード

data {
  int N;
  real Hit_Points[N];
  real Attack[N];
  real Defense[N];
  real Special_Attack[N];
  real Special_Defense[N];
  real Speed[N];
  int<lower=0, upper=1> Y[N];
}


parameters {
  real b[7];
}


model {
  for (n in 1:N)
    Y[n] ~ bernoulli_logit(b[1] + b[2]*Hit_Points[n] + b[3]*Attack[n] + b[4]*Defense[n] + b[5]*Special_Attack[n] + b[6]*Special_Defense[n] + b[7]*Speed[n]);
}

data {

int N;

real Hit_Points[N];

real Attack[N];

real Defense[N];

real Special_Attack[N];

real Special_Defense[N];

real Speed[N];

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

}

parameters {

real b[7];

}

model {

for (n in 1:N)

Y[n] ~ bernoulli_logit(b[1] + b[2]*Hit_Points[n] + b[3]*Attack[n] + b[4]*Defense[n] + b[5]*Special_Attack[n] + b[6]*Special_Defense[n] + b[7]*Speed[n]);

}

■rstanでロジスティック回帰を行い、推定結果を可視化するコード

library(rstan)

N <- nrow(pokemon_data_standardized_filtered)

data <- list(N = N,
             Hit_Points = pokemon_data_standardized_filtered$Hit_Points,
             Attack = pokemon_data_standardized_filtered$Attack,
             Defense = pokemon_data_standardized_filtered$Defense,
             Special_Attack = pokemon_data_standardized_filtered$Special_Attack,
             Special_Defense = pokemon_data_standardized_filtered$Special_Defense,
             Speed = pokemon_data_standardized_filtered$Speed,
             Y = pokemon_data_standardized_filtered$y)

fit <- stan(file = 'logistic_regression.stan',
            data = data,
            seed = 1234)

summary(fit)

traceplot(fit)

source('common.R')

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

param_names <- c('mcmc', paste0('b', 1:7))
d_est <- data.frame(1:N_mcmc, ms$b)
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(pokemon_data_standardized_filtered)

data <- list(N = N,

Hit_Points = pokemon_data_standardized_filtered$Hit_Points,

Attack = pokemon_data_standardized_filtered$Attack,

Defense = pokemon_data_standardized_filtered$Defense,

Special_Attack = pokemon_data_standardized_filtered$Special_Attack,

Special_Defense = pokemon_data_standardized_filtered$Special_Defense,

Speed = pokemon_data_standardized_filtered$Speed,

Y = pokemon_data_standardized_filtered$y)

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

data = data,

seed = 1234)

summary(fit)

traceplot(fit)

source('common.R')

ms <- rstan::extract(fit)

N_mcmc <- length(ms$lp__)

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

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

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

■MCMCのシミュレーション結果のトレースプロット

どうやら収束してそうです。

■ロジスティック回帰の推定結果

見にくいので、推定結果を松浦さんの「StanとRでベイズ統計モデリング」にあるコードを用いて可視化します。

■推定結果の可視化

どうやら、0を含まない係数について見てみると、b3（攻撃）、b5（特殊攻撃）、b6（特殊防御）が高いほど、レア度が増す傾向があるようです。珍しいポケモンは攻撃が強いという傾向があると言えるのではないでしょうか。

まとめ

rvestは簡単にスクレイピングできて便利。
ポケモンデータは色々整備されてそうで今後も分析したら面白そう。
珍しいポケモンは「攻撃」、「特殊攻撃」、「特殊防御」が高い傾向がある。
経験値が必要なのにTTVの低い、コスパの悪そうなポケモンたちがいる。

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