######## Now, do the approximation and sim for q=.02, a=2/3, various n. 
## R < approx204.s --no-save

sink("outapprox204.txt")

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

a <- 2/3
b <- 1
q <- .02
n1 <- c(1,3,4)
nsim <- 10000000



for(k in c(1:length(n1))){
n <- n1[k]
par11 <- rep(0,nsim)
for(j in c(1:n)){
par11 <- par11 + runif(nsim)^(-1/a)*b
}
par2 <- sort(par11)
c1 <- par2[round(q*nsim)]
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],c1,chat,ctild,"\n")
}

sink()


## NOTE: n = 2,5,10,20,50,100 are done!

#p1 = .136 + .235q + q^2 + .0066min(n,10) - .05max(a,1),