Bernoulli Factory Algorithms

[peteroupc]

Issue opened to seek comments or suggestions related to my page on Bernoulli Factory algorithms, which are algorithms to turn coins biased one way into coins biased another way.

https://peteroupc.github.io/bernoulli.html

You can send comments on this document on this issues page. You are welcome to suggest additional Bernoulli factory algorithms, especially specific continued fraction expansions and series expansions for the general martingale and Mendo algorithms.

[peteroupc]

I have made numerous updates to this document recently, with more algorithms and discussion.

Letting them know: @lmendo.

[peteroupc]

I have added the following open questions to the Bernoulli Factory Algorithms page. Any answers to these questions will greatly improve the article.

• Is there a simple Bernoulli factory algorithm that can simulate the probability equal to Euler's constant γ, without relying on floating-point arithmetic? This repeats an open question given in "On Buffon Machines and Numbers". (Solved.)

• See the open questions found in the section "Probabilities Arising from Certain Permutations" in the appendix.

• I request expressions of mathematical functions that can be expressed in any of the following ways:

  • Series expansions for continuous functions that equal 0 or 1 at the points 0 and 1. These are required for Mendo's algorithm for certain power series.
  • Series expansions for alternating power series whose coefficients are all in the interval [0, 1] and form a nonincreasing sequence. This is required for another class of power series.
  • Upper and lower bound approximations that converge to a given constant or function. These upper and lower bounds must be nonincreasing or nondecreasing, respectively.
  • Simple continued fractions that express useful constants.

  All these expressions should not rely on floating-point arithmetic or the direct use of irrational constants (such as π or sqrt(2)), but may rely on rational arithmetic. For example, a series expansion that directly contains the constant π is not desired; however, a series expansion that converges to a fraction of π is.

In particular, I have found it surprisingly hard to find information and research on alternating power series and lower/upper bound approximations needed to apply certain Bernoulli factory algorithms, including Algorithm 3 of "Simulating Probabilities of Unknown Events via Reverse Time Martingales". Especially helpful would be information on how to transform an arbitrary function into an alternating series (preferably one that is easy to calculate) or another construct that satisfies these algorithms.

Letting these people know: @63EA13D5 .

[lmendo]

Regarding the simulation of Euler's constant without floating-point arithmetic, I think I found a simple way to do that, as well as for other constants. Please see the attached paper (I just submitted it to arXiv, but it doesn't appear yet) arXiv

[peteroupc]

Thank you for your response. In the meantime, I found a rational interval arithmetic described in Daumas, M., Lester, D., Muñoz, C., "Verified Real Number Calculations: A Library for Interval Arithmetic", arXiv:0708.3721 [cs.MS], 2007, which computes upper and lower bounds of common functions, and I implemented this arithmetic in Python. However, as I found, even this rational interval arithmetic doesn't always lend itself to simple Bernoulli factory algorithms (and due to Python's highly slow Fraction class, I had to optimize some parts of the code), which is why I have avoided involving it in the hopes that an even simpler algorithm is possible without resorting to this kind of arithmetic.

Also, on your problem of "characteriz[ing] the minimum number of inputs that an arbitrary algorithm must use to transform 1/2-coins into a τ-coin", I am aware of the following paper by D. Kozen that relates to this: "Optimal Coin Flipping". By the results in that paper, for example, any algorithm of the kind relevant here will require at least 2 unbiased coin flips on average to simulate a coin of known bias τ, except when τ is a dyadic rational.

[lmendo]

Thanks. The paper you link seems very interesting; I'll take a look. If the lower bound is 2, then my algorithm is close to optimum.

In my algorithm, and probably in others, computations can be reused if you are going to simulate the simulation several times (page 11 of my paper). This can increase speed.

[peteroupc]

I have noticed one error in your paper: In (16), "(2j-1)2j(2j+1)" should read "2(j-1)(2(j-1)+1)(2(j-1)+2)". Otherwise, the algorithm would simulate the wrong probability.

Also, for clarity, the "1/2" and "3/2" in the algorithm's pseudocode should be in parentheses.

---

Also, in the meantime, I have another open question:

• Is there a simple and reasonably efficient algorithm that can simulate the probability z = choose(n, k)/(2*n), for n up to, say, 2³⁰, without relying on floating-point arithmetic or precalculations? This is trivial for relatively small n, but complicated for larger n even though z is a rational number; see also an appendix in Bringmann et al., 2014, which ultimately relies on floating-point arithmetic. One idea is to use approximations of ln(z), ln(choose(n, k)), or ln(n!) that converge from above and below to the true value. Then an exponential random number could be compared with ln(z), which is the same as comparing z itself with a uniform random number. This is what I have implemented (see interval.py and betadist.py), but the whole procedure was far from trivial. Is there something simpler?

This now leads to another question, which I don't yet classify as an open question, since I might check whether the algorithm I just learned today would be helpful: It would be nice to have a simple algorithm that simulates a probability of the form exp(-y), where y is a logarithm of a really huge number (or a sum of such logarithms), but whose integer part is not known in advance; one case of this was just given above. (Note that my page currently includes an algorithm for simulating exp(-z), but assumes z's integer part is known.) (I think I have an algorithm developed.)

REFERENCE:

K. Bringmann, F. Kuhn, et al., “Internal DLA: Efficient Simulation of a Physical Growth Model.” In: Proc. 41st International Colloquium on Automata, Languages, and Programming (ICALP'14), 2014.

[lmendo]

Thanks for the correction. In my code I am using equation (15), so I hadn't noticed the mistake.

As for the other problems, I don't really know, sorry

[peteroupc]

The requests and open questions for all my articles are now in a dedicated page: Requests and Open Questions.

However, since this GitHub issue is about all aspects of Bernoulli factories, not just the ones mentioned in Requests and Open Questions, I will leave this issue open for a while.

[peteroupc]

I want to draw attention to my Supplemental Notes for Bernoulli Factory algorithms:

• https://peteroupc.github.io/bernsupp.html

Especially the section "Approximation Schemes". It covers ways to build polynomial approximations** to the vast majority of functions for which the Bernoulli factory problem can be solved (also called factory functions). They include concave, convex, and twice differentiable functions. Now my goal is to find faster and more practical polynomial approximation schemes for these functions. Which is why I have the following open questions:

• Are there factory functions used in practice that are not covered by the approximation schemes in this section?
• Are there specific functions (especially those in practical use) for which there are practical and faster formulas for building polynomials that converge to those functions (besides those I list in this section or the main Bernoulli Factory Algorithms article)? One example is the scheme for min(2λ, 1−ε) found in Nacu and Peres, "Fast simulation of new coins from old", 2005.

** See my Math Stack Exchange question for a formal statement of the kind of approximations that I mean.

[peteroupc]

I want to draw attention to my Supplemental Notes for Bernoulli Factory algorithms:

* https://peteroupc.github.io/bernsupp.html

Especially the section "Approximation Schemes". It covers ways to build polynomial approximations** to the vast majority of functions for which the Bernoulli factory problem can be solved (also called factory functions).

This was part of my previous comment.

However, these polynomial approximation methods don't ensure a finite expected running time in general.

I suspect that a running time with finite expectation isn't possible in general unless the residual probabilities formed by the polynomials are of order O(1/n^(2+epsilon)) for positive epsilon, which Holtz et al. ("New coins from old, smoothly", Constructive Approximation, 2011) proved is possible only if the function to be simulated is C2 continuous. I suspect this is so because sums of residual probabilities of the form O(1/n^2) don't converge, but such sums do converge when the order is O(1/n^(2+epsilon)). (By residual probabilities, I mean the probability P(N > n) given in Result 3, condition (v) of that paper.)

Thus I have several questions:

1. Given a C2 continuous function f(x), are there practical algorithms for building polynomials that converge to that function in a manner needed for the Bernoulli factory problem, where the expected number of p-coin flips is finite? What about C^k functions for k>2?
2. The Holtz et al. paper seems to answer the previous question, except I find it impractical since it has no examples and uses constants such as "s", "θ_α", and "D" that are hard to determine and have no easy lower bounds. Thus another question is: What are practical lower bounds for those constants? Can you provide (programmer-friendly) examples of the method (in terms of Bernstein polynomial approximations) for various functions such as Lipschitz continuous, C2 continuous, concave, convex, real analytic, C-infinity, etc.?
3. Is it true that if f is a Lipschitz continuous function (or in the alternative, a C2 continuous function), then it can be simulated with a finite expected number of flips?

Other interesting questions:

1. Mossel and Peres ("New coins from old: computing with unknown bias", Combinatorica, 2005) investigated which functions can be simulated by pushdown automata (state machines with a stack) when those machines are given an infinite "tape" of flips of a biased coin. They showed that only algebraic functions can be simulated this way. I would like to better characterize this class of algebraic functions, and have written an article appendix showing my progress. Are there other results on this question?
2. I have written an article appendix on the question of which functions are strongly simulable; they admit a Bernoulli factory with nothing but the biased coin as randomness. Of course, if the domain includes neither 0 or 1, then that's all functions that admit a Bernoulli factory. The difficulty is to characterize which functions are strongly simulable even if the domain includes 0 and/or 1, and the appendix shows my progress. My question is: Are the conditions in Proposition 2 necessary and sufficient for being strongly simulable?

[peteroupc]

The following is an updated version of most of my requests and open questions on the Bernoulli factory problem:

1. What simulations exist that are "relatively simple" and succeed with an irrational probability between 0 and 1? What about "relatively simple" Bernoulli factory algorithms for factory functions? Here, "relatively simple" means that the algorithm:

   • Should use only uniform random integers (or bits) and integer arithmetic.
   • Does not use floating-point arithmetic or make direct use of square root or transcendental functions.
   • Does not calculate base-n expansions directly.
   • Should not use rational arithmetic or increasingly complex approximations, except as a last resort.

   See also Flajolet et al. (2010)[^2]. There are many ways to describe the irrational probability or factory function. I seek references to papers or books that describe irrational constants or factory functions in any of the following ways:

   • For irrational constants:
     • Simple continued fraction expansions.
     • Closed shapes inside the unit square whose area is an irrational number. (Includes algorithms that tell whether a box lies inside, outside, or partly inside or outside the shape.) Example.
     • Generate a uniform (x, y) point inside a closed shape, then return 1 with probability x. For what shapes is the expected value of x an irrational number? Example.
     • Functions that map [0, 1] to [0, 1] whose integral (area under curve) is an irrational number.
   • For Bernoulli factory functions:
     • Functions with any of the following series expansions, using rational arithmetic only:
       • Power series where f(0) is 0 and f(1) is rational or vice versa (see "Certain Power Series").
       • Series with non-negative terms that can be "tucked" under a discrete probability mass function (see "Convex Combinations").
       • Alternating power series whose coefficients are all in the interval [0, 1] and form a nonincreasing sequence (see "Certain Power Series").
       • Series with non-negative terms and bounds on the truncation error (see "Certain Converging Series").
     • A way to compute two sequences of polynomials written in Bernstein form that converge from above and below to a factory function as follows: (a) Each sequence's polynomials must have coefficients lying in [0, 1], and be of increasing degree; (b) the degree-n polynomials' coefficients must lie at or "inside" those of the previous upper polynomial and the previous lower one (once the polynomials are elevated to degree n). For a formal statement of these polynomials, see my question on MathOverflow.
       
       The supplemental notes include formulas for computing these polynomials for large classes of factory functions, but none of them ensure a finite expected number of coin flips in general, and it is suspected that a finite number of flips isn't possible unless the factory function is C² continuous (has two or more continuous "slope" functions). Thus one question is: Given a C² continuous factory function, are there practical algorithms for building polynomials described here, where the expected number of coin flips is finite (besides the algorithms in this article or the supplemental notes)? Examples worth thinking about are sin(π*λ)/2 and sin(2*π*λ)/4 + 1/2.

2. Is there a simpler or faster way to implement the base-2 or natural logarithm of binomial coefficients? See the example in the section "Certain Converging Series".

3. A pushdown automaton is a state machine that holds a stack of symbols. Mossel and Peres (2005)[^4] investigated which functions (f(λ)) can be simulated by these machines when they're given an infinite "tape" of flips of a coin that shows heads with probability λ. They showed that pushdown automata can simulate only algebraic functions, but perhaps not all of them. The question is: What is the exact class of algebraic functions a pushdown automaton can simulate? Can it simulate the functions min(λ, 1−λ) and λ^1/p where p>2 is a prime number? I have written an article appendix showing my progress, but are there other results on this question?

4. The following is an open question in Nacu and Peres 2005. Let J be a closed interval on (0, 1), such as [1/100, 99/100], and let f(λ) be a function that admits a Bernoulli factory. Suppose there is an algorithm that takes a coin with unknown probability of heads λ and produces one or more samples of the probability f(λ). When the probability λ can be any value in J, is it possible for this algorithm to have an expected number of input coin flips per sample that is arbitrarily close to the so-called entropy bound? The entropy bound is h(f(λ))/h(λ) where h(x) = −x*ln(x)−(1−x)*ln(1−x) is related to the Shannon entropy function. Does the answer change if the algorithm can also use a separate source of unbiased random bits? See my section "Multiple-Output Bernoulli Factory".

5. A factory function f(λ) is strongly simulable if there is an algorithm to toss heads with probability f(λ) using only a coin that shows heads with probability λ and no other randomness. Keane and O'Brien (1994) showed already that f(λ) is strongly simulable if neither 0 nor 1 is in f's domain. It's also easy to show that if f is strongly simulable, then f(0) and f(1) must each be 0, 1, or undefined. Is this a sufficient condition to be strongly simulable? I have written an article appendix showing my progress, but are there other results on this question?

[^1]: Keane, M. S., and O'Brien, G. L., "A Bernoulli factory", ACM Transactions on Modeling and Computer Simulation 4(2), 1994.

[^2]: Flajolet, P., Pelletier, M., Soria, M., "On Buffon machines and numbers", arXiv:0906.5560 [math.PR], 2010.

[^3]: Ł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.

[^4]: Mossel, Elchanan, and Yuval Peres. New coins from old: computing with unknown bias. Combinatorica, 25(6), pp.707-724.

[peteroupc]

I have collected all my questions on the Bernoulli factory problem in one place:

• Open Questions on the Bernoulli Factory Problem.