nreps<-10000

set.seed(84)
gam<-2.5
xold<-0
old.p<-dgamma(xold,lambda,log=T)

x<-rep(0,nreps)

png(file="Gamma%04d.png", width=600, height=600)

sig.q<-1

xv<-seq(0,10,by=0.001)
yv<-dgamma(xv,gam,1)

par(mar=c(4,4,4,4))
for(istep in 1:nreps){
	
	xnew<-abs(rnorm(1,xold,sig.q))
	new.p<-dgamma(xnew,gam,log=T)
	log.ratio<-(new.p-old.p)
	if(log(runif(1)) < log.ratio){
		xold<-xnew
		old.p<-new.p
	}
	x[istep]<-xold

	if(istep %% 50 == 0){
		hist(x[1:istep],breaks=seq(0,20,by=0.5),freq=FALSE,ylim=range(0,0.4),col='black',xlab='x',border='white',main='Gamma example')
		box()
		lines(xv,yv,col='red',lwd=2)
		text(15,0.4,istep)
	}

}
dev.off()

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

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

##############################################################################


nreps<-200

set.seed(84)
gam<-2.5
xold<-0
old.p<-dgamma(xold,lambda,log=T)

x<-rep(0,nreps)

png(file="Gamma%04d.png", width=600, height=600)

sig.q<-1

xv<-seq(0,10,by=0.001)
yv<-dgamma(xv,gam,1)

par(mar=c(4,4,4,4))
for(istep in 1:nreps){
	
	xnew<-abs(rnorm(1,xold,sig.q))
	new.p<-dgamma(xnew,gam,log=T)
	log.ratio<-(new.p-old.p)
	if(log(runif(1)) < log.ratio){
		xold<-xnew
		old.p<-new.p
	}
	x[istep]<-xold

	plot(1:istep,x[1:istep],type='s',xlim=range(0,nreps),ylim=range(0,10),xlab='Step',ylab=expression(x[t]))
}
dev.off()

#Home
#conv.cm<-'"C:\\Program Files\\ImageMagick-7.0.7-Q16\\convert.exe" -delay 25 *.png gamma_2.gif'

conv.cm<-'"C:\\Program Files\\ImageMagick-6.9.9-Q16-HDRI\\convert.exe" -delay 25 *.png gamma_2.gif'

system(conv.cm)

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

##############################################################################


nreps<-200

set.seed(84534)
gam<-2.5
xold<-6
old.p<-dgamma(xold,lambda,log=T)

x<-rep(0,nreps)

png(file="Gamma%04d.png", width=600, height=600)

sig.q<-1

xv<-seq(0,10,by=0.001)
yv<-dgamma(xv,gam,1)

par(mar=c(4,4,4,4))
for(istep in 1:nreps){
	
	xnew<-abs(rnorm(1,xold,sig.q))
	new.p<-dgamma(xnew,gam,log=T)
	log.ratio<-(new.p-old.p)
	if(log(runif(1)) < log.ratio){
		xold<-xnew
		old.p<-new.p
	}
	x[istep]<-xold

	plot(1:istep,x[1:istep],type='s',xlim=range(0,nreps),ylim=range(0,10),xlab='Step',ylab=expression(x[t]))
}
dev.off()

#Home
#conv.cm<-'"C:\\Program Files\\ImageMagick-7.0.7-Q16\\convert.exe" -delay 25 *.png gamma_3.gif'

conv.cm<-'"C:\\Program Files\\ImageMagick-6.9.9-Q16-HDRI\\convert.exe" -delay 25 *.png gamma_3.gif'

system(conv.cm)

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

##############################################################################


nreps<-200

set.seed(84534)
gam<-2.5
xold<-10
old.p<-dgamma(xold,lambda,log=T)

x<-rep(0,nreps)

png(file="Gamma%04d.png", width=600, height=600)

sig.q<-0.25

xv<-seq(0,10,by=0.001)
yv<-dgamma(xv,gam,1)

par(mar=c(4,4,4,4))
for(istep in 1:nreps){
	
	xnew<-abs(rnorm(1,xold,sig.q))
	new.p<-dgamma(xnew,gam,log=T)
	log.ratio<-(new.p-old.p)
	if(log(runif(1)) < log.ratio){
		xold<-xnew
		old.p<-new.p
	}
	x[istep]<-xold

	plot(1:istep,x[1:istep],type='s',xlim=range(0,nreps),ylim=range(0,10),xlab='Step',ylab=expression(x[t]))
}
dev.off()

#Home
#conv.cm<-'"C:\\Program Files\\ImageMagick-7.0.7-Q16\\convert.exe" -delay 25 *.png gamma_4.gif'

conv.cm<-'"C:\\Program Files\\ImageMagick-6.9.9-Q16-HDRI\\convert.exe" -delay 25 *.png gamma_4.gif'

system(conv.cm)

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