二項

• 信頼区間というのがある
• 母分布からデータセットをとって、母分布の平均が知りたいときに、昨日は点推定をしたわけだが、今日は区間推定
• 確率pで生起するときN回観測して、k回起きたという観測結果を基にpを区間推定する
• Rのbinomパッケージではそのヘルプ文書で『Fisher Information for Various Binomial Parameterizations』と題して、複数の信頼区間算出関数を書いている
• こちらでは、それを含めて、『信頼区間』の保守性、とかをコメントしている
• 実際、どうなの？ということで、pを振って、母分布からN=20でサンプリングして、その観測結果に基づき、それぞれの方法で信頼区間を出して、その区間内にpがどのくらいの割合で入るか(95%、入るとよいが…という方法のはず)をぐるぐるをループを回してやってみる
• 記事などで指摘がある通り、pが0(1でも…)に近い付近での歪みが大きいわけだが、どの方法がどういう風に、どれくらい信頼区間に影響を与えているかを図にする。

library(binom)
ps <- seq(from = 0,to=0.5,by=0.01)
ps <- ps[-1]
ps <- ps[-length(ps)]

n.iter <- 1000
N <- 20
results <- matrix(0,length(ps),11)
for(i in 1:length(ps)){
	res <- matrix(0,n.iter,11)
	for(j in 1:n.iter){
		X <- sample(0:1,N,replace=TRUE,prob=c(1-ps[i],ps[i]))
		binom.out <- binom.confint(length(which(X==1)),N,conf.int=0.95,prior.shape1=1,prior.shape2=1)
		ci <- cbind(binom.out$lower,binom.out$upper)
		res[j,] <- (ci[,1]<ps[i])*(ci[,2]>ps[i])
	}
	results[i,] <- apply(res,2,sum)
}

matplot(ps,results,type="l")
abline(h=950)
legend(0.4,600,binom.out[[1]],lty=1:11,col=1:11)