### Coverage probabilities of CIs for p

# Compute the Wald CI

alpha = 0.05

# parameters

n = 30
p = 0.5

x = rbinom(1,n,p)

# Compute the Wald CI

p.hat = x/n
p.hat

p.hat.se = sqrt((p.hat*(1-p.hat))/n)
p.hat.se

z = c(qnorm(alpha/2),qnorm(1-alpha/2))
z

p.ci.wald = p.hat + z*p.hat.se
p.ci.wald

# Compute the variance stabalized CI

p.hat.vst = asin(sqrt(p.hat))
p.hat.vst

p.hat.vst.se = sqrt(1/(4*n))
p.hat.vst.se

p.ci.vst = p.hat.vst + z*p.hat.vst.se
p.ci.vst

p.ci.vs = (sin(p.ci.vst))^2
p.ci.vs

# Compute the +2 CI, Score CI

p.hat2 = (x+2)/(n+4)
p.hat2

p.hat2.se = sqrt((p.hat2*(1-p.hat2))/n)
p.hat2.se

p.ci.2 = p.hat2 + z*p.hat2.se
p.ci.2

# Compute the Clopper and Pearson CI

p.ci.cp = c(qbeta( alpha/2, x, n-x+1),qbeta( (1-alpha/2), x+1, n-x))
p.ci.cp

# Compute the coverage probabilities

B = 100000

n = 1000
p = 0.1

alpha = 0.05
z = c(qnorm(alpha/2),qnorm(1-alpha/2))

Count.wald = 0
Count.vs = 0
Count.2 = 0
Count.cp = 0

for (i in 1:B){
	x = rbinom(1,n,p)
	p.hat = x/n
	p.hat.se = sqrt((p.hat*(1-p.hat))/n)
	p.ci.wald = p.hat + z*p.hat.se
	if (p > p.ci.wald[1] && p < p.ci.wald[2]) Count.wald = Count.wald + 1
	p.hat.vst = asin(sqrt(p.hat))
	p.hat.vst.se = sqrt(1/(4*n))
	p.ci.vs = (sin(p.hat.vst + z*p.hat.vst.se))^2
	if (p > p.ci.vs[1] && p < p.ci.vs[2]) Count.vs = Count.vs + 1
	p.hat2 = (x+2)/(n+4)
	p.hat2.se = sqrt((p.hat2*(1-p.hat2))/n)
	p.ci.2 = p.hat2 + z*p.hat2.se
	if (p > p.ci.2[1] && p < p.ci.2[2]) Count.2 = Count.2 + 1
	p.ci.cp = c(qbeta( alpha/2, x, n-x+1),qbeta( (1-alpha/2), x+1, n-x))
	if (p > p.ci.cp[1] && p < p.ci.cp[2]) Count.cp = Count.cp + 1
}

Count.wald/B
Count.vs/B
Count.2/B
Count.cp/B