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

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

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