PERFECT SIMULATION OF CONDITIONALLY SPECIFIED MODELS

May 26, 1998

Jesper Møller^{1}, MaPhySto^{2}

Abstract: We discuss how the ideas of producing perfect simulations based on coupling from the past for finite state space models naturally extend to multivariate distributions with infinite or uncountable state spaces such as autogamma, auto-Poisson and auto-negative-binomial models, using Gibbs sampling in combination with sandwiching methods originally introduced for perfect simulation of point processes.

Keywords: Coupling from the past; exact simulation; Gibbs sampling; locally specified exponential family distributions; Markov chain Monte Carlo; Metropolis-Hastings algorithm; spatial statistics.

AMS 1991 MATHEMATICS SUBJECT CLASSIFICATION: 62H11, 62M30, 60K35.

1 Introduction

Since Propp and Wilson’s (1996) seminal work on perfect simulation there has been an extensive interest in developing and applying their ideas in different contexts (see the survey in Propp and Wilson (1997)). Briefly, the main idea is to use coupling from the past (CFTP) and repeatedly use the same sampler for generating upper and lower Markov chains started further and further back in time until a pair of upper and lower chains coalesce at time 0, and then return the result as a perfect (or exact) simulation from a given target distribution (in Kendall (1996) and Kendall and Møller (1997) it is argued why the terminology ‘perfect’ is preferable). To do this Propp and Wilson (1996) assume that the state space is finite and equipped with a partial ordering so that the sampler is monotone and there exist a unique minimal element, in which the lower processes are started, and a unique maximal element, in which the upper processes are started. Then a chain produced by the sampler, started at any time n ≤ 0 in an arbitrary initial state, sandwiches between that pair of lower and upper chains which was started at the same time n. Thereby it can be established under weak (ergodicity) conditions for the sampler that coalescence will happens for all sufficiently large n, and by considering a ‘virtual simulation from time minus infinity’ it follows that the output is a simulation from the target distribution.

1Address for correspondence: Department of Mathematics, Aalborg University, Fredrik Bajers Vej 7E, 9220 Aalborg, Denmark, Email:jm@math.auc.dk

2MaPhySto – Centre for Mathematical Physics and Stochastics, funded by a grant from The Danish National Research Foundation.

These ideas have now been extended in various ways. Kendall (1996) and Häggström, van Lieshout and Møller (1996) outlined how to do perfect simulation for point processes, where the state space is uncountable. Especially, as Propp and Wilson (1996) in their examples required the target distribution to be attractive in a certain sense so that the Gibbs sampler becomes monotone, Kendall’s work showed how to handle the opposite repulsive case. This has been further generalised in Kendall and Møller (1997), where the role of the minimal element (the empty point configuration) is emphasized (there exists no maximal element in a point process setting); see also Kendall (1997) and Häggström and Nelander (1997). An even more general approach, but for simulating multivariate continuous distributions has very recently been studied in Murdoch and Green (1997). This and the other mentioned papers will be commented further on in this paper.

The purpose of this paper is to show how these ideas can be further extended to produce perfect simulation of multivariate discrete (Section 2) and continuous (Section 3) target distributions, where the target distribution is naturally specified through its conditional distributions of one component given the others so that Gibbs sampling is the obvious way of producing samples. The examples of such distributions to be discussed will mainly be locally specified exponential family distributions (Besag (1974); Cressie (1993)) with applications in spatial statistics such as auto-binomial, auto-Poisson, auto-negative-binomial, auto-gamma models and certain pairwise-difference interaction models (Sections 2.3 and 3.3). Indeed many other examples of models could be included, e.g. combinations of the local characteristics from different types of models may specify a joint distribution from which we can make perfect simulations. The auto-gamma model has also been used in other papers on Markov chain Monte Carlo methods, in particular in connection to a Bayesian analysis of a dataset on pump reliability (Gelfand and Smith (1990); Murdoch and Green (1997)). In relation to this some empirical results will be reported at the end of Section 3.

The techniques used in Section 2 are much inspired by Kendall and Møller (1997) and some of the terminology and notation will be borrowed from that work. Compared to Propp and Wilson (1996, 1997) and Häggström and Nelander (1997) the extension is mainly that infinite discrete state spaces are covered as well provided the model is repulsive. For definiteness an auto-gamma model is considered in Section 3 and it is demonstrated how we can make perfect simulations with an arbitrarily good accuracy by Gibbs sampling. Also the possibility of using the Metropolis-Hastings algorithm will be investigated in Section 3.

Though the practicality of doing perfect simulation in general is yet not so clear, the present paper and the contributions mentioned above at least demonstrate that it is feasible to tackle many complex distributions by exploring the monotonicity properties of the model and the sampler. One important model, which I couldn’t handle, is a conditional auto-regression as it seems impossible to bound the lower and upper processes by dominating chains (the construction of upper and lower processes will be similar to the coupling construction for the auto-binomial model in (ii), Section 2.3 - the problem is to realize (if possible) how to start these processes in the right manner - obviously, one needs both a ‘positive’ dominating chain and a ‘negative’ dominating chain).

2 Discrete multivariate distributions

Suppose we want to make simulations from a target distribution π = D(X) with X = (X_{1}, …, X_{k}) a k-dimensional discrete random vector, which is specified by its conditional distributions D(X_{i}|X_{−i}), i = 1, …, k, where X_{−i} = (X_{1}, …, X_{i−1}, X_{i+1}, …, X_{k}). Assume also that the support Ω of X is a subset of {0, 1, …}^{k} and that it contains the minimal element 0 = (0, …, 0) ∈ Ω. Let F_{i}(·|x_{−i}) denote the cumulative distribution function of D(X_{i}|X_{−i} = x_{−i}) when P (X_{−i} = x_{−i}) > 0. We can then generate a Markov chain X(t) = (X(t, 1), …, X(t, k)), t = 0, 1, …, started in X(0) = 0 and using cyclic Gibbs sampling by setting

X(t, i) = F_{i}^{−}(R(t, i)|X(t, i)_{←}) for i = 1, …, k, t = 1, 2, …

where

X(t, i)_{←} = (X(t, 1), …, X(t, i − 1), X(t − 1, i + 1), …, X(t − 1, k))

is the k − 1 states of the components just before the ith update at time t and the R(t, i) are iid uniform numbers between 0 and 1. Here the inverse F^{−} of a cumulative distribution function F is defined by F^{−}(r) = min{s : F (s) ≥ r}. In fact, if the Markov chain X(t) is irreducible it is also aperiodic and it converges weakly towards X (see, for example, Roberts and Smith (1994)). In the sequel irreducibility of X(t) is assumed.

In this section we show how CFTP can be used for producing a perfect simulation from π within finite time. It is assumed that the conditional distribution functions satisfy a certain monotonicity condition with respect to the natural partial ordering ≤ on R^{d} (d = 1, 2, …) given by (x_{1}, …, x_{d}) ≤ (y_{1}, …, y_{d}) if x_{i} ≤ y_{i}, i = 1, …, d: We assume that for any i, F_{i}(·|X_{−i}) is increasing in X_{−i}, that is,

F_{i}(·|x_{−i}) ≥ F_{i}(·|x́_{−i}) if x_{−i} ≥ x́_{−i}, x, x́ ∈ Ω.  (1)

In analogy with Kendall (1996) and Kendall and Møller (1997) we refer to (1) as the repulsive case since X_{i} tends to be smaller as X_{−i} increases (another terminology is used in Häggström and Nelander (1997)). The opposite attractive case has earlier been studied (in the binary case where X_{i} ∈ {0, 1}) in Propp and Wilson (1996) and (for point processes) in Kendall (1996), Häggström, van Lieshout and Møller (1996) and Kendall and Møller (1997). This case will be commented on further in (ii), Section 2.3.

2.1 Perfect simulation algorithm (discrete repulsive models)

Let the situation be as described above and suppose that the iid uniform variates R(t, i) are defined back in time t = −1, −2, …, too. CFTP is then obtained by reusing these random numbers in the following construction of lower and upper processes L_{n}(t) = (L_{n}(t, 1), …, L_{n}(t, k)) and U_{n}(t) = (U_{n}(t, 1), …, U_{n}(t, k)), which are started at times n and generated forwards in time: For each integer n ∈ Z, set

L_{n}(n) = 0, L_{n}(t, i) = F_{i}^{−}(R(t, i)|U_{n}(t, i)_{←}), i = 1, …, k, t > n (2)

and

U_{n}(n) = D(n), U_{n}(t, i) = F_{i}^{−}(R(t, i)|L_{n}(t, i)_{←}), i = 1, …, k, t > n (3)

with the dominating chain D(n) = (D(n, 1), …, D(n, k)) given by the mutually independent components

D(n, i) = F_{i}^{−}(R(n, i)|0_{−i}), i = 1, …, k, n ∈ Z.

Due to the conditioning in (2)-(3) we need to extend the definition of F_{i}(·|x_{−i}) when P (X_{−i} = x_{−i}) = 0. For i ∈ {1, …, k} and x_{−i} ∈ {0, 1, …}^{k−1}, define

F_{i}(·|x_{−i}) = max{F_{i}(·|y_{−i}) : y_{−i} ≤ x_{−i}, P (X_{−i} = y_{−i}) > 0}

which ensures that F_{i}(x_{i}|x_{−i}) is increasing in x_{−i}.

Now, in the perfect simulation algorithm we use a strictly decreasing sequence of non-positive starting times n = n_{j} (0 ≥ n_{1} > n_{2} > …) and repeat to generate lower and upper processes (L_{n}, U_{n}) until coalescence happens at time 0: For j = 1, 2, … set n = n_{j} and generate (L_{n}(t), U_{n}(t)), t = n, …, 0, until L_{n}(0) = U_{n}(0); return then Y = L_{n}(0) as a perfect simulation from the target distribution π (see Section 2.2).

In principle an arbitrary sequence of starting times {n_{j}} as above may be used, but following the nearly optimal choice of doubling n (see Propp and Wilson (1996)) we propose to set n_{j} = −2^{j} in applications. Notice that by the definition (2)-(3), when generating (L_{n}, U_{n}) for n = −2^{j} and j ≥ 2 we are reusing the random number stream R(t, i), i = 1, …, k, n/2 ≤ t ≤ 0, used in the generation of (L_{n/2}, U_{n/2}), …, (L_{0}, U_{0}) together with the new random number stream R(t, i), i = 1, …, k, n ≤ t < n/2. Further comments are given in Section 2.3.

2.2 Theoretical results (discrete repulsive models)

The perfect simulation algorithm actually works (at least in theory) according to a general result presented in Kendall and Møller (1997, Theorem 1). In this section we restate and verify this for the present setup, partly for completeness and since we refer to the proof later on and partly because the results presented below extend those in Häggström and Nelander (1997). Actually, Häggström and Nelander assume that Ω is bounded by some maximum x^{max}; they start each upper process in x^{max}, but compared to the perfect simulation algorithm in Section 2.1, this can only increase the coalescence time by at most one.

For n ∈ Z let X_{n}(t) denote the target chain defined by cyclic Gibbs sampling and started in the minimum 0 at time n:

X_{n}(n) = 0, X_{n}(t, i) = F_{i}^{−}(R(t, i)|X_{n}(t, i)_{←}), i = 1, …, k, t > n.

Proposition 1 For all times m ≤ n ≤ u ≤ t we have that the following sandwiching properties hold:

L_{n}(t) ≤ X_{n}(t) ≤ U_{n}(t) ≤ D(t),  (4)
L_{n}(t) ≤ L_{m}(t) ≤ U_{m}(t) ≤ U_{n}(t),  (5)
L_{m}(t) = U_{m}(t) if L_{n}(u) = U_{n}(u).  (6)

Proposition 1 is easily verified by induction. By (6), if coalescence happens at time u in the nth pair of lower and upper processes (when n ≤ u ≤ 0), we have also coalescence from time u to time 0 in the mth pair of lower and upper processes (whenever m ≤ n), and so Y = L_{m}(0) = U_{m}(0). Consequently, in order to verify that coalescence happens within finite time in the perfect simulation algorithm and the output Y follows the target distribution, it suffices to consider the case where n_{j} = −j and define

N = sup{n ≤ 0 : L_{n}(0) = U_{n}(0)},

so −N is the number of pairs of lower and upper processes needed for obtaining coalescence at time 0 when all the non-positive numbers are used as starting times (we set sup ∅ = −∞).

Theorem 1 With probability one, N > −∞, and Y = L_{N}(0) follows π.

Proof: Clearly with probability one, D(t) = 0 for some t ≤ 0, and so by (4), N > −∞ almost surely. Hence we can define a ‘virtual simulation’ Y = lim_{n→−∞} L_{n}(0) from time minus infinity, and using (4)–(6), we get with probability one that

Y = L_{N}(0) = X_{N}(0) = U_{N}(0) = lim_{n→−∞} X_{n}(0).

This together with the fact that X_{n}(0) converges weakly to X as n tends to minus infinity (since, by stationarity of R(t), D(X_{n}(0)) = D(X(−n))) imply that Y follows π = D(X).

Remarks: By stationarity of R(t),

D(L_{n}(t), U_{n}(t)) = D(L_{0}(t − n), U_{0}(t − n)) if t ≥ n.  (7)

Hence, as observed in Propp and Wilson (1996), D(−N) = D(M), where M = inf{n ≥ 0 : L_{0}(n) = U_{0}(n)} is the first time of coalescence when the lower and upper processes are started at time 0. But notice that L_{0}(M) is in general a biased sample from π (Propp and Wilson (1996)).

The proof of Theorem 1 is only based on the sandwiching properties and the fact that 0 is an ergodic atom for the dominating chain. In the terminology of Foss and Tweedie (1997), −N is a succesful backward coupling time, so the proof that Y ∼ π follows also from Foss and Tweedie’s Theorem 3.1. Clearly, if P (·, ·) is the transition probability matrix for a target chain, then by (1), P (x, 0) ≥ ∏_{1}^{k} F_{i}(0|0_{−i}) > 0 for any x ∈ Ω. Hence the state space Ω is a small set or equivalently P (·, ·) is uniformly ergodic (see, e.g., Meyn and Tweedie (1993, Theorem 16.0.2)). Although this gives uniform ergodicity, it is not for the same reason as in Foss and Tweedie (1997, Theorem 4.2), since −N is not ‘vertical’ in the sense of Foss and Tweedie: in fact, if the state space Ω is infinite, we cannot achieve a backward coupling by considering sample paths generated by the Gibbs sampler and started at any state in Ω and any time (≤ 0). However, due to the stochastic domination given by D, we need only to consider target chains started in 0.

The following lemma can be applied for establishing that the models considered in Section 2.3 below are attractive or repulsive.

Lemma 1 Let λ denote counting measure on {0, 1, …} (or Lebesgue measure on (0, ∞)). Suppose that F_{i}(·|x_{−i}) (whenever well-defined for x_{−i}) has density

f_{i}(x_{i}|x_{−i}) ∝ b_{i}(x_{i}) exp(θ_{i}(x_{−i})x_{i})  (8)

with respect to λ, where b_{i}(·) ≥ 0 and θ_{i}(·) are measurable real functions. Then F_{i}(·|x_{−i}) is decreasing (increasing) in x_{−i} if θ_{i}(x_{−i}) is increasing (decreasing) in x_{−i}.

Proof: Letting a(θ) = ∫_{0}^{∞} exp(θy)dν(y) with dν(y) = b_{i}(y)dλ(y), then Θ = {θ : a(θ) < ∞} is an interval. For θ ∈ Θ, let F_{θ}(·) denote the cumulative distribution function with density exp(θy)/a(θ) with respect to ν. For any r ≥ 0 with ν([0, r]) > 0, it is easily seen by differentation that g_{r}(θ) = ∫_{r}^{∞} exp(θy)dν(y)/ ∫_{0}^{r} exp(θy)dν(y) is increasing in θ ∈ Θ (if θ is an endpoint of Θ, we consider just left or right derivatives of g_{r}(θ)). Hence F_{θ}(r) = 1/(1+g_{r}(θ)) is decreasing in θ if ν([0, r]) > 0, and F_{θ}(r) = 0 otherwise, whereby the lemma is verified.

2.3 Comments and examples (discrete attractive and repulsive models)

(i) Some empirical findings for simple examples of binary repulsive models (i.e. when Ω ⊆ {0, 1}^{k}) such as the hard-core model, Ising anti-ferromagnet, and random cluster model (with the parameter for the number of connected components chosen so that the random cluster model is repulsive) are reported in Häggström and Nelander (1997).

(ii) Notice that L_{n} and U_{n} are not individually Markov chains in the repulsive case. In the attractive case where for any i, F_{i}(·|X_{−i}) is decreasing in X_{−i}, we need an upper bound x^{max} on Ω and we redefine the lower and upper chains by setting L_{n}(n) = D(n), U_{n}(n) = x^{max} and interchanging the role of L_{n}(t, i)_{←} and U_{n}(t, i)_{←} when conditioning in (2)-(3). Then both L_{n} and U_{n} become Markovian, and the results in Proposition 1 and Theorem 1 are still valid except that (4) has to be modified so that D becomes a lower dominating chain (correspondingly the role of the minimum 0 and the maximum x^{max} are now interchanged). Propp and Wilson (1996,1997) report on simulation studies of the attractive Ising model (defined on very huge lattices) and its accompanying random cluster model.

Another model is the auto-binomial model. Here D(X_{i}|X_{−i} = x_{−i}) is a binomial distribution with parameters n_{i} and p_{i}(x_{−i}), where log(p_{i}(x_{−i})/(1 − p_{i}(x_{−i}))) = β_{i} + ∑_{j:j≠i} β_{ij}x_{j} and β_{i}, β_{ij} = β_{ji} are real parameters. The model is attractive if all β_{ij} ≥ 0 and repulsive if all β_{ij} ≤ 0. For illustrative purposes, if the β_{ij} have different signs, then we can start the upper processes in x^{max} = (n_{1}, …, n_{k}) and the lower ones in 0, and set

L_{n}(t, i) = F_{i}^{−}(R(t, i)|L_{n}(t, j)1_{[β_{ij}≥0]} + U_{n}(t, j)1_{[β_{ij}≤0]}, 1 ≤ j < i;
L_{n}(t − 1, j)1_{[β_{ij}≥0]} + U_{n}(t − 1, j)1_{[β_{ij}≤0]}, i < j ≤ k),
U_{n}(t, i) = F_{i}^{−}(R(t, i)|L_{n}(t, j)1_{[β_{ij}≤0]} + U_{n}(t, j)1_{[β_{ij}≥0]}, 1 ≤ j < i;
L_{n}(t − 1, j)1_{[β_{ij}≤0]} + U_{n}(t − 1, j)1_{[β_{ij}≥0]}, i < j ≤ k),

when t > n, where 1_{[·]} denotes indicator function.

Yet another interesting model is the ‘pairwise-difference prior’ model with (in the present setup) finite state space Ω = {0, …, m}^{k} and

f_{i}(x_{i}|x_{−i}) ∝ exp(− ∑_{j:j≠i} β_{ij}(x_{i} − x_{j})^{2})

where β_{ij} = β_{ji} are real parameters (see, for example, Besag (1989) and Green (1996)). In applications of image analysis one takes all β_{ij} ≥ 0 so that the model becomes attractive. When doing conventional MCMC forwards simulations from this model, unless m is small enough, it is often more convenient to use another Metropolis-Hastings algorithm than the Gibbs sampler.

Note that if we multiply the pairwise-difference prior density (or any other density with ‘full conditionals’ of the form (8)) with a ‘likelihood term’ ∏_{i} g_{i}(x_{i}) in order to get a posterior density, the monotonicity properties of the prior and the posterior are the same.

(iii) Examples where Ω = {0, 1, …}^{k} is infinite are provided by the auto-Poisson model, where D(X_{i}|X_{−i} = x_{−i}) is a Poisson distribution with mean λ_{i}(x_{−i}) = exp(β_{i} + ∑_{j:j≠i} β_{ij}x_{j}) and parameters β_{i} ∈ R and β_{ij} = β_{ji} ≤ 0, and an auto-negative-binomial model given by that

f_{i}(x_{i}|x_{−i}) = ( \begin{matrix}α_{i} + x_{i} − 1 \\ α_{i} − 1 \end{matrix} ) λ_{i}(x_{−i})^{x_{i}}(1 − λ_{i}(x_{−i}))^{α_{i}}

with the parameter α_{i} > 0 and where now β_{i} < 0. Both models are repulsive; the joint distribution is not well-defined if we allow some β_{ij} to be positive (Besag (1974)), so it is not possible to include the attractive case.

Of course higher order interaction terms may be included in the autobinomial, auto-Poisson and auto-negative-binomial models. For example, we can extend the exponent in the definition of p_{i}(x_{−i}) and λ_{i}(x_{−i}) with the term ∑_{j<k:i∉{j,k}} β_{ijk}x_{j}x_{k}, where β_{ijk} does not depend on the ordering of i, j, k; for the auto-Poisson and auto-negative-binomial models, β_{ijk} ≤ 0.

(iv) Suppose that X_{1}, …, X_{l} are conditionally independent given X_{l+1}, …, X_{k} and that X_{l+1}, …, X_{k} are conditionally independent given X_{1}, …, X_{l}, where 1 ≤ l < k. Then in the repulsive case (1), an alternative perfect simulation algorithm is provided by replacing (2)-(3) by

L_{n}(n) = (0, …, 0, D(n, l + 1), …, D(n, k)), L_{n}(t, i) = F_{i}^{−}(R(t, i)|L_{n}(t, i)_{←}), U_{n}(n) = (D(n, 1), …, D(n, l), 0, …, 0), U_{n}(t, i) = F_{i}^{−}(R(t, i)|U_{n}(t, i)_{←}),

for t > n, whereby L_{n} and U_{n} become individually Markov chains with state space Ω. This perfect simulation algorithm shares the properties in Proposition 1 and Theorem 1 if the partial order ≤ used in (4)-(5) is replaced by another partial order ⪯ defined by

x ⪯ y ⟺ x_{i} ≤ y_{i} for i = 1, …, l, x_{j} ≥ y_{j} for j = l + 1, …, k.

A comparative study of this perfect simulation algorithm, the algorithm in Section 2.1, and certain perfect simulation algorithms based on extensions of Fill (1997) will appear in Møller and Schladitz (1998).

3 Perfect simulation of auto-gamma models

In the sequel we use notation as in Section 2, but it is now assumed that each conditional distribution of the target π = D(X) is given by D(X_{i}|X_{−i} = x_{−i}) = Γ(α_{i}, γ_{i}(x_{−i})), the gamma distribution with shape parameter α_{i} > 0 and inverse scale parameter

γ_{i}(x_{−i}) = β_{i} + ∑_{j:j≠i} β_{ij}x_{j},

where the parameters β_{i} > 0 and β_{ij} = β_{ji} ≥ 0. Thus X = (X_{1}, …, X_{k}) is stochastically dominated by k independent gamma variates with parameters (α_{1}, β_{1}), …, (α_{k}, β_{k}), respectively, and π has density

f(x) ∝ exp(− ∑_{i} β_{i}x_{i} − ∑_{i<j} β_{ij}x_{i}x_{j}) ∏_{1}^{k} x_{i}^{α_{i}−1}  (9)

with respect to Lebesgue measure on Ω = (0, ∞)^{k} (the density is only welldefined for nonnegative β_{ij}). Note that the dependence is expressed through the scale parameter γ_{i} (compare with Cressie (1993, p. 440) who instead considers a dependence through the shape parameter).

The auto-gamma model is also repulsive as F_{i}(·|x_{−i}) increases as x_{−i} increases, so like in Section 2 we can use cyclic Gibbs sampling based on ‘inversion’ when simulating from π. However, it is more obvious to utilize the fact that the conditional distributions are scale models: So let G(t, i) ∼ Γ(α_{i}, 1), i = 1, …, k, t ∈ Z, denote mutually independent gamma variates (for the simulation of the gamma distribution, see, e.g., Ripley, 1987). Then for the Markov chain X(t) = (X(t, 1), …, X(t, k)), t = 0, 1, …, we now set

X(0) = 0, X(t, i) = G(t, i)/γ_{i}(X(t, i)_{←}) for i = 1, …, k, t = 1, 2, …

whereby it is easily verified (see, for example, Roberts and Smith (1994)) that X(t) converges weakly to X.

Apart from one major point, CFTP and the perfect simulation algorithm in Section 3.1 below follow the same lines as in Section 2.1: As we are using Gibbs sampling on continuous distributions we never get coalescence at time 0 of the lower and upper processes (Proposition 2), but we shall establish in Theorem 2 that with an arbitrarily good accuracy we can deliver a perfect simulation from π within finite time. These aspects are further discussed in Section 3.3.

3.1 Perfect simulation algorithm (auto-gamma)

Define the lower and upper processes by

L_{n}(n) = 0, L_{n}(t, i) = G(t, i)/γ_{i}(U_{n}(t, i)_{←}), i = 1, …, k, t > n  (10)

and

U_{n}(n) = D(n), U_{n}(t, i) = G(t, i)/γ_{i}(L_{n}(t, i)_{←}), i = 1, …, k, t > n (11)

where the dominating chain D(t) is defined by the mutually independent gamma variates

D(t, i) = G(t, i)/β_{i}, i = 1, …, k, t ∈ Z.

Let n_{1} > n_{2} > … be a given strictly decreasing sequence of non-positive integers. Then, in the perfect simulation algorithm, for j = 1, 2, … set n = n_{j} and generate (L_{n}(t), U_{n}(t)), t = n, …, 0, until U_{n}(0, i) − L_{n}(0, i) ≤ ϵ, i = 1, …, k; return Z = (L_{n}(0) + U_{n}(0))/2 as a perfect simulation from π with accuracy provably within ϵ of π, where ϵ > 0 is a ‘user-specified parameter’.

3.2 Theoretical results (auto-gamma)

In this section we verify that for any given ϵ > 0 the stopping time

N(ϵ) = N(ϵ, {n_{j}}) = sup{n_{j} : U_{n_{j}}(0, i) − L_{n_{j}}(0, i) ≤ ϵ, i = 1, …, k}

is almost surely finite, and that there is a coupling with a random vector W so that D(W) = π and the output Z = (L_{N(ϵ)}(0) + U_{N(ϵ)}(0))/2 satisfies that

|Z(i) − W (i)| ≤ ϵ, i = 1, …, k.  (12)

This justifies that the ‘accuracy is ϵ’.

Let X_{n}(t) denote the target chain defined by cyclic Gibbs sampling and started in the minimum 0 at time n:

X_{n}(n) = 0, X_{n}(t, i) = G(t, i)/γ_{i}(X_{n}(t, i)_{←}), i = 1, …, k, t > n.

Proposition 2 For all times m ≤ n ≤ u ≤ t we have that the sandwiching properties (4)-(6) remain valid in the present situation. However, with probability one, N(0) = −∞, and if n ≤ s < t, then

[U_{n}(s, i) − L_{n}(s, i) ≤ ϵ] ⇏ [U_{n}(t, i) − L_{n}(t, i) ≤ ϵ].  (13)

Proof: The sandwiching properties are straightforwardly verified by induction. Clearly, P (U_{n}(t, i) > L_{n}(t, i) for all i and n ≤ t) = 1, so N(0) = −∞ almost surely. Since U_{n}(t, i) − L_{n}(t, i) equals

D(t, i)\frac{∑_{j=1}^{i−1} β_{ij}(U_{n}(t, j) − L_{n}(t, j)) + ∑_{j=i+1}^{k} β_{ij}(U_{n}(t − 1, j) − L_{n}(t − 1, j))}{γ_{i}(L_{n}(t, i)_{←})γ_{i}(U_{n}(t, i)_{←})}

where the gamma variate D(t, i) is independent of the fraction, it follows that (13) holds.

Remarks: As N(0) = −∞ we have to take ϵ > 0 in order to get a finite stopping time (as discussed in (iii), Section 3.3, there exist perfect Metropolis-Hastings algorithms with N(0) > −∞). Further, (6) is really not useful anymore, and because of (13) we cannot stop the backward sampling the first time we get intermediate ‘ϵ-coupling’. This is in contrast to the perfect simulation algorithm in the discrete case, but for the auto-gamma model it turns out that only the sandwiching properties (4)-(5) are needed in order to establish (12).

Note that the argument in Propp and Wilson (1996) for preferring the sequence {−2^{j}} still applies in the present situation. Clearly, if we consider a subsequence {m_{j}} of {n_{j}}, then for the corresponding stopping times we have that N(ϵ, {m_{j}}) ≤ N(ϵ, {n_{j}}), so −N(ϵ, {0, −1, …}) is the smallest number of pairs of lower and upper processes needed for obtaining coalescence at time 0. Observe also that as (7) is satisfied, D(−N(ϵ, {0, −1, …})) = D(M(ϵ)) with

M(ϵ) = inf{n ≥ 0 : L_{0}(n, i) = U_{0}(n, i) ≤ ϵ, i = 1, …, k}.

Theorem 2 Let ϵ > 0. With probability one, we have that N(ϵ) > −∞ and the limits

L_{−∞}(t) = lim_{n→−∞} L_{n}(t), U_{−∞}(t) = lim_{n→−∞} U_{n}(t)  (14)

exist and agree. Moreover,

L_{−∞}(t) follows the auto-Gamma model π  (15)

and (12) is satisfied with W = L_{−∞}(0).

Proof: By Proposition 2 and (5), the limits in (14) exist almost surely. Further,

(5) also gives that U_{n}(t) − L_{n}(t) decreases to U_{−∞}(t) − L_{−∞}(t) as n(≤ t) tends to minus infinity. To show that both L_{−∞}(t) = U_{−∞}(t) and N(ϵ) > −∞ hold almost surely, it suffices to show that E(U_{0}(t) − L_{0}(t)) tends to 0 as t → ∞, cf. (7). Hence, letting

l(tk + i) = L_{0}(t, i), u(tk + i) = U_{0}(t, i), d(tk + i) = D(t, i), t ≥ 0

we just have to verify that

lim_{t→∞} E(u(tk + i) − l(tk + i)) = 0.  (16)

Consider now any t ≥ 1 and i ∈ {1, …, k}. Set α = min β_{i} > 0, β = max β_{ij} > 0, and

A(s) = max{D(s, 1), …, D(s, k)}

B(s) = \frac{β(k − 1) max\{A(s − 1), A(s)\}}{α + β(k − 1) max\{A(s − 1), A(s)\}}.

Further, for i, j ∈ {1, …, k} with i ≠ j, let β(i, j) = β_{ij} and set β(i, i − j) = β(i, i − j + k) if i < j. We prove by induction that

(u(tk + i) − l(tk + i))/u(tk + i) ≤ B(1) ⋯ B(t)  (17)

for t ≥ 1 and i = 1, …, k: If t = i = 1, then by (10)-(11),

\frac{u(k + 1) − l(k + 1)}{u(k + 1)} = \frac{∑_{j=1}^{k−1} β(1, 1 − j)u(k + 1 − j)}{β_{1} + ∑_{j=1}^{k−1} β(1, 1 − j)u(k + 1 − j)} ≤ B(1)

where we have used that l(k + 1 − j) = 0, u(·) ≤ d(·), and the real function a → a/(a + b) is increasing for b > 0. So for any t ≥ 1 and i = 1, …, k with tk + i > k + 1 we get by similar arguments and by combining the induction hypothesis with the fact that 0 < B(t) < 1,

\frac{u(tk + i) − l(tk + i)}{u(tk + i)} = ∑_{j=1}^{k−1} \frac{β(i, i − j)u(tk + i − j)}{γ_{i}(u(tk + i)_{←})} \frac{u(tk + i − j) − l(tk + i − j)}{u(tk + i − j)}

≤ B(1) ⋯ B(t − 1)\frac{∑_{j=1}^{k−1} β(i, i − j)u(tk + i − j)}{β_{i} + ∑_{j=1}^{k−1} β(i, i − j)u(tk + i − j)} ≤ B(1) ⋯ B(t)

where we set B(1) ⋯ B(t − 1) = 1 if t = 1. Thereby (17) is verified.

Since d(tk + i) is independent of B(1), …, B(t − 1) and B(t) < 1 we get from (17) that

E(u(tk + i) − l(tk + i)) ≤ (α_{i}/β_{i})E(B(1) ⋯ B(t − 1)).

As the A(s) are iid and B(·) < 1, we have that

E(B(1)B(2) ⋯ B(t − 1)) ≤ E(B(1)B(3) ⋯ B(t − 1)) = (EB(1))^{t/2}

for t even and where EB(1) < 1. So using again that B(·) < 1, we conclude that E(B(1) ⋯ B(t − 1)) → 0 as t → ∞. This implies that (16) holds, and so we have shown that with probability one, L_{−∞}(t) = U_{−∞}(t) and N(ϵ) > −∞.

To verify (15) is similar to the last part in the proof of Theorem 1 by using (4)-(5) and observing that the Markov chain L_{−∞}(t) = U_{−∞}(t) = X_{−∞}(t) is in equilibrium for all t ∈ Z as D(X_{n}(t)) = D(X_{0}(t − n)), t ≥ n, where X_{0}(t − n) converges weakly towards X as n → ∞.

Finally, as L_{n}(0) ≤ W (0) ≤ U_{n}(0) if n ≤ 0, (12) follows from the triangle inequality.

3.3 Comments and examples (auto-gamma)

(i) In computing the state space is of course finite, and in order to avoid technical subtleties we may simply refer to the results for the discrete case as outlined in Section 2. If we take ϵ = 0, the perfect simulation algorithm (Section 3.1) terminates when the lower and upper processes become equal because of rounding error on the machine; this may possibly cause substantial numerical errors so that we are not quite ensured that the accuracy is given by the precision “ϵ_{machine}” of the numbers used in the computations. So it may be safer to keep the size of ϵ several orders of magnitude larger than ϵ_{machine}.

(ii) The main example discussed in Murdoch and Green (1997) is how to do perfect simulations from a posterior distribution for a dataset on pump reliability (Gelfand and Smith (1990)), which is actually an auto-gamma model (9) with k = 11 and pairwise interactions between just one variate and each of the remaining variates. For this model I have produced two C-programs PerfectGammaBackwards.c (the algorithm in Section 3.1 for generating a perfect simulation) and PerfectGammaForwards.c (for forward runs and with output as shown below). These programs are available by anonymous ftp (ftp://ftp.math. auc.dk/pub/jm/); the accuracy ϵ is specified in the beginning of the programs. The average of M(ϵ) based on 10000 simulations for each value of ϵ was 9.3047 (ϵ = 10^{−3}), 11.3170 (ϵ = 10^{−4}), 13.3262 (ϵ = 10^{−5}), 19.3508 (ϵ = 10^{−8}), 31.3775 (ϵ = 10^{−14}), 34.8263 (ϵ = 0); the standard error of the average was 0.0050 (ϵ = 10^{−3}), 0.0052 (ϵ = 10^{−4}), 0.0054 (ϵ = 10^{−5}), 0.0061 (ϵ = 10^{−8}), 0.0072 (ϵ = 10^{−14}), 0.0120 (ϵ = 0).

The CFTP algorithms in Murdoch and Green (1997) use a larger number of paths than just the two paths for the lower and upper processes in the perfect simulation algorithm considered in Section 3.1, so it would not make much sense to use M(ϵ) in a comparative analysis. Duncan Murdoch has suggested to me that a hybrid algorithm could possibly do better than either theirs or mine, e.g. by first using my algorithm so that the lower and upper processes become close enough and then using their multigamma coupler.

(iii) It is at least theoretically feasible to make simulations with accuracy 0 by the following perfect Metropolis-Hastings algorithm by noticing that the density (9) of the auto-gamma model (as well as the auto-binomial, auto-Poisson, auto-negative-binomial and many other models) is of a particular form:

f(x) ∝ h(x) ∏_{1}^{k} q_{i}(x_{i})

where (in the case of the auto-gamma model) q_{i} is the density of Γ(α_{i}, β_{i}) and

h(x) = exp(− ∑_{i<j} β_{ij}x_{i}x_{j}).

If we are using cyclic updates and the usual notation as introduced in Section 2, then it is now a Metropolis-Hastings algorithm which generates a Markov chain X(t), t ≥ 0: At time t, when the ith coordinate is to be updated and the current state is given by X(t − 1, i) and X(t, i)_{←}, first a proposal D(t, i) is generated from q_{i} together with a uniform number R(t, i) (between 0 and 1), and secondly we set X(t, i) = D(t, i) if R(t, i) ≤ min{1, a_{i}(D(t, i), X(t − 1, i), X(t, i)_{←})} and we retain X(t, i) = X(t − 1, i) otherwise (where all gamma variates D(t, i) and random numbers R(t, i) are assumed to be mutually independent). Here

a_{i}(d_{i}, x_{i}, x_{−i}) = \frac{f(x_{1}, …, x_{i−1}, d_{i}, x_{i+1}, …, x_{k})q_{i}(x_{i})}{f(x_{1}, …, x_{i−1}, x_{i}, x_{i+1}, …, x_{k})q_{i}(d_{i})} = exp((x_{i}−d_{i}) ∑_{j:j≠i} β_{ij}x_{j})

is the Metropolis-Hastings ratio (Hastings, (1970)). This is larger than 1 if d_{i} ≤ x_{i} (i.e. moving down is always accepted), whilst it decreases from 1 to 0 as a function of x_{−i} when d_{i} > x_{i}. It is easily seen that X(t) converges weakly to X.

CFTP should now be obvious: Set L_{n}(n) = 0, U_{n}(n) = D(n), and for t > n set

L_{n}(t, i) = D(t, i) if R(t, i) ≤ a_{i}(D(t, i), L_{n}(t − 1, i), U_{n}(t, i)_{←})

and L_{n}(t, i) = L_{n}(t − 1, i) else, and

U_{n}(t, i) = D(t, i) if R(t, i) ≤ a_{i}(D(t, i), U_{n}(t − 1, i), L_{n}(t, i)_{←})

and U_{n}(t, i) = U_{n}(t−1, i) else. Clearly, the usual sandwiching properties remain valid with D(t) as the dominating chain, and for a perfect Metropolis-Hastings algorithm similar to the algorithm in Section 2.1 it is not hard to verify that the coalescence time is almost surely finite and the output follows the target distribution. (Sketch of proof: Let ϵ > 0. Conditioning on the event that D(t − 1, 1) ≤ ϵ, …, D(t, k) ≤ ϵ, then for t > n the conditional probability of the event L_{n}(t) = D(t) is greater than δ(ϵ)^{k}, where δ(ϵ) = inf{a_{i}(d_{i}, x_{i}, x_{−i}) : 0 ≤ x_{i} ≤ d_{i} ≤ ϵ, 0 ≤ x_{j} ≤ ϵ, i ≠ j} > 0. From this we get that the coalescence time is almost surely finite; that the output follows the target distribution is verified in the same way as in the proofs of Theorems 1 and 2.)

However, in practise (even for the simple example considered in (ii) above) this perfect Metropolis-Hastings algorithm may become extremely slow as the lower processes get stuck in 0 for a very long time unless the parameters β_{ij}/(β_{i}β_{j}), i < j, are sufficiently small or equivalently if the interactions are sufficiently weak. One may object that this is clearly due to the way proposals are generated. But we need a coupling construction, where it is feasible to generate a dominating chain, which is in equilibrium. Actually we are using the smallest possible dominating chain as this has both to dominate the proposals and the Markov chain X_{−∞}(t) started ‘in equilibrium at time minus infinity’ and generated by the Metropolis-Hastings algorithm (see the proofs of Theorems 1 and 2).

(iv) The coupling construction in Section 3.1 and the proofs of Proposition 2 and Theorem 2 are essentially only based on a few properties: (a) The ‘local characteristics’ D(X_{i}|X_{−i} = x_{−i}) are scale families, where the inverse scaling factor γ_{i}(x_{−i}) is an increasing function of x_{−i} and the support is the positive half-line. Thereby lower and upper processes are naturally constructed (Section 3.1). (b) There is a (natural) dominating and stationary Markov chain D(t) (in the proof of Theorem 2 we referred to that the D(t, i) are iid, but this may of course be weakened). (c) ED(t) is finite. Apart from (a)-(c) the fact that D(X_{i}|X_{−i} = x_{−i}) is a gamma distribution was really never used.

Certain models can be transformed into a auto-gamma model such as autogeneralized-gamma and hierarchical models obtained by conditional independent normally distributed variates with the precisions (inverse variances) given by an auto-gamma model. It would be interesting to see if the ideas in Sections 3.1-3.2 become useful for other multivariate continuous distributions.

Acknowledgement: This research was supported by MaPhySto – Centre for Mathematical Physics and Stochastics, funded by a grant from the Danish National Research Foundation – and by the Danish Informatics Network in the Agricultural Sciences, funded by a grant from the Danish Research Councils. Many thanks to my colleagues in Aalborg, Martin B. Hansen and Uffe Kjærulff for programming assistance ((ii) in Section 3.3) and Bjarne Højgaard for stimulating discussions. Also the encouragement and helpful comments of Julian Besag, Olle Häggström, Duncan Murdoch, and James Propp are gratefully acknowledged.

References

[1] Besag, J.E. (1974). Spatial interaction and the statistical analysis of lattice systems. Journal of Royal Statistical Society B 36, 192–225.

[2] Besag, J.E. (1989). Towards Bayesian image analysis. Applied Statistics 16, 395–407.

[3] Cressie, N. (1993). Statistics for Spatial Data, revised edn. New York: Wiley.

[4] Fill, J.A. (1997). An interruptible algorithm algorithm for perfect sampling via Markov chains. Preprint, Department of Mathematical Sciences, The John Hopkins University. To appear in Annals of Applied Probability.

[5] Foss, S.G. and Tweedie, R.L. (1997). Perfect simulation and backward coupling. Stochastic Models. (To appear)

[6] Gelfand, A.E. and Smith, A.F.M. (1990). Sampling-based approaches to calculating marginal densities. Journal of the American Statistical Association 85, 398–409.

[7] Green, P.J. (1996). MCMC in image analysis. Markov Chain Monte Carlo in Practice (eds. W.R. Gilks, S. Richardson and D.J. Spiegelhalter), 381– 399. London: Chapman and Hall.

[8] Häggström, O., van Lieshout, M.N.M. and Møller, J. (1996). Characterisation results and Markov chain Monte Carlo algorithms including exact simulation for some spatial point processes. Research Report R-96-2040, Department of Mathematics, Aalborg University. To appear in Bernoulli.

[9] Häggström, O. and Nelander, K. (1997). Exact sampling from antimonotone systems. Research Report 1997–03, Department of Mathematics, Chalmers University of Technology and Göteborg University.

[10] Hastings, W.K. (1970). Monte Carlo sampling methods using Markov chains and their applications. Biometrika 57, 97–109.

[11] Kendall, W.S. (1996). Perfect simulation for the area-interaction point process. Research Report 292, Department of Statistics, University of Warwick. To appear in Probability Perspective (eds. C.C Heyde and L. Accardi). Singapore: World Scientific Press.

[12] Kendall, W.S. (1997). On some weighted Boolean models. To appear in Advances in Theory and Applications of Random sets (ed. D. Jeulin), 105– 120. Singapore: World Scientific Publishing Company.

[13] Kendall, W.S. and Møller, J. (1997). Perfect Metropolis-Hastings simulation of locally stable point processes. (In preparation)

[14] Meyn, S.P. and Tweedie, R.L. (1993). Markov Chains and Stochastic Stability. London: Springer-Verlag.

[15] Murdoch, D.J. and Green, P.J. (1997). Exact sampling from a continuous state space. Research Report S–97–01, Department of Mathematics, University of Bristol. To appear in Scandinavian Journal of Statistics.

[16] Møller, J. and Schladitz, K. (1998). Extensions of Fill’s algorithm for perfect simulation. (In preparation)

[17] Propp, J.G. and Wilson, D.B. (1996). Exact sampling with coupled Markov chains and applications to statistical mechanics. Random Structures and Algorithms 9, 223–252.

[18] Propp, J.G. and Wilson, D.B. (1997). How to get a perfectly random sample from a generic Markov chain and generate a random spanning tree of a directed graph. Journal of Algorithms. (To appear)

[19] Ripley, B.D. (1987). Stochastic Simulation New York: Wiley.

[20] Roberts, G.O. and Smith, A.F.M. (1994). Simple conditions for the convergence of the Gibbs sampler and the Metropolis-Hastings algorithms. Stochastic Processes and their Applications. 49, 207–216.