ベータ二項分布を用いたアイテム推薦におけるCTR推定（２）
最近、以下の書籍が出版されました。
レコメンドシステムの分析手法が書かれている和書はあまり出版されておらず*1、以前読んだことがあるとすれば、「集合知プログラミング*2」くらいでした。
教師あり学習の２値問題などは最近書籍で取り上げられることは多いのですが、この手の問題は教師あり学習手法だけでは不十分であり、ベイズモデリングやその他の知識や技術が求められることから、より専門性が高くなります。またマーケティング領域等における現実的な分析ニーズでは、２値問題を定義し、アプローチするケースが多い気がしますので*3、中々目に触れる機会がありません（私の努力不足でもありますが・・・）。
- 作者: Deepak K. Agarwal,Bee‐Chung Chen,島田直希,大浦健志
- 出版社/メーカー: 共立出版
- 発売日: 2018/04/21
- メディア: 単行本
- この商品を含むブログ (1件) を見る
さて、前置きはここまでにして、上記書籍の中でStanを用いて解ける設問(3.6演習)がありますので、早速やっていきます。*4
アイテム１とアイテム２のCTRをそれぞれ とした場合、アイテム１のCTRがアイテム２より優位となる事後確率を求める問題です。
上記のCTRはアイテム２の方が大きいですが、各アイテムの表示回数が小さいとき、アイテム１の事後確率が結果的に大きくなる可能性がある、ということが考えられます*5。
題意より、ベータ二項分布を使用して、確率 の差分
と表示回数
の倍率
を以下のように動かし、上記ケースを見てみたいと思います。
ここで、ベータ二項分布のモデリングは以下記事に書いておりますので、割愛します。
http://blog.hatena.ne.jp/FAKED-SCIENTIST/faked-scientist.hatenablog.com/edit?entry=17391345971642131739
また今回は以下の二項分布に従うサンプルデータ（表示回数の） を生成し、数値シミュレートしました。
アイテム１のCTRは30%と設定しましたが、適宜変更して挙動を確かめてもいいかと思いますが、似たような振る舞いをすると思います。
使用した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 ( ) を推定するものです。後半は推定結果をプロットしてます。
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]) }
以下は、結果のプロット①に該当するグラフです。
横軸 に対して、
ごとのEAPを折れ線グラフでプロットしました。
ご覧頂く通り、 が小さく、かつ
まではアイテム１とアイテム２のEAPは大小関係が逆転するケースがあります。これは表示回数
が小さいとき、アイテム１・２のCTRは安定せず、それが逆転する事象を生起させている、と理解できます。
結果のプロット②では、2.5%～97.5%の信用区間 (Credit Interval) を ごとに点線で表しました。実線はプロット①で示したEAPです。
表示回数の倍率 が小さい場合は信用区間の幅が広くなり、
が大きくなるに従い、徐々に狭まることが分かります。つまり、
の大きさがCTRの安定化に寄与します。
