Abstract
In this paper we discuss different Monte Carlo (MC) approaches to generate unit-rate Poisson processes and provide an analysis of Poisson bridge constructions, which form the discrete analogue of the well-known Brownian bridge construction for a Wiener process. One of the main advantages of these Poisson bridge constructions is that they, like the Brownian bridge, can be effectively combined with variance reduction techniques. In particular, we show here, in practice and proof, how we can achieve orders of magnitude efficiency improvement over standard MC approaches when generating unit-rate Poisson processes via a synthesis of antithetic sampling and Poisson bridge constructions. At the same time we provide practical guidance as to how to implement and tune Poisson bridge methods to achieve, in a mean sense, (near) optimal performance.
A Proofs
A.1 Proof of Proposition 5.1
Without loss of generality we take t∈[tr,tr+1] for some 1≤r<R . Next we use the fact that if we condition the Poisson processes on the start of the interval, i.e. on Y(1)(tr),…,Y(M)(tr) , and the end of the interval, i.e. on Y(1)(tr+1),…,Y(M)(tr+1) , the increments of the Poisson processes are binomially distributed for t∈[tr,tr+1] , i.e.
where
Then note that in general we can decompose the MSE for
and in particular this holds for
where
In addition, since the Poisson skeleton is unbiased by assumption, we know that
and thus
Using the fact that the Poisson process conditional on its endpoints is given by the value of the Poisson skeleton at
This, in turn, implies that the covariance of the conditional expectation can be factored as
Finally, summing over all
B Poisson bridge implementation details
In this appendix we provide explicit descriptions of the algorithms used in the midpoint-based and median-based Poisson bridge constructions. This is largely based on [4, Chapter 15] with minor corrections and additions to allow for more memory efficient implementations. We use matlab-index notation for vectors throughout this appendix.
B.1 Efficient sampling of uniform order statistics
We first describe in Algorithm 1 an efficient and stable implementation of the exponential spacings method to sequentially generate the order statistics of uniform random variables on a finite interval based on [4, Section 15.6],[4] which is a numerically stable version of the exponential spacings method [2, Section V.3.1]. This subroutine is used in both the median-based and midpoint-based Poisson bridge constructions.
Algorithm 1 (
UOS
(
a
,
b
,
N
)
. (This generates N ordered uniform random variates on
(
a
,
b
)
.)).
![](/document/doi/10.1515/mcma-2021-2090/asset/graphic/j_mcma-2021-2090_fig_0033.png)
B.2 Midpoint-based Poisson bridge
The full midpoint-based Poisson bridge construction is described in Algorithm 2, which uses two subroutines; Algorithm 4 to construct a Poisson skeleton, which is based on [4, Section 15.3], and Algorithm 5 to fill in such a Poisson skeleton efficiently, using Algorithm 1. The full Poisson bridge construction uses the fact that multiple steps of the Poisson bridge can be sampled independently by the independent increment property of a Poisson process. Although in Algorithm 2 we sample a Poisson process on
Note that to use variance reduction techniques in combination with the Poisson bridge method we simply propose to sample the uniform random variables,
where
Algorithm 2 (
PB
midpoint
(
τ
f
,
τ
s
,
L
,
𝐔
)
. (This generates a single Poisson process
Y
(
t
)
when
t
∈
[
0
,
⌈
τ
f
τ
s
⌉
τ
s
)
via the midpoint-based Poisson bridge construction with L levels and step size
τ
s
.)).
![](/document/doi/10.1515/mcma-2021-2090/asset/graphic/j_mcma-2021-2090_fig_0034.png)
B.3 Median-based Poisson bridge
The full median-based Poisson bridge construction is described in Algorithm 3, which uses two subroutines; Algorithm 6 to construct a Poisson skeleton and Algorithm 7 to fill in such a Poisson skeleton. Note that to use variance reduction techniques in combination with the Poisson bridge method we simply make sure that the uniform random variables,
To implement a single-step of the median-based Poisson bridge we can almost directly follow [4, Section 15.2]. Note, however, that [4, Section 15.2] describes a (hybrid) median-based Poisson bridge of fixed internal time length, i.e. for a given internal time τ one first samples
In order to use variance reduction methods the gamma and beta random variates need to be generated via inverse transform sampling. For the beta random variates this means that we have to invert the regularised incomplete beta function, i.e. solve
Algorithm 3 (
PB
median
(
N
f
,
N
s
,
L
,
𝐔
)
. (This generates a single Poisson process
Y
(
t
)
when
0
≤
Y
(
t
)
≤
⌈
N
f
N
s
⌉
N
s
via the median-based Poisson bridge construction with L levels and step size
N
s
.)).
![](/document/doi/10.1515/mcma-2021-2090/asset/graphic/j_mcma-2021-2090_fig_0035.png)
B.4 Subroutines
Algorithm 4 (
stepPB
midpoint
(
τ
s
,
L
,
𝐔
)
. (This generates a single Poisson skeleton for the Poisson process
Y
(
t
)
when
t
∈
[
0
,
τ
s
)
via the midpoint-based Poisson bridge construction with L levels.)).
![](/document/doi/10.1515/mcma-2021-2090/asset/graphic/j_mcma-2021-2090_fig_0036.png)
Algorithm 5 (
fillPB
midpoint
(
𝐌
,
𝐭
)
. (This fills in a single Poisson skeleton generated by the midpoint-based Poisson bridge construction.)).
![](/document/doi/10.1515/mcma-2021-2090/asset/graphic/j_mcma-2021-2090_fig_0037.png)
Algorithm 6 (
stepPB
median
(
N
s
,
L
,
𝐔
)
. (This generates a single Poisson skeleton for the Poisson process
Y
(
t
)
when
0
≤
Y
(
t
)
≤
N
s
via the median-based Poisson bridge construction with L levels.)).
![](/document/doi/10.1515/mcma-2021-2090/asset/graphic/j_mcma-2021-2090_fig_0038.png)
Algorithm 7 (
fillPB
median
(
𝐣
,
𝐭
)
. (This fills in a single Poisson skeleton generated by the median-based Poisson bridge construction.)).
![](/document/doi/10.1515/mcma-2021-2090/asset/graphic/j_mcma-2021-2090_fig_0039.png)
References
[1] C. H. L. Beentjes, Variance reduction techniques for chemical reaction network simulation, Ph.D. thesis, University of Oxford, 2020. Search in Google Scholar
[2] L. Devroye, Non-Uniform Random Variate Generation, Springer, New York, 1986. 10.1007/978-1-4613-8643-8Search in Google Scholar
[3] B. L. Fox, Generating Poisson processes by quasi-Monte Carlo, Lecture notes (1996). Search in Google Scholar
[4] B. L. Fox, Strategies for Quasi-Monte Carlo, Internat. Ser. Oper. Res. Management Sci. 22, Springer, New York, 1999. 10.1007/978-1-4615-5221-5Search in Google Scholar
[5] A. Gil, J. Segura and N. M. Temme, Efficient and accurate algorithms for the computation and inversion of the incomplete gamma function ratios, SIAM J. Sci. Comput. 34 (2012), no. 6, A2965–A2981. 10.1137/120872553Search in Google Scholar
[6] M. B. Giles, Algorithm 955: Approximation of the inverse Poisson cumulative distribution function, ACM Trans. Math. Software 42 (2016), no. 1, 1–22. 10.1145/2699466Search in Google Scholar
[7] P. Glasserman, Monte Carlo Methods in Financial Engineering, Springer, New York, 2003. 10.1007/978-0-387-21617-1Search in Google Scholar
[8] P. L’Ecuyer and R. Simard, Inverting the symmetrical beta distribution, ACM Trans. Math. Software 32 (2006), no. 4, 509–520. 10.1145/1186785.1186786Search in Google Scholar
[9] C. Lemieux, Monte Carlo and Quasi-Monte Carlo Sampling, Springer Ser. Statist., Springer, New York, 2009. 10.1007/978-0-387-78165-5_5Search in Google Scholar
[10] P. A. Maginnis, M. West and G. E. Dullerud, Exact variance-reduced simulation of lattice continuous-time Markov chains with applications in reaction networks, Bull. Math. Biol. 81 (2019), no. 8, 3159–3184. 10.1007/s11538-019-00576-2Search in Google Scholar
[11] B. Moskowitz and R. E. Caflisch, Smoothness and dimension reduction in quasi-Monte Carlo methods, Math. Comput. Model. 23 (1996), no. 8–9, 37–54. 10.1016/0895-7177(96)00038-6Search in Google Scholar
[12] A. B. Owen, Monte Carlo Theory, Methods and Examples, Book in progess, http://statweb.stanford.edu/~owen/mc/. Search in Google Scholar
© 2021 Walter de Gruyter GmbH, Berlin/Boston