1−2.ベイズと分布〜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)
  • 事前分布を想定する
    • 正規分布と想定する
    • 正規分布には2つのパラメタがある
      • 平均と標準偏差
      • 平均は(よくわからないけれど)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")
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)