#Brian OMeara
#http://www.brianomeara.info
#released under GPL2
mlVSBayes.ani<-function(nreps=1000) {
	par(mfrow=c(2,4))
	likelihood<-density(c(rnorm(10000,mean=20,sd=3),rnorm(8000,mean=30,sd=4),runif(4000,10,31),rnorm(5000,mean=10,sd=2)),from=0,to=50,n=1000)
	likelihood$y=likelihood$y/sum(likelihood$y)
	prior1<-density(runif(1000000,0,50),from=0,to=50,n=1000)
	prior1$y<-prior1$y/sum(prior1$y)
	prior2<-density(rnorm(10000,30,5),from=0,to=50,n=1000)
	prior2$y<-prior2$y/sum(prior2$y)
	prior3<-density(1+rexp(10000,.1),from=0,to=50,n=1000)
	prior3$y<-prior3$y/sum(prior3$y)
	posteriors1<-prior1$y*likelihood$y/sum(prior1$y*likelihood$y)
	posteriors2<-prior2$y*likelihood$y/sum(prior2$y*likelihood$y)
	posteriors3<-prior3$y*likelihood$y/sum(prior3$y*likelihood$y)

	maxY<-1.5*max(c(prior1$y,prior2$y,prior3$y))
	titles<-c("No prior (likelihood)", "Prior 1", "Prior 2", "Prior 3","Likelihood", "Posterior 1", "Posterior 2", "Posterior 3")
	for (i in 1:8) {
		plot(x=c(0,50),y=c(0,maxY),type="n",bty="n",ylab="",xlab="",yaxt="n",main=titles[i])
	}
	par(mfg=c(2,2))
	lines(x=prior1$x,y=prior1$y*likelihood$y/sum(prior1$y*likelihood$y))
	par(mfg=c(2,3))
	lines(x=prior2$x,y=prior2$y*likelihood$y/sum(prior2$y*likelihood$y))
	par(mfg=c(2,4))
	lines(x=prior3$x,y=prior3$y*likelihood$y/sum(prior3$y*likelihood$y))

	par(mfg=c(1,2))
	lines(prior1)
	par(mfg=c(1,3))
	lines(prior2)
	par(mfg=c(1,4))
	lines(prior3)
	par(mfg=c(2,1))
	lines(x=likelihood$x, y=likelihood$y/sum(likelihood$y))
	par(mfg=c(2,2))
	lines(x=prior1$x,y=prior1$y*likelihood$y/sum(prior1$y*likelihood$y))
	par(mfg=c(2,3))
	lines(x=prior2$x,y=prior2$y*likelihood$y/sum(prior2$y*likelihood$y))
	par(mfg=c(2,4))
	lines(x=prior3$x,y=prior3$y*likelihood$y/sum(prior3$y*likelihood$y))
	maxlnLValue<-0
	maxlnLPoint<-0
	sampling<-matrix(data=sample.int(length(likelihood$x),4),nrow=nreps,ncol=4)
	scores<-c(likelihood$y[sampling[1,1]],posteriors1[sampling[1,2]],posteriors2[sampling[1,3]],posteriors3[sampling[1,3]])
	for (replicate in 2:nreps) {
			for(subplot in 1:4) {
				if (subplot==1) {
					newsample<-sample.int(length(likelihood$x),1)
					while(abs(newsample-sampling[replicate-1, 1])>100) {
											newsample<-sample.int(length(likelihood$x),1)
					}
					if (runif(1)<0.01) { #random restart
							newsample<-sample.int(length(likelihood$x),1)
					}

					if (likelihood$y[newsample]>likelihood$y[sampling[replicate -1, 1] ]) {
							sampling[replicate, 1]<-newsample
							maxlnLValue<-likelihood$y[newsample]
							maxlnLPoint<-likelihood$x[newsample]
					}
					else {
							sampling[replicate, 1]<-sampling[replicate-1, 1]
					}
					par(mfg=c(2,1))

						lines(x=likelihood$x, y=likelihood$y/sum(likelihood$y),col="black")

					points(x=likelihood$x[ sampling[replicate, 1]],y=likelihood$y[ sampling[replicate, 1]],col=rgb(0,0.1,1,0.5))
				}	
				else {
					newsample<-sample.int(length(likelihood$x),1)
					while(abs(newsample-sampling[replicate-1, subplot])>100) {
											newsample<-sample.int(length(likelihood$x),1)
					}
					if (runif(1)<0.01) { #random restart
							newsample<-sample.int(length(likelihood$x),1)
					}

					posteriors<-c()
					if (subplot==2) {
							posteriors<-posteriors1
					}
					else if (subplot==3) {
						posteriors<-posteriors2
					}
					else if (subplot==4) {
						posteriors<-posteriors3
					}
					if (posteriors[newsample]> posteriors[sampling[replicate-1, subplot] ]) {
							sampling[replicate, subplot]<-newsample
					}
					else if ((posteriors[newsample] / posteriors[sampling[replicate-1, subplot] ]) > runif(1)) {
							sampling[replicate, subplot]<-newsample

					}
					else {
							sampling[replicate, subplot]<-sampling[replicate-1, subplot]
					}
					par(mfg=c(2,subplot))


					estimatedPosterior<-density(likelihood$x[sampling[1:replicate,subplot]],from=0,to=50,n=1000)
					lines(x=estimatedPosterior$x,estimatedPosterior$y/sum(estimatedPosterior$y),col=rgb(0,.1,1,0.02))
					lines(x=likelihood$x,y=posteriors/sum(posteriors),col=rgb(0,0,0,1),lwd=2)
				}
					#print(sampling[replicate, ])

			}
	}
	for (subplot in 2:4) {
					par(mfg=c(2,subplot))
 				estimatedPosterior<-density(likelihood$x[sampling[1:nreps,subplot]],from=0,to=50,n=1000)
					lines(x=estimatedPosterior$x,estimatedPosterior$y/sum(estimatedPosterior$y),col=rgb(1,0,0,1),lwd=2)

	}
	par(mfg=c(2,1))
	points(maxlnLPoint, maxlnLValue,col=rgb(1,0,0,1),bg=rgb(1,0,0,1),cex=2,pch=19)
}


 mlVSBayes.ani()