nreps<-1000

set.seed(84)

library(MASS)
library(mnormt)
n<-25
X<-cbind(1,mvrnorm(n,mu=c(0,0),Sigma=matrix(c(1,0.8,0.8,1),2,2)))
be<-c(1,2,-2)

eta<-X %*% be
p<-exp(eta)/(1+exp(eta))
Y<-rbinom(n,1,p)


lr.loglike<-function(bev,yv,Xm){
	etav<-Xm %*% bev
	p<-exp(etav)/(1+exp(etav))
	return(sum(yv*log(p)+(1-yv)*log(1-p)))
}

be.old<-c(0,0,0)

old.like<-lr.loglike(be.old,Y,X)
old.prior<-dmnorm(be.old,rep(0,3),diag(rep(5,3)),log=T)

be.samp<-matrix(0,nrow=nreps,ncol=3)

png(file="Logisticn%04d.png", width=600, height=600)
for(istep in 1:nreps){
	
	for(j in 1:3){
		be.new<-be.old
		be.new[j]<-be.old[j]+rnorm(1)
		new.like<-lr.loglike(be.new,Y,X)
		new.prior<-dmnorm(be.new,rep(0,3),diag(rep(5,3)),log=T)
		log.ratio<-(new.like+new.prior-old.like-old.prior)
		if(log(runif(1)) < log.ratio){
			be.old<-be.new
			old.like<-new.like
			old.prior<-new.prior
		}
	}

	be.samp[istep,]<-be.old

	if(istep %% 1 == 0){
		par(mfrow=c(2,2),mar=c(2,2,2,2))
		plot(be.samp[1:istep,1],be.samp[1:istep,2],pch=19,cex=0.5,col='gray',xlim=range(-5,5),ylim=range(-5,5))
		points(be.samp[istep,1],be.samp[istep,2],pch=19,cex=0.8,col='red')
		plot(be.samp[1:istep,1],be.samp[1:istep,3],pch=19,cex=0.5,col='gray',xlim=range(-5,5),ylim=range(-5,5))
		points(be.samp[istep,1],be.samp[istep,3],pch=19,cex=0.8,col='red')
		plot(be.samp[1:istep,2],be.samp[1:istep,3],pch=19,cex=0.5,col='gray',xlim=range(-5,5),ylim=range(-5,5))
		points(be.samp[istep,2],be.samp[istep,3],pch=19,cex=0.8,col='red')
	}


}
dev.off()

#conv.cm<-'"C:\\Program Files\\ImageMagick-7.0.7-Q16\\convert.exe" -delay 5 *.png logistic_1.gif'
conv.cm<-'"C:\\Program Files\\ImageMagick-6.9.9-Q16-HDRI\\convert.exe" -delay 5 *.png logistic_1.gif'
system(conv.cm)


file.remove(list.files(pattern=".png"))