######## Now, do the approximation
## R < approx201.s --no-save

sink("outapprox201.txt")
cat("Here c-hat = approx. 0.02 quantile using p*=sqrt(.02).\n\n")
### cat("And c-tilde = approx. using p*=.136+.235q+q^2+.0066min(n,10)-.05max(a,1).\n\n")

cat("n a c-hat c-tilde\n\n")

a1 <- c(.5,2/3,1,1.5)
b <- 1
q <- .02
n1 <- c(2,5,10,20,50,100)




for(k in c(1:6)){
n <- n1[k]
for(j in c(1:4)){
a <- a1[j]
p <- sqrt(q)
p2 <- q/p
y <- b/(1-p2^(1/n))^(1/a)
if(a != 1) muy <- (a*b^a*y^(1-a)-a*b)/(1-a)/(1-b^a*y^(-a))
if(a == 1) muy <- log(y) / (1-1/y)
if(a != 2) sy <- sqrt((a*b^a*y^(2-a)-a*b^2)/(2-a)/(1-b^a*y^(-a))-muy^2)
if(a == 2) sy <- sqrt(2*log(y) / (1-y^(-2))-muy^2)
chat <- sy * qnorm(p)*sqrt(n) + n * muy
p <- .136 + .235*q + q*q + .0066*min(n,10) - .05*max(a,1)
p2 <- q/p
y <- b/(1-p2^(1/n))^(1/a)
if(a != 1) muy <- (a*b^a*y^(1-a)-a*b)/(1-a)/(1-b^a*y^(-a))
if(a == 1) muy <- log(y) / (1-1/y)
if(a != 2) sy <- sqrt((a*b^a*y^(2-a)-a*b^2)/(2-a)/(1-b^a*y^(-a))-muy^2)
if(a == 2) sy <- sqrt(2*log(y) / (1-y^(-2))-muy^2)
ctild <- sy * qnorm(p)*sqrt(n) + n * muy
cat(n1[k],a1[j],chat,ctild,"\n")
}
}

sink()




#p1 = .136 + .235q + q^2 + .0066min(n,10) - .05max(a,1),
#which seemed to be a good approximation to the optimal value of p_1 for
#a=.66 and values of q between 0 and .1 or so. However, from these results
#we now see that this value of p1 doesn't do so well for p_1 greater than
#.1, and even for other values of p1 it does poorly for certain choices of
#a and n.  So perhaps we should experiment with how to improve it.  On the
#other hand, for most cases it works pretty well. 



%ERROR WITH sqrt(p):
[1]  0.0157296288  0.0119533683  0.0079062582  0.0052009496  0.0040290792
 [6]  0.0034694769  0.0018675649  0.0012404466 -0.0049784147 -0.0019160377
[11] -0.0010041897  0.0005489468 -0.0110387866 -0.0056257151 -0.0006206704
[16]  0.0019024907 -0.0148733765 -0.0075615524  0.0003039003  0.0045187175
[21] -0.0182146010 -0.0085055324  0.0017399035  0.0058682144

%ERROR WITH .136 + .235q + q^2 + .0066min(n,10) - .05max(a,1),
[1] -2.325525e-03  2.861704e-03  4.145899e-03 -1.754425e-03 -7.999977e-03
 [6] -2.121836e-03 -4.505549e-04 -2.954831e-03  4.538090e-03  2.733287e-03
[11]  1.283003e-03 -3.926102e-04  7.089704e-04  2.567154e-05  2.275015e-03
[16]  7.113661e-04 -1.526533e-03 -1.376031e-03  3.478128e-03  3.224185e-03
[21] -4.312607e-03 -2.248740e-03  4.869064e-03  4.619396e-03