What are relatively simple simulations that succeed with an irrational probability?

Question

What are relatively simple simulations that succeed with an irrational probability?

Let me break down this question.

Relatively simple simulations

By "relatively simple" I mean only algorithms that involve only—

• unbiased random bits (and thus discrete uniform variates), and
• integer arithmetic.

The algorithm may:

• Use Flajolet's "geometric bag" technique to build up continuous uniform variates and to sample a probability equal to a uniform random variate.
• Use rational arithmetic, but only as a last resort.
• Generate Poisson random variates, since there is a way to do so using only integer arithmetic and/or coins with unknown success probability.

The algorithm may not:

• Use "multiprecision interval arithmetics or [...] functions of high complexity as primitives" (Flajolet et al. 2010).
• Make direct use of square root or transcendental functions or constants.
• Rely on "sequences of approximation functions of increasing complexity" (Flajolet et al. 2010), except as a last resort.

Succeed with an irrational probability

• The algorithm's expected value must equal a constant, namely an irrational number in (0, 1).
• If the algorithm outputs $X$, then it must be that $\mathbb{E}[X]$ is the desired constant.
• The desired irrational constant is allowed to depend on one or more parameters.

There are two ways the algorithm could "succeed with an irrational probability".

• It could return either 0 or 1 with probability equal to the constant.
• Or it could produce a value $r\ge 1$; then the algorithm could return a Bernoulli($1/r$) variate.

I ask for irrational probabilities because the rational case is quite trivial; if the rational probability is known there is no need for algorithms like the one I'm asking for. Rather, it's enough to compare a uniform($d$) random variate with $n$, where $d$ and $n$ are the denominator and numerator.

Examples

• 1 divided by the golden ratio.
• 1/π.
• $e^{-1}$.
• The Euler–Mascheroni constant (even though it's not known to be irrational).
• The frog problem with negative steps. Here, the expected number of iterations turns out to be an irrational number in a far from obvious way, which is desirable for this question. (Indeed, the more implicit the algorithm's use of the irrational constant, the better.)

Remarks

1. The algorithm or method should have been published in a research article or book, and should be different from the ones I already have in my pages on Bernoulli factory algorithms or more arbitrary-precision algorithms, especially the section "Algorithms for Specific Constants". My literature searches have shown that algorithms I'm seeking are quite hard to find.
2. My interest is not to compute irrational constants to high precision, but rather to have the algorithm serve as a "black box", a coin that "shows heads" with an irrational probability and in an exact manner; e.g., to serve as input coins to other so-called "Bernoulli Factory" algorithms that turn a Bernoulli($\lambda$) probability into an exact Bernoulli($f(\lambda)$) probability.

References

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

Answer 1

I think the details may be somewhat different for each 'desired irrational constant', but here is a strategy that may work for many such constants.

Here is a simple algorithm that estimates $\pi/4,$ the area in the unit square beneath the quarter unit-circle with center at the origin.

set.seed(2021)
x = runif(10^6); y = runif(10^6)
mean(y <= sqrt(1-x^2))
[1] 0.785459
pi/4
[1] 0.7853982

Can you find a function in a suitable square or rectangle that bounds an area equal to each of your desired irrational constants (or a simple function thereof)?

Figure with only 50,000 points in the unit square for clarity.

B = 50000; x = runif(B);  y = runif(B)
plot(x, y, pch=".")
blue = (y <= sqrt(1-x^2))
points(x[blue], y[blue], col="blue", pch=".")

Answer 2

There is a universal algorithm. It doesn't matter whether the probability is irrational or not.

It suffices to implement a procedure to output either $0$ or $1$ that will (a) almost surely terminate and (b) output $1$ with a probability $\phi,$ where $\phi$ is any number in the interval $[0,1]$ (rational or irrational). The following description relies on an arbitrarily long sequence of iid Bernoulli$(1/2)$ variables $X_1,X_2,X_3,\ldots.$

Procedure f(phi):
for i in 1, 2, 3, ...
    if phi >= 1/2 then 
        if X[i]==0 return(1) else return(f(2*phi-1))
    else 
        if X[i]==0 return(0) else return(f(2*phi))

It is manifestly simple, using only (a) comparison to $1/2,$ (b) multiplication by $2,$ and (c) subtraction of $1.$

This algorithm randomly walks the binary tree determined by binary expansions of real numbers in the interval $[0,1].$ It outputs $1$ as soon as it enters a branch all of whose ultimate values will be less than $\phi$ and it outputs $0$ as soon as it enters any branch all of whose ultimate values will be $\phi$ or greater.

You can easily establish that the chance of outputting $1$ is no less than any finite binary number less than $\phi$ and is no greater than any finite binary number greater than $\phi,$ demonstrating $f$ implements a Bernoulli$(\phi)$ variable.

It is also straightforward to show that on any call, $f$ has a $1/2$ chance of terminating, whence (a) it will terminate almost surely (b) with an expected number of calls equal to $1+1/2+1/2^2+\cdots = 2.$

Here is an R implementation. sample.int(2,1) implements the sequence of $X_i:$ it returns 1 and 2 with equal probabilities.

f <- function(phi) {
  X <- sample.int(2, 1)
  if(isTRUE(phi >= 1/2)) {
    if (isTRUE(X == 1)) return(1) else return(f(2*phi - 1))
  } else {
    if (isTRUE(X == 1)) return(0) else return(f(2*phi))
  }
}

I applied this two thousand times to each of 128 randomly-generated floating point numbers in $[0,1],$ keeping track of the calls to $f$ and comparing the mean value (which estimates $\phi$) to $\phi$ itself with a Z score. This required generating a quarter million Bernoulli$(\phi)$ values (for various $\phi$). On a single core it took 2.5 seconds, showing it is practicable and reasonably efficient.

These graphics summarize the results.

Most Z scores are between $-2$ and $2,$ as expected of a correct procedure.

All averages are close to $2,$ as claimed, and do not depend on the value of $\phi,$ as indicated by the near-horizontal Loess smooth.

It is rare, in any of these simulations, for any call to $f$ to nest more deeply than $15$ in the recursion stack. In other words, there is essentially no risk that any one call to $f$ will take an inordinately long time. (This can be proven by examining the hypergeometric distribution of the number of calls to $f.$)

R code

This is the full (reproducible) simulation study.

#
# The algorithm.  It requires 0 <= phi <= 1.
#
f <- function(phi) {
  COUNT <<- COUNT+1      # For the study only--not an essential part of `f`
  X <- sample.int(2, 1)
  if(isTRUE(phi >= 1/2)) {
    if (isTRUE(X == 1)) return(1) else return(f(2*phi - 1))
  } else {
    if (isTRUE(X == 1)) return(0) else return(f(2*phi))
  }
}
#
# Simulation study.
#
set.seed(17)
replications <- 2e3
COUNT <- 0
system.time({
  results <- sapply(seq_len(128), function(i) {
    phi <- runif(1) # A (uniformly) random probability to study
    COUNT <<- 0     # Total number of calls to `f`
    MAX <- 0        # Largest number of calls for any one value of `phi`
    sim <- replicate(replications, {
      count <- COUNT
      x <- f(phi)
      if (COUNT - count > MAX) MAX <<- COUNT - count
      x
    })
    m <- mean(sim)                     # The simulation estimate of `phi`
    se <- sqrt(var(sim) / length(sim)) # Its standard error
    c(Value=phi, Estimate=m, SE=se, Z=(m-phi)/se, 
      Calls=COUNT, Max=MAX, Replications=length(sim), Expectation=COUNT/length(sim))
  })
})
#
# Plots.
#
X <- as.data.frame(t(results))
sub <- paste(replications, "replications")
ggplot(X, aes(Value, Z)) + geom_point(alpha=1/2) + ggtitle("Z Scores", sub)
ggplot(X, aes(Value, Max)) + geom_point(alpha=1/2) + ggtitle("Most Calls to f", sub)
ggplot(X, aes(Value, Expectation)) + geom_point(alpha=1/2) + geom_smooth(span=1) +
  ggtitle("Average Calls to f", sub)