与えられた条件から、推定する〜信頼区間？〜

数学いらずの医科統計学 第2版

• 関連記事など(こちら)
• 条件とデータが与えられたときに、知りたい変数の値を推定したい
• 最尤推定
  • 最も尤もらしい値を推定する(最尤推定値)
  • 推定したい変数の尤度関数
  • 観察データから尤度関数の値を求める
  • 最尤推定値は、尤度関数の最大値を与える変数値。微分できる？→微分する
• 尤もらしい範囲を推定する
  • 尤度に基づくありそうな値とありそうもない値とがわかる
  • 尤もらしい値を範囲で選ぶとしたらどうしたいか
  • Rの信頼区間(こちら)
• メモ

p<- seq(0,1,0.00001)
pp<- rep(0,length(p))
for(i in 1:length(p)){
	pp[i]<- dbinom(4,15,p[i])
}
plot(pp)


#plot(pp/sum(pp),cex=0.1,axes=FALSE,ylim=c(0,0.0005),
#	 xlab="Ratio of Red balls",ylab="Likelyhood")
#axis(1,
#	 c(0,which(pp==max(pp)),length(p)),
#	 c(0,round(which(pp==max(pp))/length(p),3),1))
#axis(2)

a<- pp/sum(pp)
a_order<- order(a,decreasing=TRUE)

a_ordered<-a[a_order]
cumsum_a_ordered<-cumsum(a_ordered)
plot(cumsum_a_ordered)

closest<-which(abs(cumsum_a_ordered-0.95)==min(abs(cumsum_a_ordered-0.95)))
abline(v=closest)
abline(h=0.95)

plot(pp,type="l")
#abline(h=a[a_order[closest]])

thres<-a[a_order[closest]]
mores<-which(a>thres)
lessANDeq<-which(a<=thres)

sum(a[mores])
sum(a[lessANDeq])

CI<-qbeta(c(1-0.95)/2,1-(1-cirange)/2),4+1,15-4+1)
abline(v=CI*length(p),col=2)
abline(v=a_order[closest])
abline(v=a_order[closest-1])


N<-15
k<-4
p<-seq(from=0,to=1,by=0.001)
v<-dbeta(p,k+1,N-k+1)
v<-v/sum(v)
plot(p,v,type="l")
abline(v=k/N)
cirange<-0.95
ci1<-qbeta(c((1-cirange)/2,1-(1-cirange)/2),k+1,N-k+1)
abline(v=ci,col=2)

lesses<-which(p<ci1[1])
mores<-which(p>ci1[2])


a_order<- order(v,decreasing=TRUE)

a_ordered<-v[a_order]
cumsum_a_ordered<-cumsum(a_ordered)

closest<-which(abs(cumsum_a_ordered-0.95)==min(abs(cumsum_a_ordered-0.95)))
thres<-v[a_order[closest]]
mores2<-which(v>thres)
lesss2<-which(v<thres)
col<-rep(0,length(p))
col[mores2]<-2
col[lesss2]<-3
par(new=TRUE)
plot(p,v,col=col,cex=0.1)
abline(h=thres,col=3)
sum(v[c(lesses,mores)])
sum(v[-c(lesses,mores)])
sum(v[c(lesss2)])
sum(v[c(mores2)])

• dbinom()はp引数がベクトル対応なのでループがはしょれるらしい(こちら)
  • for()は、好きかもしれない