simMC_pois10 <- function(n){

  x <- numeric(n)

  x[1] <- sample(c(0,1,2), 1, prob=rep(1/3,3))

  for(i in 2:n){
    tmp <- x[i-1]
    if(tmp <= 9){
      ## Caution: the zero case is no problem here!
      ## If tmp==0 then the probability that the next
      ## state is -1 is tmp/20 and hence 0.
      x[i] <- sample(c(tmp-1, tmp, tmp+1), 1, 
	      prob=c(tmp/20, (10-tmp)/20, 1/2))
    } else {
      x[i] <- sample(c(tmp-1, tmp, tmp+1), 1, 
	      prob=c(1/2, (tmp-9)/(2*(tmp+1)), 5/(tmp+1)))
    }
  }
  return(x)
}

## show traceplots
trace1 <- simMC_pois10(1000)
trace2 <- simMC_pois10(1000)
trace3 <- simMC_pois10(1000)

par(mfrow=c(3,1))
plot(trace1, type="l", xlab="Iteration", ylab="x")
plot(trace2, type="l", xlab="Iteration", ylab="x")
plot(trace3, type="l", xlab="Iteration", ylab="x")

## show histograms

plot_dist <- function(trace, until, limits){

  if(until > length(trace)){
    stop("the trace is not long enough!")
  }

  # for values 0, 1, ..., 20
  my_x = 0:20
  # get the corresponding density for Poisson(10) distribution
  my_y = dpois(my_x, lambda=10)

  # use only the chain until a certain position.
  tmp <- factor(trace[1:until], levels=0:20)
  tab <- table(tmp)
  plot(my_x, my_y, type="b", col=2, xlab="x", ylab="Density", 
    main=paste("Iteration", until), ylim=limits)
  # plot the proportion of times we had values 0:20 in the MC
  lines(0:20, tab/length(tmp), type="h", cex=1.2)
}
trace1 <- simMC_pois10(5000)


par(mfrow=c(3,1), cex.lab=1.2)
plot_dist(trace1, 1, limits=c(0,1))
plot_dist(trace1, 2, limits=c(0,1))
plot_dist(trace1, 3, limits=c(0,1))

par(mfrow=c(3,1), cex.lab=1.2)
plot_dist(trace1, 10, limits=c(0,0.5))
plot_dist(trace1, 20, limits=c(0,0.5))
plot_dist(trace1, 30, limits=c(0,0.5))

par(mfrow=c(3,1), cex.lab=1.2)
plot_dist(trace1, 50, limits=c(0,0.2))
plot_dist(trace1, 500, limits=c(0,0.2))
plot_dist(trace1, 5000, limits=c(0,0.2))