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!
 A: 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.

*

*With probability 1/2, flip the input coin and return the result.

*(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.

*(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.

*(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.

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

*

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

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


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.

A: 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)


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


