p = '205225_at'  # ESR1

# "old" way
source('bstrap.R')

bstrap.tst(X[C$er==1,10], X[C$er==0,10], 10000, tst.student)
bstrap.tst(X[C$er==1, p], X[C$er==0, p], 10000, tst.student)


# Boostrap testing: compute the test statistic in parallel
source('bstrap2.R') # on master too
p = '205225_at'
RNGkind("L'Ecuyer-CMRG")

cl = makeCluster(4, type='PSOCK')

clusterSetRNGStream(cl, 1234)  # each worker has an independent stream of RN
  
# force each worker to load the code for bootstrapping
clusterEvalQ(cl, 
  {
   load('transbig.rdata')
   source('bstrap2.R')
   p = '205225_at'  # ESR1
   B = 100000       # total number of bootstraps
   ## you can add other code here for initialization
  }
  )
    
res = clusterEvalQ(cl,
  {
    bstrap.tst2(X[C$er==1, p], X[C$er==0, p], B/4, tst.meandiff)
  }
  )

stopCluster(cl)
    
theta.star = do.call(c, res)
theta.0    = tst.meandiff(X[C$er==1, p], X[C$er==0, p])

pv = sum(theta.star >= theta.0) / length(theta.star)