Based on simulations, I've come up with a simple function to compute
the confidence interval for the variance of the binomial variance,
where the true variance is
v = rho*(1-rho)/n
where rho = true probability of success and n = # of trials.
For x = # successes observed in n trials, p = x / n as usual.
For p < 0.25 or p > 0.75, I use the proportion-based transformed
confidence interval, as suggested by Moshe Olshansky. I used the
Wilson interval, but some might prefer the Blaker interval. For p >=
0.25 or p <= 0.75, I use the standard chi-square based confidence
interval (which is very conservative for this problem in this range).
The composite method gives reliable results over a wide range of rho.
This should suit the purpose until someone comes up with a more
theoretically sound method. Bootstrap is not reliable for n < 50, at least.
#09.26.08 02.50 binomVarCI.r
#copyright 2008 by Robert A LaBudde, all rights reserved
#CI for binomial sample variance
#created: 09.26.08 by r.a. labudde
#changes:
require('binGroup')
binomVarCI<- function (n, x, conf=0.95) {
p<- x/n #proportion
if (p<0.25 | p>0.75 | x==0 | x==n) { #use proportion-based CI
pCI<- binWilson(n, x, conf.level=conf) #CI for proportion
vCI<- sort(c(pCI[1]*(1-pCI[1])/(n-1), pCI[2]*(1-pCI[2])/(n-1)))
} else { #use chi-square-based CI
phiL<- qchisq(0.025, n-1)/(n-1)
phiU<- qchisq(0.975, n-1)/(n-1)
#vest<- p*(1-p)/(n-1)) #variance estimate
vCI<- c(vest/phiU, vest/phiL) #chi-square-based
}
return (vCI)
}
Here is a test program to measure coverage:
#09.26.08 03.10 tbinomVarCI.r
#copyright 2008 by Robert A LaBudde, all rights reserved
#test of CI for binomial sample variance
#created: 09.26.08 by r.a. labudde
#changes:
nReal <- 1000
for (POD in c(0.05, 0.10, 0.15, 0.20, 0.25, 0.30, 0.35, 0.50)) {
cat('\nPOD: ', sprintf('%8.4f', POD), '\n')
for (nRepl in c(6, 12, 20, 50)) {
vtrue<- POD*(1-POD)/nRepl
pcover<- 0
for (iReal in 1:nReal) {
x<- rbinom(1, nRepl, POD)
vCI<- binomVarCI(nRepl, x)
#vest<- (x/nRepl)*(1 - x/nRepl)/(nRepl-1)
if (vtrue >= vCI[1] & vtrue<= vCI[2]) pcover<- pcover + 1
} #iReal
pcover<- pcover/nReal
cat('n: ', sprintf('%4i', nRepl), ' Coverage: ',
sprintf('%8.4f', pcover), '\n')
} #nRepl
} #POD
================================================================
Robert A. LaBudde, PhD, PAS, Dpl. ACAFS e-mail: [EMAIL PROTECTED]
Least Cost Formulations, Ltd. URL: http://lcfltd.com/
824 Timberlake Drive Tel: 757-467-0954
Virginia Beach, VA 23464-3239 Fax: 757-467-2947
"Vere scire est per causas scire"
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
