The Bernoulli factory

April 22, 2010

(This article was first published on Xi'an's Og » R, and kindly contributed to R-bloggers)

A few months ago, Latuszyński, Kosmidis, Papaspiliopoulos and Roberts arXived a paper I should have noticed earlier as its topic is very much related to our paper with Randal Douc on the vanilla Rao-Blackwellisation scheme. It is motivated by the Bernoulli factory problem, which aims at (unbiasedly) estimating f(p) from an iid sequence of Bernoulli B(p). (The paper only considers functions f valued in [0,1]. In our case, the function is f(p)=1/p.) It appears that this problem has been studied for quite a while, in particular by Yuval Peres. Being in a train to Marseille (thanks to Eyjafjallajökull!), I do not have access to those earlier papers of Peres’, but Latuszyński et al. mentions that there are conditions on f such that it is sufficient to generate a Bernoulli event with probability

f_0(p) = min (2p, 1-2epsilon)

where varepsilon>0 is arbitrary. In particular, constructing an unbiased estimator of

f_0(p) = min ( 2p, 1 )

does not seem to be achievable (Nacu and Peres, 2005). (The way it is rephrased in Latuszyński et al. does not seem correct, though, as they state that f(p)=2p cannot be estimated in an unbiased manner, missing the constraint that the estimator must belong to [0,1], I think.)

The paper by Latuszyński et al. develops an original scheme to achieve simulation from B(f(p)) through the simulation of two bounding sequences that are respectively super- and submartingales and that both converge to f(p). (But their simulation scheme does not have to wait for the two sequences to coalesce.) This idea presumably (?) stemmed from the experience of the authors, in particular Gareth Roberts, in perfect sampling, since the most advanced perfect samplers made intensive use of this sandwiching idea of Kendall and Møller (2000, Advances in Applied Probability). The whole thing is also related to the famous Series B paper of Beskos et al. (2006). The method of Latuszyński et al. then builds the upper and lower processes via a truncated series decomposion of f(p), whose availability induces constraints on f.

The first application illustrated in Latuszyński et al. is the unbiased estimation of a transform f(p) that has a known series expansion

f(p) = sum_{i=1}^infty (1-)^k a_k p^k


1le a_1le a_2 le cdots

In that case, we could use the scheme of our paper with Randal, estimating p^k by

hat p^k = X_1ldots X_k.

The probability of using at least n simulations is then p^n, while the scheme of Latuszyński et al. leads to a probability of  a_n p^n. (Note however that the direct approach above allows to handle any series decomposition, alternating or not, with no constraint on the a_i‘s.)

What I find exciting about this Bernoulli factory problem is that the well-known issue of the absence of unbiased estimators for most transforms of a parameter p (Lehmann and Casella, 1998) vanishes when an unlimited number of iid simulations with mean p is available. Here are the slides of the talk given by Omiros last week at the Big’ MC seminar:

Filed under: R, Statistics Tagged: Bernouilli factory, control variates, MCMC, Metropolis-Hastings, Monte Carlo methods, Rao-Blackwellisation, simulation, slides

To leave a comment for the author, please follow the link and comment on their blog: Xi'an's Og » R. offers daily e-mail updates about R news and tutorials on topics such as: Data science, Big Data, R jobs, visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...

If you got this far, why not subscribe for updates from the site? Choose your flavor: e-mail, twitter, RSS, or facebook...

Tags: , , , , , , , , ,

Comments are closed.


Never miss an update!
Subscribe to R-bloggers to receive
e-mails with the latest R posts.
(You will not see this message again.)

Click here to close (This popup will not appear again)