#				
## 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(20,50,100,250,500,1000)
source("../gestCov_txRand.R")

results <- NULL
n.vess <- NULL

for(n in n.sample) {
   mean.eta <- NULL
   mle <- NULL
   vess <- NULL
   pur <- NULL
   zipi80 <- NULL
   zipi90 <- NULL
   zipi95 <- NULL
   n.used <- NULL
   
   for(j in 1:n.sim) {
#if(j %in% seq(200,n.sim,200)) print(j)
      L <- c(rep(1,n/2),rep(0,n/2))
      A <- c(rep(0,n/4),rep(1,n/4),rep(0,n/4),rep(1,n/4))
      eta <- -1 + L + A
      X <- rnorm(n,eta,1)
      psi.hat <- solve(t(cbind(1,L,A))%*%cbind(1,L,A))%*%t(cbind(1,L,A))%*%X
      psi.cov <- solve(t(cbind(1,L,A))%*%cbind(1,L,A))*sum((X-cbind(1,L,A)%*%psi.hat)^2)/(n-3)

      eta.hat <- cbind(1,L,A)%*%psi.hat
      se.eta <- sqrt(psi.cov[1,1] + L^2*psi.cov[2,2] + A^2*psi.cov[3,3] + 2*L*psi.cov[2,1] + 2*A*psi.cov[3,1] + 2*L*A*psi.cov[2,3])
      mean.eta <- cbind(mean.eta,mean(eta.hat))

      mle <- cbind(mle,mean(eta.hat*(eta.hat>0)))

      cover0.80 <- (eta.hat - qnorm(.9)*se.eta < 0) & (eta.hat + qnorm(.9)*se.eta > 0)
      cover0.90 <- (eta.hat - qnorm(.95)*se.eta < 0) & (eta.hat + qnorm(.95)*se.eta > 0)
      cover0.95 <- (eta.hat - qnorm(.975)*se.eta < 0) & (eta.hat + qnorm(.975)*se.eta > 0)

      vess <- cbind(vess,mean(eta.hat[!cover0.95]*(eta.hat[!cover0.95]>0)))
      n.used <- cbind(n.used,length(eta.hat[!cover0.95]))

      w <- 1-pnorm(-abs(eta.hat/se.eta)) ### Multiply by 2 for 2-sided test
      pur <- cbind(pur,sum(w*eta.hat*(eta.hat>0))/sum(w))

      zipi80 <- cbind(zipi80,mean(ifelse(cover0.80,0,eta.hat)*(ifelse(cover0.80,0,eta.hat)>0)))
      zipi90 <- cbind(zipi90,mean(ifelse(cover0.90,0,eta.hat)*(ifelse(cover0.90,0,eta.hat)>0)))
      zipi95 <- cbind(zipi95,mean(ifelse(cover0.95,0,eta.hat)*(ifelse(cover0.95,0,eta.hat)>0)))
   }

   results <- rbind(results,c(apply(mean.eta,1,mean),apply(mean.eta,1,sd),
                        apply(mle,1,mean),apply(mle,1,sd),
                        apply(vess,1,mean,na.rm=T),apply(vess,1,sd,na.rm=T),
                        apply(pur,1,mean),apply(pur,1,sd),
                        apply(zipi80,1,mean),apply(zipi80,1,sd),
                        apply(zipi90,1,mean),apply(zipi90,1,sd),
                        apply(zipi95,1,mean),apply(zipi95,1,sd)))
   n.vess <- c(n.vess,apply(n.used,1,mean))
print(summary(as.vector(n.used)))
print(sum(as.vector(n.used)==0))
}
#print(round(results[,-(1:2)],3))
print(n.vess)

temp <- round(results[,-(1:2)],3)
paste(temp[,1]," (",temp[,2],") & ",temp[,3]," (",temp[,4],") & ",
    temp[,5]," (",temp[,6],")",sep="")
paste(temp[,7]," (",temp[,8],") & ",
    temp[,9]," (",temp[,10],") & ",temp[,11]," (",temp[,12],")",sep="")