ふんわり R-tips

ぜんぜんわからない、俺たちは雰囲気でRをやっている

コインで理解するベータ分布

表と裏、投げるとどちらかの面が出るコインを例に、ベータ分布について説明します。

ベータ分布とは

確率変数 \theta が、\alpha, \beta をパラメータとする確率密度関数


\mathrm{Be}(\alpha, \beta) = \frac{\theta^{\alpha-1} (1-\theta)^{\beta-1}}{B(\alpha, \beta)}


を持つとき、\theta はベータ分布 \mathrm{Be}(\alpha, \beta) に従う、と言います。

分母に出てきた B(\alpha, \beta) はベータ関数で、ベータ分布を積分した結果が1になるために必要な正規化項です。\alpha, \beta を与えられると、ベータ関数は定数値を返します。

ベータ分布の期待値は、E(\theta) = \frac{\alpha}{\alpha + \beta} で計算されます。

ベータ関数とガンマ関数


B(\alpha, \beta) = \int_0^1 u^{\alpha - 1} (1-u)^{\beta - 1} du\\
= \frac{\Gamma(\alpha) \Gamma(\beta)}{\Gamma(\alpha + \beta)}


ガンマ関数 \Gamma(\alpha) は次の式で計算されます。


\Gamma(\alpha) = \int_{0}^{\infty} u^{\alpha - 1} e^{-u} du


\Gamma(\alpha + 1) = \alpha \Gamma(\alpha)より、ガンマ関数のパラメータが正の整数の場合、\Gamma(\alpha + 1) = \alpha !となり階乗関数として扱うことができます。

ベータ分布の例

 \alpha \beta の値を具体的に与えると、ベータ分布が一つに決まります。

 \alpha = 1, \beta = 1の場合

 B(1,1) = 1 より、 \mathrm{Be}(1,1) = 1 。一様分布となります。

a <- 1
b <- 1
x <- seq(0.01, 1.0, len = 500)
plot(x, dbeta(x,a,b),type = "l",col=a+b,ylim=c(0,5))

f:id:phmpk:20170502165747p:plain

 \alpha = 2, \beta = 1の場合

 B(2,1) = \frac{1}{2} より、 \mathrm{Be}(2,1) = 2 \theta 。1次関数のような分布となります。。

a <- 2
b <- 1
x <- seq(0.01, 1.0, len = 500)
plot(x, dbeta(x,a,b),type = "l",col=a+b,ylim=c(0,5))

f:id:phmpk:20170502165545p:plain

 \alpha = 2, \beta = 2の場合

 B(2,2) = \frac{1}{6} より、 \mathrm{Be}(2,2) = 6 \theta (1 - \theta) 。上に凸の2次関数のような分布となります。

a <- 2
b <- 2
x <- seq(0.01, 1.0, len = 500)
plot(x, dbeta(x,a,b),type = "l",col=a+b,ylim=c(0,5))

f:id:phmpk:20170502165405p:plain

ベータ分布を使ったパラメータ推定

ここからがコインの出番です。コインには表裏があり、コインを投げたときに表の出る確率 \theta を推定するのがパラメータ推定です。パラメータ推定には、2つの方法があります。

  • 最尤推定

  • ベイズ推定

ベータ分布が活躍するのが、ベイズ推定です。ベイズ推定では事前確率と事後確率という考え方を使って、コインを投げるたびに事前確率を事後確率へ、ベイズ更新と呼ばれる計算で更新することができます。ベイズの定理を使うとき、事前分布がベータ分布なら、事後分布もベータ分布となり計算が簡単です。この関係は、共役事前分布と呼ばれます。

最尤推定とベイズ推定

最尤推定は、コインを投げたときに表の出る確率 \theta の推定値を唯一に特定する点推定です。

最尤推定における確率の考え方を頻度主義と言い、確率を無限回の試行という理想的な状況で達成できる理論値として考えており、この考え方で導かれる確率を客観確率と呼びます。

ベイズ推定では、確率を確信の度合いとして定義し、事前確率(分布)から事後確率(分布)への変化を、確信度の変化と考えます。

ベイズ推定における確率の考え方をベイズ主義と言い、ベイズ主義で導かれる確率を主観確率と呼びます。

一応、書いてみましたが、これらの主義について覚えておいても役に立つわけではありません。2つの主義が争っているわけでもなく、私は確率って言うと2通りの解釈があるんだなーくらいに理解しています。コインを例に2つの考え方で確率を求めてみます。次のような問題を考えます。

表の出る確率 \theta が0.6のコインがある。この確率を観測者は知らず、未知のパラメータ \theta として与えられているとする。以下の2通りの試行で、\theta を推定せよ。 1. コインを10回投げて、表が6回出た。 2. コインを30回投げて、表が18回出た。


最尤推定では、どちらの場合も同じ確率を算出します。\frac{6}{10} = 0.6\frac{18}{30} = 0.6 。一方、ベイズ推定では、観測回数0のときに、どの \theta も等しく確からしいという一様分布の状態から、観測回数を増やすと、真の値近くにピークを持つベータ分布に変化していきます。

コインの問題を一般化して、

表の出る確率 \theta (0 {<} \theta {<} 1) が未知のコインがある。コインを n 回投げた結果が、x^{(n)} = x_1, x_2, \cdots, x_t, \cdots, x_n であったとし、n 回のうち、r 回が表であった。この結果をもとに、\theta を推定せよ。


求めたいのは、観測した結果 x^{(n)} が得られたもとでの、\theta の条件付き確率なので、


p(\theta | x^{(n)}) = \frac{p(x^{(n)} | \theta)}{p(x^{(n)})} \cdot p(\theta)



p(x^{(n)}) = \int_{0}^{1} p(x^{(n)} | \theta) p(\theta) d \theta


コインを1回投げて表になる確率が \theta各試行は独立であり、観測結果 x^{(n)} が得られる確率  p(x^{(n)} | \theta) は、


 p(x^{(n)} | \theta) = \prod_{t=1}^{n} p(x_t | \theta) = \theta^{r} (1 - \theta)^{n-r}


事前分布に一様分布を仮定すると、


p(x^{(n)}) = \int_{0}^{1} \theta^{r} (1 - \theta)^{n-r} d \theta = B(r+1, n-r+1) \\
= \frac{r! (n-r)!}{(n+1)!} = \frac{1}{(n+1) \cdot {}_n \mathrm{C}_{r}}


最終的にパラメータ \theta の事後分布は、n, rをパラメータとするベータ分布の確率密度関数で表されます。


p(\theta | x^{(n)}) = (n + 1) \cdot {}_{n} \mathrm{C}_{r} \theta^{r} (1 - \theta)^{n - r}


太字にした「各試行は独立」「事前分布に一様分布を仮定」となっていないと、現実の問題に適用すると無理が出てきます。

最初の問題に戻り、2通りの試行で事後分布がどのように異なるのか見てみましょう。それぞれ、\mathrm{Be}(7,5), \mathrm{Be}(19,13) で表される分布となります。

x <- seq(0.01, 1.0, len = 500)
plot(x, dbeta(x,7,5),type = "l",col=1,ylim=c(0,5))
par(new=T)
plot(x, dbeta(x,19,13),type = "l",col=i,xlab="",ylab="",axes=F,ylim=c(0,5))

f:id:phmpk:20170502165220p:plain

回数が多くなることで、\theta の確信度が高くなっていることが確かめられます。この確信度は最尤推定では得られません。

表の出る確率 \theta = 0.6 のコインの出た面を二項分布にしたがって生成し、どのように収束していくのか見てみましょう。

set.seed(30)
data <- rbinom(30,1,0.6)
data
# [1] 0 1 1 0 1 0 1 1 1 1 1 0 0 1 1 0 0 1 0 1 1 1 1 1 1 0 1 1 1 0 1 1 0 1 0 0 1 1 1 0 0 1
# [43] 1 0 0 1 0 0 1 1 1 1 0 0 1 1 1 1 0 0 1 1 1 0 1 0 1 1 1 1 1 0 1 0 1 1 1 0 0 0 1 1 0 0
# [85] 1 1 1 1 0 0 0 1 0 1 0 1 1 1 0 1

二項分布に従ってコインの表と裏の観測結果を生成します。

a <- 0
b <- 0
x <- seq(0.01, 1.0, len = 500)
for (i in 1:length(data)) {
  if(data[i] == 1) { a = a + 1 } 
  else             { b = b + 1 }
  if(i>3){
  if(i==4){ plot(x, dbeta(x,a,b),type = "l",col=i,ylim=c(0,5)) }
  else    { plot(x, dbeta(x,a,b),type = "l",col=i,xlab="",ylab="",axes=F,ylim=c(0,5)) }
  par(new=T)
  }
}

f:id:phmpk:20170502164510p:plain

さらにコイン投げの回数を増やして動画にしてみましょう。

library(magick)
a <- 0
b <- 0
data <- rbinom(100,1,0.6)
png(file="example%03d.png", width=200, height=200)
for (i in 1:length(data)) {
  if(data[i] == 1) { a = a + 1 } 
  else             { b = b + 1 }
  if(i>3){
  if(i==4){ plot(x, dbeta(x,a,b),type = "l",ylim=c(0,8)) }
  else    { plot(x, dbeta(x,a,b),type = "l",ylim=c(0,8)) }
  }
}
dev.off()
lf <- list.files(full.names = T)
pngs <- lapply(lf, image_read)
animation <- image_animate(image_join(pngs))
image_write(animation, "rbinom.gif")
file.remove(list.files(pattern=".png"))

f:id:phmpk:20170502162956g:plain

回数を重ねると、観測者が知らない \theta の値付近にピークが細く収束していく様子がわかります。

参考書籍

続・わかりやすいパターン認識―教師なし学習入門―

続・わかりやすいパターン認識―教師なし学習入門―

パターン認識と機械学習 上

パターン認識と機械学習 上

  • 作者: C.M.ビショップ,元田浩,栗田多喜夫,樋口知之,松本裕治,村田昇
  • 出版社/メーカー: 丸善出版
  • 発売日: 2012/04/05
  • メディア: 単行本(ソフトカバー)
  • 購入: 6人 クリック: 33回
  • この商品を含むブログ (18件) を見る

パターン認識と機械学習 下 (ベイズ理論による統計的予測)

パターン認識と機械学習 下 (ベイズ理論による統計的予測)

  • 作者: C.M.ビショップ,元田浩,栗田多喜夫,樋口知之,松本裕治,村田昇
  • 出版社/メーカー: 丸善出版
  • 発売日: 2012/02/29
  • メディア: 単行本
  • 購入: 6人 クリック: 14回
  • この商品を含むブログを見る

統計学を拓いた異才たち(日経ビジネス人文庫)

統計学を拓いた異才たち(日経ビジネス人文庫)

  • 作者: デイヴィッド・サルツブルグ,竹内惠行、熊谷悦生
  • 出版社/メーカー: 日本経済新聞出版社
  • 発売日: 2010/04/01
  • メディア: 文庫
  • 購入: 16人 クリック: 320回
  • この商品を含むブログ (34件) を見る

ベイズ統計の理論と方法

ベイズ統計の理論と方法

参考になるリンク

ill-identified.hatenablog.com

ベイズ統計と主観確率・客観確率について、とてもわかりやすい日本語でまとめられていて勉強になります。