set.seed(12345)
# 2
# B)
# Find 10% lower tail quartile for critical value k(p) from a 
# binom(20,p)
# Set p0=(.9,.75,.5852,.5851,.25)
n=20
alpha=.1

ki=qbinom(alpha,n,.9,lower.tail=F)
kii=qbinom(alpha,n,.75,lower.tail=F)
kiii=qbinom(alpha,n,.5852,lower.tail=F)
kiv=qbinom(alpha,n,.5851,lower.tail=F)
kv=qbinom(alpha,n,.25,lower.tail=F)

# Simulate densities for population based on binom(20,p0) and add
# lines for upper 10% bound to check
N=1000000

par(mfrow=c(2,1))

Ti=rbinom(N,n,.9)
hist(Ti)
abline(v=ki,col=2)

Tii=rbinom(N,n,.75)
hist(Tii)
abline(v=kii,col=2)

Tiii=rbinom(N,n,.5852)
hist(Tiii)
abline(v=kii,col=2)

Tiv=rbinom(N,n,.5851)
hist(Tiv)
abline(v=kiv,col=2)

Tv=rbinom(N,n,.25)
hist(Tv)
abline(v=kv,col=2)

layout(1,1)

# Upper bounds for critical value
c(ki,kii,kiii,kiv,kv)
## [1] 20 17 15 14  8
# C)
# plot critical value over all possible values of p
p=(1:9999)/10000
kp=qbinom(alpha,n,p,lower.tail=F)
plot(p,kp,ylab='k(p)',main="AR constant over p")
# see section (F) for graphic

# E)
# find the max value of p when k <= 15
max(p[kp<=15])
## [1] 0.6393
# check
plot(p,kp,ylab='k(p)',main="AR constant over p")
abline(v=max(p[kp<=15]),h=15,col=2)

# F)
# invert k to find a lower conf bound for p
# find the max value of p when k strictly < 15
max(p[kp<15])
## [1] 0.5851