sliced Poisson

slissonOne 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"))

4 Responses to “sliced Poisson”

  1. […] Poisson slice sampler of a Poisson distribution, I could not tell why the code was not working! My earlier post prompted him to do so and a somewhat optimised version is given […]

  2. The version of this that showed up on R-bloggers was missing the ” mode=floor(lambda)” line, which confused me greatly until I came here to check it out …

  3. […] article was first published on Xi'an's Og » R, and kindly contributed to […]

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

Follow

Get every new post delivered to your Inbox.

Join 667 other followers