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

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

n.sim <- 10
n.sample <- 100#c(100,250,1000)
source("../gestCov_txRand.R")
source("../covJR_txRand2.R")
source("../gestCov_txR_ineff.R")
source("../gestCov_rand_iter.R")
  
psi10 <- 1.6	
psi11 <- 0	
psi20 <- -3	
psi21 <- 2
psi22 <- 1	
psi23 <- 0	

psi<-c(psi10,psi11,psi20,psi21,psi22,psi23)
results <- psi
alpha<-c(.5,.5)

d.opt.k <- function(psi.hat.k,stat1,tx0=0) {
   if(length(as.vector(psi.hat.k))==2) return(as.numeric(cbind(1,stat1) %*% psi.hat.k > 0))
   return(as.numeric(cbind(1,stat1,tx0,stat1*tx0) %*% psi.hat.k > 0))
}

for(n in n.sample) {
psi<-c(psi10,psi11,psi20,psi21,psi22,psi23)
results <- psi
   psi.est90 <- var.est90 <- NULL
   psi.est95 <- var.est95 <- NULL
   psi.est99 <- var.est99 <- NULL
   psi.estPUR <- var.estPUR <- NULL
   psi.estJR <- var.estJR <- NULL
   p.op <- NULL
   cover.score <- rep(0,6)
   cover.scoreA <- cover.scoreAA <- rep(0,4)

   U.true <- NULL
   U.est <- NULL
   var.true <- NULL
   var.est <- NULL
   var.psi.hat <- NULL

   for(j in 1:n.sim) {
#if(j %in% seq(200,n.sim,200)) print(j)
      L1 <- rnorm(n,0,3)
      A1 <- rbinom(n,1,alpha[1])  
      L2 <- rnorm(n,.75,.5)
#      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      

#jr_2ints.fcn(Y,L,A,return.U=T,true.psi=psi,true.a=c(.5,.5))
#      temp <- g.estimate.Kint(Y,L,A,return.U=T,true.psi=psi,true.a=c(.5,.5)) 
      temp <- g_estCov.fcn(Y,L,A,return.U=T,true.psi=psi,true.a=c(.5,.5))
      starter <- psi#c(temp$psi.hat[2,1:2],temp$psi.hat[1,])
      temp <- g_est_iter.fcn(Y,L,A,start.vals=starter) 
      psi.estJR <- cbind(psi.estJR,temp)   

#      temp <- g_estCov.fcn(Y,L,A,return.U=T,true.psi=psi,true.a=c(.5,.5))
#      psi.estJR <- cbind(psi.estJR,c(temp$psi.hat[2,1:2],temp$psi.hat[1,]))   
#      var.estJR <- cbind(var.estJR,c(temp$psi.var[2,1:2],temp$psi.var[1,])) 
#      cover.score <- cover.score + ((c(temp$U.hat[2,1:2],temp$U.hat[1,]) - 1.96*sqrt(c(temp$U.var.hat[2,1:2],temp$U.var.hat[1,])) < 0) & 
#                   (c(temp$U.hat[2,1:2],temp$U.hat[1,]) + 1.96*sqrt(c(temp$U.var.hat[2,1:2],temp$U.var.hat[1,])) > 0))
#      cover.scoreA <- cover.scoreA + ((temp$U.hat[1,] - 1.96*sqrt(temp$var.U.true) < 0) & 
#                   (temp$U.hat[1,] + 1.96*sqrt(temp$var.U.true) > 0))
#      cover.scoreAA <- cover.scoreAA + ((temp$U.true - 1.96*sqrt(temp$var.U.true) < 0) & 
#                   (temp$U.true + 1.96*sqrt(temp$var.U.true) > 0))

#      U.true <- cbind(U.true,temp$U.true)
#      U.est <- cbind(U.est,temp$U.hat[1,])
#      var.true <- cbind(var.true,temp$var.U.true)
#      var.est <- cbind(var.est,temp$U.var.hat[1,])
#      var.psi.hat <- cbind(var.psi.hat,temp$var.U.hat.notrue)

#      endpts <- cbind((psi.estJR[,j] - 1.96*sqrt(var.estJR[,j])),(psi.estJR[,j] + 1.96*sqrt(var.estJR[,j]))) 
#	temp <- apply(cbind(d.opt.k(endpts[3:6,1],L2,A1),d.opt.k(psi.estJR[3:6,j],L2,A1),
 #        d.opt.k(endpts[3:6,2],L2,A1)),1,mean)
 #     p.op <- cbind(p.op,sum(!(temp==0 | temp ==1))/n*100)

   }

   print(paste("n = ",n,sep=""))
   print(summary(as.vector(p.op)))
   t0g <- apply(psi.estJR,1,sd)
   t1g <- apply(psi.estJR,1,mean)
#   t1g.bias <- abs(apply(psi.estJR-psi,1,mean))
#   t1g.t <- abs(apply((psi.estJR-psi)/sqrt(var.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
#   t4g <- apply((psi.estJR - qt(.975,100)*sqrt(var.estJR) < psi) & (psi.estJR + qt(.975,100)*sqrt(var.estJR) > psi),1,sum)/n.sim*100 # Wald CI coverage

#    results <- cbind(results,cbind(t1g,t2g,t3g,t4g))
#    results <- cbind(results,cbind(t1g,t1g.bias,t1g.t,t0g,t2g,t3g,t4g))
    results <- cbind(results,cbind(t1g,t0g))
    print(round(results,3))
#    print(round(cover.score/n.sim*100,3))
 #   print(round(cover.scoreA/n.sim*100,3))
  #  print(round(cover.scoreAA/n.sim*100,3))
   # print(round(cbind(apply(U.true,1,mean),apply(U.est,1,mean),apply(var.true,1,mean),
    #                      apply(var.est,1,mean),apply(var.psi.hat,1,mean)),5))
#    print(round(cbind(sqrt(apply(var.true,1,mean)),
 #                         sqrt(apply(var.est,1,mean)),sqrt(apply(var.psi.hat,1,mean))),5))
  #  print((sqrt(apply(var.true,1,mean))-sqrt(apply(var.psi.hat,1,mean)))/sqrt(apply(var.true,1,mean))*100)
}