sliced Poisson

(This article was originally published at Xi'an's Og » R, and syndicated at StatsBlogs.)

One of my students complained that his slice sampler of a Poisson distribution was not working when following the instructions in Monte Carlo Statistical Methods (Exercise 8.5). This puzzled me during my early morning run and I checked on my way back, even before attacking the fresh baguette I had brought from the bakery… The following R code is the check. And it does work! As the comparison above shows…

slice=function(el,u){
#generate uniform over finite integer set
   mode=floor(lambda)
   sli=mode
   x=mode+1
   while (dpois(x,el)>u){
       sli=c(sli,x);x=x+1}
   x=mode-1
   while (dpois(x,el)>u){
       sli=c(sli,x);x=x-1}
   return(sample(sli,1))}

#example
T=10^4
lambda=2.414

x=rep(floor(lambda),T)
for (t in 2:T)
   x[t]=slice(lambda,runif(1)*dpois(x[t-1],lambda))

barplot(as.vector(rbind(
   table(x)/length(x),dpois(0:max(x),
   lambda))),col=c("sienna","gold"))

Please comment on the article here: Xi'an's Og » R

Tags: books, ENSAE, Kids, Monte Carlo Statistical Methods, pictures, Poisson distribution, R, Running, slice sampling, statistics, University life