set.seed(4532)
N<-10000
lam<-2.5
X<-rep(0,N)
X[1]<-0
for(istep in 2:N){
	if(X[istep-1] == 0){
		Z<-1
		qrat<-1/2
		al<-min(1,qrat*dpois(Z,lam)/dpois(X[istep-1],lam))
	}else{
		Z<-sample(c(X[istep-1]-1,X[istep-1]+1),size=1,prob=c(1/2,1/2))
		if(Z == 0){
			qrat<-2
		}else{
			qrat<-1
		}
		al<-min(1,qrat*dpois(Z,lam)/dpois(X[istep-1],lam))
	}
	u<-runif(1)

	if(u < al){
		X[istep]<-Z
	}else{
		X[istep]<-X[istep-1]
	}
}


plot(X[1:200],type="s",xlab="t",ylab="X",lwd=2)

hist(X,br=c(-1:12),xlab="t",main=" ",col="black",border="white")
points(c(0:12)-0.5,dpois(c(0:12),lam)*N,col="red",pch=3)

plot(c(20:N),tvn,type="l",xlab="N",ylab="TV")

###################################################################
set.seed(4532)
N<-100000
lam<-2.5
X<-rep(0,N)
X[1]<-0
for(istep in 2:N){
	if(X[istep-1] == 0){
		Z<-1
		qrat<-1/2
		al<-min(1,qrat*dpois(Z,lam)/dpois(X[istep-1],lam))
	}else{
		Z<-sample(c(X[istep-1]-1,X[istep-1]+1),size=1,prob=c(1/2,1/2))
		if(Z == 0){
			qrat<-2
		}else{
			qrat<-1
		}
		al<-min(1,qrat*dpois(Z,lam)/dpois(X[istep-1],lam))
	}
	u<-runif(1)

	if(u < al){
		X[istep]<-Z
	}else{
		X[istep]<-X[istep-1]
	}
}


round(rbind(table(X)/N,dpois(c(0:11),lam)),4)

iburn<-1000
tvn<-c(iburn:N)*0
Xlim<-100
for(i in iburn:N){
	
	Xtab<-table(X[1:i])
	pihat<-Xtab/i
	pivals<-as.numeric(names(Xtab))
	pihatvec<-rep(0,Xlim);pihatvec[pivals+1]<-pihat
	pivec<-dpois(c(0:(Xlim-1)),lam)

	tvn[i-iburn+1]<-sum(abs(pivec-pihatvec))/2

}

plot(c(iburn:N),tvn,type="l",xlab="N",ylab="TV")