#				
## Simulated data to test variance  - random treatment, 2 intervals (comparing at my covariance with (3.14))
#				

rm(list=ls())
library(MASS)

n.sim <- 1000
n.sample <- c(250,500,1000)#seq(250,1000,50)#
n.boot <- 500
source("AAOgest_txS.R")
source("../gestCov_txRand.R")
source("gestCov_txRand_ZIPI.R")
 
psi20 <- -3	
psi21 <- 2
psi22 <- 1	
psi23 <- 0	
psi10 <- psi20	
psi11 <- psi21

psi<-c(psi20,psi21,psi22,psi23)
results <- psi
alpha<-c(.5,.5)
get.zipi <- 1
#bias.mat <- NULL

#for(psi21 in seq(1,3,.1)) {
#bias.vec <- NULL
for(n in n.sample) {
psi<-c(psi20,psi21,psi22,psi23)
results <- psi
   psi.est80 <- var.est80 <- NULL
   psi.est90 <- var.est90 <- NULL
   psi.est95 <- var.est95 <- NULL
   psi.est99 <- var.est99 <- NULL
   psi.estJR <- var.estJR <- NULL
   psi.estAAO <- var.estAAO <- NULL
   for(j in 1:n.sim) {
if(j %in% seq(100,n.sim,100)) print(j)
      L1 <- rnorm(n,0,3)
      A1 <- rbinom(n,1,alpha[1])  
      L2 <- round(rnorm(n,.75,.5))
      A2 <- rbinom(n,1,alpha[2])  
      L <- cbind(L1,L2)
      A <- cbind(A1,A2)
      Y.opt <- rnorm(n,1,1) 
      Y <- Y.opt - abs(psi10 + psi11*L1)*(A1 - as.numeric(psi10 + psi11*L1>0))^2 - abs(psi20 + psi21*L2 + psi22*A1 + psi23*L2*A1)*(A2 - as.numeric(psi20 + psi21*L2 + psi22*A1 + psi23*L2*A1>0))^2      

      if(get.zipi) {
       temp <- ZIPI_rand.fcn(Y,L,A,return.U=F,z.val=qnorm(0.9)) # 80% CI zeroing
       icw.est <- (temp$psi.hat[2,1:2]/temp$psi.var[2,1:2] + temp$psi.hat[1,1:2]/temp$psi.var[1,1:2])/(1/temp$psi.var[2,1:2] + 1/temp$psi.var[1,1:2]) 
       psi.est80 <- cbind(psi.est80,c(icw.est,temp$psi.hat[1,3:4]))   
       temp <- ZIPI_rand.fcn(Y,L,A,return.U=F,z.val=qnorm(0.975)) # 95% CI zeroing
       icw.est <- (temp$psi.hat[2,1:2]/temp$psi.var[2,1:2] + temp$psi.hat[1,1:2]/temp$psi.var[1,1:2])/(1/temp$psi.var[2,1:2] + 1/temp$psi.var[1,1:2])
       psi.est95 <- cbind(psi.est95,c(icw.est,temp$psi.hat[1,3:4]))   
      }
#      temp <- ZIPI_rand.fcn(Y,L,A,return.U=F,z.val=qnorm(0.9)) # 80%CI
#      temp <- ZIPI_rand.fcn(Y,L,A,return.U=F,z.val=qnorm(0.975)) # 95%CI
      temp <- g_estCov.fcn(Y,L,A,return.U=F)
      icw.est <- (temp$psi.hat[2,1:2]/temp$psi.var[2,1:2] + temp$psi.hat[1,1:2]/temp$psi.var[1,1:2])/(1/temp$psi.var[2,1:2] + 1/temp$psi.var[1,1:2])
      psi.estJR <- cbind(psi.estJR,c(icw.est,temp$psi.hat[1,3:4]))

      start.search <- c(icw.est,temp$psi.hat[1,3:4])
      temp <- nlm(AAOg_estSeq.fcn,p=start.search,Y=Y,action=A,status=L)
      psi.estAAO <- cbind(psi.estAAO,temp$estimate)

## Bootstrap the variance!
      boot.psi <- boot.psiAAO <- boot.psi80 <- boot.psi95 <- NULL
      for(i in 1:n.boot) {
       boot.samp <- sample(1:n,replace=T)
       Y.boot <- Y[boot.samp]
       L.boot <- L[boot.samp,]
       A.boot <- A[boot.samp,]
       temp <- g_estCov.fcn(Y.boot,L.boot,A.boot,return.U=F)
       icw.est <- (temp$psi.hat[2,1:2]/temp$psi.var[2,1:2] + temp$psi.hat[1,1:2]/temp$psi.var[1,1:2])/(1/temp$psi.var[2,1:2] + 1/temp$psi.var[1,1:2])
       boot.psi <- cbind(boot.psi,c(icw.est,temp$psi.hat[1,3:4]))

       start.search <- c(icw.est,temp$psi.hat[1,3:4])
       temp <- nlm(AAOg_estSeq.fcn,p=start.search,Y=Y.boot,action=A.boot,status=L.boot)
       boot.psiAAO <- cbind(boot.psiAAO,temp$estimate)

       temp <- ZIPI_rand.fcn(Y.boot,L.boot,A.boot,return.U=F,z.val=qnorm(0.9)) # 80% CI zeroing
       icw.est <- (temp$psi.hat[2,1:2]/temp$psi.var[2,1:2] + temp$psi.hat[1,1:2]/temp$psi.var[1,1:2])/(1/temp$psi.var[2,1:2] + 1/temp$psi.var[1,1:2]) 
       boot.psi80 <- cbind(boot.psi80,c(icw.est,temp$psi.hat[1,3:4]))

       temp <- ZIPI_rand.fcn(Y.boot,L.boot,A.boot,return.U=F,z.val=qnorm(0.975)) # 95% CI zeroing
       icw.est <- (temp$psi.hat[2,1:2]/temp$psi.var[2,1:2] + temp$psi.hat[1,1:2]/temp$psi.var[1,1:2])/(1/temp$psi.var[2,1:2] + 1/temp$psi.var[1,1:2])
       boot.psi95 <- cbind(boot.psi95,c(icw.est,temp$psi.hat[1,3:4]))
      }
      var.estJR <- cbind(var.estJR,apply(boot.psi,1,var)) 
      var.estAAO <- cbind(var.estAAO,apply(boot.psiAAO,1,var)) 
      var.est80 <- cbind(var.est80,apply(boot.psi80,1,var)) 
      var.est95 <- cbind(var.est95,apply(boot.psi95,1,var)) 
   }

   print(paste("n = ",n,sep=""))
   if(get.zipi) {
    t1.80 <- apply(psi.est80,1,mean)
    t2.80 <- apply(sqrt(var.est80),1,mean) # calculated SE
    t3.80 <- apply(sqrt(var.est80 + (psi.est80 - psi)^2),1,mean)  # rMSE
    t4.80 <- apply((psi.est80 - 1.96*sqrt(var.est80) < psi) & (psi.est80 + 1.96*sqrt(var.est80) > psi),1,sum)/n.sim*100 # Wald CI coverage
    t1.95 <- apply(psi.est95,1,mean)
    t2.95 <- apply(sqrt(var.est95),1,mean) # calculated SE
    t3.95 <- apply(sqrt(var.est95 + (psi.est95 - psi)^2),1,mean)  # rMSE
    t4.95 <- apply((psi.est95 - 1.96*sqrt(var.est95) < psi) & (psi.est95 + 1.96*sqrt(var.est95) > psi),1,sum)/n.sim*100 # Wald CI coverage
    res.zipi <- cbind(t1.95,t2.95,t3.95,t4.95)
   }

   t1g <- apply(psi.estJR,1,mean)
   t2g <- apply(sqrt(var.estJR),1,mean) # calculated SE
   t3g <- apply(sqrt(var.estJR + (psi.estJR - psi)^2),1,mean)  # rMSE
   t4g <- apply((psi.estJR - 1.96*sqrt(var.estJR) < psi) & (psi.estJR + 1.96*sqrt(var.estJR) > psi),1,sum)/n.sim*100 # Wald CI coverage

   t1a <- apply(psi.estAAO,1,mean)
   t2a <- apply(sqrt(var.estAAO),1,mean) # calculated SE
   t3a <- apply(sqrt(var.estAAO + (psi.estAAO - psi)^2),1,mean)  # rMSE
   t4a <- apply((psi.estAAO - 1.96*sqrt(var.estAAO) < psi) & (psi.estAAO + 1.96*sqrt(var.estAAO) > psi),1,sum)/n.sim*100 # Wald CI coverage

   if(get.zipi) {
    results <- cbind(results,cbind(t1.80,t2.80,t3.80,t4.80,t1g,t2g,t3g,t4g,t1a,t2a,t3a,t4a))
    print(round(results,3))
    print(round(res.zipi,3))
   }

  if(!get.zipi) {
    results <- cbind(results,cbind(t1g,t2g,t3g,t4g))
    print(round(results,3))
    bias.vec <- c(bias.vec,t1g[1] - psi10)
  }
}
#bias.mat <- cbind(bias.mat,bias.vec)
#}
#bias.mat <- abs(bias.mat)

#postscript("bias_by_psi22VESSa.eps")
#filled.contour(seq(250,1000,50),seq(1,3,.1),bias.mat,xlab="n",ylab=expression(psi[21]),col=rainbow(35),zlim=c(-.07,.07))
#dev.off()
#postscript("bias_by_psi21VESSb.eps")
#filled.contour(seq(250,1000,50),seq(1,3,.1),bias.mat,xlab="n",ylab=expression(psi[21]),col=rainbow(35),levels=seq(0,.07,.0025))
#dev.off()
#filled.contour(seq(250,1000,50),seq(1,3,.1),bias.mat,xlab="n",ylab=expression(psi[21]),col=rainbow(35),levels=seq(0,.07,.0025))

#filled.contour(seq(250,1000,50),seq(1,3,.1),bias.mat,xlab="n",ylab=expression(psi[21]),col=rainbow(35),zlim=c(-.07,.07))
#filled.contour(seq(250,1000,50),seq(1,3,.1),bias.mat,xlab="n",ylab=expression(psi[21]),col=rainbow(35),levels=seq(-.07,.07,.005))

#round(cbind(t1p,t2p,t3p,t4p,t1g,t2g,t3g,t4g),3)