表と裏、投げるとどちらかの面が出るコインを例に、ベータ分布について説明します。
ベータ分布とは
確率変数 が、 をパラメータとする確率密度関数
を持つとき、 はベータ分布 に従う、と言います。
分母に出てきた はベータ関数で、ベータ分布を積分した結果が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))
の場合
より、 。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))
の場合
より、 。上に凸の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))
ベータ分布を使ったパラメータ推定
ここからがコインの出番です。コインには表裏があり、コインを投げたときに表の出る確率 を推定するのがパラメータ推定です。パラメータ推定には、2つの方法があります。
最尤推定
ベイズ推定
ベータ分布が活躍するのが、ベイズ推定です。ベイズ推定では事前確率と事後確率という考え方を使って、コインを投げるたびに事前確率を事後確率へ、ベイズ更新と呼ばれる計算で更新することができます。ベイズの定理を使うとき、事前分布がベータ分布なら、事後分布もベータ分布となり計算が簡単です。この関係は、共役事前分布と呼ばれます。
最尤推定とベイズ推定
最尤推定は、コインを投げたときに表の出る確率 の推定値を唯一に特定する点推定です。
最尤推定における確率の考え方を頻度主義と言い、確率を無限回の試行という理想的な状況で達成できる理論値として考えており、この考え方で導かれる確率を客観確率と呼びます。
ベイズ推定では、確率を確信の度合いとして定義し、事前確率(分布)から事後確率(分布)への変化を、確信度の変化と考えます。
ベイズ推定における確率の考え方をベイズ主義と言い、ベイズ主義で導かれる確率を主観確率と呼びます。
一応、書いてみましたが、これらの主義について覚えておいても役に立つわけではありません。2つの主義が争っているわけでもなく、私は確率って言うと2通りの解釈があるんだなーくらいに理解しています。コインを例に2つの考え方で確率を求めてみます。次のような問題を考えます。
最尤推定では、どちらの場合も同じ確率を算出します。 と 。一方、ベイズ推定では、観測回数0のときに、どの も等しく確からしいという一様分布の状態から、観測回数を増やすと、真の値近くにピークを持つベータ分布に変化していきます。
コインの問題を一般化して、
求めたいのは、観測した結果 が得られたもとでの、 の条件付き確率なので、
コインを1回投げて表になる確率が 。各試行は独立であり、観測結果 が得られる確率 は、
事前分布に一様分布を仮定すると、
最終的にパラメータ の事後分布は、をパラメータとするベータ分布の確率密度関数で表されます。
太字にした「各試行は独立」「事前分布に一様分布を仮定」となっていないと、現実の問題に適用すると無理が出てきます。
最初の問題に戻り、2通りの試行で事後分布がどのように異なるのか見てみましょう。それぞれ、 で表される分布となります。
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))
回数が多くなることで、 の確信度が高くなっていることが確かめられます。この確信度は最尤推定では得られません。
表の出る確率 のコインの出た面を二項分布にしたがって生成し、どのように収束していくのか見てみましょう。
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) } }
さらにコイン投げの回数を増やして動画にしてみましょう。
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"))
回数を重ねると、観測者が知らない の値付近にピークが細く収束していく様子がわかります。
参考書籍
- 作者: 石井健一郎,上田修功
- 出版社/メーカー: オーム社
- 発売日: 2014/08/26
- メディア: 単行本(ソフトカバー)
- この商品を含むブログ (2件) を見る
- 作者: C.M.ビショップ,元田浩,栗田多喜夫,樋口知之,松本裕治,村田昇
- 出版社/メーカー: 丸善出版
- 発売日: 2012/04/05
- メディア: 単行本(ソフトカバー)
- 購入: 6人 クリック: 33回
- この商品を含むブログ (18件) を見る
- 作者: C.M.ビショップ,元田浩,栗田多喜夫,樋口知之,松本裕治,村田昇
- 出版社/メーカー: 丸善出版
- 発売日: 2012/02/29
- メディア: 単行本
- 購入: 6人 クリック: 14回
- この商品を含むブログを見る
- 作者: デイヴィッド・サルツブルグ,竹内惠行、熊谷悦生
- 出版社/メーカー: 日本経済新聞出版社
- 発売日: 2010/04/01
- メディア: 文庫
- 購入: 16人 クリック: 320回
- この商品を含むブログ (34件) を見る
- 作者: 渡辺澄夫
- 出版社/メーカー: コロナ社
- 発売日: 2012/03
- メディア: 単行本
- 購入: 1人 クリック: 4回
- この商品を含むブログ (5件) を見る
参考になるリンク
ベイズ統計と主観確率・客観確率について、とてもわかりやすい日本語でまとめられていて勉強になります。