作者qqpolt (忘記)
看板Statistics
標題[問題] 可以幫我看一下程式嗎(R)?關於coverage probabilty的
時間Wed Nov 12 22:00:25 2008
我是看一篇paper的數據寫的
信賴區間算出來是跟paper上的一樣
但是再算它的coverage prob的時候卻有落差
可以幫我看看嗎?
我是用R語言寫的
以下是算信賴區間的程式
N=467
ai=c(56,72,73,59,62,87,58)
X=7.235
c3=matrix(0,1,7)
c4=matrix(0,1,7)
c5=matrix(0,1,7)
c6=matrix(0,1,7)
c7=matrix(0,1,7)
c8=matrix(0,1,7)
c9=matrix(0,1,7)
c10=matrix(0,1,7)
c11=matrix(0,1,7)
for(s in 1:7)
{
c3[s]=4*ai[s]*(N-ai[s])
c4[s]=c3[s]/N
c5[s]=X*(X+c4[s])
c6[s]=sqrt(c5[s])
c7[s]=X+2*ai[s]
c8[s]=c7[s]-c6[s]
c9[s]=c7[s]+c6[s]
c10[s]=c8[s]/(2*(N+X))
c11[s]=c9[s]/(2*(N+X))
}
c10
c11
以下是我算coverage prob的程式(我是跑10000次)
N=467
ai=c(56,72,73,59,62,87,58)
X=7.235
c3=matrix(0,1,7)
c4=matrix(0,1,7)
c5=matrix(0,1,7)
c6=matrix(0,1,7)
c7=matrix(0,1,7)
c8=matrix(0,1,7)
c9=matrix(0,1,7)
c10=matrix(0,1,7)
c11=matrix(0,1,7)
p=ai/N
ccc=rep(0,10000)
for(i in 1:10000){
w=rmultinom(1,N,p)
for(j in 1:7)
{ c3[j]=4*w[j]*(N-w[j])
c4[j]=c3[j]/N
c5[j]=X*(X+c4[j])
c6[j]=sqrt(c5[j])
c7[j]=X+2*w[j]
c8[j]=c7[j]-c6[j]
c9[j]=c7[j]+c6[j]
c10[j]=c8[j]/(2*(N+X))
c11[j]=c9[j]/(2*(N+X))
}
aa=ifelse(c10[1]<=p[1]&p[1]<=c11[1]&c10[2]<=p[2]&p[2]<=c11[2]&c10[3]<=p[3]&p[3]<=c11[3]&c10[4]<=p[4]&p[4]<=c11[4]&c10[5]<=p[5]&p[5]<=c11[5]&c10[6]<=p[6]&p[6]<=c11[6]&c10[7]<=p[7]&p[7]<=c11[7],0,1)
bb=sum(aa)
ccc[i]=ifelse(bb==0,1,0)
}
sum(ccc)
跑出的sum(ccc)都是9550左右
但實際上paper寫的是9350左右~~有點落差
我已經完全不知道哪裡有問題了~想了一下午
希望高手救救我
--
※ 發信站: 批踢踢實業坊(ptt.cc)
◆ From: 60.249.21.60