[Stan]項目反応理論(IRT)の段階反応モデルでbaysemのアンケートデータの分析をしてみる

はじめに

stanのユーザーガイドを見ていて、項目反応理論（IRT）についての章があり気になりました。勉強会のLTなどで手法の名前をちらっと聞いたことはあったのですが、使い道について調べていませんでした。ビジネスにおける実活用もしやすそうだと思ったので、カジュアルに分析して備忘録として残したいと思います。

目次
・項目反応理論（Item Response Theory：IRT）とは
・ビジネスでの適用可能性について
・データ
・モデルの推定
・結果の解釈
・おわりに

項目反応理論（Item Response Theory：IRT）とは

関西学院大学の教授のブログによると、

項目反応理論とは、テストについての計量モデルで、問題に対する正解・不正解のデータから、問題の特性や、回答者の学力を推定するためのモデルです。

とあります。また、Wikipediaによると、TOEFLの問題の評価のために使われているそうです。

主に、バイナリーと順序変数のモデルがあるようで、以下の母数がモデルに想定されています。どちらもほぼ同じです。

回答が2値変数のモデル

2母数のロジスティックモデル
- 特性値（例えば、広告配信の満足度とか）
- 識別度母数（項目特性曲線の傾き）
- 困難度母数（項目特性曲線の切片）
- 定数

回答が順序変数のモデル（まずい < まぁまぁ < おいしい）

段階反応モデル
- 特性値（例えば、広告配信の満足度とか）
- 識別度母数（項目特性曲線の傾き）
- 困難度母数（項目特性曲線の切片）

※項目特性曲線は横軸に特性値、縦軸に質問の正答率を取ったものです。

ビジネスでの適用可能性について

顧客のアンケート結果の解釈
- 異質な集団間の得点を比較可能
- 異なる尺度間の得点を比較可能（昔のアンケートだと5段階、今のアンケートは7段階などの状況はビジネスデータでありうる。）
人事評価のバイアスの統制
- 採用面接時の個人特性の正当な評価
アンケート項目の項目削減によるアンケートコストの低減
- 各アンケート項目が理解されたかどうかを分析し、一つ一つのアンケート項目の精度を高める

データ

今回扱うデータはbaysemパッケージに入っているデータセットです。Yellow Pagesの広告プロダクトにおける満足度サーベイの回答データで、全ての回答は1から10のスケールで点数が付けられています（1がPoorで10がExcellent）。質問数は10個で、回答数は1811件です。

各質問の内容（baysemパッケージのドキュメントに載っていました。）
q1：全体の満足度

価格について
q2：競争的な価格設定
q3：昨年と同じ広告の最小値に対しての価格の引き上げ
q4：消費者の数に対しての適切な価格設定

効果について
q5：広告の購入の潜在的な影響
q6：広告を通じて自身のビジネスへの集客ができたか
q7：多くの消費者にリーチしたかどうか
q8：年間を通じて消費者に対する長期での露出があったか
q9：多くの家計やビジネスを必要としている人に届いたかどうか
q10：ビジネスを必要としている地理上のエリアに届いたかどうか

今回のIRT適用における特性値は、「広告プロダクトに関する満足度の傾向」としてみたいと思います。

モデルの推定

今回は教科書にならって以下の段階反応モデルを用います。

ここでaは識別力（広告の満足度が高まりやすいかどうか）、bは境界パラメータ（回答カテゴリ間の境界値）、θは特性（回答者がどれだけ広告に満足しているか）を表しています。Dは定数項で、以下では1とおいています。cはアンケートの回答のカテゴリ番号です。今回の例では10段階の評価が入ることになります。最後に、uは反応を、jは質問の番号を表しています。

実践ベイズモデリング -解析技法と認知モデル-

こちらの本のサポートサイトからダウンロードできるzipファイルにstanのコードやRコードがありますので、そちらを利用しています。

モデルですが、以下のような設定となっています。

data{
    int ni; // the number of record
    int nj; // the number of item
    int nc; // the number of grade
    real D; // Constant
    int<lower=1,upper=10> y[ni,nj]; // data
}

parameters{
  vector<lower=0,upper=5>[nj] a;
    ordered[nc-1] ba[nj];
    vector<lower=-4,upper=4>[ni] theta;
}

transformed parameters{
    real b[nj,nc];
    vector<lower=0,upper=1>[nc-1] pa[ni,nj];
    simplex[nc] p[ni,nj];
    for (j in 1:nj){
        for (c in 1:nc){
            if (c ==1){
                b[j,c] = ba[j,c];
            }else if (c ==nc){
                b[j,c] = ba[j,c-1];
            }else{
                b[j,c] = (ba[j,c-1]+ba[j,c])/2;
            }
        }
    }
    for (i in 1:ni){
        for (j in 1:nj){
            for (c in 1:nc-1){
                pa[i,j,c] = 1/(1+exp(-D*a[j]*(theta[i] - ba[j,c])));
            }       
        }
    }
    for (i in 1:ni){
        for (j in 1:nj){
            for(c in 1:nc){
                if (c==1){
                    p[i,j,c] = 1-pa[i,j,c];
                }else if(c==nc){
                    p[i,j,c] = pa[i,j,c-1];
                }else{
                    p[i,j,c] = pa[i,j,c-1] - pa[i,j,c];
                }
            }
        }
    }
}

model{
    for (i in 1:ni){
        theta[i] ~ normal(0,1);
        for (j in 1:nj){
            y[i,j] ~ categorical(p[i,j]);
        }
    }
    for (j in 1:nj){
        a[j] ~ lognormal(0,sqrt(0.5));
        for (c in 1:nc-1){
            ba[j,c] ~ normal(0,2);
        }
    }
}

generated quantities{
    real bg[nj,nc];
    bg = b;
}

data{

int ni; // the number of record

int nj; // the number of item

int nc; // the number of grade

real D; // Constant

int<lower=1,upper=10> y[ni,nj]; // data

}

parameters{

vector<lower=0,upper=5>[nj] a;

ordered[nc-1] ba[nj];

vector<lower=-4,upper=4>[ni] theta;

}

transformed parameters{

real b[nj,nc];

vector<lower=0,upper=1>[nc-1] pa[ni,nj];

simplex[nc] p[ni,nj];

for (j in 1:nj){

for (c in 1:nc){

if (c ==1){

b[j,c] = ba[j,c];

}else if (c ==nc){

b[j,c] = ba[j,c-1];

}else{

b[j,c] = (ba[j,c-1]+ba[j,c])/2;

}

for (i in 1:ni){

for (j in 1:nj){

for (c in 1:nc-1){

pa[i,j,c] = 1/(1+exp(-D*a[j]*(theta[i] - ba[j,c])));

}

for (i in 1:ni){

for (j in 1:nj){

for(c in 1:nc){

if (c==1){

p[i,j,c] = 1-pa[i,j,c];

}else if(c==nc){

p[i,j,c] = pa[i,j,c-1];

}else{

p[i,j,c] = pa[i,j,c-1] - pa[i,j,c];

}

model{

for (i in 1:ni){

theta[i] ~ normal(0,1);

for (j in 1:nj){

y[i,j] ~ categorical(p[i,j]);

}

for (j in 1:nj){

a[j] ~ lognormal(0,sqrt(0.5));

for (c in 1:nc-1){

ba[j,c] ~ normal(0,2);

}

generated quantities{

real bg[nj,nc];

bg = b;

}

こちらをキックするためのRコードです。

library(bayesm)
library(rstan)
library(shinystan)
library(bayesplot)
library(tidyverse)
library(gridExtra)

# Data Import -------------------------------------------------------------
data("customerSat")
dataset <- customerSat

# Kick Stan model ---------------------------------------------------------
ni <- nrow(dataset) # 分析対象者の数
nj <- ncol(dataset) # 項目数
nc <- length(table(as.factor(dataset$q1))) # 回答の種類
D <- 1 # 定数項

stan_data <- list(y = dataset,
                  nj = nj,
                  ni = ni,
                  nc= nc,
                  D = D)

par <- c("theta","ba","a","b")
war <- 2500
ite <- 5000
see <- 1234
dig <- 2
cha <- 4

fit <- stan(file = "model/graded_response_model.stan",
                data = stan_data,
                pars = par,
                verbose = F,
                seed = see,
                chains = cha,
                warmup = war,
                iter = ite)

# Diagnose ----------------------------------------------------------------

traceplot(fit)

print(fit, pars = par, digits_summary = dig)

summary_table <- data.frame(summary(fit)$summary)
ggplot(data = data.frame(Rhat = summary_table$Rhat), aes(Rhat)) + geom_histogram()

library(bayesm)

library(rstan)

library(shinystan)

library(bayesplot)

library(tidyverse)

library(gridExtra)

# Data Import -------------------------------------------------------------

data("customerSat")

dataset <- customerSat

# Kick Stan model ---------------------------------------------------------

ni <- nrow(dataset) # 分析対象者の数

nj <- ncol(dataset) # 項目数

nc <- length(table(as.factor(dataset$q1))) # 回答の種類

D <- 1 # 定数項

stan_data <- list(y = dataset,

nj = nj,

ni = ni,

nc= nc,

D = D)

par <- c("theta","ba","a","b")

war <- 2500

ite <- 5000

see <- 1234

dig <- 2

cha <- 4

fit <- stan(file = "model/graded_response_model.stan",

data = stan_data,

pars = par,

verbose = F,

seed = see,

chains = cha,

warmup = war,

iter = ite)

# Diagnose ----------------------------------------------------------------

traceplot(fit)

print(fit, pars = par, digits_summary = dig)

summary_table <- data.frame(summary(fit)$summary)

ggplot(data = data.frame(Rhat = summary_table$Rhat), aes(Rhat)) + geom_histogram()

処理時間としては、2014年末モデルのMacbook Proのcorei5、メモリ8GBで数時間程度でした。（正確な時間はわかりませんが、寝て起きたら計算が終わっていました。）

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

Rhatも1.1未満におさまっています。

結果の解釈

# histogram of theta
theta <- rstan::extract(fit)$theta %>% apply(2,mean)
ggplot(data = data.frame(theta_mean = theta), aes(theta_mean)) + geom_histogram()

# histogram of theta

theta <- rstan::extract(fit)$theta %>% apply(2,mean)

ggplot(data = data.frame(theta_mean = theta), aes(theta_mean)) + geom_histogram()

まず、推定した特性値の値のユーザーごとの平均値を求めて、ヒストグラムを描いてみます。どうやら、上限周辺にやたらと高い評価をしてそうなユーザーがいるようです。

最後に、項目特性曲線を質問ごとに、そして回答ごとに描いてみようと思います。

# Visualization -----------------------------------------------------------
alpha <- rstan::extract(fit)$a %>% apply(2,mean)

for (i in 1:nj){
  eval(parse(text=paste0("beta_",i," <- rstan::extract(fit)$b[1:10000,",i,",1:nc] %>% apply(2,mean)")))
}

for (i in 1:nj){
  eval(parse(text=paste0("ggdf_",i," <- data.frame(matrix(ncol = nc,nrow = length(theta))) ; colnames(ggdf_",i,") <- 1:nc")))
}

## probability
for(i in 1:nj){
  for (j in 1:nc){
    eval(parse(text=paste0("ggdf_",j,"[,",i,"] <- 1/(1+exp(-alpha[",i,"]*(theta-beta_",j,"[",i,"])))")))
    eval(parse(text=paste0("ggdf_",j,"$theta <- theta")))
  }
}

## gather
for (i in 1:nj){
  eval(parse(text=paste0("ggdf_gt_",i," <- ggdf_",i," %>% tidyr::gather(key=var,value,-theta,factor_key=TRUE)")))
}

## ggplot
for (i in 1:nj){
  eval(parse(text=paste0("p",i," <- ggplot(data = ggdf_gt_",i,", aes(x = theta, y = value, colour = var)) + geom_line() + ggtitle(\"Q",i,"\")")))
}

# 2×3でグラフを描画
grid.arrange(p1, p2, p3, p4, p5, nrow = 3)
grid.arrange(p6, p7, p8, p9, p10,  nrow = 3)

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

alpha <- rstan::extract(fit)$a %>% apply(2,mean)

for (i in 1:nj){

eval(parse(text=paste0("beta_",i," <- rstan::extract(fit)$b[1:10000,",i,",1:nc] %>% apply(2,mean)")))

}

for (i in 1:nj){

eval(parse(text=paste0("ggdf_",i," <- data.frame(matrix(ncol = nc,nrow = length(theta))) ; colnames(ggdf_",i,") <- 1:nc")))

}

## probability

for(i in 1:nj){

for (j in 1:nc){

eval(parse(text=paste0("ggdf_",j,"[,",i,"] <- 1/(1+exp(-alpha[",i,"]*(theta-beta_",j,"[",i,"])))")))

eval(parse(text=paste0("ggdf_",j,"$theta <- theta")))

}

## gather

for (i in 1:nj){

eval(parse(text=paste0("ggdf_gt_",i," <- ggdf_",i," %>% tidyr::gather(key=var,value,-theta,factor_key=TRUE)")))

}

## ggplot

for (i in 1:nj){

eval(parse(text=paste0("p",i," <- ggplot(data = ggdf_gt_",i,", aes(x = theta, y = value, colour = var)) + geom_line() + ggtitle(\"Q",i,"\")")))

}

# 2×3でグラフを描画

grid.arrange(p1, p2, p3, p4, p5, nrow = 3)

grid.arrange(p6, p7, p8, p9, p10, nrow = 3)

質問1~10に関して、10段階の回答ごとの項目反応曲線を以下に描いています。上まで戻るのが面倒なので、質問内容を再掲します。

q1：全体の満足度
q2：競争的な価格設定
q3：昨年と同じ広告の最小値に対しての価格の引き上げ
q4：消費者の数に対しての適切な価格設定
q5：広告の購入の潜在的な影響
q6：広告を通じて自身のビジネスへの集客ができたか
q7：多くの消費者にリーチしたかどうか
q8：年間を通じて消費者に対する長期での露出があったか
q9：多くの家計やビジネスを必要としている人に届いたかどうか
q10：ビジネスを必要としている地理上のエリアに届いたかどうか

これらの傾向から、9〜10点を獲得するにはある程度は特性値が高まる必要がある質問としては、q1〜q6のように見えます。価格や購買など自身のビジネスに直結しそうな質問が多い印象です。逆にふわっとした質問であるq7~q10は特性値が低くても9〜10点を取れる可能性が高い傾向があります。

おわりに

Stanのユーザーガイドを読むことで、普段自分が業務で扱っているアプローチなどが如何に限定的であることが実感できました。今回はIRTのアンケートデータへの適用事例を知れ、そこから様々な文献や便利なコードに至ることができました。社内のアンケートデータへの適用は面白そうだと思いますので業務で使ってみようと思います。

参考情報

[1] 豊田秀樹 (2017) 『実践ベイズモデリング -解析技法と認知モデル-』朝倉書店
[2] Yoshitake Takebayashi (2015) 「項目反応理論による尺度運用」 SlideShare
[3] 持主弓子・今城志保 (2011) 「IRTの組織サーベイへの応用」
[4] 清水裕士 (2017) 「項目反応理論をStanで実行するときのあれこれ」 Sunny side up!
[5] 清水裕士 (2016) 「Stanで多次元項目反応理論」 Sunny side up!
[6] 小杉考司 (2013) 「項目反応理論について」
[7] Daniel C. Furr et al. (2016) “Two-Parameter Logistic Item Response Model”
[8] daiki hojo (2018) “Bayesian Sushistical Modeling” Tokyo.R#70
[9] abrahamcow (2017) 「［RStan］項目反応理論の応用でフリースタイルダンジョン登場ラッパーの強さをランキングしてみた」

Python/Rもくもく会をプライベートで開催するための参考図書・資料をまとめる

はじめに

社内で定時後に有志で勉強会というか、その場に集まってPythonやRをもくもくと勉強をするもくもく会を開きたいと考えています。目的としては分析スキルの向上や機械学習ができるようになりたいとかいう個々人の願いを叶えることです。
色々なスキルレベルのメンバーが参加することが予想されるので、皆を幸せにするためにもレベルに応じた良い教材が必要だと思いました。
ここでは、レベルに応じて適切な教材などを忘備録として残していきたいと思います。
（私自身、全てのレベルの対象者に適切な教材を網羅しているわけではないので、随時更新していこうと思います。）

受講対象について

受講対象（PythonやRをまともに触ったことがない人）は2軸で分けるとすると以下のようになると思います。

・プログラミング経験あり/経験なし
・数学の心得あり/心得なし

プログラミング経験なし&数学の心得あり（アルキメデス）
理系出身の人がメインだと思います。学部・学科によっては全然扱わないですよね。数的な思考は得意だが、それを活かすスキルが不足しているような人でしょう。眼の前におかれた数学の問題を紙とペンで解くことはできるが、仕事で使えないという感じ。私も偉そうなことは言えないですが、コードが荒れがちなので周りに良い先生がいたほうが良いと思います。
プログラミング経験なし&数学の心得なし（葉っぱ隊）
一番習得に時間がかかると思います。野球やったことないのに、野球選手になりたいという人に皆さんは違和感を感じるでしょう。イメージはそんな感じです。一番時間がかかるからこそ、挫折しないための教材選びが重要かもしれません。スキル的に全裸なので、葉っぱ隊と名付けましょう。
プログラミング経験あり&数学の心得あり（デーサイ候補）
最も頼もしい存在です。教科書をお渡ししておけば勝手に成長すると思います。ある程度経験を積めば分析業務を任せても良いと思います。
プログラミング経験あり&数学の心得なし（進捗ありマン）
各種手法の原理を知るまではそれなりに時間がかかると思いますが、手を動かして何ができるかをすぐに味わえるので、モチベーションを維持しながら学んでいきやすいと思います。コード自体は実行できるので進捗ありマンと名付けてみましょう。

この2軸でPythonとRに関する便利な資料を探したいと思います。
ただし、どの本に関してもどのレベルの人が買っても良いとは思います。ただ、数学の心得がない中で、テイラー展開とか平均値の定理とかラグランジュ未定乗数法などの表現を目にした際に、挫折してしまう可能性があるので、適した書籍から順次広げていくのが良いと思います。なお、今回はPCでもくもくと進めれそうな書籍を選んでいます。紙とペンで進める本も重要なのですが、そのようなかた向けの書籍は取り上げていません。

アルキメデス向けの教材

Python

機械学習のエッセンス -実装しながら学ぶPython,数学,アルゴリズム- (Machine Learning)
Pythonの入門的な記述もあり、数学の話もきちんと扱っているので、アルキメデス系の人にはぴったりの本だと思います。コードが荒れがちなアルキメデス系もこの本を写経しながら学べば進捗が良いはずです。
Pythonで学ぶあたらしい統計学の教科書 (AI & TECHNOLOGY)
数式はそれほど出てきませんが、Pythonが少しわかれば統計学に関してもあらかたこの本で上達すると思います。手を動かすことで前には進めると思います。アルキメデス系であれば、これを学習した上で、さらに深く理解するために他の統計学の専門書に手を出すこともできるでしょう。

R

みんなのR 第2版
『Rによるデータサイエンス』と迷ったのですが、プログラムの実行結果がそのまま載っている印象だったので、こちらの本がプログラミング初心者には優しいと判断しました。ほとんど数式は出てこないのですが、一般化線形モデルや時系列解析などもカバーしてくれています。また、データの前処理に関する記述もこちらの本の方が手厚いです。

葉っぱ隊向けの教材

Python

Pythonスタートブック [増補改訂版]
本当にプログラミングがはじめての人向けの本です。まずはプログラミング自体に慣れたほうが良いと思います。
プロゲートのPython入門講座
妻におすすめされた講座です。無料枠でもある程度学びがあるようです。環境を構築しなくても良いという点が非常に葉っぱ隊に適しているとのことです。

R

Rによるやさしい統計学
Rのインストールあるいは統計学の初歩のところから、応用まで幅広く説明している本です。数式はあまり出てきませんがコードが載っているので、手を動かすことができると思います。

読み物

マンガでわかる統計学回帰分析編
初めて統計学に触れた人で、挫折したくない場合の強力な助っ人と言えそうな本です。オーム社なので内容はしっかりしています。

デーサイ候補向けの教材

Python

みんなのPython 第4版
他のプログラミング言語の経験があれば、どのような記述の仕方なのかを掴み取るために効率的な本だと思います。記法やメソッドなどを一通り学べば、それに続くデータ分析の勉強も捗るはずです。私は一番始めに『初めてのPython 第3版』という700ページを超える本を読んでいたのですが、いきなりあれは辛いです。
[第2版]Python 機械学習プログラミング達人データサイエンティストによる理論と実践 (impress top gear)
数式による解説が適宜なされ、Pythonのコードも多く載せられており、ほとんどの手法に関して学ぶことができる本です。

R

データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC (確率と情報の科学)
実ビジネスではロジスティック回帰などを多用する場面が多いので、数学とプログラミングに問題がなければこの本も平行して学べると思います。
データマイニング入門
業務で使いそうなデータマイニングの手法が幅広く紹介されています。一部ではありますが数式もあるので納得もしやすいと思います。
Rによるデータサイエンスデータ解析の基礎から最新手法まで
数式はそこまで多くないですが、様々な手法を扱っている上に、Rのコードも載せられているので、原理を理解しながら、コードも回せる良い本だと思います。ただ、書かれている多くは金先生のWebサイトに掲載されています。
データサイエンスのための統計学入門 ―予測、分類、統計モデリング、統計的機械学習とRプログラミング
そんなに数式は出てこないですが、統計学から機械学習まで幅広く扱われています。Rを用いてもくもくと進めることができると思います。あと、特徴量エンジニアリングに関するヒントも書かれていたりするので、長く使える一冊かもしれません。

進捗ありマン向けの教材

R

RStudioではじめるRプログラミング入門
プログラミング経験のある進捗ありマンであれば、R言語の扱い方をまずは知りたいだろうと思います。関数の書き方やヘルプページの使い方、オブジェクトの説明、S3の話などが詳しく書かれています。
新米探偵、データ分析に挑む
R Studioのインストール方法なども載っているので、進捗ありマンなら最初から最後まで実践できると思います。数式もほとんど出てきません。色んな分析事例をRで取り組むことで分析業務のイメージも付いてくると思います。
RユーザのためのRStudio[実践]入門−tidyverseによるモダンな分析フローの世界−
R言語について何となくつかめた進捗ありマンがモダンな記法であるtidyverseを効率よく学べる良い本です。データ整形・クロス集計・可視化がモダンな記法で書けるようになると結構楽しいと思います。

Python

Pythonによるデータ分析入門第2版 ―NumPy、pandasを使ったデータ処理
数式などは出てきませんが、Pythonでデータ分析を行う上で重要なモジュールであるpandasを詳しく扱っています。いきなり機械学習などをするよりかはpandasをマスターしてデータの前処理スキルを高めるのが良いかもしれません。進捗もあると思いますし。

今後について

そもそもPythonやRに触れたことがない人にとって、Tokyo.Rの初心者セッションは少し適していないのかなと思ったので、今回は取り上げていないですが、一通り使い方をわかってもらえたら初心者セッションの資料を使ったもくもく会も開きたいと思います。最終的にはKaggle部とかを作るとかになるのかもしれませんが、そこまで行けるか行けないか。

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

はじめに

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

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

Uplift Modelingとはなにか

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

Uplift Modelingの卑近な例え話

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

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

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

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

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

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

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

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

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

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

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

Uplift Modelingのサンプルデータ

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

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

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

tools4upliftについて

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

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

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

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

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

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

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

head(x,5)

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

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

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

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

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

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

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

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

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

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

library(tools4uplift)

library(tidyverse)

library(data.table)

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

head(x,5)

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

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

x_sub <- x[sampleNum,]

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

summary(x_sub)

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

set.seed(123);

split.data1 <- SplitUplift(data = x_sub,

p = 0.7,

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

train <- split.data1[[1]]

valid <- split.data1[[2]]

base.tm <- DualUplift(data = train,

treat = "treatment",

outcome = "conversion",

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

# baseline model for control group

base.tm[[1]]

# baseline model for treatment group

base.tm[[2]]

# predict the uplift on the validation set

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

treat = "treatment",

outcome = "conversion",

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

# evaluate the model performance

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

treat = "treatment",

outcome = "conversion",

prediction = "uplift_prediction",

nb.group = 5)

# Qini曲線の描画

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

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

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

# Qini係数の算出

QiniArea(base.tm.perf)

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

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

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

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

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

%matplotlib inline

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

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

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

source_df.describe()

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

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

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

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

train_sample_num = len(train_is_cv_list)

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

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

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

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

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

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

stat_data = []

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

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

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


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

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

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

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


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

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

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

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

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

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

%matplotlib inline

import pandas as pd

import matplotlib.pyplot as plt

from operator import itemgetter

plt.style.use("ggplot")

from sklearn.model_selection import train_test_split

from sklearn.linear_model import LogisticRegression

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

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

source_df.head(10)

source_df.describe()

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

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

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

train_is_cv_list, test_is_cv_list, train_is_treat_list, \

test_is_treat_list, train_feature_vector_df,\

test_feature_vector_df = train_test_split(is_cv_list, is_treat_list,

feature_vector_df,

train_size=0.5,

test_size=0.5,

random_state=42)

treat_model = LogisticRegression(C=0.01)

control_model = LogisticRegression(C=0.01)

train_sample_num = len(train_is_cv_list)

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

treat_feature_vector_list = train_feature_vector_df[train_is_treat_list]

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

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

treat_model.fit(treat_feature_vector_list, treat_is_cv_list)

control_model.fit(control_feature_vector_list, control_is_cv_list)

treat_score = treat_model.predict_proba(test_feature_vector_df)

control_score = control_model.predict_proba(test_feature_vector_df)

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

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

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

treat_uu = 0

control_uu = 0

treat_cv = 0

control_cv = 0

treat_cvr = 0.0

control_cvr = 0.0

lift = 0.0

stat_data = []

for is_cv, is_treat, score in result:

if is_treat:

treat_uu += 1

if is_cv:

treat_cv += 1

treat_cvr = treat_cv / treat_uu

else:

control_uu += 1

if is_cv:

control_cv += 1

control_cvr = control_cv / control_uu

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

lift = (treat_cvr - control_cvr) * treat_uu

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

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

quantile_data = []

for n in range(10):

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

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

quantiled_result = result[start:end]

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

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

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

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

treat_cvr = treat_cv / treat_uu

control_cvr = control_cv / control_uu

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

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

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

qdf.plot.bar()

plt.xlabel("percentile")

plt.ylabel("conversion rate")

df = pd.DataFrame(stat_data)

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

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

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

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

plt.xlabel("uplift score rank")

plt.ylabel("conversion count")

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

plt.xlabel("uplift score rank")

plt.ylabel("conversion rate")

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

plt.xlabel("uplift score rank")

plt.ylabel("lift count")

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

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

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

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

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

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

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

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

おわりに

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

参考文献

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

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

はじめに

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

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

バンディット問題とは

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

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

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

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

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

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

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

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

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

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

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

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

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

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

ε-貪欲法

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

UCB(Upper Confidence Bound)方策

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

KL-UCB

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

MED(Minimum Empirical Divergence)方策

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

トンプソン抽出

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

Contextual Bandit問題とは

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

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

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

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

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

サンプル実行

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

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

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

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

library(contextual)

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

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

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

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

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

bandit             <- BasicBernoulliBandit$new(prob_per_arm)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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


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

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

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

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

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

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

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

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

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

library(contextual)

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

## Simulation of the multi-armed Bandit examples in

## of "Bandit algorithms for website optimization"

## by John Miles White.

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

# Contextuals policies correctly picks one of the max arms.

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

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

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

horizon <- 250

simulations <- 5000

bandit <- BasicBernoulliBandit$new(prob_per_arm)

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

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

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

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

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

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

history <- simulation$run()

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

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

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

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

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

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

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

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

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

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

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

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

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

history <- simulation$run()

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

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

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

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

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

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

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

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

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

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

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

history <- simulation$run()

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

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

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

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

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

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

ε-貪欲法

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

・平均報酬額

　・累積報酬額

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

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

・平均報酬額

・累積報酬額

UCB方策

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

・平均報酬額

・累積報酬額

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

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

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

library(contextual)
library(data.table)

# Import personalization data-set

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

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

# Convert data

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

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

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

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

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

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

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

library(contextual)

library(data.table)

# Import personalization data-set

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

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

data <- fread(url, stringsAsFactors=TRUE)

# Convert data

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

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

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

# Set simulation parameters.

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

horizon <- nrow(data)

# Initiate Replay bandit with 10 arms and 100 context dimensions

log_S <- data

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

Companion_Alone + Companion_Family + Companion_Partner")

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

# Define agents.

agents <-

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

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

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

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

# Initialize the simulation.

simulation <-

Simulator$new(

agents = agents,

simulations = simulations,

horizon = horizon

)

# Run the simulation.

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

sim <- simulation$run()

# plot the results

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

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

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

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

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

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

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

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

# Import and merge files

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

rm(movies_dat,ratings_dat)

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

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

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

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

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

rm(user_features, not_50, top_50, count_movies)

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

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

simulations <- 1
horizon     <- nrow(top_50_movies)

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

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

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

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

sim  <- simulation$run()

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

library(contextual)

library(data.table)

library(splitstackshape)

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

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

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

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

# Import and merge files

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

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

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

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

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

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

rm(movies_dat,ratings_dat)

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

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

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

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

top_50_movies <- all_movies[MovieID %in% top_50]

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

user_features <- all_movies[MovieID %in% not_50]

rm(all_movies)

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

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

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

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

user_features$total <- NULL

# Add user features to top50

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

top_50_movies <- na.omit(top_50_movies)

rm(user_features, not_50, top_50, count_movies)

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

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

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

simulations <- 1

horizon <- nrow(top_50_movies)

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

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

i.V23 + i.V24")

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

agents <-

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

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

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

simulation <-

Simulator$new(

agents = agents,

simulations = simulations,

horizon = horizon

)

sim <- simulation$run()

plot(sim,

type = "cumulative",

regret = FALSE,

rate = TRUE,

legend_position = "bottomright")

おわりに

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

参考情報

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

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

はじめに

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

ORとは

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

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

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

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

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

ORの強みとしては、

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

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

ORの手法

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

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

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

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

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

Rでの実践例

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

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

数理最適化

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

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

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

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

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

library(lpSolve)

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

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

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

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

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

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

library(lpSolve)

# 目的関数の係数

f.obj <- c(18, 11)

# 制約式の左辺の係数

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

# 制約式の等号・不等号

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

# 制約式の右辺

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

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

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

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

#解を求める

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

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

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

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

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

Success: the objective function is 615

> # 解を求める

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

[1] 25 15

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

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

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

# 感度分析

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

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

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

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

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

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

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

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

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

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

library(lpSolve)

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

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

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

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

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

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

library(lpSolve)

# 目的関数の係数

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

# 制約式の左辺の係数

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

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

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

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

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

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

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

# 制約式の等号・不等号

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

"<=",

"==",

"==")

# 制約式の右辺

f.rhs <- c(25,

35,

40,

15,

20,

35,

30)

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

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

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

# 解を求める

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

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

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

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

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

Success: the objective function is 3610

> # 解を求める

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

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

ナップサック問題

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

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

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

library(lpSolve)

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

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

library(lpSolve)

#便益

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

#費用

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

#予算

budget <- 20

ans <- lp(direction = "max",

objective.in = benefit,

const.mat = cost,

const.dir = "<=",

const.rhs = budget,

all.bin = TRUE)

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

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

> print(ans)

Success: the objective function is 31

> print(ans$solution)

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

混合整数計画

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

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

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

library(Rglpk)

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

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

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

max <- TRUE

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

library(Rglpk)

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

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

1, 0, 1, 0,

1, 0, 0, 1,

0, -8, 0, 0,

0, 0, -6, 0,

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

dir <- c("<=",

"<=",

"<=")

rhs <- c(15,

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

max <- TRUE

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

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

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

$solution
[1] 8 6 0 1 1 0

$status
[1] 0

$solution_dual
[1] NA

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

$auxiliary$dual
[1] NA

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

$optimum

[1] 100

$solution

[1] 8 6 0 1 1 0

$status

[1] 0

$solution_dual

[1] NA

$auxiliary

$auxiliary$primal

[1] 14 0 0 0

$auxiliary$dual

[1] NA

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

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

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

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

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

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

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

library(orloca)

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

# Compute the minimum
sol <- distsummin(loca)

# Show the result
sol

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

library(orloca)

# A new unweighted loca.p object

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

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

# Compute the minimum

sol <- distsummin(loca)

# Show the result

sol

# Evaluation of the objective function at solution point

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

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

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

> # Show the result

> sol

[1] 4.468528 3.843755

> # Evaluation of the objective function at solution point

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

[1] 22.80259

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

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

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

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

#図にプロット

library(ggplot2)

library(png)

image <- readPNG("map_open.png")

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

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

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

g <- g + geom_point()

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

aes(x = x,y = y),

colour = "red",

size = 3)

待ち行列

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

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

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

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

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

library(tidyverse)

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

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

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

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

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

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

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

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

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

library(tidyverse)

#乱数の生成

set.seed(101)

rand_uni <- runif(500)

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

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

set.seed(5)

rand_uni_demand <- runif(500)

simulation <- data.frame(interval = itm,

x=rand_uni_demand) %>%

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

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

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

time_required=if_else(x <= 0.25, 1.2,

if_else(x <= 0.5,1.5,

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

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

start = NA)

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

if (i == 1) {

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

}

else {

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

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

}

#終了時間

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

#待ち時間

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

#開始待ち人数

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

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

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

}

#空き時間

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

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

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

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

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

最後に、空き時間です。

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

最短路問題

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

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

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

library(igraph)

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

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

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

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

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

library(igraph)

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

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

E(net)$weight <- shortpath_dataset$distance

E(net)$label <- shortpath_dataset$distance

E(net)$width <- 13

V(net)$size <- 35

plot(net)

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

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

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

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

a <- list("test")

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

b <- 0

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

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

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

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

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

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

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

}

print(a)

print(b)

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

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

> print(b)
[1] 13

> print(a)

[[1]]

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

> print(b)

[1] 13

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

巡回セールスマン問題

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

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

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

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

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

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


library(TSP)
library(proxy)

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

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

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

#図にプロット

library(ggplot2)

library(png)

library(ggrepel)

image <- readPNG("map_open.png")

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

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

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

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

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

g <- g + geom_point(size = 3)

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

aes(x = x,y = y),

colour = "red",

size = 3) + geom_text_repel()

library(TSP)

library(proxy)

#距離行列の作成

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

tsp <- TSP(data)

tsp

## use some methods

n_of_cities(tsp)

labels(tsp)

## calculate a tour

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

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

tour[1:5]

tour_length(tour)

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

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

> tour[1:5]

1 2 3 4 5

> tour_length(tour)

[1] 12.64099

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

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

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

NPV、IRR

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

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

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

[1] 0.3348214

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

[1] 0.1551911

データ包絡分析法

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

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

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

library(Benchmarking)

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

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

library(Benchmarking)

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

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

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

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

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

ポートフォリオ選択問題

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

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

編集を終えて

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

参考情報

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

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

はじめに

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

分析工程

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

データの収集

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

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

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

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

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

データの整形

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

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

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

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

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

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

library(tidyverse)

library(directlabels)

library(TSclust)

library(gghighlight)

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

trend_dataset <- data.frame()

for (i in 1:21) {

#データの読込

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

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

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

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

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

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

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

trend_dataset <- trend_dataset %>% rbind(trend_data)

}

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

trend_dataset <- trend_dataset %>%

group_by(keyword) %>%

mutate(period = 1:n())

可視化

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

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

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

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

ggplot(data = trend_dataset,

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

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

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

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

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

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

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

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

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

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

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

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

h <- hclust(d)

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

par(cex=0.6)

par(family = "HiraKakuProN-W3")

plot(h, hang = -1)

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

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

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

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

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

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


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

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

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

# クラスタ数は3とする

data.frame(cutree(h, 3))

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

clusters$keyword <- rownames(clusters)

rownames(clusters ) <- NULL

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

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

rename_keywordlist <- read_csv("rename_keywordlist.csv")

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

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

predicate = max(cluster_number) == 1)

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

predicate = max(cluster_number) == 2)

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

predicate = max(cluster_number) == 3)

まずはクラスター1

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

続いてクラスター2

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

そしてクラスター3

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

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

トレンドの推定

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

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

まずは、stanのコード

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

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

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

data {

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

int performer_num; // 芸人さんの数

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

}

parameters {

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

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

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

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

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

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

}

model {

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

for(i in 2:T) {

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

}

// ランダム効果

r ~ normal(0, s_r);

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

for(i in 1:T) {

for(j in 1:performer_num) {

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

}

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

library(rstan)
library(bayesplot)

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

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

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

mcmc_rhat(rhat(fit))

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

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

library(rstan)

library(bayesplot)

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

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

data <- list(T = T,

performer_num = performer_num ,

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

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

data = data,

seed = 1,

iter = 30000,

warmup = 10000,

thin = 10

)

mcmc_rhat(rhat(fit))

# データの整形

stan_df_1 <- fit %>%

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

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

t() %>%

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

data.frame

# 列名の変更

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

# 結果

head(stan_df_1, n = 3)

ggplot(data = trend_dataset_withcluster) +

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

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

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

geom_line(data = stan_df_1,

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

geom_ribbon(data = stan_df_1,

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

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

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

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

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

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

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

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

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

rot = bayesdfa::rotate_trends(mod_3)

names(rot)

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

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

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

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

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

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

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

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

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

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

rot = bayesdfa::rotate_trends(mod_5)

names(rot)

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

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

今後の課題

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

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

参考情報

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

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

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

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

初心者セッション

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

資料なし

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

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

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

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

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

資料なし

応用セッション

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

資料なし

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

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

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

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

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

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

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

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

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

tidyeval入門以前

Yet Another Introduction to tidyeval from yutannihilation

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

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

LT

「うまい飯が作りたい」

introduction_of_recipes from Yutaka Kuroki

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

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

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

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

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

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

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

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

資料なし

CiNii API その2

資料なし

Rで健康体

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

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

Soccer × Attribution Analysis

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

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

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

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

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

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

資料なし

結婚式について

なぜ開くのか

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

工程

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

工夫した点

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

結果としてよかったこと

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

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

反省点

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

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

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

rstanarmパッケージを使って簡単にベイズモデリングを実行する

はじめに

今回は、rstanarmというパッケージを用いて赤ワインデータを色々といじってみようと思います。
マーケティングの意思決定のための分析などでベイズ統計を使う場面が多々あるのですが、似たような属性のデータがあるのであれば、
一つ一つstanコードを書くのではなく、Rの関数でサクッと実行して試行錯誤していくという形に持っていけたらいいなぁと感じていました。
本気を出すところではstanを、ルーティンワーク的なタスクではrstanarmをみたいな形で使い分けれると良いのではないでしょうか。

rstanarmとは

バックエンドの計算をStanに実行させて、統計モデルの推定を行うためのパッケージ。R上でlm関数のように簡単にベイズ推定を行うことができる。対象ユーザーはベイズ推定に慣れ親しんでいない頻度主義系のソフトウェアユーザー。
詳しくはこちら。

インストールする

まずはrstanarmのインストールするのですが、コケまくりました。そのため、バージョンを下げてみることにします。

devtools::install_version("rstanarm", version = "2.17.3", repos = "http://cran.us.r-project.org")

1	devtools::install_version("rstanarm", version = "2.17.3", repos = "http://cran.us.r-project.org")

ここに過去のバージョンがありますが、2.17.4だと動かなかったものの、2.17.3なら動きました。

rstanarmのサンプルを回してみる

今回は、以下の文献を参考にして、大人のIrisとも言える、ワインデータを扱い、質の高いワインかどうかを決める要素を探ります。
How to Use the rstanarm Package | Jonah Gabry and Ben Goodrich

こちらの文献には、ベイズ分析の4つのステップとして以下があげられています。

1.同時分布の特定（同時分布は事前分布と条件付きの尤度をかけ合わせたもの。）
2.MCMCで事後分布を描く
3.モデルがフィットしているか評価する
4.事後予測分布を描き、結果に影響を与える予測項を確認する。

これらのステップをできるだけ素早くできると良いですね。

まずはデータを読み込んで、スケーリングしておきます。（可視化結果は前回と同じなので、載せません。）
加えて、6点以上の評価であれば1を取る二項変数を作成しておきます。

library(tidyverse)
library(GGally)
library(rstanarm)
library(shinystan)
library(loo)

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

#可視化
wine_dataset_vis <- wine_dataset
ggpairs(wine_dataset_vis)

wine_dataset <- wine_dataset %>% mutate( y = if_else(condition = quality > 6, 1, 0 ))
wine_dataset <- wine_dataset %>% mutate_at(funs(scale(.)), .vars = c(1:11))
wine_dataset <- wine_dataset %>% select(-quality)

library(tidyverse)

library(GGally)

library(rstanarm)

library(shinystan)

library(loo)

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

#可視化

wine_dataset_vis <- wine_dataset

ggpairs(wine_dataset_vis)

wine_dataset <- wine_dataset %>% mutate( y = if_else(condition = quality > 6, 1, 0 ))

wine_dataset <- wine_dataset %>% mutate_at(funs(scale(.)), .vars = c(1:11))

wine_dataset <- wine_dataset %>% select(-quality)

GLMでロジスティックモデルを推定し、rstanarmで推定した結果と比較します。rstanarmでは傾きや切片の事前分布にスチューデントのt分布を、尤度にロジスティック分布を設定しています。

# GLM Estimation ----------------------------------------------------------
wine_glm_1 <- glm(y ~ fixed.acidity + volatile.acidity + citric.acid + residual.sugar + chlorides + free.sulfur.dioxide + 
                    total.sulfur.dioxide + density + pH + sulphates + alcohol,
                        data = wine_dataset, family = binomial(link = "logit"))


# Bayesian Estimation with rstanarm ---------------------------------------
wine_bglm_1 <- stan_glm(y ~ fixed.acidity + volatile.acidity + citric.acid + residual.sugar + chlorides + free.sulfur.dioxide + 
                          total.sulfur.dioxide + density + pH + sulphates + alcohol,
                              data = wine_dataset,
                              family = binomial(link = "logit"), 
                              prior = student_t(df = 7), 
                              prior_intercept = student_t(df = 7),
                              chains = 4, cores = 4, seed = 123)

#glmの結果との確認
cbind("model_glm" = coef(wine_glm_1),
      "model_bayesianglm" = coef(wine_bglm_1 ))

                       model_glm model_bayesianglm
(Intercept)          -2.81452789       -2.85848089
fixed.acidity         0.47871946        0.48115910
volatile.acidity     -0.46215347       -0.46939630
citric.acid           0.11060698        0.10547044
residual.sugar        0.33762729        0.32948671
chlorides            -0.41494490       -0.43932815
free.sulfur.dioxide   0.11318519        0.11910296
total.sulfur.dioxide -0.54377987       -0.56297397
density              -0.48655012       -0.47686448
pH                    0.03461116        0.03067029
sulphates             0.63563064        0.64074793
alcohol               0.80280901        0.82144931

# GLM Estimation ----------------------------------------------------------

wine_glm_1 <- glm(y ~ fixed.acidity + volatile.acidity + citric.acid + residual.sugar + chlorides + free.sulfur.dioxide +

total.sulfur.dioxide + density + pH + sulphates + alcohol,

data = wine_dataset, family = binomial(link = "logit"))

# Bayesian Estimation with rstanarm ---------------------------------------

wine_bglm_1 <- stan_glm(y ~ fixed.acidity + volatile.acidity + citric.acid + residual.sugar + chlorides + free.sulfur.dioxide +

total.sulfur.dioxide + density + pH + sulphates + alcohol,

data = wine_dataset,

family = binomial(link = "logit"),

prior = student_t(df = 7),

prior_intercept = student_t(df = 7),

chains = 4, cores = 4, seed = 123)

#glmの結果との確認

cbind("model_glm" = coef(wine_glm_1),

"model_bayesianglm" = coef(wine_bglm_1 ))

model_glm model_bayesianglm

(Intercept) -2.81452789 -2.85848089

fixed.acidity 0.47871946 0.48115910

volatile.acidity -0.46215347 -0.46939630

citric.acid 0.11060698 0.10547044

residual.sugar 0.33762729 0.32948671

chlorides -0.41494490 -0.43932815

free.sulfur.dioxide 0.11318519 0.11910296

total.sulfur.dioxide -0.54377987 -0.56297397

density -0.48655012 -0.47686448

pH 0.03461116 0.03067029

sulphates 0.63563064 0.64074793

alcohol 0.80280901 0.82144931

ほとんど係数の大きさが同じであることが確認できます。

ベイズ推定の良いところは事後分布から関心のある係数に関しての取りうる値などをシミュレーションできるところですが、
posterior_interval関数で簡単に計算することができます。

> round(posterior_interval(wine_bglm_1, prob = 0.95, pars = "fixed.acidity"), 2)
              2.5% 97.5%
fixed.acidity 0.03   0.9
> round(posterior_interval(wine_bglm_1, prob = 0.95, pars = "volatile.acidity"), 2)
                  2.5% 97.5%
volatile.acidity -0.76  -0.2

> round(posterior_interval(wine_bglm_1, prob = 0.95, pars = "fixed.acidity"), 2)

2.5% 97.5%

fixed.acidity 0.03 0.9

> round(posterior_interval(wine_bglm_1, prob = 0.95, pars = "volatile.acidity"), 2)

2.5% 97.5%

volatile.acidity -0.76 -0.2

肝心のMCMCの収束診断ですが、shinystanを使います。

やや余談ですが、他のデータセットでshinystanを用いた際に、予測結果にNAsが含まれている場合に、
shinystanが起動しないという問題があり、以下のようなエラー文が吐かれます。

Error in validate_y(y) : NAs not allowed in 'y'.

1	Error in validate_y(y) : NAs not allowed in 'y'.

調べたところ、こちらのgithubにあるように、

launch_shinystan(womensrole_bglm_1,ppd=FALSE)

1	launch_shinystan(womensrole_bglm_1,ppd=FALSE)

のように引数でppd=FALSEのように設定することで、立ち上げることができました。

3つの基準をクリアしているため、収束しています。

係数の分布についても可視化します。

rstanarmの良い点の一つとして、モデルのアップデートが容易に行える点があげられると思いますが、実際、以下のように先程のモデルに変数を追加して推定することができます。
今回は、alcoholを二乗したものを新しい変数として加えます。

#モデルのアップデート
(wine_bglm_2 <- update(wine_bglm_1, formula. = . ~ . + I(alcohol^2)))

1 2	#モデルのアップデート (wine_bglm_2 <- update(wine_bglm_1, formula. = . ~ . + I(alcohol^2)))

次に、looパッケージを用いて、更新したモデルと元のモデルの性能の比較を行います。
looパッケージは統計モデルの予測精度の指標として扱われる、WAIC(Widely Applicable Information Criterion)を計算するためのパッケージで、WICは事後分布から得られる対数尤度の平均や分散からなる値として表されます。looはleave-one-out cross-validationのleave-one-outの頭文字。

さっそく入れようと思ったところ、

Error: is.data.frame(data) || is.matrix(data) is not TRUE

1	Error: is.data.frame(data) \|\| is.matrix(data) is not TRUE

というエラーが出ました。
こちらでも議論されていましたが、

remove.packages("loo")
devtools::install_github("stan-dev/loo", ref = "v1.1.0")

1 2	remove.packages("loo") devtools::install_github("stan-dev/loo", ref = "v1.1.0")

でバージョンを2.0.0から1.1.0に落としたら動きました。

ここで、事後分布が特定のサンプルデータに対して敏感であるかどうかをlooパッケージを用いて可視化します。

loo_bglm_1 <- loo(wine_bglm_1)
loo_bglm_2 <- loo(wine_bglm_2)

par(mfrow = 1:2, mar = c(5,3.8,1,0) + 0.1, las = 3)
plot(loo_bglm_1, label_points = TRUE)
plot(loo_bglm_2, label_points = TRUE)

loo_bglm_1 <- loo(wine_bglm_1)

loo_bglm_2 <- loo(wine_bglm_2)

par(mfrow = 1:2, mar = c(5,3.8,1,0) + 0.1, las = 3)

plot(loo_bglm_1, label_points = TRUE)

plot(loo_bglm_2, label_points = TRUE)

縦軸のshape parameter kは推定の信頼性の指標とされ、大きければ大きいほど信頼できないと見なし、横軸は今回推定したワインデータのデータの番号で、左が元のモデル、右が変数を追加したモデルのものです。
どうやらどちらも0.4未満のkに収まっているようです。参考情報の事例では0.5を超えていましたが、moderate outliersと説明されていたので、今回の推定は問題ないと思われます。

続いてモデルの比較を行います。

> #モデルの比較
> compare_models(loo_bglm_1, loo_bglm_2)
elpd_diff        se 
      5.4       2.8

> #モデルの比較

> compare_models(loo_bglm_1, loo_bglm_2)

elpd_diff se

5.4 2.8

elpd_diffに関しては右のモデルの精度が高ければ正の値を、低ければ負の値を取るようになっています。標準誤差も返されます。
どうやら変数を追加したモデルの方が、ちょっとだけ良さそうです。

続いて、事後予測分布から、どの変数がどのように予測に影響を与えるのかを確かめます。
比較のためにデータを２つほど作成し、両者において一つだけ変数が違うという状況下での、予測される確率の比較を行います。

#パラメータを比較して変数が与える影響を確認する。
newdata <- data.frame(fixed.acidity = c(8.319637,8.319637) ,
                      volatile.acidity = c(0.5278205,0.5278205) ,
                      citric.acid = c(0.2709756,0.2709756),
                      residual.sugar = c(2.538806,2.538806),
                      chlorides = c(0.08746654,0.08746654) ,
                      free.sulfur.dioxide = c(15.87492,15.87492),
                      total.sulfur.dioxide = c(46.46779,46.46779),
                      density = c(0.9967467,0.9967467),
                      pH = c(3.311113,3.311113),
                      sulphates = c(0.6581488,0.6581488),
                      alcohol = c(15.42298,14.42298),
                      'I(alcohol^2)' = c(15.42298^2,14.42298^2) )

y_rep <- posterior_predict(wine_bglm_2, newdata)
summary(y_rep)
summary(apply(y_rep, 1, diff))


> summary(y_rep)
       1                2         
 Min.   :0.0000   Min.   :0.0000  
 1st Qu.:0.0000   1st Qu.:0.0000  
 Median :0.0000   Median :0.0000  
 Mean   :0.2402   Mean   :0.2717  
 3rd Qu.:0.0000   3rd Qu.:1.0000  
 Max.   :1.0000   Max.   :1.0000  
> 
> summary(apply(y_rep, 1, diff))
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
-1.0000  0.0000  0.0000  0.0315  0.0000  1.0000

#パラメータを比較して変数が与える影響を確認する。

newdata <- data.frame(fixed.acidity = c(8.319637,8.319637) ,

volatile.acidity = c(0.5278205,0.5278205) ,

citric.acid = c(0.2709756,0.2709756),

residual.sugar = c(2.538806,2.538806),

chlorides = c(0.08746654,0.08746654) ,

free.sulfur.dioxide = c(15.87492,15.87492),

total.sulfur.dioxide = c(46.46779,46.46779),

density = c(0.9967467,0.9967467),

pH = c(3.311113,3.311113),

sulphates = c(0.6581488,0.6581488),

alcohol = c(15.42298,14.42298),

'I(alcohol^2)' = c(15.42298^2,14.42298^2) )

y_rep <- posterior_predict(wine_bglm_2, newdata)

summary(y_rep)

summary(apply(y_rep, 1, diff))

> summary(y_rep)

1 2

Min. :0.0000 Min. :0.0000

1st Qu.:0.0000 1st Qu.:0.0000

Median :0.0000 Median :0.0000

Mean :0.2402 Mean :0.2717

3rd Qu.:0.0000 3rd Qu.:1.0000

Max. :1.0000 Max. :1.0000

> summary(apply(y_rep, 1, diff))

Min. 1st Qu. Median Mean 3rd Qu. Max.

-1.0000 0.0000 0.0000 0.0315 0.0000 1.0000

他の要素をある一定水準で保った際に、alcoholだけ1度下げることで、平均3%ほど高い評価が得られる確率が高まるという考察となります。

以上で、rstanarmの一連の使い方となるのですが、
一部の関数においては、階層ベイズモデルも行えるので、試してみようと思います。

ただ、階層ベイズにするにも、赤ワインのデータしかないので、グループ変数をどうにか作らないといけません。
あまりやりたくはありませんが、データがないので、説明変数を元にk-means(K=3)によるクラスタリングを行い、それをグループ変数とします。

stan_glmer関数を使えば、以下のような簡単な記述で定数項や係数がグループごとに異なるパラメータの分布に従うとする階層ベイズモデルを推定できます。

stan_glmer(非説明変数 ~ 変数 + ( 変数 - 1 | グループ変数 ),
                         data = wine_dataset_hc,
                         family = binomial(link = "logit"),
                         prior = student_t(df = 7), 
                         prior_intercept = student_t(df = 7),
                         chains = 4,cores = 4, seed = 123,iter = 500)

stan_glmer(非説明変数 ~ 変数 + ( 変数 - 1 | グループ変数 ),

data = wine_dataset_hc,

family = binomial(link = "logit"),

prior = student_t(df = 7),

prior_intercept = student_t(df = 7),

chains = 4,cores = 4, seed = 123,iter = 500)

今回は、切片だけがグループごとに異なるモデル、傾きだけがグループごとに異なるモデル、切片も傾きも異なるモデルを作成します。

# Hierachical_Models ------------------------------------------------------
wine_dataset_hc <- wine_dataset
wine_dataset_hc$wine_cluster <- as.numeric(kmeans(wine_dataset_hc %>% select(1:11) ,3)$cluster)
wine_dataset_hc %>% group_by(wine_cluster) %>% summarize(mean(as.numeric(y)),n())

# Hierarchical Models(Varying Intercept) -------------------------------------------------------
wine_bglm_3 <- stan_glmer(y ~ fixed.acidity + volatile.acidity + citric.acid + residual.sugar + chlorides + free.sulfur.dioxide + 
                           total.sulfur.dioxide + density + pH + sulphates + alcohol + ( 1 | wine_cluster ),
                         data = wine_dataset_hc,
                   family = binomial(link = "logit"),
                   prior = student_t(df = 7), 
                   prior_intercept = student_t(df = 7),
                   chains = 4,cores = 4, seed = 123,iter = 500)
coef(wine_bglm_3)

# Hierarchical Models(Varying Slope) ---------------------------------------------------------
wine_bglm_4 <- stan_glmer(y ~ fixed.acidity + volatile.acidity + citric.acid + residual.sugar + chlorides + free.sulfur.dioxide + 
                           total.sulfur.dioxide + density + pH + sulphates + alcohol +
                           ( fixed.acidity + volatile.acidity + citric.acid + residual.sugar + chlorides + free.sulfur.dioxide + 
                               total.sulfur.dioxide + density + pH + sulphates + alcohol - 1 | wine_cluster ),
                         data = wine_dataset_hc,
                         family = binomial(link = "logit"),
                         prior = student_t(df = 7), 
                         prior_intercept = student_t(df = 7),
                         chains = 4,cores = 4, seed = 123,iter = 500)
coef(wine_bglm_4)

# Hierarchical Models(Varying Intercept and Slope) -------------------------------------------
wine_bglm_5 <- stan_glmer(y ~ fixed.acidity + volatile.acidity + citric.acid + residual.sugar + chlorides + free.sulfur.dioxide + 
                           total.sulfur.dioxide + density + pH + sulphates + alcohol +
                           ( fixed.acidity + volatile.acidity + citric.acid + residual.sugar + chlorides + free.sulfur.dioxide + 
                               total.sulfur.dioxide + density + pH + sulphates + alcohol | wine_cluster ),
                         data = wine_dataset_hc,
                         family = binomial(link = "logit"),
                         prior = student_t(df = 7), 
                         prior_intercept = student_t(df = 7),
                         chains = 4,cores = 4, seed = 123,iter = 500)
coef(wine_bglm_5)

loo_bglm_1 <- loo(wine_bglm_1)
loo_bglm_3 <- loo(wine_bglm_3)
loo_bglm_4 <- loo(wine_bglm_4)
loo_bglm_5 <- loo(wine_bglm_5)

# Hierachical_Models ------------------------------------------------------

wine_dataset_hc <- wine_dataset

wine_dataset_hc$wine_cluster <- as.numeric(kmeans(wine_dataset_hc %>% select(1:11) ,3)$cluster)

wine_dataset_hc %>% group_by(wine_cluster) %>% summarize(mean(as.numeric(y)),n())

# Hierarchical Models(Varying Intercept) -------------------------------------------------------

wine_bglm_3 <- stan_glmer(y ~ fixed.acidity + volatile.acidity + citric.acid + residual.sugar + chlorides + free.sulfur.dioxide +

total.sulfur.dioxide + density + pH + sulphates + alcohol + ( 1 | wine_cluster ),

data = wine_dataset_hc,

family = binomial(link = "logit"),

prior = student_t(df = 7),

prior_intercept = student_t(df = 7),

chains = 4,cores = 4, seed = 123,iter = 500)

coef(wine_bglm_3)

# Hierarchical Models(Varying Slope) ---------------------------------------------------------

wine_bglm_4 <- stan_glmer(y ~ fixed.acidity + volatile.acidity + citric.acid + residual.sugar + chlorides + free.sulfur.dioxide +

total.sulfur.dioxide + density + pH + sulphates + alcohol +

( fixed.acidity + volatile.acidity + citric.acid + residual.sugar + chlorides + free.sulfur.dioxide +

total.sulfur.dioxide + density + pH + sulphates + alcohol - 1 | wine_cluster ),

data = wine_dataset_hc,

family = binomial(link = "logit"),

prior = student_t(df = 7),

prior_intercept = student_t(df = 7),

chains = 4,cores = 4, seed = 123,iter = 500)

coef(wine_bglm_4)

# Hierarchical Models(Varying Intercept and Slope) -------------------------------------------

wine_bglm_5 <- stan_glmer(y ~ fixed.acidity + volatile.acidity + citric.acid + residual.sugar + chlorides + free.sulfur.dioxide +

total.sulfur.dioxide + density + pH + sulphates + alcohol +

( fixed.acidity + volatile.acidity + citric.acid + residual.sugar + chlorides + free.sulfur.dioxide +

total.sulfur.dioxide + density + pH + sulphates + alcohol | wine_cluster ),

data = wine_dataset_hc,

family = binomial(link = "logit"),

prior = student_t(df = 7),

prior_intercept = student_t(df = 7),

chains = 4,cores = 4, seed = 123,iter = 500)

coef(wine_bglm_5)

loo_bglm_1 <- loo(wine_bglm_1)

loo_bglm_3 <- loo(wine_bglm_3)

loo_bglm_4 <- loo(wine_bglm_4)

loo_bglm_5 <- loo(wine_bglm_5)

先程紹介した、looパッケージを使って、ベースとなるモデルとの比較を行います。

> #モデルの比較
> compare_models(loo_bglm_1, loo_bglm_3)
elpd_diff        se 
     -0.8       0.5 
> compare_models(loo_bglm_1, loo_bglm_4)
elpd_diff        se 
      1.2       2.5 
> compare_models(loo_bglm_1, loo_bglm_5)
elpd_diff        se 
      1.0       2.5

> #モデルの比較

> compare_models(loo_bglm_1, loo_bglm_3)

elpd_diff se

-0.8 0.5

> compare_models(loo_bglm_1, loo_bglm_4)

elpd_diff se

1.2 2.5

> compare_models(loo_bglm_1, loo_bglm_5)

elpd_diff se

1.0 2.5

うーん、残念ながらどのモデルもベースモデルよりも圧倒的に強いものは無さそうです。

感想

まだまだrstanarmの関数やら機能やら定義を全て把握しきれていないですが、そこらへんがクリアーになれば、これまでのstanでの推定業務において生産性が高まる可能性を感じました。
簡単な階層ベイズモデルくらいなら、非常に直感的に書ける点や、変数の追加によるモデルのアップデートが容易な点などはポイント高いなぁと思います。
とはいえ、実務としてマニュアルでstanコードを作成していくのは必須なので、このパッケージを使うことによって、stanコードの改善に時間をより一層割けるようになるなら、それが一番だと思いました。
あと、「ベイズ初めてです！」という新入りの方とかには慣れ親しんでもらうには良さそうですね。lm関数レベルで実行できてしまうので。
今回、mc-stan.orgの配下にあるページなどを漁る過程で、ベイズ推定結果の可視化などで知らないことにも色々と出会えたので、今後も読み進めていきます。

追記

2018-09-10: Stanを使って変数選択したいにprojpredというパッケージが紹介されており、これを使えば、情報量基準に従った変数選択を簡単に行えるそうです。こうなると、「ベイズ推定に慣れ親しんでいない頻度主義系のソフトウェアユーザー」に限らず多くの人が幸せになれるパッケージなのかもしれませんね。

参考情報

Using the loo package (version >= 2.0.0)
Leave-one-outクロスバリデーションの２つのデメリット、からの解決方法
 stan_glmer | Bayesian Generalized Linear Models With Group-Specific Terms Via Stan
WAICを計算してみる
 Package ‘bayesm’
Hierarchical Partial Pooling for Repeated Binary Trials
Leave-one-out cross-validation for non-factorizable models
priors | Prior Distributions And Options
StanとRでベイズ統計モデリング (Wonderful R)

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でグラフィカルモデル