IBM SPSS Statisticsにベイズ推論が実装されるらしい。
これ以外にも最近は身の回りでやけにBayesian t-testに関する話題を耳にすることが多いので、t-検定をベイズ化することのどこにニーズがあるのかを考えていたのだけれど、どうやらA/B testにおいて結果を途中で確認（Peeking）したり、サンプルサイズを変更したいという状況が結構多いそうだ。
古典的な、そして計画立てられたt-検定では分析の事前に必要なサンプルサイズを求めておく。通常は実験の開始以降にデータを追加したりなどしないのだが、ビジネスではA/Bテストの結果の如何によってデータを追加したり、あるいは途中経過を見て早めに施策の実施判断を行いたいと状況がよくあるため、そういったニーズにマッチしたのだろう（実際のところ伝統的なt-検定を使う場合であってもデータの追加や結果の覗き見は行われていると思うが）。
さて、Rにおけるベイジアンモデリングのためのパッケージとしては{MCMCpack}があり、より複雑なモデルを構築したいのであればStanを呼び出す方法があげられる。しかしt-検定ぐらい自分で書いてみようと思ったので、以下にスクリプトを示す。
データの準備
# パッケージの読み込み library(tidyverse) library(ggplot2)
サンプルデータとして以下のような正規分布に従う２標本を用意する。数値に意味はないので関心があれば標本サイズなどをそれぞれ変更してみると良いが、目的は標本の平均の差の分布について確認することである。
N_01 <- 100 mu_01 <- 100 sig_01 <- 10 N_02 <- 100 mu_02 <- 103 sig_02 <- 10 set.seed(123) Pop_01 <- rnorm(N_01, mu_01, sig_01) Pop_02 <- rnorm(N_02, mu_02, sig_02)
以下のようなサンプルが得られるだろう：
> Pop_01 %>% as_data_frame() %>% mutate("Group" = "01") %>% bind_rows(as_data_frame(Pop_02) %>% mutate("Group" = "02")) %>% ggplot(aes(value, fill = Group)) + geom_histogram(alpha = 0.5, position = "identity") + theme_bw()
このデータに対して古典的なt-検定を実施すると、平均値の差が0であるという帰無仮説はp = 0.44で棄却されない。ちなみにt.testはデフォルトでは等分散性を仮定しないWelchのt-検定となる。
> t.test(Pop_01, Pop_02) Welch Two Sample t-test data: Pop_01 and Pop_02 t = -0.7674, df = 197.35, p-value = 0.4438 alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: -3.642862 1.601916 sample estimates: mean of x mean of y 100.9041 101.9245
関数の定義
このデータに対して何らかのパラメータ（ここでは平均と標準偏差）を与えた時の尤度を求める関数を以下のように定義する：
likelihood <- function(pars) { mu_01 <- pars[1] sig_01 <- pars[2] mu_02 <- pars[3] sig_02 <- pars[4] sum( dnorm(Pop_01, mu_01, sig_01, log = TRUE), dnorm(Pop_02, mu_02, sig_02, log = TRUE) ) }
データの尤度は0 ~ 1の数字の乗算となるため標本サイズが大きい場合には極めて小さな数値となり、桁あふれが発生してしまう。そのためdnormにおいてlog = TRUEとし、対数尤度の総和を求めることとする。
この関数を用いて最尤法で平均と分散を求めてみる。optimを使用するが、デフォルトは最小化なのでfnscaleを負の数として最大化を行う：
> pars <- c(100, 10, 100, 20) > optim(par = optim(par = pars, fn = likelihood, control=list(fnscale = -1))$par, + fn = likelihood, control=list(fnscale = -1)) $par [1] 100.904769 9.082697 101.923277 9.620816 $value [1] -730.8205 $counts function gradient 101 NA $convergence [1] 0 $message NULL
続いて事前分布を定義する。事前分布の設定は色々と工夫すべきポイントではあるが、簡単のため-10000 ~ 10000の一様分布とした。尤度はデータと同じく対数としてある。
prior <- function(pars) { mu_01 <- pars[1] sig_01 <- pars[2] mu_02 <- pars[3] sig_02 <- pars[4] sum( dunif(mu_01, -10000, 10000, log = TRUE), dunif(mu_02, -10000, 10000, log = TRUE), dunif(sig_01, -10000, 10000, log = TRUE), dunif(sig_02, -10000, 10000, log = TRUE) ) }
事後分布は事前分布と尤度の積であるが、対数を取ってあるためここでは和となる。
posterior <- function(pars) { likelihood(pars) + prior(pars) }
MCMCの実行
以上で準備ができたので以下のようにMCMCを実行する。なおサンプラーとしてGibbs Samplerを実装したつもりだが、Gibbs Samplingにおいてパラメータの採択確率があったかは記憶が不確かなので、この辺りは後日確認しておくとしてひとまず進めることとする。サンプルの数は30,000としたが、実際には割と早く収束しているようなのでこんなには必要ない。
N_iter <- 30000 pars <- c(50, 30, 20, 50) results <- matrix(0, nrow = N_iter, ncol = 8) results[1, ] <- c(pars, 0, 0, 0, 0) for (i in 2:N_iter) { for (j in 1:4) { candidate <- pars candidate[j] <- candidate[j] + rnorm(1, sd = 0.5) ratio <- exp(posterior(candidate) - posterior(pars)) t <- runif(1) if (t < ratio) { pars <- candidate } } results[i, ] <- c(pars, ratio, t, posterior(pars), posterior(candidate)) }
以下のような結果となる。
> par(mfrow = c(2, 2)) > plot(results[, 1]) > plot(results[, 2]) > plot(results[, 3]) > plot(results[, 4])
要約統計量も見ておこう。burn-inとして最初の15,000サンプルを切り捨て、自己回帰性を小さくするために10刻みでサンプルを取得する。
> Res_MCMC <- results[seq(15000, N_iter, 10), ] > summary(Res_MCMC[, c(1:4)]) V1 V2 V3 V4 Min. : 97.79 Min. : 7.198 Min. : 98.41 Min. : 7.997 1st Qu.:100.24 1st Qu.: 8.732 1st Qu.:101.34 1st Qu.: 9.266 Median :100.86 Median : 9.191 Median :101.98 Median : 9.753 Mean :100.86 Mean : 9.233 Mean :101.98 Mean : 9.774 3rd Qu.:101.45 3rd Qu.: 9.675 3rd Qu.:102.63 3rd Qu.:10.218 Max. :104.16 Max. :11.959 Max. :105.20 Max. :12.632
単純な推定なので精度よく推定できているようである。
平均値の差の分布
では最後に、この分析の目的であった２標本の平均の差の分布について確認してみよう。
par(mfrow = c(1, 1)) hist(Res_MCMC[, 1] - Res_MCMC[, 3])
これを見ると、平均値の差の分布が比較的中心に近い位置で0を跨いでおり、伝統的なt-検定で言う所の有意であるとは言えなさそうだ。実際に平均値の差の経験分布を確認すると、95%確信区間は-3.76 ~ 1.62である。
> quantile(Res_MCMC[, 1] - Res_MCMC[, 3], seq(0, 1, 0.025)) 0% 2.5% 5% 7.5% 10% 12.5% 15% 17.5% -5.62647895 -3.76661174 -3.35462735 -3.11314731 -2.89345678 -2.70960554 -2.56386382 -2.43985754 20% 22.5% 25% 27.5% 30% 32.5% 35% 37.5% -2.30487729 -2.18744684 -2.07270632 -1.98450427 -1.89198288 -1.77273198 -1.67356108 -1.57205147 40% 42.5% 45% 47.5% 50% 52.5% 55% 57.5% -1.47598439 -1.40161516 -1.30787797 -1.21558957 -1.12424331 -1.03598088 -0.95253546 -0.85994024 60% 62.5% 65% 67.5% 70% 72.5% 75% 77.5% -0.78004294 -0.68704363 -0.59508922 -0.49850858 -0.42349299 -0.30716666 -0.19718374 -0.10822330 80% 82.5% 85% 87.5% 90% 92.5% 95% 97.5% 0.01189818 0.14609443 0.30742495 0.45296435 0.69891710 0.92152549 1.18895282 1.62042273 100% 3.58972078
改めてt-検定の結果を見てみると、95%信頼区間として同様の結果が得られていることがわかるだろう。
> t.test(Pop_01, Pop_02) Welch Two Sample t-test data: Pop_01 and Pop_02 t = -0.7674, df = 197.35, p-value = 0.4438 alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: -3.642862 1.601916 sample estimates: mean of x mean of y 100.9041 101.9245
おまけ
以上でt-検定のベイズ化は終わりだが、ベイズ化した分析のメリットとして以下のような疑問に答えることが可能となる：１つ目の標本の平均が２つ目の標本の平均よりも大きい確率はどの程度だろうか？
> mean(Res_MCMC[, 3] - Res_MCMC[, 1] < 0) [1] 0.2031979
上記のような質問は、古典的な仮説検定（頻度論的アプローチ）では答えることができない。仮説検定は常に「帰無仮説が棄却されるか否か」（例えば平均値と平均値に差が0であるか否か）であり、「何%の確率で棄却される」とは表現できないのだ。しかしこれをベイズ化することでより直感的な表現が可能となる。これはビジネスへの応用において確かにメリットになると感じた。