Fake Data-Scientistの学習日記

似非データサイエンティストの学習記録です。

ベータ二項分布を用いたアイテム推薦におけるCTR推定（２）

R Stan 推薦システム ベイズ統計学

最近、以下の書籍が出版されました。
レコメンドシステムの分析手法が書かれている和書はあまり出版されておらず*1、以前読んだことがあるとすれば、「集合知プログラミング*2」くらいでした。

教師あり学習の２値問題などは最近書籍で取り上げられることは多いのですが、この手の問題は教師あり学習手法だけでは不十分であり、ベイズモデリングやその他の知識や技術が求められることから、より専門性が高くなります。またマーケティング領域等における現実的な分析ニーズでは、２値問題を定義し、アプローチするケースが多い気がしますので*3、中々目に触れる機会がありません（私の努力不足でもありますが・・・）。


推薦システム: 統計的機械学習の理論と実践

推薦システム: 統計的機械学習の理論と実践


さて、前置きはここまでにして、上記書籍の中でStanを用いて解ける設問(3.6演習)がありますので、早速やっていきます。*4

アイテム１とアイテム２のCTRをそれぞれ p, p+\delta(>0)とした場合、アイテム１のCTRがアイテム２より優位となる事後確率を求める問題です。
上記のCTRはアイテム２の方が大きいですが、各アイテムの表示回数が小さいとき、アイテム１の事後確率が結果的に大きくなる可能性がある、ということが考えられます*5。

題意より、ベータ二項分布を使用して、確率 p の差分 \delta と表示回数 n の倍率 k を以下のように動かし、上記ケースを見てみたいと思います。

\begin{equation} \delta = \begin{cases} 0 & (アイテム１の場合) \\ 0.005, 0.0075, 0.01, 0.025,0.05 & (アイテム２の場合) \\ \end{cases} \end{equation}
\begin{equation} k = 1, 2, 3, \cdots, 100 \end{equation}

ここで、ベータ二項分布のモデリングは以下記事に書いておりますので、割愛します。
http://blog.hatena.ne.jp/FAKED-SCIENTIST/faked-scientist.hatenablog.com/edit?entry=17391345971642131739


また今回は以下の二項分布に従うサンプルデータ（表示回数の） X を生成し、数値シミュレートしました。
アイテム１のCTRは30%と設定しましたが、適宜変更して挙動を確かめてもいいかと思いますが、似たような振る舞いをすると思います。

\begin{eqnarray} X &\sim& \text{Binomial} \hspace{1.0mm}(kn, p+\delta) \\ p &=& 0.3 \\ n &=& 5, 10, 15, \cdots, 100 \end{eqnarray}

使用したStanコードです。

#stanコード：CTR_pred.stan
data{
  int<lower=0> n;
  int<lower=0> X[n];
  int<lower=0> N[n];
}

parameters{
  real<lower=0, upper=1> p;
  real<lower=0> alpha;
  real<lower=0> beta;
}

model{
  #事前分布
  p ~ beta(alpha, beta);
  #尤度
  X ~ binomial(N, p);
  #弱情報事前分布
  alpha ~ uniform(0,1);
  beta ~ uniform(0,1);
}

使用したRコードです。
前半はサンプルデータを生成し、それに対してMCMCでCTR ( p, p + \delta ) を推定するものです。後半は推定結果をプロットしてます。

library(rstan)
# Stanモデルのコンパイル
stanmodel <- stan_model(file='~/CTR_pred.stan')
#Set p値（事後分布のパラメータの真値）
p <- 0.3

# 差分δの値
delta <- c(0, 0.005, 0.0075, 0.01, 0.025,0.05);len.delta <- length(delta)

# 出力パラメータの設定
pars <- c("p")

# kの最大値を設定
K <- 100
#結果を格納するデータフレーム
res <- data.frame(matrix(NA, ncol=9, nrow=len.delta*K));
colnames(res) <- c('k', 'delta.num', 'p+delta', 'mean', '2.5%', '25%', '50%', '75%', '97.5%')
# カウンター変数
cnt <- 0

#########################
#サンプリング START
#########################
for(delta.cnt in 1:len.delta){
  for(k in 1:K){
    #表示回数の設定
    N <- seq(5, 100, by=5) * k; t_max <- length(N)
    
    # p + δの設定
    p2 <- p + delta[delta.cnt]

    # 二項分布に従うサンプルデータの生成
    X_t <- rep(NA, t_max)
    for(i in 1:t_max)
      X_t[i] <- rbinom(n=1, size=N[i], prob=p2 )
    
    # 読込データの設定
    data <- list(n=t_max, X=X_t, N=N)
    
    # サンプリング
    fit_stan <- sampling(
      stanmodel,
      data = data,
      par=pars,
      seed = 1234,
      chains=1, iter=10000, warmup=2500
    )
    
    # 結果格納
    cnt <- cnt + 1
    res[cnt, 1] <- k; res[cnt, 2] <- delta.cnt
    res[cnt, 3] <- p2; res[cnt, 4] <- mean(as.data.frame(fit_stan)$p)
    res[cnt, 5:9] <- quantile(as.data.frame(fit_stan)$p, probs = c(0.025, 0.25, 0.5, 0.75, 0.975))
  }
}
#########################
#サンプリング END
#########################

# 結果のプロット①
for(d in 1:len.delta){
  res.sub <- subset(res, delta.num==d)
  plot(res.sub$k, res.sub$mean, type="l", 
       xlim=c(0,K), ylim=c(0.28, 0.38), xlab="k", ylab="EAP of p+δ", col=d, lty=1, lwd=2)
  if(d < len.delta) par(new=T)
}
labels <- c("0.10", "0.105", "0.1075", "0.11", "0.125","0.15")
legend("topleft", legend = labels, col = 1:7,lty = 1, lwd=2)

# 結果のプロット②
xlims <- c(0, K); ylims<- c(0.28, 0.34)
mains <- c("δ = 0", "δ = 0.005", "δ = 0.0075", "δ = 0.01");
par(mfrow=c(2, 2)) 
for(d in 1:4){
  res.sub <- subset(res, delta.num==d)
  plot(x=res.sub$k, res.sub[,'50%'], type="l", 
       xlim=xlims, ylim=ylims, xlab="", ylab="", col=d, lty=1, lwd=2)
  par(new=T)
  plot(res.sub$k, res.sub[,'2.5%'], type="l", 
       xlim=xlims, ylim=ylims, xlab="", ylab="", col=d, lty=2, lwd=1)
  par(new=T)
  plot(res.sub$k, res.sub[,'97.5%'], type="l", 
       xlim=xlims, ylim=ylims, xlab="k", ylab="", col=d, lty=2, lwd=1,main=mains[d])
}


以下は、結果のプロット①に該当するグラフです。
横軸 k に対して、 \delta ごとのEAPを折れ線グラフでプロットしました。
ご覧頂く通り、 k が小さく、かつ \delta = 0.01 まではアイテム１とアイテム２のEAPは大小関係が逆転するケースがあります。これは表示回数 kn が小さいとき、アイテム１・２のCTRは安定せず、それが逆転する事象を生起させている、と理解できます。
f:id:FAKED-SCIENTIST:20180509230219j:plain:w450

結果のプロット②では、2.5%～97.5%の信用区間 (Credit Interval) を \delta = 0, 0.005, 0.0075, 0.01 ごとに点線で表しました。実線はプロット①で示したEAPです。
表示回数の倍率 k が小さい場合は信用区間の幅が広くなり、k が大きくなるに従い、徐々に狭まることが分かります。つまり、k の大きさがCTRの安定化に寄与します。
f:id:FAKED-SCIENTIST:20180509232537j:plain

*1:洋書の状況は分かりません。

*2:www.oreilly.co.jp

*3:一部のデータサイエンティストや機械学習エンジニアには使う場面があるかと思いますが。

*4:あくまで個人の主観に基づくものですので、著者が意図しない解答になっている可能性があります。ご了承ください。

*5:本の同じ章にヒントは記述されてました

ベータ二項分布を用いたアイテム推薦におけるCTR推定

R Stan ベイズ統計

ベータ二項分布を用いた分析例に取り上げます。

ECサイトに訪問したユーザに対して、あるアイテムを推薦した場合、どの程度の反応があるかを測定したいケースが考えられます。
効果測定の指標をCTR(Click Through Rate)とすれば、アイテムの推薦(表示)回数に対して、ユーザがどのくらい応答(クリック)したかを測ります。

このとき、推薦回数を N 、CTRを p とし、応答数 X は以下の二項分布に従うと仮定します。

X \sim \text{Binomial}(N, p) \hspace{5.0mm} p \in [0, 1]


繰り返しにはなりますが、運用者が知りたいことは、アイテムの推薦回数 N に対して、どのくらい応答があるか、ということです。すなわち、CTRの確率 p を推測し、アイテムの推薦を効率的に遂行できる環境を構築していきたい、と言えます。


さて、ベイズの定理より、 p の事後分布のカーネルを求めると、

\begin{equation} f(p=\text{CTR} \mid X= \text{応答数}) \propto \hspace{2mm}f(X \mid p) \hspace{1.0mm} f(p) \end{equation}

となります。このとき、カーネルの尤度 f(X \mid p) が二項分布ですので、自然共役事前分布はベータ分布です。
この構造をグラフィカルモデリングで表すと、以下のようになります。

f:id:FAKED-SCIENTIST:20180508122918j:plain:w300
CTR推定のプレート図
\begin{eqnarray} p & \sim & \text{Beta} \hspace{1.0mm} (\alpha, \beta) \\ X_t & \sim & \text{Binomial} \hspace{1.0mm}(N_t, p) \\ \alpha, \beta & \sim & \text{Uniform} \hspace{1.0mm}(0, 1) \end{eqnarray}

ここで、 N_t, X_t は時刻 t の表示回数とそれに対する応答数を表します。
また \alpha, \beta もデータから推定するパラメータとなっており、その事前分布に弱情報事前分布を設定しています。このケースでは、凡庸な手段もしれませんが、一様分布としました。

Stanコードは以下です。

#stanコード：CTR_pred.stan
data{
  int<lower=0> n;
  int<lower=0> X[n];
  int<lower=0> N[n];
}

parameters{
  real<lower=0, upper=1> p;
  real<lower=0> alpha;
  real<lower=0> beta;
}

model{
  #事前分布
  p ~ beta(alpha, beta);
  #尤度
  X ~ binomial(N, p);
  
  #弱情報事前分布
  alpha ~ uniform(0,1);
  beta ~ uniform(0,1);
}


Rコードは以下となります。
今回はCTRを p = 15\% と設定し、表示回数を50～500回まで10ずつ動かし、各々の応答数を生成しました。

library(rstan)
#乱数の設定
set.seed(1234)

#パラメータ設定
N <- seq(50, 500, by=10)
t_max <- length(N)
p <- 0.15

#応答数 を生成
X_t <- rep(NA, t_max)
for(i in 1:t_max)
  X_t[i] <- rbinom(n=1, size=N[i], prob=p)

#読込データを設定
data <- list(n=t_max, X=X_t, N=N)

#モデルのコンパイル
stanmodel <- stan_model(file='~/CTR_pred.stan')

#出力パラメータの設定
pars <- c("p","alpha", "beta")

#サンプリング
fit_stan <- sampling(
  stanmodel,
  data = data,
  par=pars,
  seed = 1234,
  chains=4, iter=10000, warmup=2500
)

#結果の表示
dig <- 3
print(fit_stan, pars=pars, digits=dig)

#トレースプロットの出力
traceplot(fit_stan, pars=pars, alpha=0.5)


EAP推定値はこんな感じでした。

mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
p 0.149 0.000 0.003 0.142 0.146 0.148 0.151 0.155 2282 1.001
\alpha 0.517 0.006 0.260 0.072 0.303 0.509 0.730 0.971 1637 1.006
\beta 0.622 0.006 0.249 0.120 0.433 0.651 0.832 0.984 1730 1.003


トレースプロットで推定値が収束していることを確認できました。
f:id:FAKED-SCIENTIST:20180508151014j:plain:w450

ベイズモデリングの実践編_#1-2

R Stan 実践ベイズモデリング

「基礎からのベイズ統計学」と並行して、「ベイズモデリングの実践編」も学習中です。
徐々に、Stanの使い方やReferenceの見方をキャッチアップできてきました。

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

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


早速ですが、実装した教科書の例題(分析例: 1-2)を載せていきます。
ガンベル分布を利用して、極限統計を解く問題です。
使用したStanコードとRコードは以下となります。

Stanコード

#ana1-2.stan
data{
  int<lower=0> N;
  real<lower=0> Y[N];
}

parameters{
 #RQ.1 最頻値と確信区間
  real mu;
  real<lower=0> sigma;
}

model{
  for(i in 1:N){
    Y[i] ~ gumbel(mu, sigma);
  }
}

generated quantities{
  real<lower=0> s;
  real p;
  real x;
  real<lower=0, upper=1> U;
  real inv_r;
  real y;
  real<lower=0, upper=1> Uy;
  
 #RQ.2 各都市の最長記録の散らばり(標準偏差)
  s = pi() * sigma * inv_sqrt(6);

 #RQ.3 実現値の確率
  p = exp(-exp((mu-8.95)/sigma));

 #Rq.4 ガンベル分布の99%点の再現レベル
  x = mu - sigma * log(-log(0.99));

 #Rq.4 世界記録の再現期間が100年より長い確率
  U = x < 8.95 ? 1 : 0;

 #RQ.5 再現期間の必要年数
  inv_r = inv(1 - p);

 #RQ.6 推定パタメータを用いたガンベル分布の予測値
  y = gumbel_rng(mu, sigma);
 #翌年に世界記録が更新される確率
  Uy = y > 8.95 ? 1 : 0;
}

Rコード

library(rstan)

#各年の最長記録データ
max_record <- c(8.95, 8.68, 8.70, 8.74, 8.71, 8.58, 8.63, 8.48, 8.60, 8.65,
                8.41, 8.52, 8.53, 8.60, 8.60, 8.56, 8.66, 8.73, 8.74, 8.47,
                8.54, 8.35, 8.56, 8.51, 8.52)

#読込データの設定
data <- list(N=length(max_record), Y=max_record)

#モデルのコンパイル
stanmodel <- stan_model(file='./ana1-2.stan')

#表示パラメータの設定
pars <- c("mu","sigma", "s", "p", "x", "U", "inv_r", "y", "Uy")

#サンプリング
fit_stan <- sampling(
  stanmodel,
  data = data,
  par=pars,
  seed = 1234,
  chains=4, iter=20000, warmup=5000
)

#結果確認
dig <- 3
print(fit_stan,pars=pars, digits_summary=dig)

#トレースプロット
traceplot(fit_stan, pars=pars, alpha=0.5)

結果はこんな感じでした。

Inference for Stan model: ana1-2.
4 chains, each with iter=20000; warmup=5000; thin=1; 
post-warmup draws per chain=15000, total post-warmup draws=60000.

        mean se_mean     sd   2.5%    25%    50%    75%   97.5% n_eff Rhat
mu     8.542   0.000  0.026  8.492  8.525  8.542  8.559   8.594 30623    1
sigma  0.121   0.000  0.019  0.090  0.107  0.119  0.133   0.166 28628    1
s      0.156   0.000  0.025  0.115  0.138  0.153  0.170   0.213 28628    1
p      0.963   0.000  0.023  0.905  0.952  0.968  0.979   0.991 24923    1
x      9.100   0.001  0.100  8.935  9.030  9.090  9.159   9.326 26889    1
U      0.041   0.001  0.198  0.000  0.000  0.000  0.000   1.000 36026    1
inv_r 39.402   0.171 30.153 10.497 20.848 31.157 48.012 117.175 31026    1
y      8.613   0.001  0.161  8.378  8.502  8.586  8.694   9.007 55916    1
Uy     0.038   0.001  0.191  0.000  0.000  0.000  0.000   1.000 60000    1


最後に \muと \sigmaのEAP推定値を用いて、ガンベル分布の密度関数をプロットします。

library(evd)
x <- 8:10
dgumbel(x, loc=8.542, scale=0.121)
curve(dgumbel(x, loc=8.542, scale=0.121), 8, 10, type="l", ylab="f(x)")

f:id:FAKED-SCIENTIST:20180504210928p:plain

ベイズ統計の学習記録（演習#6-3-3）

R Stan 基礎からのベイズ統計学

章末問題6.3-3)
題意より、昨年大会のボーダーラインを0.90分位数のEAP推定値を以下で求めます。

\begin{align} \displaystyle \xi_{0.90}^{(t)} = g(\mu^{(t)}, \sigma^{(t)}) = \mu^{(t)} + z_{0.90}\sigma^{(t)} = \mu^{(t)} + 1.282\sigma^{(t)} \end{align}

続いて、B子さんの実績(平均87点、標準偏差5点)より、B子さんが本選に進める確率を定式化すると、

\displaystyle \begin{align} p_{1-\Psi(\xi_{0.90})}^{(t)} = g(\xi_{0.90}^{(t)}) = 1 - \Psi(\xi_{0.90}^{(t)}) \end{align}
のようになります。

ここで、normal_cdf関数を用いれば、 1-\Psi(z)に標準化( z = \frac{\xi_{0.90}^{(t)}-87}{5} )された正規累積分布関数の確率を返してくれます。


実際に使用したStanコードとRコードは以下の通りです。

Stanコード

#ex_6-3-3.stan

#データサイズと確率変数Yを定義
data{
  int<lower=0> N;
  real<lower=0> Y[N];
}

#正規分布の平均と分散を設定
parameters{
  real mu;
  real<lower=0> sigma;
}

#パラメータの数値変換
#0.90分位数を平均と標準偏差より算出
transformed parameters {
  real xi;
  xi = mu + 1.282*sigma;
}

#確率変数Yを平均mu, 標準偏差sigmaの正規分布に従う確率変数Yを生成
model{
  for(i in 1:N){
    Y[i] ~ normal(mu, sigma);
  }
}

#生成量の計算
#0.90分位数を超える確率を計算
generated quantities{
  real<lower=0, upper=1> prob_over;
  prob_over = 1- normal_cdf(xi , 87, 5);
}

Rコード

library(rstan)

#データセットの作成
karaoke_score <- c(75, 82, 77, 86, 98, 91, 85, 84, 82, 79,
                88, 79, 84, 87, 69, 93, 84, 82, 89, 78,
                90, 74, 75, 84, 89, 81, 74, 79, 82, 84)

#Stanで読み込むデータに整形
data <- list(N=length(karaoke_score), Y=karaoke_score)

#Stanのmodelファイルをコンパイル
stanmodel <- stan_model(file='./ex_6-3-3.stan')

#表示するパラメータの設定
pars <- c("xi", "prob_over")

#Stanのmodelファイルをコンパイル
fit_stan <- sampling(
  stanmodel,
  data = data,
  par=pars,
  seed = 1234,
  chains=1, iter=11000, warmup=1000
)

#結果の表示
#有効桁数は3と設定
dig <- 3
print(fit_stan,pars=pars, digits_summary=dig)

#traceplotでMCMCサンプルが収束しているどうかを確認
traceplot(fit_stan, pars=pars, alpha=0.5)

結果はこんな感じでした。

Inference for Stan model: ex_6-3-3.
1 chains, each with iter=11000; warmup=1000; thin=1; 
post-warmup draws per chain=10000, total post-warmup draws=10000.

            mean se_mean    sd   2.5%    25%    50%    75%  97.5% n_eff Rhat
xi        91.297   0.019 1.714 88.348 90.100 91.155 92.337 95.024  7962    1
prob_over  0.208   0.001 0.089  0.054  0.143  0.203  0.268  0.394  8772    1

ベイズ統計の学習記録（演習#6-3-1）

Stan R 基礎からのベイズ統計学

ベイズ統計の学習にあたり、以下の書籍で解いたプログラミング演習を備忘録の意味も込め、掲載しています。*1


基礎からのベイズ統計学: ハミルトニアンモンテカルロ法による実践的入門

基礎からのベイズ統計学: ハミルトニアンモンテカルロ法による実践的入門


章末問題6.3-1)
正規分布を用いて、12週の売上平均を推測する問題となっています。
実際に使用したStanコードとRコードは以下の通りです。

Stanコード

#ex_6-3-1.stan

#データサイズと確率変数Yを定義
data{
  int N;
  real Y[N];
}

#正規分布の平均と分散を設定
parameters{
  real mu;
  real<lower=0> sigma;
}

#確率変数Yを平均mu, 標準偏差sigmaの正規分布に従う確率変数Yを生成
model{
  for(i in 1:N){
    Y[i] ~ normal(mu, sigma);
  }
}

#生成量の計算
#delta_over: 12週の平均売上が7万円を超える確率を計算
generated quantities{
  real delta_over;
  real delta_over2;
  delta_over = step(mu - 70000);
  delta_over2 = step(mu - 75000);
}

Rコード

library(rstan)

#データセットの作成
sales <- c(76230, 73550, 80750, 71500, 75420, 74840,
           71580, 76920, 68450, 76990, 64070, 76200)

#Stanで読み込むデータに整形
data <- list(N=length(sales), Y=sales)

#Stanのmodelファイルをコンパイル
stanmodel <- stan_model(file='./ex_6-3-1.stan')

#表示するパラメータの設定
pars <- c("mu", "delta_over", "delta_over2")

#Stanのmodelファイルをコンパイル
fit_stan <- sampling(
  stanmodel,
  data = data,
  par=pars,
  init=function(){
    list(mu=runif(1, 0,1000), sigma=runif(1, 0, 5))
  },
  seed = 1234,
  #11000回のサンプリング、ワームアップは1000回
  chains=1, iter=11000, warmup=1000
)

#結果の表示
#有効桁数は3と設定
dig <- 3
print(fit_stan,pars=pars, digits_summary=dig)

#traceplotでMCMCサンプルが収束しているどうかを確認
traceplot(fit_stan, pars=pars, alpha=0.5)

結果はこんな感じでした。

Inference for Stan model: ex_6-3-1.
1 chains, each with iter=11000; warmup=1000; thin=1; 
post-warmup draws per chain=10000, total post-warmup draws=10000.

                 mean se_mean       sd     2.5%      25%      50%      75%    97.5% n_eff Rhat
mu          73873.476  21.176 1517.223 70825.24 72927.85 73877.61 74795.38 76959.06  5133    1
delta_over      0.992   0.001    0.091     1.00     1.00     1.00     1.00     1.00  4970    1
delta_over2     0.207   0.005    0.405     0.00     0.00     0.00     0.00     1.00  5862    1

*1:あくまでの個人の見解ですので、間違いのある可能性もございます。その場合はご了承ください。