Generating a min{p, 0.5} coin from a p-coin - Bernoulli factory type problem

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!

• 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.$ – whuber Oct 27 at 16:02
• 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 – Bravo Oct 27 at 20:23
• 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. – whuber Oct 27 at 20:28
• good point, just edited! – Bravo Oct 27 at 21:46

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:

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

### 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
pbinom(cutoff,n,0.5)
}

### some function that we want to convert the coin to
conv_p <- function(p) {
0.5+0.5*sin(p*pi*4)
}

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
set.seed(1)
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

$$\frac{p^2(1-p)+p(1-p)^2}{p^2(1-p)+(1-p)^3}=\frac{p}{1-2p+2p^2}$$

### function in numerator and denominator
simulated <- function(x,par){
k <- length(par)-1
xp <- x^(k:0)*(1-x)^(0:k)
return(sum(par*xp))
}
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") { int1value = 10^6 } if (int2$$message != "OK") {
int2$$value = 10^6 } return(int1$$value+int2$value) } ### compute best option for 3 tosses ### by just trying every option choose(3,0:3) 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)  • 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. – Bravo Nov 16 at 17:45 • @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. – Sextus Empiricus Nov 23 at 9:58 • 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 – Luis Mendo 2 days ago • @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. – Sextus Empiricus 2 days ago
• @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 – Luis Mendo 2 days ago

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).

The following code in Python simulates a function of the form min(λ, c) using the reverse-time martingale framework (Algorithm 4) of Łatuszyński et al. 2009/2011. (I also tried using the Thomas and Blanchet method, but all my attempts at implementing it doesn't lead to the correct results.) Note that currently, this code is inefficient, largely because it has to calculate an expectation that requires a growing amount of work to calculate as the polynomial's degree gets large. Note also that usually the algorithm terminates quickly but occasionally the algorithm can take a very long time, especially if g0 is close to 0. This is in part because the function is not differentiable at c, resulting in a low rate of convergence.

import math
import random
import statistics
from fractions import Fraction

def binco(n, k):
# Binomial coefficient
ret=1
for i in range(n-k+1, n+1):
ret*=Fraction(i,(n-i+1))
return int(ret)

def hypergeom(f, m, n, h):
# Hypergeometric expectation
ret=0
nh=binco(n, h)
for i in range(h+1):
fr=Fraction(f(m,i)*binco(n-m,h-i)*binco(m,i))
ret+=fr/nh
return ret

def funcabove(n, k, c):
# Calculates the upper nth-degree Bernstein
# coefficient for min(k/n, c).
# A Cross Validated answer https://math.stackexchange.com/a/3889525
# contains references on convergence using polynomials
ret=min(k/n, c)
kk=(4306+837*math.sqrt(6))/5832
above=kk*n**(-0/2)*1*(n**(-1/2))
return Fraction(ret+above)

def funcbelow(n, k, c):
# Calculates the lower nth-degree Bernstein
# coefficient for min(k/n, c)
return Fraction(min(k/n, c))

def simulate(p, c=0.5): # Simulates min(p, c)
l=0
lt=0
u=1
ut=1
ones=0
count=0
prevcount=0
g0=random.random()
fb=lambda a, b: funcbelow(a, b, c)
fa=lambda a, b: funcabove(a, b, c)
i=0
while True:
olddegree=0 if i==0 else 1<<(i-1)
degree=1<<i
for j in range(prevcount, degree):
ones+=(1 if random.random()<p else 0)
if ones>degree: raise ValueError
l=funcbelow(degree, ones, c)
u=funcabove(degree, ones, c)
if olddegree>0:
ls=hypergeom(fb,olddegree,degree,ones)
us=hypergeom(fa,olddegree,degree,ones)
else:
ls=0
us=1
mult=(ut-lt)/(us-ls)
ltnew=lt+(l-ls)*mult
utnew=ut-(us-u)*mult
lt=ltnew
ut=utnew
if g0<=lt: return 1
if g0>ut: return 0
prevcount=degree
i+=1
return 0

# Mean should be close to 0.2
print(statistics.mean(simulate(0.2) for i in range(2000)))


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.

EDIT:

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). An algorithm that finds an unbiased estimate of the coin's bias p, then calculates f(p), then generates a uniform(0,1) random number, then compares that number with f(p) will simulate that probability in theory (but will not be exact in practice because it can cause errors due to the use of fixed precision, such as rounding and cancellations, and thus might not be considered a Bernoulli factory). However, 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, including the algorithm just given, 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.

REFERENCES:

• 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 an 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)


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);
end
plot(p_test, result, 'o')


• thanks! It seems that it is intuitive to find a `consistent' coin which has bias going to zero, but an unbiased coin seems difficult – Bravo Nov 10 at 4:58
• 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. – Peter O. Nov 12 at 4:23