Suppose we are given a coin with arbitrary (unknown) head probability $p$, I am wondering if there is an easy-to-implement algorithm for generating a $\min\{p, 0.5\}$ coin for any $p\in [0,1]$. Function $f(p) = \min\{p, 0.5\} \geq \min\{p, 1-p\}$ so the famous paper by Keane and O’Brien guarantees such an algorithm exists. Meanwhile, the function looks slightly simpler than the function discussed in Nacu and Peres, so I am wondering if there is any existing simple algorithm for simulating this. Similarly, can we simulate $f(p) = \min\{p,a\}$ for fixed $a\in (0,1)$ efficiently? Thanks!

  • 1
    $\begingroup$ The result by Nacu and Peres shows the answer is in the negative, because the function $p\to\min(p,1/2)$ is not real analytic on any neighborhood of $1/2.$ $\endgroup$
    – whuber
    Oct 27, 2020 at 16:02
  • $\begingroup$ Thanks! I agree it's not analytic, but here I am wondering if there are simple algorithms which may not be fast (exponential on the number of coins used) but easy to implement $\endgroup$
    – Bravo
    Oct 27, 2020 at 20:23
  • $\begingroup$ So by "efficient algorithm" do you mean "easy to implement"? If so, it would be better to change your terminology to reflect what you mean. $\endgroup$
    – whuber
    Oct 27, 2020 at 20:28
  • $\begingroup$ good point, just edited! $\endgroup$
    – Bravo
    Oct 27, 2020 at 21:46

3 Answers 3


Note: This answer is wrong/incomplete and will need an update according to the comments

Y̶o̶u̶ ̶c̶a̶n̶ ̶a̶p̶p̶r̶o̶a̶c̶h̶ ̶t̶h̶e̶ ̶f̶u̶n̶c̶t̶i̶o̶n̶ ̶b̶u̶t̶ ̶n̶o̶t̶ ̶g̶e̶t̶ ̶e̶q̶u̶a̶l̶

Suppose we are given a coin with arbitrary (unknown) head probability $p$, I am wondering if there is an easy-to-implement algorithm for generating a $\min\{p, 0.5\}$ coin for any $p\in [0,1]$. Function $f(p) = \min\{p, 0.5\} \geq \min\{p, 1-p\}$ so the famous paper by Keane and O’Brien guarantees such an algorithm exists

W̶e̶ ̶s̶h̶o̶u̶l̶d̶ ̶b̶e̶ ̶m̶o̶r̶e̶ ̶n̶u̶a̶n̶c̶e̶d̶.̶ ̶T̶h̶e̶r̶e̶ ̶i̶s̶ ̶̶n̶o̶̶ ̶s̶u̶c̶h̶ ̶a̶l̶g̶o̶r̶i̶t̶h̶m̶ ̶t̶h̶a̶t̶ ̶a̶l̶l̶o̶w̶s̶ ̶y̶o̶u̶ ̶t̶o̶ ̶g̶e̶n̶e̶r̶a̶t̶e̶ ̶a̶ $min\lbrace 0,0.5 \rbrace$ c̶o̶i̶n̶.̶ ̶W̶h̶a̶t̶ ̶K̶e̶a̶n̶e̶ ̶a̶n̶d̶ ̶O̶'̶B̶r̶i̶e̶n̶ ̶s̶t̶a̶t̶s̶ ̶i̶s̶ ̶t̶h̶a̶t̶ ̶y̶o̶u̶ ̶c̶a̶n̶ ̶̶a̶p̶p̶r̶o̶a̶c̶h̶̶ ̶t̶h̶e̶ ̶d̶e̶s̶i̶r̶e̶d̶ ̶f̶u̶n̶c̶t̶i̶o̶n̶ ̶a̶s̶ ̶c̶l̶o̶s̶e̶ ̶a̶s̶ ̶y̶o̶u̶ ̶l̶i̶k̶e̶.̶

Intuitive example

We can generalize the method from Luis Mendo to get a procedure which allows you to generate a Bernoulli variable with probability $f(p)$ given a Bernoulli variable with probability $p$.

The steps are as following:

  • Use $m$ coin flips to estimate $\hat p$, the bias of the coin.

  • Based on the bias of the coin, use $n$ unbiased coin flips (for which you can use John von Neumann's algorithm) to generate a coin with approximately $\hat{f}(\hat{p})$ probability.

    In the computational example below the $\hat{f}(\hat{p})$ is equal to one of the quantiles of a Bernoulli variable with size $n$ and probability $0.5$

By increasing $m$ you can make $\hat p$ get closer to the true value $p$. By increasing $n$ you can make $\hat{f}(\hat p)$, get closer to ${f}(\hat p)$. In the end we can make $\hat{f}(\hat p)$ as close to ${f}(p)$ as we like by increasing $m$ and $n$.

The code below demonstrates this construction for some function $f(p) = \frac{1}{2} + \frac{1}{2} \sin(4\pi p)$ with a simulation.

We see that with $m = n = 1000$ you get to the following approximation:

enter image description here

### estimate p by **tossing m coins** and use average number of 1's
estimate_p <- function(m,p) {

### create a coin toss with probability p **based on n tosses** with a fair coin
toss_p <- function(n,p,k=10^3) {
  rng <- 0:n
  cutoff <- qbinom(p,n,0.5) 
  #sum(rbinom(k,n,0.5) <= cutoff)/k

### some function that we want to convert the coin to 
conv_p <- function(p) {

ps <- seq(0.001,0.999,0.001)
plot(ps, conv_p(ps), type = "l",
     xlab = "p coin", ylab = "p process")

### estimate f(p) with coin tosses of coin p
m <- 1000 ### number of coin tosses to estimate p
n <- 1000 ### number of fair coin tosses to approach f(p)

k_trials = 1:10^2
p_trials <- seq(0.01,0.99,0.01)

pt = 0.1
k_trials = 1:100

for (pt in p_trials) {
  for (k in k_trials) {
    p_est <- estimate_p(m,conv_p(pt))
    p_out <- toss_p(n,p_est)
    points(pt,p_out, pch = 21 , col = rgb(0,0,0,0.1), bg = rgb(0,0,0,0.1), cex = 0.3)

Trick for a more efficient example

The above example is simple and helps to see intuititively how an algorithm can approach as close as we like. However, we could try to see how to make a faster method.

The trick is that we can tabulate the results like

toss     probability
HH       p²
HT       p(1-p)
TH       (1-p)p
TT       (1-p)²

then select a few of the tosses that create a desirable ratio. For instance in the case of John von Neumann's trick (generating a fair coin out of a biased coin) we could decide

  • 'heads' if we observe HT with the unfair coin

  • 'tails' if we observe TH with the unfair coin

  • 'toss again' if we did not observe either.

    (Yes, this approach throws away results which is not efficient. There's a lot of ways to improve, and complicate, the strategies to solve these kind of problems).

Then the ratio's of heads and tails is $\frac{p(1-p)}{(1-p)p} = 1$ so you get the fair coin independent of $p$.

Generalizing the trick

What we will try to do now is find out the best combination of events for tails and heads given $n$ coin flips to get a ratio like

$$\frac{a_1 p^n + a_2 p^{n-1}(1-p) + \dots + a_n (1-p)^n }{b_1 p^n + b_2 p^{n-1}(1-p) + \dots + b_n (1-p)^n}$$

that approximates our desired function the best. (the coefficients $a_i$ and $b_i$ will need to be integers and are bounded by binomial coefficients)

The r-code below computes the optimal for 3 coin tosses (I am computing just 3 tosses because the least-squares problem with integer constraints is for the moment computed in a dumb way)

We will be selecting 'heads' for HHT and HTT and 'tails' for HTH and TTT which will give an odds ratio of



### function in numerator and denominator
simulated <- function(x,par){
  k <- length(par)-1
  xp <- x^(k:0)*(1-x)^(0:k)
simulated <- Vectorize(simulated, vectorize.args = "x")

### function to compare estimate of f(p) with given desired f(p)
fn <- function(x) {
  l <- length(x) ### number of pars
  a <- x[1:(l/2)]
  b <- x[(l/2+1):l]
  ### integrate difference in first half
  int1 <- integrate(f = function(x) {(simulated(x,a)/(simulated(x,a)+simulated(x,b))-x)^2},
                    lower = 0.001, upper = 0.5, stop.on.error = F)
  int2 <- integrate(f = function(x) {(simulated(x,a)/(simulated(x,a)+simulated(x,b))-0.5)^2},
                    lower = 0.5, upper = 0.999, stop.on.error = F)
  if (int1$message != "OK") {
    int1$value = 10^6
  if (int2$message != "OK") {
    int2$value = 10^6

### compute best option for 3 tosses
### by just trying every option

RSS = 10^6
solution = c(0,0,0,0,0,0,0,0)
for (a1 in 0:1) {
  for (a2 in 0:3) {
    for (a3 in 0:3) { 
      for (a4 in 0:1) {
        for (b1 in 0:(1-a1)) {
          for (b2 in 0:(3-a2)) {
            for (b3 in 0:(3-a3)) {
              for (b4 in 0:(1-a4)) {
                if ((a1+a2+a3+a4+b1+b2+b3+b4)>0) {
                  test_RSS <- fn(c(a1,a2,a3,a4,b1,b2,b3,b4))
                  if (test_RSS < RSS) {
                    RSS = test_RSS
                    solution <- c(a1,a2,a3,a4,b1,b2,b3,b4)
                a <- c(a1,a2,a3,a4)
                b <- c(b1,b2,b3,b4)
                #lines(xs, simulated(xs,a)/(simulated(xs,a)+simulated(xs,b)),col = rgb(0,0,0,0.1))

l <- length(solution)
a <- solution[1:(l/2)]
b <- solution[(l/2+1):l]

xs <- seq(0.001,0.999,0.001)
plot(-10,-10, xlim = c(0,1), ylim = c(0,1), 
     xlab = "p coin", ylab = "p simulated")
lines(xs, simulated(xs,a)/(simulated(xs,a)+simulated(xs,b)))
lines(c(0,0.5,1), c(0,0.5,0.5), col = 2)
  • 2
    $\begingroup$ Just to make sure, if I am not mistaken, Keane and O'Brien actually says we can make exactly the new coin (not approximate) through a random number of draws. Meanwhile, Nacu and Peres says, Though the number of required draws is a random variable, we may get controls for its tail given different smoothness assumptions. Please correct me if I am wrong. $\endgroup$
    – Bravo
    Nov 16, 2020 at 17:45
  • 1
    $\begingroup$ @Bravo I do not see it stated in their article that you can get exactly the new coin. The proof also makes use of the fact that the function is polynomial bounded, and we can get within any of the bounds. But that does not mean that there is an algorithm that can get exactly the new coin for some given finite number of flips or steps. $\endgroup$ Nov 23, 2020 at 9:58
  • $\begingroup$ If I understand it correctly, the theorem in Keane and O'Brien's paper does say that theres is an (exact) Bernoulli factory algorithm for any continuous function that is polynomially bounded away from 0 and 1. The OP's function satisfies those conditions, so there is one such algorithm. The number of observations it requires, according to the referred paper, is finite with probability 1 $\endgroup$
    – Luis Mendo
    Nov 26, 2020 at 17:03
  • $\begingroup$ @LuisMendo I find their theorem not so easy to understand. They speak about 'simulable' if some random variable defined by a function of the coin flips and an auxiliary variable $\varphi(a,x_1,x_2,\dots)$ follows the distribution of the $Q$ that is to be simulated. But does that mean that it is exactly equal or that it is as close as we like? I will have to dig again through their proof I guess. I had the idea that it was a proof based on bounds/limits and not on being exactly equal. $\endgroup$ Nov 26, 2020 at 17:31
  • $\begingroup$ @SextusEmpiricus To me, "follows the distribution of the Q that is to be simulated" means the exact distribution, not an arbitrarily close approximation thereof. Section 1 of Nacu and Peres' paper is clearer, I think $\endgroup$
    – Luis Mendo
    Nov 26, 2020 at 17:35

Since the function min(λ, c) meets the Keane–O'Brien theorem, this implies that there are polynomials that converge from above and below to the function. This is discussed, for example, in Thomas and Blanchet (2012) and in Łatuszyński et al. (2009/2011). However, neither paper shows how to automate the task of finding polynomials required for the method to work, so that finding such a sequence for an arbitrary function (that satisfies the Keane–O'Brien theorem) remains far from trivial (and the question interests me to some extent, too).

But fortunately, there is an alternative way to simulate min(λ, 1/2) without having to build a sequence of polynomials explicitly. This algorithm I found is given below, and I have a page describing the derivation of this algorithm.

  1. With probability 1/2, flip the input coin and return the result.
  2. (Random walk.) Generate unbiased random bits until more zeros than ones are generated this way for the first time. Then set m to (n −1)/2+1, where n is the number of bits generated this way.
  3. (Build a degree-m*2 polynomial equivalent to (4*λ*(1− λ ))m/2.) Let z be (4m/2)/choose(m*2, m). Define a polynomial of degree m*2 whose (m*2)+1 Bernstein coefficients are all zero except the mth coefficient (starting at 0), whose value is z. Elevate the degree of this polynomial enough times so that all its coefficients are 1 or less (degree elevation increases the polynomial's degree without changing its shape or position; see the derivation in the appendix). Let d be the new polynomial's degree.
  4. (Simulate the polynomial, whose degree is d (Goyal and Sigman 2012).) Flip the input coin d times and set h to the number of ones generated this way. Let a be the hth Bernstein coefficient (starting at 0) of the new polynomial. With probability a, return 1. Otherwise, return 0.

I suspected that the required degree d would be floor(m*2/3)+1. With help from the MathOverflow community, steps 3 and 4 of the algorithm can be described more efficiently as follows:

  • (3.) Let r be floor(m*2/3)+1, and let d be m*2+r.
  • (4.) (Simulate the polynomial, whose degree is d.) Flip the input coin d times and set h to the number of ones generated this way. Let a be (1/2) * 2m*2*choose(r, h_−_m)/choose(d, h) (the polynomial's hth Bernstein coefficient starting at 0; the first term is 1/2 because the polynomial being simulated has the value 1/2 at the point 1/2). With probability a, return 1. Otherwise, return 0. (Here, choose(n, k) is a binomial coefficient.)

In addition, there is an approximate way to sample min(λ, c) and most other continuous functions f that map (0, 1) to (0, 1). Specifically, it's trivial to simulate an individual polynomial with Bernstein coefficients in [0, 1], even if the polynomial has high degree and follows the desired function closely (Goyal and Sigman 2012):

  • Flip the input coin n times (where n is the polynomial's degree), and let j be the number of ones.
  • With probability a[j], that is, the j-th control point, starting at 0, for the polynomial's corresponding Bézier curve, return 1. Otherwise, return 0.

To use this algorithm, simply calculate a[j] = f(j/n), where n is the desired degree of the polynomial (such as 100). Each a[j] is one of the Bernstein coefficients of a polynomial that closely approximates the function; the higher n is, the better the approximation.


Let me clarify two things: Generating fair bits from a biased coin, and simulation vs. Estimation.

  1. First, generating fair bits. You can generate unbiased bits either by generating them separately from the coin, or by using biased coin tosses and applying a randomness extraction procedure to turn them into unbiased bits. Ways to do so include not just the von Neumann algorithm itself, but also randomness extractors that assume no knowledge of the coin's bias, including Yuval Peres's (1992) iterated von Neumann extractor as well as an "extractor tree" by Zhou and Bruck (2012). See also my note on randomness extraction.
  2. Second, the difference between simulating and estimating probabilities. Essentially, "simulation" means generating the same distribution, and "estimation" means generating the same expected value (Glynn 2016). However, a Bernoulli factory for simulating f(p) acts as an unbiased estimator for f(p) (Łatuszyński et al. 2009/2011). But a function that doesn't meet the Keane–O'Brien theorem, such as min(2 p, 1 − (2 p)), can't be simulated by any algorithm without further knowledge of p, because the estimate will not be unbiased (Łatuszyński et al. 2009/2011). (However, it is possible to simulate min(2 p, 1 − (2 p), 1−ε) this way.) See also my note.


  • Goyal, V. And Sigman, K., 2012. On simulating a class of Bernstein polynomials. ACM Transactions on Modeling and Computer Simulation (TOMACS), 22(2), pp.1-5.
  • Łatuszyński, K., Kosmidis, I., Papaspiliopoulos, O., Roberts, G.O., "Simulating events of unknown probabilities via reverse time martingales", arXiv:0907.4018v2 [stat.CO], 2009/2011.
  • Thomas, A.C., Blanchet, J., "A Practical Implementation of the Bernoulli Factory", arXiv:1106.2508v3 [stat.AP], 2012.
  • Glynn, P.W., "Exact simulation vs exact estimation", Proceedings of the 2016 Winter Simulation Conference, 2016.
  • Zhou, H. And Bruck, J., "Streaming algorithms for optimal generation of random bits", arXiv:1209.0730 [cs.IT], 2012.
  • Peres, Y., "Iterating von Neumann's procedure for extracting random bits", Annals of Statistics 1992,20,1, p. 590-597.

I don't know how to generate that function, but there is a simple way to generate increasingly good approximations of it by consuming an increasingly large of inputs. Namely, you observe $N$ inputs, where $N$ is an odd integer. Then:

  • If the majority of those inputs is $0$ you "decide" that $p<1/2$. So you take an addtional input an output that.
  • If the majority of he $N$ inputs is $1$ you "decide" that $p>1/2$. So you output a Bernoulli auxiliary random variable with parameter $1/2$ (if needed you can generate that variable from additional inputs using the well-known von Neumann procedure).

Here is some Matlab code that plots the resulting $f_N(p)$ as a function of $p$ for several values of $N$. It's easy to prove that $f_N(p)\rightarrow f(p)$ as $N \rightarrow \infty$.

N = 101;
p_axis = 0:.01:1;
y = 1/2 + (p_axis-1/2).*binocdf((N-1)/2, N, p_axis);
plot(p_axis, y)

enter image description here

As a check, here's an experiment to estimate $f_{101}(p)$ using $10^5$ realizations for each value of $p$. The markers represent the proportion of $1$ in the output for each set of $10^5$ realizations.

N = 101;
R = 1e5;
p_test = 0:.05:1;
result = NaN(size(p_test));
for k = 1:numel(p_test)
    p = p_test(k);
    t = mean(rand(N,R)<=p, 1)<=1/2;
    out_k = NaN(1,R);
    out_k(t) = rand(1,sum(t))<=p;
    out_k(~t) = rand(1,sum(~t))<=1/2;
    result(k) = mean(out_k);
plot(p_test, result, 'o')

enter image description here

  • $\begingroup$ thanks! It seems that it is intuitive to find a `consistent' coin which has bias going to zero, but an unbiased coin seems difficult $\endgroup$
    – Bravo
    Nov 10, 2020 at 4:58
  • $\begingroup$ Your results together with the method of Thomas and Blanchet 2012 suggest that an exact simulation method for this kind of factory function is possible, as long as we have an automated way to compute monotonic sequences of Bernstein polynomials that converge from above and below to that function. Unfortunately, the Thomas and Blanchet paper is not systematic enough in the problem of finding such polynomials, and it leaves out a number of details needed for a programmer to check their implementation's correctness. $\endgroup$
    – Peter O.
    Nov 12, 2020 at 4:23

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.