#				
## 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 <- seq(250,1000,50)
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)
bias0.matJR <- bias0.matAAO <- bias0.mat80 <- bias0.mat95 <- NULL
bias1.matJR <- bias1.matAAO <- bias1.mat80 <- bias1.mat95 <- NULL

for(psi22 in seq(0,2,.1)) {
psi10 <- psi20	
psi11 <- psi21
bias0.vecJR <- bias1.vecJR <- NULL
bias0.vecAAO <- bias1.vecAAO <- NULL
bias0.vec80 <- bias1.vec80 <- NULL
bias0.vec95 <- bias1.vec95 <- NULL
for(n in n.sample) {
psi<-c(psi20,psi21,psi22,psi23)
results <- psi
   psi.est80 <- NULL
   psi.est95 <- NULL
   psi.estJR <- NULL
   psi.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      

      temp <- ZIPI_rand.fcn(Y,L,A,return.U=F,z.val=qnorm(0.9)) # 80%CI
      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
#      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 <- 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)
   }

   t1.80 <- apply(psi.est80,1,mean)
#   t1.95 <- apply(psi.est95,1,mean)
   t1g <- apply(psi.estJR,1,mean)
   t1a <- apply(psi.estAAO,1,mean)

   bias0.vecJR <- c(bias0.vecJR,t1g[1] - psi20)
   bias0.vecAAO <- c(bias0.vecAAO,t1a[1] - psi20)
   bias0.vec80 <- c(bias0.vec80,t1.80[1] - psi20)
#   bias0.vec95 <- c(bias0.vec95,t1.95[1] - psi20)
#   bias1.vecJR <- c(bias1.vecJR,t1g[2] - psi21)
#   bias1.vecAAO <- c(bias1.vecAAO,t1a[2] - psi21)
#   bias1.vec80 <- c(bias1.vec80,t1.80[2] - psi21)
#   bias1.vec95 <- c(bias1.vec95,t1.95[2] - psi21)

   print(paste("n = ",n,sep=""))
   print(paste("psi = ",psi20,sep=""))

}
bias0.matJR <- cbind(bias0.matJR,abs(bias0.vecJR))
#bias1.matJR <- cbind(bias1.matJR,abs(bias1.vecJR))
bias0.matAAO <- cbind(bias0.matAAO,abs(bias0.vecAAO))
#bias1.matAAO <- cbind(bias1.matAAO,abs(bias1.vecAAO))
bias0.mat80 <- cbind(bias0.mat80,abs(bias0.vec80))
#bias1.mat80 <- cbind(bias1.mat80,abs(bias1.vec80))
#bias0.mat95 <- cbind(bias0.mat95,abs(bias0.vec95))
#bias1.mat95 <- cbind(bias1.mat95,abs(bias1.vec95))
}

postscript("bias0_psi22r.eps")
filled.contour(seq(250,1000,50),seq(0,2,.1),bias0.matJR,xlab="n",
    ylab=expression(psi[22]),col=rainbow(50),levels=seq(0,.18,.005))
dev.off()
postscript("bias0_psi22s.eps")
filled.contour(seq(250,1000,50),seq(0,2,.1),bias0.matAAO,xlab="n",
    ylab=expression(psi[22]),col=rainbow(50),levels=seq(0,.18,.005))
dev.off()
postscript("bias0_psi22ZIPI80.eps")
filled.contour(seq(250,1000,50),seq(0,2,.1),bias0.mat80,xlab="n",
    ylab=expression(psi[22]),col=rainbow(50),levels=seq(0,.18,.005))
dev.off()
postscript("bias0_psi22ZIPI95.eps")
filled.contour(seq(250,1000,50),seq(0,2,.1),bias0.mat95,xlab="n",
    ylab=expression(psi[22]),col=rainbow(50),levels=seq(0,.18,.005))
dev.off()

#postscript("bias1_psi20r.eps")
#filled.contour(seq(250,1000,50),seq(-4,-2,.1),bias1.matJR,xlab="n",
#    ylab=expression(psi[20]),col=rainbow(35),levels=seq(0,.07,.0025))
#dev.off()
#postscript("bias1_psi20s.eps")
#filled.contour(seq(250,1000,50),seq(-4,-2,.1),bias1.matAAO,xlab="n",
#    ylab=expression(psi[20]),col=rainbow(35),levels=seq(0,.07,.0025))
#dev.off()
#postscript("bias1_psi20ZIPI80.eps")
#filled.contour(seq(250,1000,50),seq(-4,-2,.1),bias1.mat80,xlab="n",
#    ylab=expression(psi[20]),col=rainbow(35),levels=seq(0,.07,.0025))
#dev.off()
#postscript("bias1_psi20ZIPI95.eps")
#filled.contour(seq(250,1000,50),seq(-4,-2,.1),bias1.mat95,xlab="n",
#    ylab=expression(psi[20]),col=rainbow(35),levels=seq(0,.07,.0025))
#dev.off()