#install.packages(c('foreach','doMC','multicore'))

###########################################################
require(foreach)
require(multicore)
require(doMC)
library(MASS)

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

set.seed(1234)
n<-50
X<-runif(n,-2,2)
be0<-1
be1<-0.5
lpY<-be0+be1*X
Y<-rbinom(n,1,expit(lpY))
fit.mle<-glm(Y ~ X, family = binomial(logit))

beta.hat<-coef(fit.mle)
beta.Sig<-1.5*summary(fit.mle)$cov.unscaled

logistic.loglike<-function(Xv,Yv,bhat){
  lp<-bhat[1]+bhat[2]*Xv
  pv<-expit(lp)
  return(sum(dbinom(Yv,1,pv,log=T)))
}

run.mcmc.metropolis<-function(nr,nt,Xv,Yv,bhat,bSig){

  Prec<-solve(bSig)
  old.be<-bhat
  old.like<-logistic.loglike(Xv,Yv,bhat)
  old.prior<--0.5*sum(bhat^2)/10
  old.prop<--0.5*t(old.be-bhat) %*% Prec %*% (old.be-bhat)
  val.vec<-matrix(0,nrow=floor(nr/nt),ncol=2)
  ico<-0
  for(irep in 1:nr){
      new.be<-mvrnorm(1,bhat,bSig)
      new.like<-logistic.loglike(Xv,Yv,new.be)
      new.prior<--0.5*sum(new.be^2)/10
      new.prop<--0.5*t(new.be-bhat) %*% Prec %*% (new.be-bhat)
      if(log(runif(1)) < (new.like+new.prior-new.prop)-(old.like+old.prior-old.prop)){
        old.be<-new.be
        old.prior<-new.prior
        old.prop<-new.prop
      }
      if(irep %% nt==0){
        ico<-ico+1
        val.vec[ico,]<-old.be
      }      
  }
  return(val.vec)
}

#Serial
trials <- 10
chunkSize <- ceiling(trials/getDoParWorkers())
smpopts <- list(chunkSize = chunkSize)
system.time(mcmc.met.samples <- foreach(icount(trials), .combine = cbind, .options.smp = smpopts) %do% {
  
  nreps<-10000
  nthin<-10
  vals<-run.mcmc.metropolis(nreps,nthin,X,Y,beta.hat,beta.Sig)
  
})

#Parallel
registerDoMC(cores=2)
getDoParWorkers()
trials <- 10
chunkSize <- ceiling(trials/getDoParWorkers())
smpopts <- list(chunkSize = chunkSize)
system.time(mcmc.met.samples <- foreach(icount(trials), .combine = cbind, .options.smp = smpopts) %dopar% {

    nreps<-10000
    nthin<-10
    vals<-run.mcmc.metropolis(nreps,nthin,X,Y,beta.hat,beta.Sig)
    
})

be0.mat<-mcmc.met.samples[,c(1:trials)*2-1]
be1.mat<-mcmc.met.samples[,c(1:trials)*2]
boxplot(be0.mat)
boxplot(be1.mat)
boxplot(cbind(be0.mat,be1.mat))

hist(be0.mat)
hist(be1.mat)

range(be1.mat)

par(pty="s")
plot(be0.mat[,1],be1.mat[,1],ylim=range(-2,2),xlim=range(0,4),pch=19,cex=0.25,xlab=expression(beta[0]),ylab=expression(beta[1]))
for(j in 2:trials){
    points(be0.mat[,j],be1.mat[,j],pch=19,col=heat.colors(trials)[j],cex=0.25)
}

ones<-rep(1,n)
be0v<-seq(0,4,length=101)
be1v<-seq(-2,2,length=101)
be0m<-outer(ones,be0v,"*")
be1xm<-outer(X,be1v,"*")

like.compute<-array(0,c(n,101,101))
for(i in 1:n){
  like.compute[i,,]<-outer(be0v,be1v*X[i],function(x,y){pv<-expit(x+y);return(dbinom(Y[i],1,pv,log=T))})
  
}
like.compute.surf<-apply(like.compute,2:3,sum)
contour(be0v,be1v,like.compute.surf,add=T,nlevels=50)