Sys.time()
rm(list=ls())
seed <- 1000
set.seed(seed)  
library(MASS)
library(Matrix)
library(mvtnorm)
library(ada)
library(qvalue)
library(pi0)

lambda <- seq(0.2,0.5,by=0.05)
num.gene <- 1000
num.sample <- 10        ### num.sample=5 or 10
rho <- 0                ### rho=0 or 0.4 or 0.8
nu <- num.sample-1 

#generate the AR(1) block covariance matrix          
b <- 50
unit <- rep(1,b)
times <- 1:b
H <- abs(outer(times, times, "-"))
Block.Cov <- rho^H
#Block.Cov <- diag(rep(1-rho,b)) + rho*unit%*%t(unit)
Cov <- diag(rep(0,num.gene))
adjust.const <- gamma((num.sample-1)/2)/gamma((num.sample-2)/2)*sqrt(2/(num.sample-1))

pi0 <- 1:9/10
pi0.Storey <- pi0.Jiang <- pi0.Langaas <- pi0.New <- NULL
repeats <- 1000
for (s in 1:length(pi0)){
 pi0.storey <- pi0.jiang <- pi0.langaas <- pi0.new <- NULL
 for (k in 1:repeats){ 
  mu <- c(rep(0,trunc(pi0[s]*num.gene)),runif(num.gene-trunc(pi0[s]*num.gene),0.5,1.5))
  mu <- sample(mu) #resample the mean values
  #block covariance matrix with 20 blocks
  Sig <- rchisq(num.gene/b,df=10)/10
  for (j in 1:(num.gene/b))
    Cov[((j-1)*b+1):(j*b),((j-1)*b+1):(j*b)] <- Sig[j]*Block.Cov
  X <- t(rmvnorm(num.sample, mu, Cov))

  #compute p-values
  avg.X <- apply(X,1,mean)
  s2.X <- apply(X,1,var)
  tstat <- (avg.X/sqrt(s2.X))*sqrt(num.sample) 
  pval <- 2*pt(-abs(tstat),num.sample-1) 

  #Storey method
  pi0.storey <- c(pi0.storey, qvalue(pval)$pi0)

  #Jiang method
  W <- NULL
  for (i in 1:length(lambda))
     W <- c(W, sum(pval>lambda[i]))
  pi0.tmp <- W/(num.gene*(1-lambda)) 
  pi0.jiang <- c(pi0.jiang, min(1,max(0,mean(pi0.tmp))))

  #Langaas method
  pi0.langaas <- c(pi0.langaas, convest(pval))

  ###New method### 
  #step (a)   
  pi0.ini <- pi0.storey[k]    
#  pi0.ini <- pi0.langaas[k]  
  #step (b)
  delta.est <- avg.X/sqrt(s2.X)*adjust.const
  #step (cd)
  z.lambda <- qnorm(1-lambda/2) 
  t.lambda <- qt(1-lambda/2,df=nu) 
  Q.lambda <- NULL
  for (i in 1:length(lambda)){
    Q.tmp <- pt(t.lambda[i],df=nu,ncp=sqrt(num.sample)*delta.est) - pt(-t.lambda[i],df=nu,ncp=sqrt(num.sample)*delta.est)
    Q.lambda <- c(Q.lambda, mean(sort(Q.tmp)[1:trunc(num.gene*(1-pi0.ini))]))
  }
  #step (e)   
  pi0.new.tmp <- (W-num.gene*Q.lambda)/(num.gene*(1-lambda)-num.gene*Q.lambda)
  pi0.new <- c(pi0.new, min(1,max(0,mean(pi0.new.tmp))))    
 }
 pi0.Storey <- cbind(pi0.Storey, pi0.storey)
 pi0.Jiang <- cbind(pi0.Jiang, pi0.jiang)
 pi0.Langaas <- cbind(pi0.Langaas, pi0.langaas)
 pi0.New <- cbind(pi0.New, pi0.new)
}

result <- NULL
result$seed <- seed
result$pi0.Storey <- pi0.Storey
result$pi0.Jiang <- pi0.Jiang
result$pi0.Langaas <- pi0.Langaas
result$pi0.New <- pi0.New
result$num.gene <- num.gene
result$num.sample <- num.sample
result$rho <- rho
result$repeats <- repeats
result$pi0 <- pi0 
result$lambda <- lambda
result$b <- b
save("result", file="result.n10.rho00")


Sys.time()