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

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.

## 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

• Welcome to CV, Peter O.! I quite like this question. :) Dec 25, 2021 at 1:47

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

• Indeed, a similar approach is used to estimate $4/(3\pi)$, except using either x or y coordinate. In the specific case of your example, it's possible to implement the specific hit test using only integer arithmetic and without resorting to square root or transcendental functions, but this might not be possible for hit tests involving other shapes. Also, to be clear, my question is neither a homework nor a self-study assignment, nor is this coursework. Aug 25, 2021 at 19:11
• I think you made it clear from the start that your main interest is in simulation, but I wanted to be sure you know that other methods are possible. Aug 25, 2021 at 23:43
• For your information, here is an example of the sort of algorithms I had in mind; simulation experiments whose expected value is an irrational number. This example, however, like many others, is suspected to have appeared before in the probability literature. Aug 26, 2021 at 12:04

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)

• A problem here might be the requirement that the algorithm may not 'Calculate an irrational number's base-n expansion directly.' (but maybe this is a very strict, too strict, requirement if you want to get to a universal algorithm) Dec 24, 2021 at 17:11
• @Sextus The expansion is never fully calculated, though (obviously). Ruling out the simple transformations needed for partial identification of rational (dyadic) approximants would seem to be overly artificial. One could circumvent such an arbitrary rule by using another base (even nonintegral or irrational bases!) and, if that isn't satisfactory, the base could even be changed arbitrarily from one step to the next. At some point, ruling out all such solutions just becomes ridiculous, IMHO.
– whuber
Dec 24, 2021 at 17:18
• An algorithm like yours does not meet the requirements stated in my question, since it obviously works off an irrational number's $n$-adic expansion. It doesn't matter whether $n$ is 2, 10, or another real number. To get a sense of the kind of algorithms I have in mind, see the other answer or the introduction of "On Buffon machines and numbers". However, a "less than universal" algorithm can be interesting if it uses a number's series expansion but only indirectly, such as Flajolet's algorithm for $1/\pi$. Dec 24, 2021 at 18:48
• As I pointed out in an earlier comment, Peter, at that point your requirements take on an appearance of being ludicrous. What say you to varying $n,$ for instance? That is, where is the dividing line between a series of simple arithmetic operations and some kind of numerical expansion?? Moreover, it's more than a little disingenuous to alter your question after it has been answered specifically to bar that kind of answer!
– whuber
Dec 24, 2021 at 21:45
• In any case, your method requires the direct use of the constant in question, namely phi (which in R is a floating-point number which is ultimately a rational number, not an irrational one). When I wrote (before you posted your answer) that the simulation must involve only "unbiased random bits" and "integer arithmetic", that does not include arithmetic that involves the direct use of the constant in question (phi) or any other non-integer constant (rather, the only arithmetic involved is between integers). Dec 24, 2021 at 22:29