Designing Perfect Simulation Algorithms using Local Correctness

TL;DR

This paper introduces the Fundamental Theorem of Perfect Simulation (FTPS), providing a unified framework with necessary and sufficient conditions for recursive algorithms to generate exact samples from a distribution, and demonstrates its application to improve Bernoulli Factory efficiency.

Contribution

The paper presents the FTPS, a unifying theoretical foundation for perfect simulation algorithms, and applies it to develop a faster Bernoulli Factory for linear functions.

Findings

FTPS characterizes correctness of recursive perfect simulation algorithms.

The framework applies to multiple existing methods like Acceptance/Rejection and CFTP.

New Bernoulli Factory is 41% faster than previous methods.

Abstract

Consider a randomized algorithm that draws samples exactly from a distribution using recursion. Such an algorithm is called a perfect simulation, and here a variety of methods for building this type of algorithm are shown to derive from the same result: the Fundamental Theorem of Perfect Simulation (FTPS). The FTPS gives two necessary and sufficient conditions for the output of a recursive probabilistic algorithm to come exactly from the desired distribution. First, the algorithm must terminate with probability 1. Second, the algorithm must be locally correct, which means that if the recursive calls in the original algorithm are replaced by oracles that draw from the desired distribution, then this new algorithm can be proven to be correct. While it is usually straightforward to verify these conditions, they are surprisingly powerful, giving the correctness of Acceptance/Rejection, Coupling from the Past, the Randomness Recycler, Read-once CFTP, Partial Rejection Sampling, Partially Recursive Acceptance Rejection, and various Bernoulli Factories. We illustrate the use of this algorithm by building a new Bernoulli Factory for linear functions that is 41\% faster than the previous method.

Equations

$$0<Z_h={\int }_{x}h(x)d\mu <\infty .$$

$$\mathbb{P}(X\in A)=\frac{{\int }_{x\in A}h(x)d\mu }{Z_h}.$$

$$\mathbb{P}(X=i)$$

$$\mathbb{P}(Y=i)=\frac{1}{{\alpha }_0}+\left(1-\frac{5}{{\alpha }_0}\right)\frac{1}{5}=\frac{1}{5}.$$

$$\mathbb{P}(X\in A)$$

$$=\mathbb{P}(Y_n\in A)-\mathbb{P}(Y_n\in A,T>n)+\mathbb{P}(X\in A,T>n).$$

$$|\mathbb{P}(X\in A)-{\pi }_{\alpha }(A)|$$

$$=|\mathbb{P}(X\in A,T>n)-\mathbb{P}(Y_n\in A,T>n)|$$

$$\le \mathbb{P}(T>n).$$

$$\mathbb{P}(U\le h(X)/g(X))$$

$$=\mathbb{E}[\mathbb{E}[\text{mathds}1(U\le h(X)/g(X))|X]]$$

$$=\mathbb{E}[h(X)/g(X)]$$

$$={\int }_x(h(x)/g(x))(g(x)/Z_g)d\mu$$

$$={\int }_{x}h(x)/Z_{g}d\mu =Z_h/Z_g>0$$

$$\mathbb{P}(W\in A)$$

$$=\frac{{\int }_{x\in A}g(x)\frac{h(x)}{g(x)}d\mu }{Z_g}+\left(1-\frac{Z_h}{Z_g}\right)\mathbb{P}(Y\in A)$$

$$=\frac{{\int }_{x\in A}h(x)d\mu }{Z_h}\cdot \frac{Z_h}{Z_g}+\left(1-\frac{Z_h}{Z_g}\right)\mathbb{P}(Y\in A)$$

$$=\mathbb{P}(Y\in A)\frac{Z_h}{Z_g}+\left(1-\frac{Z_h}{Z_g}\right)\mathbb{P}(Y\in A)$$

$$=\mathbb{P}(Y\in A),$$

$$(\forall r\in A_{\alpha })(\forall x_0,x_1\in \Omega )({\varphi }_{\alpha }(x_0,r)={\varphi }_{\alpha }(x_1,r)).$$

$$\mathbb{P}(W\in A)$$

$$=\mathbb{P}({\varphi }_{\alpha }(Y,R)\in A,R\in A_{\alpha })+\mathbb{P}({\varphi }_{\alpha }(Y,R)\in A,R\notin A_{\alpha })$$

$$=\mathbb{P}({\varphi }_{\alpha }(Y,R)\in A)$$

$$=\pi (A),$$

$$\mathbb{P}(W=1)$$

$$=\mathbb{P}(B=1)[\mathbb{P}(X=1)+\mathbb{P}(X=0)\mathbb{P}(Y=1)]$$

$$=(C/(1+C))[p+(1-p)Cp/(1+Cp)]=Cp/(1+Cp),$$

$$\mathbb{P}(W=1)$$

$$=(Cp{)}^{i-1}\left[\frac{Cp}{1+Cp}+\frac{1}{1+Cp}(Cp{)}^2\right]=(Cp{)}^i.$$

$$\mathbb{E}[T]=\frac{C}{1+C}\left[1+(1-p)\mathbb{E}[T]\right],$$

$$\mathbb{E}[T]=\frac{n}{2r-1}\cdot \frac{1-\left((1-r)/r\right)}{1-{\left((1-r)/r\right)}^n}-\frac{1}{2r-1}.$$

$$\mathbb{E}[T]\le \frac{n-i}{2r-1}.$$

$$\mathbb{E}[T]\le \frac{C(1+{ϵ}^{-1})}{(1-\exp (-3.55))(1-2\exp (-3.55/2))}$$

$$\mathbb{E}[T_0]\le \lceil \frac{3.55}{ϵ}\rceil \cdot \frac{1+Cp}{1-Cp}\cdot \frac{1-Cp}{1-(Cp{)}^{3.55/ϵ}}\cdot \frac{C}{1+Cp}=\lceil \frac{3.55}{ϵ}\rceil \cdot \frac{C}{1-(Cp{)}^{3.55/ϵ}}.$$

$$\mathbb{E}[T_0]\le \lceil \frac{3.55}{ϵ}\rceil \frac{C}{1-\exp (-3.55)}.$$

$$\mathbb{E}[T_k]\le \lceil 3.55\cdot 2^k/ϵ\rceil \frac{1-Cp}{1+Cp}\cdot \frac{C}{1-Cp}=\frac{\lceil 3.55\cdot 2^k/ϵ\rceil C}{1-Cp}.$$

$$Cp=(1-q)\cdot 1+q\cdot (Cp{)}^i\Rightarrow q=(1-Cp)/(1-(Cp{)}^i)\le (1-Cp)/(1-\exp (-3.55)).$$

$${\beta }^{-i}\le {\left(\frac{1-ϵ/2}{1-ϵ}\right)}^{-3.55/ϵ}\le \exp (-3.55/2)$$

Full text

1 Introduction

The ability to sample from complicated unnormalized densities is very beneficial in building $(\epsilon,\delta)$-randomized approximation schemes for many #P complete problems. Examples include approximation algorithms for the volume of a convex body [4], the permanent of 0-1 matrices [14, 9], the normalizing constant of the Ising model [15], the number of solutions to a Disjunctive Normal Form logical statement [16], the number of $s-t$ paths in a network [20] and many more.

Definition 1**.**

For a nonnegative measurable function $h$, say that $X$ has unnormalized density $h$ with respect to measure $\mu$ (write $X\sim h$) if

[TABLE]

and for all measurable sets $A$,

[TABLE]

For most such problems, exact computation of $Z_{h}$ is a #P complete problem, as it is for all the examples listed in the first paragraph. Hence randomized approximation algorithms are widely used to give estimates of the solution rather than the exact answer.

Given a target distribution $\pi$, a valuable tool in the construction of algorithms for drawing a random variate $X$ exactly from $\pi$ is recursion.

Definition 2**.**

Say that Alg$(\alpha)$ is a probabilistic recursive scheme (PRS) if the algorithm makes random choices, and in the course of running is allowed to call itself recursively, possibly with different values of the parameter $\alpha$.

As an example, consider the following acceptance/rejection (AR) style algorithm that draws uniformly from $\{1,2,3,4,5\}$.

[TABLE]

AR1 uses a repeat loop, but this is equivalent to the following recursive form.

[TABLE]

In practice, code is typically faster using the repeat loop formulation rather than the recursive version. However, from an execution point of view AR1 and AR2 are the same algorithm. Hence we will say that both AR1 and AR2 are examples of a PRS.

Definition 3**.**

Given a parameter set $\mathcal{P}$ and family of target distributions $\{\pi_{\alpha}:\alpha\in\mathcal{P}\}$, a PRS Alg$(\alpha)$ is a perfect simulation or perfect sampling algorithm if for all $\alpha$, the algorithm terminates in finite time with probability 1 and returns a draw with distribution equal to $\pi_{\alpha}$.

Lemma 1**.**

AR2 is a perfect simulation algorithm.

Proof.

Here $\pi_{\alpha}$ is the uniform distribution over $\{1,2,3,4,5\}$ for all $\alpha$. At each step there is a $5/6$ chance the algorithm will terminate without a recursive call, so the algorithm terminates in finite time with probability 1. Also, for $i\in\{1,\ldots,5\}$, the recursive call is exactly the same as the original call, so the probability the output is $i$ equals the chance that the initial draw of $X$ is $i$ plus the chance that the initial draw of $X$ is 6 times the chances that the recursive call returns $i$. That is,

[TABLE]

which yields $\mathbb{P}(X=i)=1/5$ for all $i\in\{1,\ldots,5\}$. ∎

Here the proof was easy because the recursive call was exactly the same as the original call. But now consider an adaptive AR algorithm that modifies the random choice based on the input parameter. Suppose that the input parameter $\alpha$ is a positive integer.

[TABLE]

Intuitively, this should still return output uniformly distributed over $\{1,2,3,4,5\}$. However, this minor change has made the analysis much more difficult, as the technique used in the last proof no longer applies since there is no reason initially to believe (for instance) that $\textbf{Alg}(100)$ and $\textbf{Alg}(10)$ have the same output distribution.

The purpose of this paper is to present a new version of the Fundamental Theorem of Perfect Simulation (FTPS) that gives two necessary and sufficient conditions for when a PRS is actually a perfect simulation algorithm. These conditions are typically both easy to verify. The first condition is that the algorithm must terminate in finite time with probability 1. The second condition is called local correctness.

Definition 4**.**

In a call $\textbf{Alg}(\alpha)$ to a PRS, suppose that any recursive calls of the form $\textbf{Alg}(\beta)$ within the algorithm are replaced with oracles that draw exactly from their target distribution $\pi_{\beta}$. If the output of the resulting algorithm with oracles can be shown to come from $\pi_{\alpha}$ for all $\alpha$, then the algorithm is locally correct with respect to $\{\pi_{\alpha}:\alpha\in\mathcal{P}\}$.

When these two conditions hold for all input parameters, the output of the PRS will be exactly from the desired result.

Theorem 1** (Fundamental Theorem of Perfect Simulation).**

Suppose $\textbf{Alg}(\alpha)$ is a PRS that satisfies the following.

It terminates in finite time with probability 1 for all $\alpha$. 2. 2.

It is locally correct with respect to $\{\pi_{\alpha}:\alpha\in\mathcal{P}\}$.

Then the output of the algorithm is an exact draw from $\pi_{\alpha}$.

This simple condition is surprisingly powerful. Consider the example from earlier.

Lemma 2**.**

$\textbf{AR3}(\alpha)$* is a perfect simulation algorithm whose output is uniform over $\{1,2,3,4,5\}$ for all positive integers $\alpha$.*

Proof.

Here $\pi_{\alpha}$ is uniform over $\{1,2,3,4,5\}$ for all positive integers $\alpha$. Since the parameter in any recursive call is at most the initial value $\alpha_{0}$ passed to the algorithm, the probability of termination is at least $5/\alpha_{0}$ and the algorithm will terminate in finite time with probability 1.

Then to show correctness using the FTPS, we assume that the recursive call at line 5 actually is a draw from the correct target distribution that is uniform over $\{1,\ldots,5\}$. Let $\alpha_{0}$ be the initial parameter passed to the algorithm. If $\alpha_{0}\leq 5$ then line 1 just draws uniformly from $1$ to $5$, so assume $\alpha_{0}>5$. Then for $i\in\{1,\ldots,5\}$, the chance that $Y=i$ is the chance that $X=i$ plus the chance that $X>5$ times the chance that the oracle returns $Y=i$. That is,

[TABLE]

Hence the algorithm is locally correct, and the FTPS gives that the overall algorithm is correct. ∎

The value of the FTPS is similar to the use of the Ergodic Theorem in Markov chain Monte Carlo. Without going into too much detail, that theorem states that a Markov chain that is irreducible, aperiodic, and stationary with respect to $\pi$ will have $\pi$ as the limiting distribution regardless of the starting state. Irreducibility and aperiodicity is similar to our first requirement that the algorithm must terminate with probability 1. Stationarity is similar to local correctness in that it is usually not too difficult to verify (although there are exceptions.)

The simplicity of the Ergodic Theorem has led to a multitude of Markov chain Monte Carlo algorithms. In the same way, the simplicity of the local correctness criterion encourages variety of design in perfect simulation algorithms. In particular, what the FTPS does for us is three-fold.

It gives easy proofs of the most common perfect simulation protocols, such as Acceptance/Rejection (AR), Coupling from the Past (CFTP), and the Randomness Recycler (RR). 2. 2.

It allows us to build variants of these algorithms, such as adaptive AR, partially recursive AR, time-heterogeneous CFTP, and fractal time CFTP, without having to start over from scratch on the proofs. We can make these minor changes without worrying that a minor modification of the algorithm will change the resulting output. 3. 3.

It allows the building of more complicated algorithms that use recursion in more unusual ways, such as the first polynomial expected time Bernoulli Factory [12].

The remainder of this paper is organized as follows. In the next section we give the proof of the FTPS. In Section 2 we employ the result to give new proofs of many of the most important protocols in the construction of perfect simulation algorithms. In particular we present a new Bernoulli Factory algorithm whose proof of correctness is greatly eased by using the FTPS. Finally, in Section 4 we consider ramifications of the model of computation that we are using.

2 Proof of the FTPS

In order to prove the FTPS, we will consider a coupled sequence of algorithms. Our initial call to $\textbf{Alg}(\alpha)$ is allowed to use recursion an unbounded number of times. Suppose that for each $n$ we construct a new algorithm $\textbf{Alg}_{n}(\alpha)$ in the following way. First we consider what level of recursion we are at in our algorithm.

Definition 5**.**

An initial call to a PRS Alg is said to be at level of recursion 0. For $i$ a positive integer, say the call to Alg occurs at the level of recursion $i$ if it was called by a call that was at level of recursion $i-1$.

Then set up $\textbf{Alg}_{n}(\alpha)$ as follows. If the level of recursion $\ell$ of a call to $\textbf{Alg}_{n}(\alpha)$ is less than $n$, then the algorithm behaves exactly like $\textbf{Alg}(\alpha)$. However, if the the call is at level of recursion $\ell=n$, the algorithm replaces the recursive calls with oracles that generate draws exactly from the target distribution.

So for example, we might call $\textbf{AR3}(100)$ and have the following outcome:

[TABLE]

Throughout this section, we will use $X$ to denote the output from $\textbf{Alg}(\alpha)$, and $Y_{n}$ to denote the output of $\textbf{Alg}_{n}(\alpha)$. For an original call, $\textbf{Alg}(\alpha)$, let $T$ denote the supremum of the levels of recursion employed by the algorithm in the call. By definition then $Y_{n}=X$ for any $n\geq T$.

In the sample run, since $\textbf{Alg}(11)$ was the last call at level of recursion 2, $T=2$, and $X=4=Y_{2}=Y_{3}=\cdots$. The outputs $Y_{0}$ and $Y_{1}$ used oracles at level of recursion 0 and 1 respectively to determine their value and so might have different values than $X$.

The first step in the proof of the FTPS is to use local correctness to show that for every $n$, $Y_{n}$ has the correct distribution.

Lemma 3**.**

For all nonnegative integers $n$, the output of $\textbf{Alg}_{n}(\alpha)$ (where Alg is locally correct) is an exact draw from $\pi_{\alpha}$.

Proof.

When $n=0$, this is just the definition of local correctness. Suppose it holds for $\textbf{Alg}_{n}(\alpha)$, and consider $\textbf{Alg}_{n+1}(\alpha)$.

For the output of any call at level of recursion $n+1$, any further recursive calls use oracles. So local correctness guarantees that the output of level of recursion $n+1$ comes exactly from the desired distribution. So it is the same as if that output came from an oracle with the correct distribution. But that means $Y_{n+1}$ has the same distribution as $Y_{n}$. That allows us to use the induction hypothesis to say that $Y_{n+1}$ must have the correct distribution, which completes the induction. ∎

We are now ready to prove the FTPS.

Proof of the FTPS.

Let $A$ be any measurable set, and $n$ be any positive integer. Then

[TABLE]

If $T\leq n$, then $X=Y_{n}$, so

[TABLE]

The probability that $Y_{n}\in A$ is the probability that a draw from $\pi_{\alpha}$ falls in $A$. Denote this probability $\pi_{\alpha}(A)$. Then

[TABLE]

Since we assumed that $T$ was finite with probability 1, as $n\rightarrow\infty$, $\mathbb{P}(T>n)\rightarrow 0$, hence $\mathbb{P}(X\in A)=\pi_{\alpha}(A)$, and we are done. ∎

3 Perfect simulation protocols

The two most common perfect simulation protocols are acceptance/rejection (sometimes called rejection sampling) and Coupling from the Past, so let us begin with those.

3.1 Acceptance Rejection

Consider the following acceptance/rejection algorithm for unnormalized densities. Suppose that $g\geq h$ are two unnormalized densities such that it is easy to draw random variates from $g$ and we wish to obtain draws from $h$. The following algorithm goes back to [22] but has been extended and used for applications ranging from counting DNF satisfying assignments [16] to drawing exact instances of solutions to Stochastic Differential Equations [2].

[TABLE]

Lemma 4**.**

For $h$ and $g$ unnormalized densities with $g\geq h$, the output of AR4 is $X\sim h$.

Proof.

By the FTPS we need termination in finite time and local correctness. The chance of terminating at each step is

[TABLE]

Here $\mathds{1}(p)$ is the usual indicator function that is 1 if its argument is true and 0 otherwise. Since the chance of terminating at each step is a fixed positive number, the chance of running infinitely often is 0.

Now for local correctness. Let $A$ be a set measurable with respect to $\mu$. Then the chance that the output lies in $A$ is the chance that the initial draw of $X\in A$ times the chance we accept plus the chance that we reject the initial draw and the recursive call $Y$ falls in $A$ once we substitute in the true distribution for the recursive call. In notation, if we call the output $W$,

[TABLE]

and so the output has the correct distribution. ∎

We call $g$ the envelope density. In adaptive AR, the envelope density is refined based on the value of $X$. For example, the approach for Log-concave density functions of Wild and Gilks [23] works like this. In general, the envelope function is itself modified as a function of the rejection sample.

[TABLE]

Lemma 5**.**

Suppose for each $g$ and rejected $X$, $g_{a}$ is still an envelope density of $h$. Then the output of AR4 is $X\sim h$.

Proof.

At each recursive call the chance of terminating can only stay the same or increase, and so the chance that we never terminate is still bounded above by 0.

For local correctness, we substitute the recursive call with a direct oracle, and so the change in the argument for the recursive call does not matter. The same proof as in the previous lemma still works. ∎

3.2 Coupling from the Past

Coupling from the Past (CFTP) was invented by Propp and Wilson [18] as a way to turn approximate sampling by a Markov chain into perfect simulation. It has since become a mainstay of perfect simulation algorithms. There have been many variants and extensions, for instance Read-Once CFTP [24], clan of ancestors [5], Bounding chain CFTP [8], time-heterogeneous CFTP [10], and fractal time CFTP [11].

Rather than repeat the entire evolution of the idea, here we move straight to the most general version, which revolves around the notion of a stationary update function.

Definition 6**.**

Say that $\phi:\Omega\times\mathcal{R}\rightarrow\Omega$ is a stationary update function with respect to $\pi$ if there is a probability distribution $\mathbb{P}$ over $\mathcal{R}$ such that if $R\sim\mathbb{P}$ and $X\sim\pi$ then $\phi(X,R)\sim\pi$ as well.

We refer to $X$ as the state and $R$ as the random choices that update the state. For instance, $\phi$ might encode taking one or more steps in a Metropolis-Hastings Markov chain whose stationary distribution is $\pi$.

For some random choices, the current state might be immaterial in determining the next state. That is, suppose for all $x\in\Omega$, $\phi(x,R)=y$. Then in this case we say the state has completely coupled.

With this in mind, we present here a very general version of CFTP. For a given parameter $\alpha$, suppose we use stationary update function $\phi_{\alpha}$. Moreover, let $A_{\alpha}$ be any set such that

[TABLE]

[TABLE]

Lemma 6**.**

Suppose $\phi_{\alpha}$ is stationary with respect to $\pi$ for all parameters $\alpha$. If a call to $\textbf{CFTP}(\alpha)$ terminates with probability 1, then the output comes from $\pi$.

Proof.

Since we are assuming finite termination, we need only verify local correctness. Let $A$ be a $\pi$-measurable set. Let $W$ be the output of $\textbf{CFTP}(\alpha)$. Let $x$ be any element of $\Omega$.

[TABLE]

since $\phi_{\alpha}$ is a stationary update function. ∎

3.3 Bernoulli Factories

Bernoulli factories were introduced by Asmussen et. al. [1] as part of an algorithm for generating samples exactly from the stationary distribution of a regenerative Markov process. The idea is as follows.

Suppose you have access to a stream of independent, identically distributed (iid) Bernoulli random variables $X_{1},X_{2},\ldots$ with unknown parameter $p$. So $\mathbb{P}(X_{i}=1)=p$, $\mathbb{P}(X_{i}=0)=1-p$. Write $X_{i}\sim\textsf{Bern}(p)$. Now suppose that we wish to construct a Bernoulli random variable with a parameter that is a function of $p$ using a random number of $\{X_{i}\}$. That is, we want $Y\sim\textsf{Bern}(f(p))$ for a known function $p$ where $Y=f_{\text{BF}}(X_{1},\ldots,X_{T})$ for some function $f_{\text{BF}}$ and stopping time $T$. Such an algorithm is a Bernoulli factory.

For instance, to draw $Y\sim\textsf{Bern}(p(1-p))$, just let $Y=X_{1}(1-X_{2})$. For the application in [1], they needed to be able to draw $W\sim\textsf{Bern}(Cp)$ where $C$ is a known constant. This problem, although simple to state, turns out to be surprisingly difficult. Nacu and Peres [17] showed that for any function that is analytic and bounded away from 1, it suffices to have a Bernoulli factory for $2p$.

The first polynomial expected time algorithm for this problem was given in [12]. Here we present a new version of the algorithm that is slightly simpler to implement and analyze and takes advantage of the FTPS.

To construct this linear factory, it helps to have a factory for $Cp/(1+Cp)$ available.

[TABLE]

Lemma 7**.**

The output of BF1$(C)$ is $\textsf{Bern}(Cp/(1+Cp))$.

Proof.

If $B=0$ the algorithm does not call itself recursively, therefore the chance that the algorithm does not terminate in finite time is 0.

Let $W$ be the output of the algorithm. To show local correctness, we assume the final call in line 8 returns $Y\sim\textsf{Bern}(Cp/(1+Cp))$. Then $W$ must be in $\{0,1\}$, and

[TABLE]

giving local correctness. ∎

Rather than just build a $Cp$ Bernoulli factory, we will build a $(Cp)^{i}$ Bernoulli factory where $C\geq 1$ and $i\in\{0,1,\ldots\}$. For this factory, we suppose that we know $\epsilon>0$ such that $Cp\leq 1-\epsilon$. If no such $\epsilon$ exists, then it is impossible to build such an algorithm, see [12] for the proof.

[TABLE]

Lemma 8**.**

For $Cp<1-\epsilon$, $i\in\{0,1,2,\ldots\}$, the output of $\textbf{BF2}(C,i,\epsilon)$ is $\textsf{Bern}([Cp]^{i})$.

Proof.

Let $W$ denote the output of the call to the algorithm. When $i=0$ the result is true since $W=1$. If $i>3.55/\epsilon$, then to show local correctness, consider $\mathbb{P}(W=1)$.

If $Cp<1-\epsilon$, then $C\beta p=C(1-\epsilon/2)p/(1-\epsilon)\leq 1-\epsilon/2$. So we assume the recursive call in line 6 comes from an oracle, and hence returns $\textsf{Bern}([\beta Cp]^{i})\sim\textsf{Bern}(\beta^{i}[Cp]^{i})$. Then the probability $W=1$ is $(\beta^{-i})(\beta^{i}[Cp]^{i})=[Cp]^{i}$.

Now suppose $i>0$ and $i<3.55/\epsilon$. Then $W$ is the result of line 10, which depends on $B_{2}$. So

[TABLE]

Hence the algorithm is locally correct.

If $i>3.55/\epsilon$, then $\beta^{-i}\leq\exp(-3.55/2)$. So every time that $i$ hits 0 or exceeds $3.55/\epsilon$ there is a $1-\exp(-3.55/2)$ chance of stopping. When $i$ falls inside those extremes, there is at least a $1/2$ chance that $i$ increases by 1, hence there will be a finite number of steps with probability 1 until $i$ reaches 0 or exceeds $3.55/\epsilon$, giving a finite number of steps with probability 1 until termination. ∎

See the appendix for an evaluation of the expected running time.

3.4 Other protocols

There are many other protocols for building perfect simulation algorithms that can employ the FTPS. For instance, the design property of the Randomness Recycler [6, 13] gives local correctness for that approach. Sink popping algorithms [19, 3] can also be written recursively, which makes their proof of correctness amenable to use of the FTPS.

Similarly, Partial Rejection Sampling [7] can be implemented recursively, allowing the FTPS to be used to show correctness.

4 Models of Computation

Until now we have been looking at these algorithms through the lens of floating point operations that can be carried out exactly. What changes when we move to a Turing machine model of computation?

When dealing with randomized algorithms, we can use a Probabilistic Turing Machine (PTM), which is essentially a Turing Machine that can access a tape with an iid stream of $\textsf{Bern}(1/2)$ random variates. In this framework, we cannot even simulate exactly a single $U$ uniform over $[0,1]$. Instead, for any possible error $\delta$, we can compute a result that falls within distance $\delta$ of the actual random variable.

When we call the algorithm for a PTM, we must specify the error tolerance that we are willing to accept. Presumably we are only using Turing computable functions. Therefore, for a potential recursive call, it is possible to determine what the error for the recursive call should be in order to obtain the desired error for the original call. The FTPS guarantees that the exact algorithm has the correct distribution, and so the actual computed output will come within the target error of a truly exact draw.

Of course, none of this is done in practice since usually machine epsilon is typically small enough to give the desired error for the original call.

5 Appendix: Expected running time of

BF2

We first consider BF1.

Lemma 9**.**

The expected number of Bernoulli draws needed for BF1 is $C/(1+Cp)$.

Proof.

Let $T$ be the number of Bernoulli draws used by the algorithm. Each call to the algorithm uses either 0 or 1 draws, and then either calls itself recursively or does not. At each step there is at least a $1/(1+C)$ chance of not recursively calling the algorithm, and so $\mathbb{E}[T]$ is finite. This gives rise to the following recursion

[TABLE]

which solves to $\mathbb{E}[T]=C/(1+Cp)$. ∎

When we call BF1 from BF2 we use the result to change the number of $Cp$ coins we need to flip, either increasing by one with probability $1/(1+Cp)$ or decreasing it by one with probability $Cp/(1+Cp)$. This is known as the gambler’s ruin problem. It has been well studied (see, for instance, [21].)

The gambler’s ruin problem is a Markov chain where, the state moves from $i$ to $i+1$ with probability $r$, and to $i-1$ with probability $1-r$. The following result is well-known.

Lemma 10**.**

Let the initial state be $1$ and suppose the gambler’s ruin Markov chain stops when the state reaches [math] or $n$. Let $T$ denote the number of steps that are taken until this is reached. Then for $r>1/2$,

[TABLE]

Also, for any starting state $i\in\{1,2,\ldots,n\}$,

[TABLE]

Using this, we can upper bound the running time of $\textbf{BF2}(C,1,\epsilon)$.

Lemma 11**.**

Let $T$ be the number of Bernoulli draws in $\textbf{BF2}(C,1,\epsilon)$. Then

[TABLE]

Proof.

Here $r=1/(1+Cp)$, so $1/(2r-1)=(1-Cp)/(1+Cp)$ and $(1-r)/r=Cp$. Initially we begin with $(Cp)^{i}$ and stop when $i=\lceil 3.55/\epsilon\rceil$. So from the previous lemma the number of draws $T_{0}$ obeys

[TABLE]

Using $1+x\leq\exp(x)$ and $Cp\leq 1-\epsilon$ gives

[TABLE]

If we fail to exit then we change from $i=\lceil 3.55/\epsilon\rceil$ to $\lceil 3.55\cdot 2/\epsilon\rceil$, and so on. Consider the chance that we need to make draws that move our problem from $i=\lceil 3.55\cdot 2^{k}/\epsilon\rceil$ to $i=\lceil 3.55\cdot 2^{k+1}/\epsilon\rceil$. The difference between these two numbers is at most $\lceil 3.55\cdot 2^{k}/\epsilon\rceil$. So the expected number of steps needed by the previous lemma is

[TABLE]

So now consider the question of the probability that we need to make these moves in the first place. Let $q$ denote the probability that in the initial run, $i>2/\epsilon$ is reached before $i=0$. Then by local correctness,

[TABLE]

Furthermore, there is only a $\beta^{-i}$ chance that we do not return a 0 and enter the $k$th level of recursion given that we are already at the $k-1$ level of recursion. Note

[TABLE]

Combining, the chance of entering level of recursion $k$ is at most

[TABLE]

This makes (by the Monotone Convergence Theorem) the upper bound on the expected number of steps

[TABLE]

Using $\lceil x\rceil\leq x+1$ together with standard facts about geometric series then completes the proof. ∎

The constant of $3.55$ was found by numerically optimizing the constant in front of the $C(1+\epsilon^{-1})$ term to minimize it. In [12], it was shown that any algorithm of this type must take $\Omega(\epsilon^{-1}C)$ steps on average and so the best that we can hope for is to improve the constant on the bound on the running time. Note that $[(1-\exp(-3.55))(1-2\exp(-3.55/2))]\leq 5.53$. The algorithm presented in [12] was shown to take at most $9.5\epsilon^{-1}C$ so this represents a marked improvement in the constant factor.

6 Acknowledgements

This work supported by U.S. National Science Foundation grant DMS-1418495.