probability of coverage

[cro77]

hi,

yesterday i started to mess around with this book : Bayesian computation with R and i stopped at the first chapter (on exercises) . in this exercise i have to estimate the true probability of coverage (P(p ele. C(y)) where C is the confidence). so what i have to do is to simulate for sample size n=20 and p = 0.5 using the rbinom function, the value y. then write a function to calculate binomial confidence interval

bin.conf = function (y,n){
  z = qnorm(0.95);
  phat = y/n;
  se = sqrt(phat*(1-phat)/n);
  return(c(phat-z*se,phat+z*se))
}

and then repeat this procedure 20 times to estimate the true probability of coverage.

so the problem is how do i do that ?

first i iterate :

for (i in 1:20){
  y = rbinom(1,0.5,20);
  bin = bin.conf(y,20);
  # and now what how do i estimate the true  probability of coverage
}

ok, i just realized this looks like homework but i assure you it is not (bit to old for school )

so if anyone has some spare time to help me i would appreciate it ! thank you

baxy

 

[Masteras]

bin=matrix(ncol=2,nrow=1000)
for (i in 1:1000){
y = rbinom(1,20,0.5)
bin[i,] = bin.conf(y,20)
}

the bin matrix contains the lower 90% conf limmit (first column) and the upper 90% conf limmit (2nd column)

I say 90% confidence interval because z=qnorm(0.95) in the previous function. If you want a 95% conf interval you must use z=qnorm(0.975).

To estimate the coverage now if it is indeed 90% as you want you must see how many rows (confidence intervals) contain 0.5 form the 1000 rows. Or how many do not include it and subtract them from 1000 and then divide by 1000.

sum(bin[,1]>0.5)
intervals who are above 0.5
sum(bin[,2]<0.5)
intervals who are below 0.5.
the sum of these two will be close to 100.
Thus, the coverage is 1-100/1000=0.9.

 

[Masteras]

By the way, you were very close. too old for school but still young to rock