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

library(MASS)

n.sim <- 2000
n.sample <- c(100)#,200,300,500,1000)
source("gestCov_txRand.R")
source("gestCov_txRand_iterate.R")
  
psi10 <- 1.6	
psi11 <- 0	
psi20 <- -3	
psi21 <- 2
psi22 <- 1.1	
psi23 <- 0	

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

se.sim <- NULL
var.mySE <- NULL
mean.mySE <- NULL
median.mySE <- NULL

mean.mySEJR <- NULL
median.mySEJR <- NULL


for(n in n.sample) {
results <- psi
   psi.est <- var.est <- NULL
   wrong.rule <- NULL
   U.est <- var.U <- var.U.hat <- NULL
   psi.estJR <- var.estJR <- NULL
   U.estJR <- var.UJR <- var.U.hatJR <- NULL
   for(j in 1:n.sim) {
if(j %in% seq(25,n.sim,25)) print(j) #print(sum(A1==1 & L2==1))
      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 <- g_estCov.fcn(Y,L,A,return.U=F) 
      psi.est <- cbind(psi.est,c(temp$psi.hat[2,1:2],temp$psi.hat[1,]))   
      var.est <- cbind(var.est,c(temp$psi.var[2,1:2],temp$psi.var[1,])) 
   }

## NB: MEDIAN sd looks like Monte Carlo SD (mean too small), though still a bit small
   print(paste("n = ",n,sep=""))
   t1 <- apply(psi.est,1,mean)
   t2a <- apply(psi.est,1,sd)    # simulated SE
#   t2 <- apply(sqrt(var.est),1,median) # calculated SE
   t2 <- apply(sqrt(var.est),1,mean) # calculated SE
   t3 <- apply(sqrt(var.est + (psi.est - psi)^2),1,mean)
   t4 <- apply((psi.est - 1.96*sqrt(var.est) < psi) & (psi.est + 1.96*sqrt(var.est) > psi),1,sum)/n.sim*100 # Wald CI coverage
#   results <- cbind(results,cbind(t1,t2a,t2,t3,t4))
   results <- cbind(results,cbind(t1,t2,t3,t4))

#print("Wrong rule:")
#print(t(apply(wrong.rule,1,summary)))
#print(t(apply(wrong.rule,1,summary)/n*100))
#se.sim <- cbind(se.sim,apply(psi.est,1,sd))    
#mean.mySE <- cbind(mean.mySE,apply(sqrt(var.est),1,mean)[1:2])
#median.mySE <- cbind(median.mySE,apply(sqrt(var.est),1,median)[1:2])

   print(round(results,3))
}



#postscript("psi10bias_100.eps")
truehist(psi.est[1,],main="",xlab=expression(psi[10]),col=0)
lines(density(psi.est[1,]),lwd=2)
abline(v=psi10,lwd=2)
abline(v=mean(psi.est[1,]),lty=2,lwd=2)
#dev.off()
#postscript("psi11bias_100.eps")
truehist(psi.est[2,],main="",xlab=expression(psi[11]),col=0)
lines(density(psi.est[2,]),lwd=2)
abline(v=psi11,lwd=2)
abline(v=mean(psi.est[2,]),lty=2,lwd=2)
#dev.off()

par(mfrow=c(2,2))
hist(psi.est[3,])
hist(psi.est[4,])
hist(psi.est[5,])
hist(psi.est[6,])
plot(density(psi.est[3,]))
abline(v=psi20)
plot(density(psi.est[4,]))
abline(v=psi21)
plot(density(psi.est[5,]))
abline(v=psi22)
plot(density(psi.est[6,]))
abline(v=psi23)

#postscript("SEpsi10.eps")
plot(n.sample,mean.mySE[1,],type="l",ylab="SE(psi10)",xlab="n",col=6)
lines(n.sample,median.mySE[1,],lty=2,col=6)
lines(n.sample,se.sim[1,])
lines(n.sample,mean.mySEJR[1,],col=4)
lines(n.sample,median.mySEJR[1,],lty=2,col=4)
#dev.off()
#postscript("SEpsi11.eps")
plot(n.sample,mean.mySE[2,],type="l",ylab="SE(psi11)",xlab="n",col=6)
lines(n.sample,median.mySE[2,],lty=2,col=6)
lines(n.sample,se.sim[2,])
lines(n.sample,mean.mySEJR[2,],col=4)
lines(n.sample,median.mySEJR[2,],lty=2,col=4)
#dev.off()

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))
}

d.opt.k(psi.est[3:6],L2,A1)
d.opt.k(psi[3:6],L2,A1)
table(d.opt.k(psi.est[3:6],L2,A1),d.opt.k(psi[3:6],L2,A1))