統計コンサルの議事メモ

統計や機械学習の話題を中心に、思うがままに

Bayesian t-testを書いてみたい

IBM SPSS Statisticsにベイズ推論が実装されるらしい。

news.mynavi.jp

これ以外にも最近は身の回りでやけにBayesian t-testに関する話題を耳にすることが多いので、t-検定をベイズ化することのどこにニーズがあるのかを考えていたのだけれど、どうやらA/B testにおいて結果を途中で確認(Peeking)したり、サンプルサイズを変更したいという状況が結構多いそうだ。

古典的な、そして計画立てられたt-検定では分析の事前に必要なサンプルサイズを求めておく。通常は実験の開始以降にデータを追加したりなどしないのだが、ビジネスではA/Bテストの結果の如何によってデータを追加したり、あるいは途中経過を見て早めに施策の実施判断を行いたいと状況がよくあるため、そういったニーズにマッチしたのだろう(実際のところ伝統的なt-検定を使う場合であってもデータの追加や結果の覗き見は行われていると思うが)。

さて、Rにおけるベイジアンモデリングのためのパッケージとしては{MCMCpack}があり、より複雑なモデルを構築したいのであればStanを呼び出す方法があげられる。しかしt-検定ぐらい自分で書いてみようと思ったので、以下にスクリプトを示す。

データの準備

# パッケージの読み込み
library(tidyverse)
library(ggplot2)

サンプルデータとして以下のような正規分布に従う2標本を用意する。数値に意味はないので関心があれば標本サイズなどをそれぞれ変更してみると良いが、目的は標本の平均の差の分布について確認することである。

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()

f:id:ushi-goroshi:20170830224339p:plain

このデータに対して古典的な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])

f:id:ushi-goroshi:20170830231541p:plain

要約統計量も見ておこう。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  

単純な推定なので精度よく推定できているようである。

平均値の差の分布

では最後に、この分析の目的であった2標本の平均の差の分布について確認してみよう。

par(mfrow = c(1, 1))
hist(Res_MCMC[, 1] - Res_MCMC[, 3])

f:id:ushi-goroshi:20170830232942p:plain

これを見ると、平均値の差の分布が比較的中心に近い位置で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-検定のベイズ化は終わりだが、ベイズ化した分析のメリットとして以下のような疑問に答えることが可能となる:1つ目の標本の平均が2つ目の標本の平均よりも大きい確率はどの程度だろうか?

> mean(Res_MCMC[, 3] - Res_MCMC[, 1] < 0)
[1] 0.2031979

上記のような質問は、古典的な仮説検定(頻度論的アプローチ)では答えることができない。仮説検定は常に「帰無仮説が棄却されるか否か」(例えば平均値と平均値に差が0であるか否か)であり、「何%の確率で棄却される」とは表現できないのだ。しかしこれをベイズ化することでより直感的な表現が可能となる。これはビジネスへの応用において確かにメリットになると感じた。