nsteps<-1000

rho<-0.6
x<-c(0,0)
xmat<-matrix(0,nrow=nsteps,ncol=2)

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

par(pty='s',mar=c(4,4,4,4))
plot(x[1],x[2],xlim=range(-3.5,3.5),ylim=range(-3.5,3.5),type='n',
           xlab=expression(x[1]),ylab=expression(x[2]))

for(istep in 1:nsteps){

	x[1]<-rnorm(1,rho*x[2],sqrt(1-rho^2))
	x[2]<-rnorm(1,rho*x[1],sqrt(1-rho^2))
	xmat[istep,]<-x

	plot(xmat[1:istep,1],xmat[1:istep,2],xlim=range(-3.5,3.5),ylim=range(-3.5,3.5),col='gray',pch=19,
           xlab=expression(x[1]),ylab=expression(x[2]))
	points(xmat[istep,1],xmat[istep,2],col='red',pch=19)

}
dev.off()

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

#system('mv *.png .\\Gibbs01\\')


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

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

nsteps<-1000

rho<--0.95
x<-c(0,0)
xmat<-matrix(0,nrow=nsteps,ncol=2)

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

par(pty='s',mar=c(4,4,4,4))
plot(x[1],x[2],xlim=range(-3.5,3.5),ylim=range(-3.5,3.5),type='n',
           xlab=expression(x[1]),ylab=expression(x[2]))

for(istep in 1:nsteps){

	x[1]<-rnorm(1,rho*x[2],sqrt(1-rho^2))
	x[2]<-rnorm(1,rho*x[1],sqrt(1-rho^2))
	xmat[istep,]<-x

	plot(xmat[1:istep,1],xmat[1:istep,2],xlim=range(-3.5,3.5),ylim=range(-3.5,3.5),col='gray',pch=19,
           xlab=expression(x[1]),ylab=expression(x[2]))
	points(xmat[istep,1],xmat[istep,2],col='red',pch=19)

}
dev.off()

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



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