Another Bernoulli factory

The paper “Exact sampling for intractable probability distributions via a Bernoulli factory” by James Flegal and Radu Herbei got posted on arXiv without me noticing, presumably because it came out just between Larry Brown’s conference in Philadelphia and my skiing vacations! I became aware of it only yesterday and find it quite interesting in that it links the Bernoulli factory method I discussed a while ago and my ultimate perfect sampling paper with Jim Hobert. In this 2004 paper in Annals of Applied Probability, we got a representation of the stationary distribution of a Markov chain as

\sum_{n=1}^{\infty} p_n Q_n(dx)


p_n = \mathbb{P}(\tau\ge n)\qquad\text{and}\qquad Q_n(A)=\mathbb{P}(X_n\in A|\tau\ge n),

the stopping time τ being the first occurrence of a renewal event in the split chain. While Q_n is reasonably easy to simulate by rejection (even tohugh it may prove lengthy when n is large, simulating from the tail distribution of the stopping time is much harder.

“Obtaining a single draw from π would require an obscene number of draws.” J. Flegal and R. Herbei

The solution of Flegal and Herbei is to use (a) accept-reject, (b) a bound obtained in our paper, and (c) the Bernoulli factory. That is,

\mathbb{P}(\tau\le n) \le\dfrac{d(n) M}{\mathbb{E}[\tau]} \quad\text{with}\quad d(n)=(\beta-1)^{-1}\beta^{-n}

where β and M depend on the minorisation conditions for the Markov chain. The Bernoulli factory is used to generate Bernoulli variables with probability ap where

a=1/[Md(n)]\qquad\text{and}\qquad p=\mathbb{P}(\tau\ge n)

The exact contribution of the paper consists in improving upon Nacu and Peres‘ original algorithm by turning the original function min{2p,1-w} into a twice differentiable function. This makes the bound in the sandwiching algorithm of Latuszyński et al. moving from n^{-1/2} to n^{-1}. As Krzysztof Latuszyński pointed out to me, it is even possible to improve the convergence to an exponential rate for the Bernoulli factory, but this does not bring an improvement in the [our] perfect sampling algorithm that motivates the study. The later remains with an infinite expected running time (Blanchet and Meng). This difficulty is reflected in the experiments where an implementation for the normal-gamma posterior required 35 hours of computation, one single call to the Bernoulli factory accounting for 70% of the computing time! The realistic random effect model is rather anticlimactic in that the final outcome of the three pages 19-22 is a single realisation from the exact sampler. Noticeably, no output from the Bernoulli factory was used in any of the examples implemented in the paper. Nonetheless, I find the paper quite interesting from the Bernoulli factory point of view and prone to open new research directions on this “cute” problem.

2 Responses to “Another Bernoulli factory”

  1. […] Another Bernoulli factory « Xi'an's Og […]

  2. James Flegal had a very nice poster about this at MCMSki 3. He mentioned that the calls to the Bernoulli factory that take the longest are virtually certain to be rejected anyway, so he’s thinking about ways to detect and avert unnecessary computation.

Leave a Reply

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

You are commenting using your 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 )

Connecting to %s

This site uses Akismet to reduce spam. Learn how your comment data is processed.

%d bloggers like this: