１−２．ベイズと分布〜Rを使う ぱらぱらめくる『はじめての統計データ分析』

• 真の分布は理論分布ではない(かもしれない)
  • 真の分布として、平均100、標準偏差10の正規分布と平均120、標準偏差30の正規分布の0.7:0.3混合分布を仮定する
  • そのような真の分布からの500人分の観察データをシミュレート作成する

n <- 500
pr <- c(0.7,0.3)
AorB <- sample(c(1,2),n,replace=TRUE,prob=pr)
x <- rep(0,n)
for(i in 1:n){
	if(AorB[i]==1){
		tmp <- rnorm(1,100,10)
	}else{
		tmp <- rnorm(1,120,30)
	}
	x[i] <- tmp
}
hist(x)

• 事前分布を想定する
  • 正規分布と想定する
  • 正規分布には２つのパラメタがある
    • 平均と標準偏差
    • 平均は(よくわからないけれど)50から150なら、どれもありそうで、特にどの値が、というのはないけれど、50より小、50より大はあり得ないと思う
    • 標準偏差は(よくわからないけれど)5から10が強くありそうで、10から50はまあ、ありそうだけれど、その外側はあり得ないと思うとする

n <- 20
pr <- c(0.7,0.3)
AorB <- sample(c(1,2),n,replace=TRUE,prob=pr)
x <- rep(0,n)
for(i in 1:n){
	if(AorB[i]==1){
		tmp <- rnorm(1,100,10)
	}else{
		tmp <- rnorm(1,120,30)
	}
	x[i] <- tmp
}
hist(x)

mm <- seq(from=0,to=200)
mm.pre <- rep(0,length(mm))
mm.pre[which(mm>50 & mm<150)] <- 1
plot(mm,mm.pre,type="l")

ss <- seq(from=0,to=100)
ss.pre <- rep(0,length(ss))
ss.pre[which(ss>5 & ss<=10)] <- 1
ss.pre[which(ss>10 & ss<50)] <- 0.5
plot(ss,ss.pre,type="l")

• • 面積は1に

mm.pre2 <- mm.pre/(150-100)
ss.pre2 <- ss.pre/(5+20)
plot(mm,mm.pre2,type="l")
plot(ss,ss.pre2,type="l")

• 事前確率
  • (m,sd) = (70,9)の正規分布と、(75,11)の正規分布との事前確率は(95,7)の方が２倍

m1 <- 100; sd1 <- 9
m2 <- 105; sd2 <- 11
preL1 <- mm.pre2[which(mm==m1)] * ss.pre2[which(ss==sd1)]
preL2 <- mm.pre2[which(mm==m2)] * ss.pre2[which(ss==sd2)]
preL1
preL2

• 尤度

L1 <- dnorm(x,m1,sd1)
L2 <- dnorm(x,m2,sd2)
plot(x,L1)
plot(x,L2)
hist(L1)
hist(L2)

• • 同時確率・同時尤度

prod(L1) # ??
prod(L2) # ??

• • 小さすぎる→対数尤度

sum(log(L1))
sum(log(L2))

• • 尤度比

exp(sum(log(L1))-sum(log(L2)))

• 事後確率

logL1.post <- sum(log(L1)) + log(preL1)
logL2.post <- sum(log(L2)) + log(preL2)
exp(logL1.post-logL2.post)