共役分布　ベルヌーイ分布の場合 ( ソフトウェア ) - OneSegの長寿と健康を願うひとりがき

ベルヌーイな世界の有様

　ベルヌーイな世界では　事象には　0　OR 　1　の札が付けられている。

例えば　A党支持派　なら 1　そうでなければ 0 という具合である。

　統計ソフトRでは　

　　ｘ1 <- rbinom(30,1,0.3)

　と書いて　3割の支持率がある世界から　30人のサンプルを抽出したことになる。　

　例えば　x1　は
[1] 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 1 0 1 0
である。　サンプルからの支持率は　6/30=0.2…　となるが　0.3からはだいぶ遠い。

世界の実情（支持率=0.3という現実）はなかなか測りがたいが

Bern<-function(x,lambda){
b<-1
for(xi in x) b<- b* lambda^xi*(1-lambda)^(1-xi)
b}　

　から支持率＝lambda　の場合に、上のようなデータが抽出される確率を推定できる。尤度（ゆうど）　likelifoodと呼んでいる。

データのみからの推計

f<-function(l) -Bern(x1,l)

optimize(f,c(0,1))$minimum

[1] 0.2000167

　lambda=0.2 のとき　データはもっともらしいということになる。

もっともサンプル数を増やせばBernを最大にするｌは0.3に近づいてくる。

xk<-rbinom(1000,1,.3)
f<-function(l)-(Bern(xk,l))
optimize(f,c(0,1))$minimum

[1] 0.3039811

尤度のグラフを描けば以下のようになる。

Likelifood.Bern<-function(x){
Lf<-NULL
for(lambda in seq(0,1,0.01)){
l<-(Bern(x,lambda))
Lf<-rbind(Lf,c(lambda,l)) }
Lf<-data.frame(lambda=Lf[,1],Likelifood=Lf[,2])
s<-sum(Lf$Likelifood)
Lf$Likelifood<-Lf$Likelifood/s
Lf}
#
plot(Likelifood.Bern(xk),type="l",col=3,main="Likelihood distribution")
lines(Likelifood.Bern(x1),col=4)

　青が n=30　緑が n=1000 の場合である。サンプル数が増えると隠された世界の値(0.3)に近づいてゆく。

思い込みと事後分布

　あらかじめ事前分布（経験　言い伝え　思い込み　占い　信念　お告げ　など）が知られていれば　パラメータまたは母数（ここでは平均値）の推定（事後分布）は効率的に（あまり手間と金をかけずに）できるだろう(できたらいいな)。

　ベーズによれば

　Posterior = Likelihood * Prior / Evidence

である。

　Likelihood = Π lambda^xi*(1-lambda)^(1-xi)

であったから　Priorに数学的に形の似たBetaを使えば計算が楽になるのである（理論的に正しいと言っているわけではない）。

　Prior = beta(a,b)*renge^(a-1)*(1-renge)^(b-1)

もっともらしくConjugate Distributionなどと呼ぶ。

Beta<-function(a,b,renge=seq(0,1,0.01)){
prob<-beta(a,b)*renge^(a-1)*(1-renge)^(b-1)
data.frame(lambda=renge,prb=prob)}

plot(Beta(a=20,b=40),type="l",col=2,main="Distribution of the prior a=20 b=40")

このとき（a=20,b=40 から得られる値はもっともらしく悪くない経験値であり）Posteriorの分布は以下のようである。

Posterior.Bern<-function(x,a,b,renge){
prob<-NULL
for(l in renge){
    p<-Bern(x=x,lambda=l)*Beta(a=a,b=b,renge=l)
    prob<-c(prob,p$prb)}
sp<-sum(prob)      #1/k(x,c,d) const
prob<- prob/sp     #Beta(c,d) prob
data.frame(lambda=renge,prob=prob)}
#

Post1<-Posterior.Bern(x=x1,a=20,b=40,renge=seq(0,1,0.01))
Postk<-Posterior.Bern(x=xk,a=20,b=40,renge=seq(0,1,0.01))

plot(Postk$lambda,Postk$prob,type="l",col=3,main="Posterior distribution a=20, b=40")

lines(Post1$lambda,Post1$prob,col=4)

　n = 30　でもLikelihoodよりはるかにましになっている。

　ところが予測がいい加減だととんでもないことになる。

Post1<-Posterior.Bern(x=x1,a=20,b=4,renge=seq(0,1,0.01))
Postk<-Posterior.Bern(x=xk,a=20,b=4,renge=seq(0,1,0.01))
plot(Postk$lambda,Postk$prob,type="l",col=3,main="Posterior distribution a=20, b=4")
lines(Post1$lambda,Post1$prob,col=4)