# R < sim202.s --no-save &
## 25952
#########################################################################
## Find optimal p1, for different alpha (fixed n, q): ##  
#########################################################################


sink("out202b.txt")
cat("Here c = the .02 quantile for the sum of n Paretos:\n\n")
cat("And we're doing 10 million simulations of each sum.\n")
cat("Really 10 million; i.e. when n=5, we're actually \n")
cat("generating 50 million Paretos altogether, 5 for each sum.\n\n\n")

a <- c(.5,2/3,1,1.5)
b <- 1
q <- .02
nsim <- 10000000 ## 10 million
n1 <- c(2,5,10,20,50,100)


### fast simulation!  ONLY NEED n SIMs for this!
for(k in c(1:6)){
n <- n1[k]
par11 <- rep(0,nsim)
par12 <- rep(0,nsim)
par13 <- rep(0,nsim)
par14 <- rep(0,nsim)
par3 <- rep(0,length(a))
for(j in c(1:n)){
par11 <- par11 + runif(nsim)^(-1/a[1])*b
par12 <- par12 + runif(nsim)^(-1/a[2])*b
par13 <- par13 + runif(nsim)^(-1/a[3])*b
par14 <- par14 + runif(nsim)^(-1/a[4])*b
# cat(j)
}
par2 <- sort(par11)
par3[1] <- par2[round(q*nsim)]
par2 <- sort(par12)
par3[2] <- par2[round(q*nsim)]
par2 <- sort(par13)
par3[3] <- par2[round(q*nsim)]
par2 <- sort(par14)
par3[4] <- par2[round(q*nsim)]
for(j in c(1:4)){
cat("n = ",n1[k]," a = ",a[j]," c = ",par3[j],"\n")
}
}

sink()