ジョンとヨーコのイマジン日記

多くのABテストではクリック数やコンバージョン数などのカウントデータを比較します。

ABテストで問題になるのが、

どのくらい差があれば、十分「Bのほうが良さそう」と判断できるのか
どのくらいのデータ（クリック数やコンバージョン数）がたまるまでまてばいいのか

などです。

そこで仮説検定の枠組みを用いて、例数（サンプルサイズ）設計をきっちり行ってABテストを実施しようと主張する人もいます。

が、たぶん無理だと思います。

その理由は次の通り：
伝統的な仮説検定の枠組みで例数設計を行うには以下を想定（設定）する必要がある

有意水準
検出力
A の期待CV（コンバージョン数）
B の期待CV

しかし AB テストの多くは探索的なものなので、A のCV、B のCVを想定するのは難しいし、有意水準とか検出力とかの設定も難しい。
また「検定して有意にならなかった場合、どうすればいいのか」という指針が与えられません。
（「確実な結果が欲しいので有意になるまで検定を続けます」とやるなら検定とはいえません）

ぼくは、ABテストで興味があるのは「有意差」ではないと思い、問題設定をこのように置きかえます。

『本当は A のほうがCVしやすいのに、誤って B の方を選んでしまう確率を $\alpha$ くらいに抑えたい場合、どのくらい試行回数が必要か』

今回はCVの比を使ってABテストの試行回数を見積もる方法を考えます。

用意するものは以下の2つ。

「AのCVに対してのBのCVの比がこれくらいあれば実質的に意味があるとみなそう」と決めた値（ $\beta$ と呼ぶことにする）
「本当は A のほうがCVしやすいのに、誤って B の方を選んでしまう確率を $\alpha$ くらいに抑えたい」という $\alpha$

AとBのCV数のモデルとしてそれぞれ次のポアソン分布を仮定します。

$a \sim \operatorname{Poisson}(\lambda)$
$b \sim \operatorname{Poisson}(\beta \lambda)$

B の CV数の期待値は A のほうの $\beta$ 倍です。

「本当は A のほうがCVしやすいのに、誤って B の方を選んでしまう確率」は $b/(a+b)$ が0.5以下になる確率と同値です。

ポアソン分布のモデルが正しいとして、 $r= b/(a+b)$ の分布はベータ分布で近似できます。

$s \sim \operatorname{Beta}(\beta \lambda, \lambda)$

なので、ベータ分布の分布関数より、 $s$ が0.5以下になる確率が $\alpha$ を下回るように $\lambda$ を定めることができます。

求めた $\lambda$ 以上のCV数がたまったら、コンバージョンが少しでも多い方を選べば、「本当は A のほうがCVしやすいのに、誤って B の方を選んでしまう確率」は $\alpha$ 以下に抑えられます。

$\beta$ や $\alpha$ を小さくとるほど必要な $\lambda$ は大きくなります。

求めた $\lambda$ が、予算や時間の制約から見て大きすぎる場合は、より大きな $\alpha$ を許容するか、 $\beta$ をより大きくして、「この程度の差なら微差だから気にしないことにする」という範囲を広げます。

実装例

R です。

target <- function(lambda,beta,alpha){
  pbeta(0.5,lambda*beta,lambda)-alpha  
}

sol1 <- uniroot(target,beta=1.03,alpha=0.05,interval = c(0,1e+05))
print(sol1$root)
#6101.871

3%のCV上昇（ $\beta=1.03$ ）を、誤りの確率5%（ $\alpha=0.05$ ）で検知したい場合、必要なCV数は約6102でした。

こういう風にすればシミュレーションができます。

simfunc <- function(lambda,beta){
  a <- rpois(100000,lambda)
  b <- rpois(100000,lambda*beta)
  data.frame(a=a,b=b)
}

sol1 <- uniroot(target,beta=1.03,alpha=0.05,interval = c(0,1e+05))
print(sol1$root)
#6101.871
set.seed(1);dat <- simfunc(round(sol1$root),1.03)
print(round(mean(dat$a>dat$b),2))
#0.05

sol2 <- uniroot(target,beta=1.02,alpha=0.05,interval = c(0,100000))
print(sol2$root)
#13662.36
set.seed(2);dat2 <- simfunc(round(sol2$root),1.03)
print(round(mean(dat$a>dat$b),2))
#0.05

なんパターンかためしてみれば、誤りの確率が設定した $\alpha$ とほぼ同じになることがわかると思います。

$\alpha$ を 0.05 に固定して、 $\beta$ を変化させて、必要な $\lambda$ の値をプロットしてみました。

f:id:abrahamcow:20210104062401p:plain

おしまい。

ここでは
[math/0604410] Discrete Component Analysis
の Gamma-Poisson モデル（GaP）の ELBO (evidence lower bound) を導出する。

まず、行列の分解がトピックモデルの一種として解釈できること説明する。

次に、モデルのパラメータ推定方法について述べる。

そしてELBOを導出し、最後にRによるシミュレーションを行う。

行列の復習からはじめてトピックモデルの結果だけ理解する

トピックモデルと総称されるモデルの中には様々なものがありますが、ここでは
[math/0604410] Discrete Component Analysis
の Gamma-Poisson モデル（GaP）を紹介します。

トピックモデルはカウントデータ（なにかの数を数えたデータ）の行列が与えられたとき、それをより小さい行列で近似して理解しやすくするためのモデルです。

行列というのはこんな感じの表のことです。

f:id:abrahamcow:20190211093407p:plain

多くのトピックモデルでは、単語の並び順は無視して、単語の出現頻度にのみ注目し、文書ごとに集計した表を利用します。
ここでは簡単のため3単語のみを集計したことにしていますが、もっと数が増えても大丈夫です。
その場合行列が横に長くなっていきます。

トピックモデルではこの行列がいくつかのトピックが混ざった結果生成されたと考えます。

トピックは単語の出現確率の違いで特徴づけられます。

例えば「刀」とか「印籠」という単語が出てくる確率は時代劇だったら高そうですが、現代劇にはめったに出てこなさそうですよね。

このような単語の出現確率の違いが文書のトピックを特徴づけていると考えるのがトピックモデルです。

トピックも行列で表しておくと便利です。

f:id:abrahamcow:20190211095029p:plain

今回の行列はトピックごとの単語の出現確率を表しているので、横合計が1になります。

トピック1は単語3が出る確率が高め、トピック2は単語2が出る確率が高めです。

ここでは簡単のため2種類のトピックのみを考えていますが、もっと数が増えても大丈夫です。

その場合行列が縦に長くなっていきます。

さて、各文書はトピックが混ざった結果出現すると書きました。

たとえば文書1はトピック1が40、トピック2が10出現したとしましょう。

表と見比べながら慎重に計算していくと、この文書の単語の出現頻度は次のようになります。

f:id:abrahamcow:20190211100755p:plain

トピック1に40をかけて、トピック2に10をかけています。単語はどちらのトピックにおいても出現確率を持っているので、トピック1から出現した数とトピック2から出現した数を足しています。

このような計算を行列では以下のように省略して書きます。

$\begin{pmatrix} 40 & 10\\ \end{pmatrix} \begin{pmatrix} 0.1 & 0.2 & 0.7\\ 0.3 & 0.5 & 0.2 \end{pmatrix}=\begin{pmatrix} 7 & 13 & 30\\ \end{pmatrix}$

意味は

$40 \times 0.1 + 10 \times 0.3 =7\\ 40 \times 0.2 + 10 \times 0.5 = 13\\ 40 \times 0.7 + 10 \times 0.2 = 30$

とまったく一緒です。

次に文書2はトピック1が20、トピック2が90出現したとしましょう。

行列を使って省略して書くと

$\begin{pmatrix} 20 & 90\\ \end{pmatrix} \begin{pmatrix} 0.1 & 0.2 & 0.7\\ 0.3 & 0.5 & 0.2 \end{pmatrix}=\begin{pmatrix} 29 & 49 & 32\\ \end{pmatrix}$

となります。

行列ではこの2回の計算をさらに省略して、

$\begin{pmatrix} 40 & 10\\ 20 & 90\\ \end{pmatrix} \begin{pmatrix} 0.1 & 0.2 & 0.7\\ 0.3 & 0.5 & 0.2 \end{pmatrix} = \begin{pmatrix} 7 & 13 & 30 \\ 29 & 49 & 32 \end{pmatrix}$

と書きます。右辺は1行目に文書1についての計算結果、2行目に文書2について計算結果が入ります。

文書の数が増えても左辺の一個目の行列が長くなるだけで、まったく同じように書くことができます。

$\begin{pmatrix} 40 & 10\\ 20 & 90\\ \vdots & \vdots\\ 40 & 70 \end{pmatrix} \begin{pmatrix} 0.1 & 0.2 & 0.7\\ 0.3 & 0.5 & 0.2 \end{pmatrix} = \begin{pmatrix} 7 & 13 & 30 \\ 29 & 49 & 32 \\ \vdots & \vdots & \vdots\\ 25 & 43 & 42 \\ \end{pmatrix}$

以上がトピックモデルのデータ生成過程です。トピックモデルではこれの逆問題をやります。つまり文書と単語の行列がつぎのように与えられたとき、

$\begin{pmatrix} 7 & 13 & 30 \\ 29 & 49 & 32 \\ \vdots & \vdots & \vdots\\ 25 & 43 & 42 \\ \end{pmatrix}$

この行列の分解

$\begin{pmatrix} 40 & 10\\ 20 & 90\\ \vdots & \vdots\\ 40 & 70 \end{pmatrix} \begin{pmatrix} 0.1 & 0.2 & 0.7\\ 0.3 & 0.5 & 0.2 \end{pmatrix}$

を探すのです。

トピックモデルは確率モデルですので本当は

$\begin{pmatrix} 7 & 13 & 30 \\ 29 & 49 & 32 \\ \vdots & \vdots & \vdots\\ 25 & 43 & 42 \\ \end{pmatrix} \approx \begin{pmatrix} 40 & 10\\ 20 & 90\\ \vdots & \vdots\\ 40 & 70 \end{pmatrix} \begin{pmatrix} 0.1 & 0.2 & 0.7\\ 0.3 & 0.5 & 0.2 \end{pmatrix}$

となる行列を探します。 $\approx$ はニアリーイコールです。トピックモデルで得られる行列のかけ算の結果は、もとの行列と厳密には一致しませんが、一致しない分は確率的な誤差とみなすのです。

モデルと変分推論

カウントデータの行列 $Y = (y_{n,k})$ があるとする。これを2つの正の実数の行列 $W=(w_{n,l})$ と $H=(h_{l,k})$ で、
$Y \approx WH$
と近似することを考える。

そのために、以下のような生成モデルを考え、 $w_{n,l}$ と $h_{l,k}$ を推定する。

$y_{n,k} \sim \mathrm{Poisson}(\sum_{l=1}^{L}w_{n,l}h_{l,k})$
$w_{n,l} \sim \mathrm{Gamma}(a,b)$
$h_{l,1:K} \sim \mathrm{Dirichlet}(\alpha_{1:K})$

このままだと対数尤度にしたとき対数の中に足し算が入る形になり、計算しづらい。なので、ポアソン分布の部分は以下のように書き換える。

$u_{n,l} \sim \mathrm{Poisson}(w_{n,l})$
$s_{n,l,1:K} \sim \mathrm{Multinomial}(u_{n,l}, h_{l,1:K})$
$y_{n,k}= \sum_{l=1}^{L}s_{n,l,k}$

ここで $u_{n,l} = \sum_{k=1}^{K} s_{n,l,k}$ である。

こうやって潜在変数 $s_{n,l,k}$ をかますことで平均場近似による近似事後分布は以下のように導ける。

$\log q(w_{n,l}) =( \sum_{k=1}^{K}\hat a-1) \log w_{n,l} -\hat bw_{n,l} + \mathrm{const.}$
$\log q(h_{l,1:K}) =\sum_{k=1}^{K}( \hat \alpha_{l,k}-1) \log h_{l,k} + \mathrm{const.}$

これらはそれぞれガンマ分布、ディリクレ分布である。

ここで
$\hat a =s_{n,l,k} +a$
$\hat b =b+1$
$\hat \alpha_{l,k} = \sum_{n=1}^{N} s_{n,l,k} +\alpha_k$
とおいた。

また、ポアソン分布を和で条件付けてやると多項分布になるので、潜在変数については以下の更新式を用いる。

$E[s_{n,l,k}] = y_{n,k} \exp(E[\log w_{n,l}]E[\log h_{l,k}])/ Z_n$
$Z_n = \sum_{l=1}^{L}\exp(E[\log w_{n,l}]E[\log h_{l,k}])$

以上をそのまま実装してもアルゴリズムは動くが、更新式をもう少し簡単にすることができる。

潜在変数 $s_{n,l,k}$ についての結果を $\sum_{k=1}^{K}s_{n,l,k}$ に代入すると、

$\sum_{k=1}^{K}s_{n,l,k} \\= \sum_{k=1}^{K} y_{n,k} \exp(E[\log w_{n,l}]E[\log h_{l,k}])/\{ \sum_{l=1}^{L}\exp(E[\log w_{n,l}]E[\log h_{l,k}]) \}$

$\sum_{n=1}^{N}s_{n,l,k}$ に代入すると、
$\sum_{n=1}^{N}s_{n,l,k} \\= \sum_{n=1}^{N} y_{n,k} \exp(E[\log w_{n,l}]E[\log h_{l,k}])/\{ \sum_{l=1}^{L}\exp(E[\log w_{n,l}]E[\log h_{l,k}]) \}$

これにより潜在変数 $s_{n,l,k}$ （3次元配列）を保存することなく、パラメータ推定を行うことができる。

添字に注意して整理すると、 $\sum_{k=1}^{K}s_{n,l,k}$ の行列は、

$E[\exp(\log W)] \circ \{(Y/(\exp(E[\log W])\exp(E[\log H])))(E[\log H]^\top)\}$

$\sum_{n=1}^{N}s_{n,l,k}$ の行列は、

$E[\exp(\log H)] \circ \{(E[\log W]^\top)(Y/(\exp(E[\log W])\exp(E[\log H])))\}$

となる。ここで $\circ$ は要素ごとのかけ算、/ は要素ごとのわり算を表す。

ELBO

潜在変数を $Z$ 、観測されるデータを $Y$ としたとき、ELBO は次の式で与えられる.

$\operatorname{ELBO} = E[ \log p (Y|Z) + \log p (Z) -\log q(Z) ]$

ここでの期待値は変分事後分布 $q(Z)$ による期待値。

GaP に出てくる確率（密度）関数の対数をすべて書き下す。

$\log p(u_{n,l}) = u_{n,l} \log w_{n,l} - w_{n,l} - \log \Gamma(u_{n,l}+1)$

$\log p(s_{n,l,1:K}) = \sum_k s_{n,l,k} \log h_{l,k} + \log \Gamma(u_{n,l}+1)- \sum_k \log \Gamma(s_{n,l,k}+1)$

$\log q(s_{n,l,1:K}) = s_{n,l,1:K} \log (\exp(E[\log w_{n,l}]E[\log h_{l,k}])/ Z_n) + \log \Gamma (y_{n,k}+1) \\- \sum_l \log \Gamma (s_{n,l,k}+1)$

$\log p(w_{n,l}) = (a-1) \log w_{n,l} -b w_{n,l} + a\log b -\log \Gamma(a)$

$\log q(w_{n,l}) = (\hat a-1) \log w_{n,l} -\hat b w_{n,l} + a\log b -\log \Gamma(a)$

$\log p (h_{l,1:K}) = \sum_{k} (\alpha_k-1) \log h_{l,k} + \log \Gamma (\sum_{k} \alpha_k) - \sum_{k} \log \Gamma(\alpha_k)$

$\log q (h_{l,1:K}) = \sum_{k} (\hat \alpha_{l,k}-1) \log h_{l,k} + \log \Gamma (\sum_{k} \hat \alpha_{l,k}) \\- \sum_{k} \log \Gamma(\hat \alpha_{l,k})$

期待値の計算の難しい $\log \Gamma(s_{n,l,k}+1)$ と $\log \Gamma(u_{n,l}+1)$ の項はキャンセルされて消える。

また、
$E[u_{n,l}] + a = \hat a$
$b+1 = \hat b$
$E[\sum_{n=1}^{N} s_{n,l,k} +\alpha_k] = \hat \alpha_{l,k}$
であるから、
$(\hat a-1) \log w_{n,l}$ と、
$\hat b w_{n,l}$ と、
$\sum_{k} (\hat \alpha_{l,k}-1) \log h_{l,k}$
の項もキャンセルされて消える。

また、
$E[s_{n,l,k} ]\log (\exp(E[\log w_{n,l}]E[\log h_{l,k}])/ Z_n)\\ = y_{n,k} \exp(E[\log w_{n,l}]E[\log h_{l,k}])/ Z_n ) \log (\exp(E[\log w_{n,l}]E[\log h_{l,k}])/ Z_n)\\ = y_{n,k} [ (\{ \exp( E[\log w_{n,l}] ) E[\log w_{n,l}] \} \exp(E[\log h_{l,k}]) +\\ \exp( E[\log w_{n,l}] ) \{\exp(E[\log h_{n,l}]) E[\log h_{l,k}]\})/ Z_n -\log Z_n ]$

R による実装例

GaPVB <-function(Y,L=2,a=0.5,b=1e-8,alpha=rep(1,ncol(Y)),maxit=500){
  N <- nrow(Y)
  K <- ncol(Y)
  EelW <- matrix(rgamma(N*L,shape=a,rate=1),N,L)
  unnorm <- matrix(alpha+colSums(Y),L,K,byrow = TRUE)
  EelH <- unnorm/rowSums(unnorm)
  log1pb <- log1p(b)
  for(iter in 1:maxit){
    Sw <- EelW * (((Y)/(EelW %*% EelH)) %*% t(EelH))
    Sh <- EelH * (t(EelW) %*% (Y/(EelW %*% EelH)))
    alpha_W <- a + Sw
    alpha_H <- alpha + Sh
    den <- rowSums(alpha_H)
    EelW <-exp(digamma(alpha_W)-log1pb)
    EelH <-exp(digamma(alpha_H)-digamma(den))
  }
  EW <- alpha_W/(b+1)
  EH <- alpha_H/den
  ELBO <- -sum(lfactorial(Y))+
    L*sum(lgamma(sum(alpha))-sum(lgamma(alpha))) +
    sum(-lgamma(rowSums(alpha_H))+rowSums(lgamma(alpha_H))) + 
    sum(-alpha_W*log1pb+a*log(b)+lgamma(alpha_W)-lgamma(a)) + 
    sum(-Y*(((EelW*log(EelW))%*%EelH + EelW%*%(EelH*log(EelH)))/(EelW%*%EelH)-log(EelW%*%EelH))) 
  return(list(W=EW,H=EH,ELBO=ELBO))
}

簡単なシミュレーションを行ってみる。

N <- 100
K <- 10
L <- 3
W <- matrix(rgamma(N*L,shape=1,scale=1000),N,L)
W <- W[,order(W[1,])]
unnormH <- matrix(rgamma(L*K,shape=1,rate=1),L,K)
H <- unnormH/rowSums(unnormH)

Y <- matrix(rpois(N*K,W%*%H),N,K)

out <- NMFVB(Y,L=L,a=0.5,b=0)


plot(out$W[,order(out$W[1,])],W, main="W", xlab = "estimates", ylab = "true value")
abline(0,1,lty=2)
plot(out$H[order(out$W[1,]),],H, main="H", xlab = "estimates", ylab = "true value")
abline(0,1,lty=2)

f:id:abrahamcow:20190323162911p:plain

f:id:abrahamcow:20190323162922p:plain

Wのスケールが大きい場合、シミュレーションで設定した真値と推定値がそこそこ一致することがわかる。

ただしラベルスイッチングの問題があるので、 $w_{11} < w_{12} < w_{13}$ となるようソートし直している。

シミュレーション

トピック数を5にしてデータを生成した。

library(parallel)
simfunc <- function(i){
  set.seed(i)
  W <- matrix(rgamma(100*5,1,0.001),100,5)
  H <- matrix(rgamma(5*20,1,1),5,20)
  H <- H/rowSums(H)
  Y <- matrix(rpois(100*20,W%*%H),100,20)
  sapply(1:10, function(L)GaPVB(Y,L)$ELBO)
}
out <- mclapply(1:100,simfunc,mc.cores = detectCores())
print(table(apply(simplify2array(out),2,which.max)))
boxplot(t(simplify2array(out)))

f:id:abrahamcow:20190830221703p:plain

100回のシミュレーション中100回とも正しいトピック数を選んでいる。

> print(table(apply(simplify2array(out),2,which.max)))

  5 
100

<br />

ジョンとヨーコのイマジン日記

キョウとアンナのラヴラブダイアイリー改め、ジョンとヨーコのイマジン日記です。

検定いらずのABテスト：ポアソン分布とベータ分布によるサンプルサイズ設計

実装例

トピックモデル（GaP; Gamma-Poisson Model）の ELBO の導出

行列の復習からはじめてトピックモデルの結果だけ理解する

モデルと変分推論

ELBO

R による実装例

シミュレーション