Perfect simulation for correlated Poisson random variables conditioned to be positive.
In this paper we present a perfect simulation method for obtaining perfect samples from collections of correlated Poisson random variables conditioned to be positive. We show how to use this method to produce a perfect sample from a Boolean model conditioned to cover a set of points: in W.S. Kendall and E. Thnnes (Pattern Recognition 32(9): 15691586, 1999), this special case was treated in a more complicated way. The method is applied to several simple examples where exact calculations can be made, so as to check correctness of the program using &khgr;2-tests, and some small-scale experiments are carried out to explore the behaviour of the conditioned Boolean model.
Introduction
Markov chain Monte Carlo (McMC) provides the opportunity to carry out simulation-based
statistical inference, by allowing one to sample approximately from a target distribution
realised as the equilibrium of a Markov chain. The approximate nature of
McMC arises because the output is subject to initialization bias. Recent work on perfect
simulation [4, 13] shows that in favourable circumstances McMC simulation algorithms
can be modified so as to replace the approximation by an exact simulation
procedure.
In particular the Coupling from the Past (CFTP) method due to Propp & Wilson
[13] has been developed intensively since its discovery in 1994. For example, dominated
CFTP is a variant of CFTP developed in [8], which permits perfect simulation
of a general class of point processes (see also [10]). Dominated CFTP has an importance
going beyond point processes, since it allows application of CFTP to a wide
class of non-uniformly ergodic Markov chains which need not be monotonic. Kendall
Thonnes [11] used dominated CFTP to carry out perfect simulation in stochastic ge-
ometry, showing how to get perfect samples from Boolean models conditioned to cover
a finite set of points.
The purpose of this paper is to explain how to get a perfect sample from collections
of correlated Poisson random variables conditioned to be positive (briefly, correlated
Poisson random variables arise as sums of a basic set of independent Poisson random
variables). Following Kendall & Thonnes [11], the motivating application is to obtain
a perfect sample from a Boolean model conditioned to cover a set of points. However
the construction in [11] is very involved. Conditional Boolean model can be viewed
as arising from correlated conditioned Poisson random variables, and this leads to a
simpler simulation algorithm; the more general random variables set-up draws attention
away from the geometry (which here is a distracting complexity). As well as being of
interest to stochastic geometers wishing to investigate the conditional Boolean model,
the results reported on here add to a sequence of case-studies [6, 7, 8, 11, 10] aimed
at developing CFTP by using it in challenge-problems from stochastic geometry. The
particular challenge of perfect simulation of conditional Boolean models has resulted
in the formulation of a rather general context for CFTP described as "extended state-space
CFTP" below.
The construction of the paper is as follows. In x2 we introduce the notion of correlated
Poisson random variables, describe a McMC method for obtaining an approximate
sample from the conditional correlated Poisson distribution, and then specify
an algorithmic scheme for obtaining a perfect sample, in the context of extended state-space
CFTP. In common with most other practical CFTP algorithms this is based on the
construction of interleaved sequences of upper- and lower-sandwich processes, bounding
intermediate Markov processes which actually have the required equilibrium distri-
bution. Fundamental properties of our CFTP algorithm are established by Theorem 1
in x3: namely that the algorithm will terminate in finite time and that the output has the
desired distribution. In x4 we describe the relationship between conditional Boolean
models and collections of correlated Poisson random variables conditioned to be positive
and thus show how to use our algorithm to conduct CFTP for the conditional
Boolean model. Several examples are described in x5. For these examples, we derive
explicit conditional distribution functions of the correlated Poisson random variables,
enabling us to carry out  2 -tests in order to check the correctness of the CFTP algo-
rithm. Relationships between the intensity of the germ process, radius of the typical
grain, number of conditioning points and the coalescence time are investigated in x6.
Finally, some comments and conclusions are presented in x7.
Correlated Poisson random variables
conditioned to be positive
We begin this section by introducing the concept of collections of correlated Poisson
random variables. We then describe a McMC method for obtaining an approximated
sample from the collection when the individual random variables are all conditioned
to be positive. Finally we describe the CFTP modification which allows us to obtain a
perfect sample from the desired distribution.
We fix a given index set and assign a rate A to each non-empty
subset A  S. Consider a collection of independent Poisson random variables fYA :
A  Sg indexed by these subsets, such that YA has mean A , and define
is a collection of correlated Poisson random variables. We also
define
and set
abbreviating
As described in x4, such collections arise naturally in the context of the Boolean
model in stochastic geometry, when the index set replaced by a finite
set of conditioning points and YB counts the number of grains covering
exactly the subset B of conditioning points and no others, while X fx i g counts
the number of grains covering the conditioning point x i . In this context, motivated
by problems from the oil industry, Lantuejoul [12] has investigated simulation of such
Boolean models conditioned to cover the conditioning points: the question of converting
this to perfect simulation, as discussed in [11], leads directly to the question of how
to conduct perfect simulation of correlated Poisson random variables conditioned to be
positive.
2.1 McMC for collection of correlated Poisson
random variables conditioned to be positive
To carry out McMC for fX we replace random variables by immigration-
death processes. For each non-empty subset A  S we define an immigration-death
process YA of immigration rate A and unit death rate, requiring the YA to be indepen-
dent. By detailed balance YA has Poisson equilibrium of mean A . Now we define XA
for any non-empty subset A  S by
Abbreviating that the equilibrium distribution for fX
is a collection of correlated Poisson random variables. The McMC procedure is now
straightforward: simulate the YA processes for a long time, then construct the X i using
Equation (1). Of course this is using a sledgehammer to crack a nut: there is no reason
not to sample directly from the YA equilibria. However the McMC approach can easily
be modified to allow for conditioning.
Considered alone, the process XA is an immigration-death process of immigration
rate A and unit death rate. Moreover
where XA e
XB means that
if there is an immigration in XA at time t then there is an immigration in XB at
that time,
and if there is a death in XB at time t then there is a death in XA at that time.
It is possible to establish a converse to the above as follows: given a collection of
immigration-death processes fXA : A  Sg satisfying the relation (2), and co-adapted
in the sense that the rate-computations are based on conditioning on a common past,
then independent immigration-death processes YA of immigration rate A and unit
death rate can be constructed such that Equation (1) applies, using the Mobius formula
[1, chapter 5, page 53] to deduce
A:BA
We now consider how to modify the above McMC procedure so as to produce
correlated Poisson random variables conditioned to be positive. Since the underlying
Markov chain is reversible, all that has to be done is to forbid transitions which lead
to one or more of the X i becoming zero. To do this in a way which maintains the
coupling required for CFTP methods, we need to prescribe replacement death times for
the potential death times which are thus prevented from becoming actual; accordingly
we introduce independent unit-rate Poisson processes ZA , one for each non-empty
subset A  S.
We define the immigration-death process e
YA to be based on the same immigration-
death structure as YA , but forbidding deaths which would lead to one of the target
random variables becoming zero. Namely, if t is a death time of YA such
that then the death time t of YA will not be accepted by e
YA ,
and we set e
Furthermore, we use the next
incident ^ t of ZA as a new potential death time of e
YA . We call this perpetuation of e
Y
from time t to ^ t. The interval [t; ^ t) is a perpetuation interval, and ^ t is a perpetuation
time of e
YA . Summarizing, if e
experiences a death at t, and if there
exists A such that
e
A such that
then perpetuation takes place. In particular, perpetuation increases e
YA by at most 1
compared with the original YA .
Since we use independent unit rate Poisson processes ZA (t) to provide
perpetuation times, e
YA inherits the unit death rate of YA subject to censorship
of those transitions which would lead to X i
Now if we define
e
A:BA
e
e
then detailed balance shows that the equilibrium distribution of
e
is that of a collection of correlated Poisson random variables conditioned to be positive.
The pseudo-code McMC-CCP displayed in Table 1 gives an explicit description of
how we can derive a simulation of e
X over a time interval [a; b] from a realization of
over the same time
interval.
initialize e
repeat
e
A:i2A
e
YA
A index of next incident
if next incident is immigration then
e
if next incident is death then
if e
e
else
label next ZA-incident as death
until no more incidents
return
The algorithm McMC-CCP performs McMC to obtain an approximate sample of conditioned
correlated Poisson random variables, using unconditioned immigration-death
processes fYA g and associated perpetuation Poisson processes fZA g.

Table

1: McMC-CCP algorithm
2.2 Extended state-space CFTP in a general context
Dominated CFTP is presented in [10] in a general context. In this paper we require
a further modest extension, extended state-space CFTP. Here we describe this in a
general context, then describe how to use it to get perfect samples from a correlated
Poisson random variables conditioned to be positive.
Our aim is to obtain a perfect sample from the equilibrium distribution  of a discrete
or continuous time Markov chain taking values in a state space
X. The difficulties here are that the Markov chain R may not be uniformly ergodic, and
even if a relevant partial order exists we do not wish to assume either monotonicity or
the existence of maximal and minimal elements. The main idea to overcome this problem
is: insert the original unordered state-space X as a subset of a bigger ordered space
e
X, with an "upper/lower-sandwich" process living on e
X and eventually being absorbed
in X. The evolution of the stochastic process e
R is chosen so that once it visits X it
follows the stochastic dynamics of R.
Specifically, we identify a sequence of majorizing processes
e
R
; such that
which are identically distributed up to a shift in time, which live on the superset e
X of
X and satisfying the following conditions:
1. e
X possesses a partial order relation , and if a 2 e
X are such that
a  b then a = b (so X is at the base of the partially ordered space e
2. e
R (n) is eventually absorbed in the subspace X, so that
for any fixed u as n !1;
3. e
R (n) evolves according to the stochastic dynamics of R once it hits X;
4. For all m  n  t  0
e
R (m) (t)  e
Then we have the following CFTP result.
Theorem 1 Suppose that the Markov chain R takes values in a state space X, and the
R (n) live on the partially ordered superset e X of X and satisfy conditions 1-4.
Suppose further that R(t) converges weakly to a unique equilibrium distribution  as
R (n)
almost surely and e
the equilibrium distribution .
Proof : It follows from Equation (4) in condition 2 above that T < 1 almost surely.
By Equation (5) of condition 4, if n > T then
e
R (n) (0)  e
Consequently by condition 1 above e
R (n)
R (T ) (0) for such n. Accordingly
e
e
R (n) (0)
exists almost surely. By condition 3, and the fact that the e
R (n) are identically distributed
up to time shifts, e
R (n) (0) has the same equilibrium distribution as that of
R(n), so
lim
R (n)
where R is started at time 0 with the common hitting distribution of the e
R (n) on X. But
converges weakly to . Therefore, e
R (T ) (0) has distribution  as required. 2
In the following subsection we describe how to construct e
R (n) as an ordered triple
of processes, and in the next section we show that conditions 1-4 hold, so that we can
use e
R (n) to get a perfect sample from the desired distribution.
2.3 CFTP for correlated Poisson random variables
conditioned to be positive
We need to introduce some further notation in order to explain howwe can use extended
state-space CFTP to produce perfect samples of correlated Poisson random variables
conditioned to be positive.
X as follows
are n-tuples of nonnegative integersg ;
e
are n-tuples of nonnegative integers
with v min  v max
where the n-tuples are all of length 2 k 1 indexed by the non-empty subsets of S, the
inequality v min  v max is interpreted element-wise, and in fact v, v min , v max only
contain entries which are 0 or 1. The embedding X  e
X is given by
The partial ordering on e
X is given by
exactly when y interpreting
equalities and inequalities element-wise.
This ordering on e
corresponds exactly to the inclusion ordering for order intervals
which was used in set-valued CFTP and, implicitly, in dominated CFTP in earlier formulations
[6, 7, 8].
It is apparent that this partial ordering and embedding of X in e
1 preceding Theorem 1.
We now describe how to use the unconditional immigration-death processes (Y A :
non-empty A  S) (assumed to be in statistical equilibrium) and the perpetuating
incident processes (ZA : non-empty A  S) to construct realizations of e
R (n) over the
time interval [ n; 0]. Note that we can use time-reversibility of the YA to extend their
trajectories as far back into the past as may be required, fixing their values at time 0 to
be samples from the unconditioned Poisson equilibrium of these processes. Extension
of the ZA is even easier!
For
e
R (n)
where YA is the unconditional immigration-death process assigned to the subset A,
and the V (n);min , V (n);max processes are defined together recursively in the following
way.
Initially, at time n, we set
for all non-empty A  S. Updates of these processes occur only at death times of
corresponding Y processes and perpetuation times of corresponding Z processes.
For convenience of notation we define
e
e
Suppose YA experiences a death at time
A , respectively
A , may be set equal to 1 when there is a danger of causing zeros in some of the
opposite processes e
respectively e
[min-death] If e
so that after the YA -death there
would be a zero in some upper process e
then we set the lower component
[max-death] If e
so that after the YA -death there
would be a zero in some lower process e
then we set the upper component
On the other hand suppose that ZA experiences a perpetuation incident at time t.
Then we can set V (n);min
A , respectively V (n);max
A , equal to 0 when there is no danger
of causing zeros in any of the opposite processes e
respectively e
[min-perpetuation] If e
so that there is no danger
of zeros in any upper process e
then we set the lower component
[max-perpetuation] We need to consider the following two cases.
1. In this case, if e
we set the upper component V (n);max
1. In this case, if e
1 for all i 2 A, then we set the upper component V (n);max
It remains to establish the conditions 1-4 required by Theorem 1, and we shall
do that in the next section. First, we give a brief description of our CFTP algo-
rithm. It begins by testing for whether e
R (n) (0) lies in X, which is to say, whether
A (0) for all non-empty A  S. If so then the vector of common
values X (n);min
is returned as a sample
of the corresponding vector of correlated Poisson random variables conditioned to be
positive. If there exist A  S such that V (n);min
A (0) 6= V (n);max
A (0) then we must
choose a new starting time m < n and repeat the simulation replacing n by m (as
suggested in [13], in this case a sensible choice is to double n). Note that we must
insist on extending backwards in such a way that we re-use the previous realization of
as we noted above, the reversibility
arising from detailed balance is of technical assistance here. We proceed to extend
backwards in time until minimal and maximal processes coalesce, and then return the
common value of the coalesced processes at time 0. Table 2 gives pseudo-code summarizing
the resulting algorithm CFTP-CCP.
3 Properties of the CFTP algorithm
To prove that our CFTP algorithm CFTP-CCP will indeed provide us with perfect
samples from the desired distribution, we must establish a sequence of lemmas relating
to the conditions 1-4.
extend Y , Z to [ n; 0]
if e
R (n)
e
A (0))
return
else
The CFTP-CCP algorithm performs extended state-space CFTP to obtain an exact
sample of conditioned correlated Poisson random variables.

Table

2: CFTP-CCP algorithm
First we must make sure that condition 1 applies, which is to say, the process e
R (n)
genuinely stays in the state-space e
X.
defined as above we have
A (t)  V (n);max
A (t) (8)
for all non-empty A  S and for n  t  0. In particular, the state-space for e
R (n)
can be taken to be e
X.
First note that by construction Equation (8) must hold at time n.
Now suppose that Equation (8) holds up to, but possibly not including, time t.
Therefore we have
e
for all i 2 S.
It suffices to consider t such that either YA experiences a death at that time, or ZA
delivers a perpetuation incident. In the case of a death the relationship (9) and the
crossover between [min-death] and [max-death] assure us that V (n);min can be set to
1 only if V (n);max is, so Equation (8) must also hold at time t.
The case of a perpetuation incident is entirely similar: If the perpetuation happens
in both upper and lower processes then V (n);max can be set to 0 only if V (n);min is,
so Equation (8) must also hold at time t. If the perpetuation happens only in the upper
process, we already have V (n);min holds no matter what is the value of
(t). The result follows by induction. 2
Setting aside condition 2 for the moment, we turn to condition 3.
persists over all
[t; 0] and over that time interval the coupled process Y
evolves as a copy of e
Y .
Proof : This is simply a matter of observing that if V (n);min , V (n);max coincide
at time t then so do e
and so there is no longer any distinction
between min and max in the conditions [min-death], [max-death] for death and [min-
perpetuation], [max-perpetuation] for perpetuation above, note that the case of perpetuation
in the upper process only will not happen any more, while birth decisions are
already identical. Consequently from time t onwards the death and perpetuation decisions
are the same for V (n);min , V (n);max so they must continue to coincide. Furthermore
the effect of the death and perpetuation decisions on Y +V (n);min , Y +V (n);max
agrees with the prescription for the evolution of e
Y in our description of the McMC-CCP
Algorithm in Table 1. 2
Now consider condition 4. This follows if we can establish the following "fun-
nelling" result:
Lemma 3 For all n  m  t  0 and all non-empty A  S we have
A (t)  V (n);min
A (t)  V (n);max
A (t)  V (m);max
Proof : The initial values V (m);min ( m), V (m);max ( m) are respectively minimal
and maximal, so this funnelling certainly holds at time m. Consider therefore
an incident at time t > m, and suppose the funnelling relationship holds at time t .
This carries through into a relationship
e
If the incident is an immigration then all of V (m);min , V (m);max , V (n);min and V (n);max
are not affected, and the relationship persists.
If the incident is a death then examination of [min-death] shows that V (m);min can
be set to 1 only if V (n);min is set to 1, since e
similar arguments using [max-death] apply to show V (n);max can be set to 1 only if
V (m);max is set to 1.
If the incident is a perpetuation in both upper and lower processes then similar
arguments concerning [min-perpetuation] and [max-perpetuation] show that V (n);min
can be set to 0 only if V (m);min is set to 0, V (m);max can be set to 0 only if V (n);max
is set to 0.
Now we consider the case, where the incident is a perpetuation only in the upper
process. In this case, we have
and
Then it follows from Lemma 1 that
A (t)  V (n);min
A (t)  V (n);max
So we only need to prove
A (t)  V (m);max
If there is i 2 A such that e
according to [max-perpetuation],
A (t) can not be set to 0. So V (m);max
holds.
If e
then we set V (m);max
and we have
e
for all i 2 A. There are two cases to be considered:
In this case, since V (m);min
actually know
e
and hence also e
A. Then it follows from
[max-perpetuation] and [min-perpetuation] conditions that we have
In this case, we have either V (n);max
e
it follows from [max-perpetuation] that V (n);max
is automatic since t is a
perpetuation time not a death time.
So in any case, (10) holds.
This means the funnelling relationship persists to time t if it has been maintained
up to time t . The result follows by induction along the sequence of incidents over the
Finally we turn to condition 2. This is a consequence of the following lemma.
Lemma 4 For all sufficiently large n, the processes V (n);min , V (n);max coalesce at
some time in the interval [ n; 0].
Proof : The lemma follows if there are times at which the upper processes V (n);max
A
all vanish. For by construction it then follows coalescence occurs at that time.
Fix T < t < 0,  > 0 and for
experiences an immigration in the time interval [t;
Then the D i are increasing events based on independent events of the form that YA
has an immigration in [t; t the FKG inequality [5, Theorem 2.4] applies and
consequently
Y
On the other hand there is  > 0 such that P [D i ] >  for all i. It follows that
alternative line of reasoning constructs an estimate based on all
of the YA experiencing immigration, thus dispensing with the FKG inequality.)
Now suppose that at time t +  the event
occurs. Let S t be the family
of subsets A  S such that YA experiences at least one immigration in [t; t
occurs then the union
A must cover S, moreover clearly S t has no
more than 2 k 1 members.
Now there are at most 2 k 1 perpetuated values among V (T );max
A , for all A  S.
Coalescence of V (T );min
A and V (T );max
A will occur if, for all A yielding perpetuated val-
ues, both the Poisson process ZA experiences incidents in the time interval [t+; t+2]
and also the YA does not hit zero in the same time interval, since then all the V (T );max
A
will be zero. But a simple stochastic comparison, using the Strong Markov property
and the representation of immigration-death processes in terms of Poisson inputs and
exponential death times, shows that this happens with a conditional probability (condi-
tioned on
Let C t be the event that
occurs together with the event whose probability is
given in Equation (11) above. Then
However C t is actually independent of events in the past of time t. We may therefore
apply the second Borel-Cantelli lemma to show that, almost surely, infinitely many of
the C 2r occur (for It follows that coalescence is almost sure to happen.The following result is an immediate consequence of the above lemmas and Theorem
1 and the existence of a unique limiting equilibrium for processes e
X in x2.
Theorem 2 The algorithm described in x2 almost surely terminates in finite time,
moreover the output of the algorithm has the distribution of the corresponding array of
correlated Poisson random variables conditioned to be positive.
The above work shows that our CFTP algorithm does indeed produce a perfect sample
of the corresponding collection of correlated Poisson random variables conditioned
to be positive.
4 Application to simulation of the conditional Boolean
model
In this section, we show how our CFTP algorithm can be used to obtain a perfect
sample from a conditional Boolean model.
4.1 Relationship between the conditional Boolean model and
correlated Poisson random variables conditioned to be positive
Recall that in section 2, we introduced the distribution of correlated Poisson random
variables showed it could be obtained as the equilibrium distribution
of the process (X 1
A:i2A
YA (t); for
for independent immigration-death processes YA (t)
Now if we consider a Boolean model as an equilibrium distribution for a spatial
immigration-death process then with appropriate assignation of conditioning points
rates A , the process XB (t) can be interpreted as the time evolving
number of germ-grain pairs whose grains cover all the conditioning points in B, while
the process YA (t) describes the number of germ-grain pairs whose grains cover all
conditioning points in A, but no conditioning point which is not in A. For example,
counts the number of germ-grain
pairs whose grains cover the conditioning point x 1 . We divide this process XB (t) into
counts the number of germ-grain
pairs whose grains cover x 1 , but not x 2 , and Y fx1 ;x2 g (t) counts the number of germ-
grain pairs whose grains cover both x 1 and x 2 . If we convert the Boolean model into a
Poisson process on a space marked using the grains, then the space of germs leading to
grains covering x 1 is divided into two disjoint subregions; Y fx1g (t) counts the number
of germs which fall in one of the subregions, and Y fx1 ;x2 g (t) counts the other. Thus
the number of germ-grain pairs in the subregion associated with A is an immigration-
death process YA (t), independent of other such immigration-death processes until the
conditioning is applied by means of forbidding certain transitions.
We have noted that some of the processes YA (t) might have zero immigration rate.
In the conditional Boolean model, this corresponds to a set A of conditioning points
which can not simultaneously be covered by a single germ-grain pair. In the case of
grains which are disks of constant radius, there will then be at least two points in A
which are separated by a distance of greater than 2r, where r is the radius of a typical
grain.
From the above it can be seen how we can use perfect samples from conditional
correlated Poisson random variables to produce perfect samples from a conditional
Boolean model.
4.2 Perfect simulation for the conditional Boolean model
For the sake of simplicity we consider only the case of disks of constant radius, though
the method can be extended to cover the case when the typical grain is random but of
bounded diameter.
Perfect simulation follows once one establishes how to simulate sorted incident
times of all YA (t) and ZA (t) over [ In essence all one has to do is to simulate
a stream of candidate incident times of time-varying intensity, then for each incident
time choose whether it is a Z-type or a Y -type incident using appropriate probabilities.
The time-varying rate of the underlying stream of candidate incidents should be
taken to be
A
YA
where N is the total number of processes ZA with A 6= 0, and assignations of type
should be made so as to deliver corresponding rates
A
For Z-type incidents, one assigns the type (indicating to which of the various ZA
processes the incident belongs) uniformly at random. For Y -type death incidents one
assigns the type (indicating to which of the various YA processes the incident belongs)
randomly but weighted by the sizes of the corresponding YA . For Y -type birth inci-
dents, one chooses a location uniformly at random from the bounded window W and
then assigns the incident to YA if this location falls in the region
x
G
y
GA (12)
where
G is a typical grain (since then the corresponding
disk will cover exactly those conditioning points constituting the indexing set A).
Birth incidents can be discarded if their locations do not fall in any of these regions.
Perfect simulation is then conducted according to the algorithm CFTP-CCP as
given in Table 2. At the conclusion of the algorithm the perfectly sampled conditional
Boolean model is built up using the locations of surviving birth incidents, supplemented
by a sample from the unconditional Boolean model on the region
y
which by construction covers no conditioning points. (This can be done easily by
thinning the unconditional Boolean model on the full region W .)
In the following, we consider plots of conditional Boolean models obtained from
the above perfect simulation method. Clearly, given the configuration of conditioning
points and the radius of a typical grain, the results will be affected by the value of the
intensity. We give plots for simulations using varying patterns of conditioning points
(random pattern, oblique lattice, square lattice) and medium and high intensities of
the underlying Poisson germ process. For low intensity (less than germs per unit
area) and closely-spaced conditioning points (such that a grain can cover more than one
conditioning point) the algorithm failed to coalesce within the memory constraints of
our computer (Sun Ultra 1, 64 MBytes RAM); the maximum and minimum processes
then take a long time to come to agreement, since initially the perpetuated maximum
grains disappear only when minimum grains cover conditioning points (itself a rare
event), and moreover perpetuated minimum grains tend to disappear quickly back to
zero since conditioning points are typically covered by several perpetuated maximum
grains.
Our plots are based on intensities of germs per unit
area. The same general features are displayed in each case: there are many extra grains
induced by the requirement to cover all conditioning points. For high intensity the
resulting figure differs less from the unconditional case, since much of the window
is covered anyway, but there is still an excess of grains covering conditioning points.
In each case we compute the number of grains covering the conditioning points and
the excess of this over the number to be expected for a comparable simulation of an
unconditional Boolean model. The excess is much larger in the case of a random
pattern of conditioning points; this is reasonable since a randomly scattered pattern
affects a much greater area of the window. The excess does not appear to depend
strongly on the underlying germ intensity. This comes as no surprise. The excess,
when generated by our algorithm, is merely the sum of all perpetuated grains:
in the notation of section 2. Since the VA are only over 0 or 1, the scope for variation
of the total excess is limited.
Note that in all the figures conditioning points are marked as "+".
20.20.611.41.8Perfect sample from a conditional Boolean model with
placed randomly in a window W of width units.
- Intensity of underlying germ process
- Radius of typical grain
Simulation yields grains, of which 103 cover at least one conditioning point.
The unconditional model would yield a Poisson number of grains, mean number 120,
of which we would expect 47:41 to cover at least one of the
points. So the excess number of grains from the conditional Boolean model covering
the conditioning points is 55:59.

Figure

1: Conditional Boolean model, random pattern
placed randomly in a window W of width units.
- Intensity of underlying germ process
- Radius of typical grain
Simulation yields grains, of which 131 cover at least one conditioning point.
The unconditional model would yield a Poisson number of grains, mean number 180,
of which we would expect 71:11 to cover at least one of the
points. So the excess number of grains from the conditional Boolean model covering
the conditioning points is 59:89.

Figure

2: Conditional Boolean model, random pattern
points placed on a lattice based on a parallelogram of horizontal
side length 2=19, height 2=19 and angle 63:43 o in a window W of width
units.
- Intensity of underlying germ process
- Radius of the typical grain
Simulation yields grains, of which 71 cover at least one conditioning point.
The unconditional model would yield a Poisson number of grains, mean number 120,
of which we would expect 47:41 to cover at least one of the
points. So the excess number of grains from the conditional Boolean model covering
the conditioning points is 23:59.

Figure

3: Conditional Boolean model, parallelogram grid
points placed on a lattice based on a parallelogram of horizontal
side length 2=19, height 2=19 and angle 63:43 o in a window W of width
units.
- Intensity of underlying germ process
- Radius of the typical grain
Simulation yields grains, of which 97 cover at least one conditioning point.
The unconditional model would yield a Poisson number of grains, mean number 180,
of which we would expect 71:11 to cover at least one of the
points. So the excess number of grains from the conditional Boolean model covering
the conditioning points is 25:89.

Figure

4: Conditional Boolean model, parallelogram grid
points placed on a square lattice of side length 2=19 in a
window W of width units.
- Intensity of underlying germ process
- Radius of the typical grain
Simulation yields grains, of which 74 cover at least one conditioning point.
The unconditional model would yield a Poisson number of grains, mean number 120,
of which we would expect 47:41 to cover at least one of the
points. So the excess number of grains from the conditional Boolean model covering
the conditioning points is 26:59.

Figure

5: Conditional Boolean model, square grid
20.20.611.41.8Perfect sample from a conditional Boolean model with
points placed on a square lattice of side length 2=19 in a
window W of width units.
- Intensity of underlying germ process
- Radius of the typical grain
Simulation yields grains, of which 98 cover at least one conditioning point.
The unconditional model would yield a Poisson number of grains, mean number 180,
of which we would expect 71:11 to cover at least one of the
points. So the excess number of grains from the conditional Boolean model covering
the conditioning points is 26:89.

Figure

Conditional Boolean model, square grid
5 Checking for correctness
Implementations of simulation methods should always be checked out carefully, to
protect against mistakes in coding as well as to confirm the theory. In this section we
summarize several simple examples which provided checks for the correctness of the
method. The programs are implemented in C. The simulation time interval is initialized
at [ 2; 0] in each case, though of course the algorithm can extend this further back in
time if required. Window W is a unit square. G is set to be a disk of radius r,  is
the intensity of the germ process. To check correctness, in each case we collect 10000
samples and then carry out a  2 test using appropriate grouping. In the following we
only report summarized results. Full details of implementation, results, and derivations
of the exact equilibrium distribution functions for the examples can be found in the
associated research report [3], available as Preprint 350 at
http://www.warwick.ac.uk/statsdept/Staff/WSK/ppt.html.
1. Single conditioning points.
5. The resulting  2
5 value was 4.804.
2. Two separated conditioning points.
5. The
resulting  2
8 value was 9.368.
3. Two close conditioning points.
5. The
resulting  2
value was 12.344.
4. Three close conditioning points in 3-fold symmetry.
g. Place the three points in a configuration of 3-fold symmetry.
8. The resulting  2
26 value was 19.993.
None of the above  2 values are significant.
6 Investigating algorithm performance
In this section we carry out some experiments to explore the way in which the run-time
of the perfect simulation algorithm depends on various parameters.
It follows from the analysis in Section 4.2 that each unconditioned process YA is
a spatial immigration-death process of constant immigration rate  and unit death rate
on EA , where A  S and EA is given by Equation (12).
Recall that k is the number of conditioning points and r is the radius of a typical
grain.
In this section we will investigate the dependence of T c on , r, k. Note that
throughout this section, we use the actual coalescence time T c , not the time produced
by doubling during the performance of Algorithm CFTP-CCP.
Experiment 1: dependence of coalescence time on grain radius
Holding fixed the intensity of the germ process and the number of conditioning points,
we investigate the relation between coalescence time T c and the radius r of a typical
grain.
We place a set conditioning points on the
vertices of a square of side-length 0:2 in a unit-width window W , and obtain independent
samples of size 200 each for 99 equally spaced values of r ranging from 0:01 to
0:5.

Figure

7 plots log(sample mean) of coalescence time versus r (the dotted line),
and includes two theoretically derived asymptotics (the curved solid line for small ra-
dius, and the horizontal solid line for large radius). A discontinuity in the sample curve
signals the effect of the possibility of single grains covering several conditional points.
The theoretical asymptotics are derived as follows.
First consider the case of grains so small that any one grain can cover at most one
conditioning point (r < 0:1). Then the only relevant YA are the k processes for which
fig is a singleton, and they are independent of each other. Coalescence occurs
exactly when each of the Y fig has experienced a death at a time when there is no
perpetuation in V max
fig . For an individual Y fig the first such event occurs approximately
at the first instance of a Poisson process of rate r 2 which is independently thinned
with deletion probability
Using this approximation, the thinned process is also a Poisson process of rate r 2 =2.
Therefore
At the alternative extreme of large r (r >> 0:1), the conditioning points are likely
to be covered by just one disk of the unconditional Boolean model: coverage probability
for just one conditioning point on its own at time zero is
and the coverage events are positively correlated. Accordingly the coalescence time T c
will be approximated by the maximum of the first perpetuation times for each of the N
YA processes, so
not depending on grain radius r.
We deduce the following theoretical asymptotics for the mean coalescence time:
for small r, and
for large r.
Radius
log(mean)
of
the
Coalescence
Time
-226

Figure

7: Graph of the log(sample mean) of T c versus grain radius r for fixed conditioning
points value of 20. The reason for the big jump at 0:1 is as
follows. When 0:1, it is possible for one grain to cover more than one conditioning
point. The jump shows the cost we pay for the consequent interaction between e
Y
processes.
Experiment 2: dependence of coalescence time on intensity
Holding fixed the radius of a typical grain and the number of conditioning points, we
now investigate the relation between coalescence time T c and the intensity  of the
germ process.
We place a set
on the vertices of a square of side-length 0:2 in a unit-width window W , and obtain
independent samples of size 200 each for 75 equally spaced values of  ranging from 2
to 150.

Figure

8 plots the log(sample mean) of coalescence time versus  (dotted line)
and again includes two theoretically derived asymptotics (solid lines).
The theoretical asymptotics are derived essentially as for Experiment 1, as follows.
First consider the case where  is very large. Then all the conditioning points
are likely to be covered by a single disk of the unconditional Boolean model so, with
the same arguments as those in Experiment 1, the coalescence time T c can also be
approximated by the maximum of the N first perpetuation times for the YA processes.
Hence the mean is approximated by asymptotic (16), which is independent of the value
of .
On the other hand, for small values of , the coalescence time T c will depend on the
the value of r. In this experiment, we take r = 0:08, so it follows from the arguments
in Experiment 1 that when  is not very large, the distribution of T c is approximately
given by asymptotic (13). Hence the approximated theoretical mean is given by (15),
but viewed as a function of .
Intensity
log(mean)
of
the
Coalescence
Time
Figure

8: Graph of the log(sample mean) of T c versus intensity  for fixed conditioning
points radius of the grain
Experiment 3: dependence of coalescence time on the number of conditioning
points
Holding fixed the intensity of the germ process and the radius r of a typical grain, we
investigate the relation between coalescence time T c and the number of conditioning
points.
On a unit-width window W , when k is small, it is likely that a single grain can
cover at most one conditioning point as long as r and  are not very large. Therefore,
for small k, the distribution of T c can also be approximated by equation (13), hence the
mean is given by asymptotic (15) but now viewed as a function of k.
On the other hand, when k is very large, any disk of the unconditional Boolean
model are likely to cover at least one conditioning point. So the distribution of T c in
this case can be approximated by equation (14), where again, N is the number of YA
processes with A  S, A 6= 0. Hence the mean is given by asymptotic (16), but now
viewed as a function of k.
We illustrate the relation between T c and k when r is small. The experiment is
designed as follows.
We take conditioning points on
the vertices of rectangles in a unit-width window W , each rectangular has width 1=11
and length 1=6. We obtain independent samples of size 200 for each k. By above
settings, any grain can cover at most one conditioning point. Therefore, in this case,
we would expect that the sample mean versus k should agree with the approximated
theoretical mean curve given by asymptotic (15), as is confirmed by the plot of sample
mean (dotted lines) and the two asymptotics of the theoretical mean versus number of
conditioning points k in Figure 9.
Mean
of
the
Coalescence
Time
Figure

9: Graph of the sample mean of T c versus k for fixed value of
of the grain
Experiment 4: dependence of actual duration of algorithm on intensity
Finally we investigate the relation between actual duration ("wall-clock time") of the
algorithm and , r and k. In each experiment, we recorded 200 independent actual
duration times for collecting one sample for each value of r,  or k and plotted the
average duration time and 95% confidence limits about the mean versus r,  or k
respectively in the following figures.
We first consider the duration time as in Experiment 1. Theoretical arguments
suggest that mean duration time should be a convex function of grain radius r, tending
to infinity as r ! 0 and to a positive constant as r ! 1. Simulations support this
conclusion, but also show that the significant features of dependence occur over only a
relatively small range of r.
In

Figure

we plot the log(sample mean) over 200 realizations and the 95% confidence
limits about the mean over the range r 2 [0:01; 0:5], from which we see that
there are two peaks corresponding to two values of 0:145. This
corresponds to the following observations: as r increases through
possible for one grain to cover at most 2 conditioning points, while as r increases
through becomes possible for one grain to cover at most 4 conditioning
points.
Now consider the dependence of mean duration time on intensity , as in Experiment
2. We gain a similar theoretical picture, except that for large intensity the mean
duration time should increase from a minimum (since the algorithm needs more time
to deal with more and more incidents of YA and ZA processes as  increases). Figure
11 plots the log(sample mean) over 200 realizations and 95% confidence limits about
the mean, and confirms this argument from theory.
The sharp increase in duration time as intensity is reduced is mirrored by our ex-
Radius
log
Actual
Duation
(in
seconds)
of
The
Algorithm
-22

Figure

10: Graph of the log(mean actual duration) (in seconds) and 95% confidence
limits over 200 realizations of the algorithm versus r for fixed value of
4. The two peaks correspond to the two values of r when one grain can cover at
most two and four conditioning points respectively.
perience of larger-scale simulations as displayed in Figures 1-6, particularly by the
failure of the algorithm to produce coalescence for low intensity in such cases.
Things are much simpler in the case of dependence on number of conditioning
points k as in Experiment 3. As in Experiment 3, fix 0:04. In this
case, all the e
Y processes are independent with each other. Therefore, if we need h
amount of time for coalescence for each conditioning point, then the total amount of
time for coalescence for k conditioning points will be approximately kh. Figure 12
plots the log(sample mean) over 200 realizations and 95% confidence limits about the
mean. This figure shows that the actual time increase as number of conditioning points
increases, in agreement with the above.
7 Conclusion
In this paper we have presented a perfect simulation method for correlated Poisson
random variables conditioned to be positive. Working in a general context (as opposed
to the geometric context of [11]) has made the algorithm simpler to implement and to
present. Moreover, although the method in this paper does not fall in the category of
dominated CFTP defined in Kendall &Moller [10] nevertheless it is a modest variation
of it. In this paper, "domination" does not occur with respect to the partial order defined
in the state space, but the lower- and upper sandwich processes are "dominated" by a
process in extended state space formed from YA and ZA for all A  S.
While we cannot give theoretical bounds on the range of parameters involved in
the particular kind of perfect simulation, it may be seen from the results of Sections
Intensity
log
Actual
Duation
(in
seconds)
of
The
Algorithm

Figure

11: Graph of the log(mean actual duration) (in seconds) and 95% confidence
limits over 200 realizations of the algorithm versus  for fixed value of
Number of k
log
Actual
Duation
(in
seconds)
of
The
Algorithm
-22
Figure

12: Graph of the log(mean actual duration) (in seconds) and 95% confidence
limits over 200 realizations of the algorithm versus k for fixed value of
4.2 and 6 that low-intensity conditional Boolean models present challenges which this
particular perfect simulation algorithm has not overcome. Of course there is a contrast
here with conventional McMC methods: perfect simulation clearly indicates when it
is failing to give a correct answer, while conventional McMC may fail without overt
warning. This should not however be overstated. Firstly, there is a price to pay for
CFTP [2], so that CFTP methods may fail when judiciously chosen McMC algorithms
still give good answers. Secondly, CFTP exhibits bias if the user is liable to terminate
runs early: this is overcome by Fill's rejection method [4] and this has been generalized
to the point process context in [14] though, as Wilson has pointed out to one of us, a
problem remains because of time-variability due to simulation of point processes as
part of the point process algorithm. Despite these caveats, perfect simulation remains
an attractive member of the statistician's tool-box, and the above algorithm delivers on
its promise for medium and high intensity conditional Boolean models.
Further development of the particular problem of perfect simulation of conditional
Boolean models leads to the consideration of conditioning constraints based on connectivity
[12]. Perfect simulation in such situations is likely to prove very difficult,
since the conditioning appears to involve arbitrarily large numbers of extra grains.
As with work on perfect simulation of Gibbs point processes [10, 9] we have taken
care to test our simulation program carefully. This is a particular issue for perfect simulation
implementations, which are inevitably more complicated than straightforward
McMC algorithms, and raises interesting issues of how best to test correctness in such
cases. In the above we have concentrated on straightforward  2 tests and compatibility
with simple theoretical asymptotics: [9] gives details of a pseudo-likelihood approach
which deals directly with the coupled random processes lying at the heart of the sim-
ulation. We hope to return to this issue in later work, as it has ramifications beyond
perfect simulation in general issues of programming complicated simulations.

Acknowledgements

This research was supported by EPSRC grant GR/L56831 and by the EU TMR network
ERB-FMRX-CT96-0095 on "Computational and Statistical methods for the analysis of
spatial data".



--R

Asymptotic Techniques for use in Statistics.
Efficient Markovian couplings: examples and counterexamples.
Perfect implementation of simulation for conditioned Boolean model via correlated Poisson random variables.
An interruptible algorithm for exact sampling via Markov Chains.

On some weighted Boolean models.
Perfect simulation for spatial point processes.
Perfect simulation for the area-interaction point process
Perfect implementation of a Metropolis-Hastings simulation of Markov point processes
Perfect simulation of point processes using Metropolis-Hastings algorithms
Perfect simulation in stochastic geometry.

Exact sampling with coupled Markov chains and applications to statistical mechanics.
Perfect simulation of some point processes for the impatient user.
--TR

--CTR
James D. Brown , Gerard B. M. Heuvelink, The Data Uncertainty Engine (DUE): A software tool for assessing and simulating uncertain environmental variables, Computers & Geosciences, v.33 n.2, p.172-190, February, 2007
Markov chain Monte Carlo;correlated Poisson random variables conditioned to be positive;conditional Boolean model;perfect simulation;extended state-space CFTP;dominated CFTP