## fine-sliced Poisson [a.k.a. sashimi]

**A**s my student Kévin Guimard had not mailed me his own 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 below:

nsim = 10^4 lambda = 6 max.factorial = function(x,u){ k = x parf=1 while (parf*u<1){ k = k + 1 parf = parf * k } k = k - (parf*u>1) return (k) } x = rep(floor(lambda), nsim) for (t in 2:nsim){ v1 = ceiling((log(runif(1))/log(lambda))+x[t-1]) ranj=max(0,v1):max.factorial(x[t-1],runif(1)) x[t]=sample(ranj,size=1) } barplot(as.vector(rbind( table(x)/length(x),dpois(min(x):max(x), lambda))),col=c("sienna","gold"))

**A**s you can easily check by running the code, it does not work. My student actually majored my MCMC class and he spent quite a while pondering why the code was not working. I did ponder as well for a part of a morning in Warwick, removing causes for exponential or factorial overflows (hence the shape of the code), but not eliciting the issue… (This now sounds like lethal fugu sashimi! ) Before reading any further, can you spot the problem?!

**T**he corrected R code is as follows:

x = rep(lambda, nsim) for (t in 2:nsim){ v1=ceiling((log(runif(1))/log(lambda))+x[t-1]) ranj=max(0,v1):max.factorial(x[t-1],runif(1)) if (length(ranj)>1){ x[t] = sample(ranj, size = 1) }else{ x[t]=ranj} }

**T**he culprit is thus the R function * sample* which simply does not understand Dirac masses and the basics of probability! When running

> sample(150:150,1) [1] 23

you can clearly see where the problem stands…! Well-documented issue with * sample* that already caused me woes… Another interesting thing about this slice sampler is that it is awfully slow in exploring the tails. And to converge to the centre from the tails. This is not very pronounced in the above graph with a mean of 6. Moving to 50 makes it more apparent:

**T**his is due to the poor mixing of the chain, as shown by the raw sequence below, which strives to achieve a single cycle out of 10⁵ iterations! In any case, thanks to Kévin for an interesting morning!

March 20, 2014 at 9:47 pm

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

March 20, 2014 at 12:26 pm

The diag() function has a similar problem. I’m afraid that we are locked into these mis-features because changing them would mean breaking old code.

The slice sampler is normally very efficient. If you are seeing poor mixing then it means that your slice is too small. In fact I can see three problems with the code:

1) The upper limit should depend on lambda

2) You are calculating the upper and lower bounds of distinct slices, i.e. you should draw u = runif(1) only once and use u for both limits.

3) I don’t understand the short-cut used to calculate the lower bound. If this works at all then it must be only when the current value of x is in the lower tail. Otherwise you have to do a similar recursive calculation to the upper limit.

March 20, 2014 at 10:48 pm

Thanks Martyn, I do not get 2). If I have two slices, I have two uniforms, otherwise the marginal gets wrong. As for 3), this is a log-transform v of the uniform for the power of λ to remain under control. Incorporating the Jacobian since v is then a truncated exponential distribution.