MCbias.study.par <- function(nr,nv,be0,be1){
  require(foreach)
  est.coef <- function(nr,nv,be0,be1) {
    X<-runif(nv,-10,10)
    Y<-be0+be1*X+rnorm(nv)
    return(coef(lm(Y~X))[2])
  }
  nodes <- getDoParWorkers()
  print(paste('Nodes being used: ',nodes))
  p <- foreach(mr = rep(nr/nodes, nodes),.combine='c') %dopar% est.coef(mr,nv,be0,be1)
  #print(p)
  return(bias <- mean(p)-be1)
}
print(b)
set.seed(123)
M<-1e7
system.time(b<-MCbias.study.par(M,1000000,1,2))



require(snowfall)
sfInit(parallel = TRUE, cpus = 2)
cl <- sfGetCluster()
clusterSetupRNG(cl, seed = rep(123, 2))
require(foreach)
require(doSNOW)
registerDoSNOW(cl)

getDoParWorkers()

set.seed(123)
M<-1e7
system.time(b<-MCbias.study.par(M,1000000,1,2))
sfStop()
registerDoSEQ()