Simulating an Exponential Random Variable, given Bernoulli, Uniform

Question

Suppose I am able to perform the following actions as many times as I would like:

• Simulate a Bernoulli random variable, with success of probability $p > 0$, and
• Simulate a uniform random variable on $[0, 1]$.

Suppose furthermore that $p$ is unknown to me.

Using these resources, is it possible to construct an algorithm which exactly simulates an exponential random variable of rate $p$?

The algorithm can have a random run-time, but I insist that it have an almost-surely finite run-time.

Answer 1

Here's one approach. You can use the uniforms to simulate a unit-rate Poisson process. Simulate say $n$ i.i.d. uniforms, denoted $\{U_i\}_{i=1}^n$, and let $$ E_i = -\log(U_i), $$ then $E_i \stackrel{\text{i.i.d.}}{\sim} \text{Exp}(1)$.

Think of these $E_i$ as the interarrival times of your Poisson process. It's a well-known fact that if you thin out the Poisson process by accepting the jumps with probability $p$ (i.e. reject with probability $1-p$), then you get a Poisson process with rate $p$.

Thus you can simulate $n$ Bernoulli(p) random variables, denote them $\{B_i\}_{i=1}^n$, and for each $i$ ignore the jump if the $B_i=0$. Then you the new interarrival times correspond to a Poisson process with parameter $p$, whose interarrivals times follow an Exp($p$) distribution.

For example, if you first 6 $E_i$ are 1.2, 1, 0.5, 2, 3, 0.6. And your first 6 $B_i$ are $0,0,1,1,0,1$. Then the Exp($p$) samples are $1.2+1+0.5$, $2$, $3+0.6$.

Edit: I went and implemented the approach in Python to make sure I described it correctly. It seems to work, and here's the code if you're interested.

import numpy as np

n = 10**6
p=0.4

x = np.random.exponential(1, size=n)
ber = np.random.binomial(1, p, size=n)

indices = np.argwhere(ber).flatten() + 1

a = np.split(x, indices)
interarrivals = [np.sum(arr) for arr in a if arr.size > 0]

print(np.mean(interarrivals))
print(np.var(interarrivals))
print(1./p)
print(1./p**2)

which gave the output

2.499064603937182
6.230236009875876
2.5
6.249999999999999

Thus the first and second moments seem to be correct, so I have confidence in the implementation.