cutoff = .001
p<-.0341

#adopt a jeffreys prior
alpha<-.5
beta<-.5
flag = 0
N_vec = rep(0,100)
for (j in 1:100){
N = 100000
while (flag == 0){
y=rbinom(1,N,p)
lci=qbeta(0.05,shape1=y+alpha,shape2=N-y+beta)
uci=qbeta(0.95,shape1=y+alpha,shape2=N-y+beta)
if (uci - lci < cutoff) {
  flag = 1
  N_vec[j] = N
}
N = N + 1000
}
flag = 0
}
hist(N_vec)

uci=1
lci=0
flag = 0
N_vec2 = rep(0,100)
for (j in 1:100){
  N = 350000
  while (flag == 0){
    y=rbinom(1,N,p)
    lci=qbeta(0.025,shape1=y+alpha,shape2=N-y+beta)
    uci=qbeta(0.975,shape1=y+alpha,shape2=N-y+beta)
    if (uci - lci < cutoff) {
      flag = 1
      N_vec2[j] = N
    }
    N = N + 1000
  }
  flag = 0
}
hist(N_vec2)

uci=1
lci=0
flag = 0
N_vec3 = rep(0,100)
for (j in 1:100){
  N = 500000
  while (flag == 0){
    y=rbinom(1,N,p)
    lci=qbeta(0.005,shape1=y+alpha,shape2=N-y+beta)
    uci=qbeta(0.995,shape1=y+alpha,shape2=N-y+beta)
    if (uci - lci < cutoff) {
      flag = 1
      N_vec3[j] = N
    }
    N = N + 1000
  }
  flag = 0
}
hist(N_vec3)

binomial.conf.interval=function (y,conf){
  z1=abs(qnorm((1 - conf)/2))
  z2=abs(qnorm((1 + conf)/2))
  nn=length(y)
  phat=sum(y)/nn
  se=sqrt(phat∗(1−phat)/nn)
  return (c(phat−z1∗se, phat+z2∗se))
}

flag = 0
N_vec_clt = rep(0,100)
for (j in 1:100){
  N = 100000
  while (flag == 0){
    y=rbinom(N,1,p)
    I = binomial.conf.interval(y,.9)
    if (I[2] - I[1] < cutoff) {
      flag = 1
      N_vec_clt[j] = N
    }
    N = N + 1000
  }
  flag = 0
}
hist(N_vec)

sample_est = c(median(N_vec), median(N_vec2),median(N_vec3))
sample_est