Simulating an Exponential Random Variable, given Bernoulli, Uniform 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.
 A: 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.
