###########################################################
require(foreach)
require(multicore)
require(doMC)

expit<-function(x){1/(1+exp(-x))}

n<-50
X<-runif(n,-2,2)
be0<-1
be1<-0.5
lpY<-be0+be1*X
Y<-rbinom(n,1,expit(lpY))
summary(glm(Y ~ X, family = binomial(logit)))
  
trials <- 10000
chunkSize <- ceiling(trials/getDoParWorkers())
smpopts <- list(chunkSize = chunkSize)
system.time(r <- foreach(icount(trials), .combine = cbind, .options.smp = smpopts) %dopar% {
  ind <- sample(n, n, replace = TRUE)
  result1 <- glm(Y[ind] ~ X[ind], family = binomial(logit))
  coefficients(result1)
})

hist(r[2,],nclass=50);abline(v=be1,lwd=2,col="red")