A Review on the Exact Monte Carlo Simulation

Abstract

Perfect Monte Carlo sampling refers to sampling random realizations exactly from the target distributions (without any statistical error). Although many different methods have been developed and various applications have been implemented in the area of perfect Monte Carlo sampling, it is mostly referred by researchers to coupling from the past (CFTP) which can correct the statistical errors for the Monte Carlo samples generated by Markov chain Monte Carlo (MCMC) algorithms. This paper provides a brief review on the recent developments and applications in CFTP and other perfect Monte Carlo sampling methods.

1. Introduction

In the past 30 years, substantial progress has been made in popularizing Bayesian methods for the statistical analysis of complex data sets. An important driving force has been the availability of different types of Bayesian computational methods, such as Markov chain Monte Carlo (MCMC), sequential Monte Carlo (SMC), approximate Bayesian computation (ABC) and so on. For many practical examples, these computational methods can provide rapid and reliable approximations to the posterior distributions for unknown parameters.

The basic idea that lies behind these methods is to obtain Monte Carlo samples from the distribution of interest, in particular the posterior distribution. In Bayesian analysis of complex statistical models, the calculation of posterior normalizing constants and the evaluation of posterior estimates are typically infeasible either analytically or by numerical quadrature. Monte Carlo simulation provides an alternative. One of the most popular Bayesian computational methods is MCMC, which is based on the idea of constructing a Markov chain with the desired distribution as its stationary distribution.

By running a Markov chain, MCMC methods generate statistically dependent and approximate realizations from the limiting (target) distribution. A potential weakness of these methods is that the simulated trajectory of a Markov chain will depend on its initial state. A common practical recommendation is to ignore the early stages, the so-called burn-in phase, before collecting realizations of the state of the chain. How to choose the length of the burn-in phase is an active research area. Many methods have been proposed for convergence diagnostics; [1] gave a comparative review. Rigorous application of diagnostic methods requires either substantial empirical analysis of the Markov chain or complex mathematical analysis. In practice, judgments about convergence are often made by visual inspection of the realized chain or the application of simple rules of thumb.

Concerns about the quality of the sampled realizations of the simulated Markov chains have motivated the search for Monte Carlo methods that can be guaranteed to provide samples from the target distribution. This is usually referred to as perfect sampling or exact sampling. In some cases, for example, the multivariate normal, perfect samples are readily available. For more complicated distributions, perfect sampling can be achieved, in principle, by the rejection method. This involves sampling from a density that bounds a suitable multiple of the target density, followed by acceptance or rejection of the sampled value. The difficulty here is in finding a bounding density that is amenable to rapid sampling while at the same time providing sample values that are accepted with high probability. In general this is a challenging task, although efficient rejection sampling methods have been developed for the special class of log-concave densities; see, for example, [2, 3].

A surprising breakthrough in the search for perfect sampling methods was made by [4]. The method is known as coupling from the past (CFTP). This is a sophisticated MCMC-based algorithm that produces realizations exactly from the target distribution. CFTP transfers the difficulty of running the Markov chain for extensive periods (to ensure convergence) to the difficulty of establishing whether a potentially large number of coupled Markov chains have coalesced. An efficient CFTP algorithm depends on finding an appropriate Markov chain construction that will ensure fast coalescence. There have been a few further novel theoretical developments following the breakthrough of CFTP, including [2, 5, 6]. Since then, perfect sampling methods have attracted great attention in various Bayesian computational problems and applied probability areas.

Apart from coupling from the past, many other perfect sampling methods were proposed for specific problems, for example, perfect sampling for random spanning trees [7, 8] and path-space rejection sampler for diffusion processes [9, 10, 11]. Very recently, a type of divide-and-conquer method has been developed in [12, 13]. These methods use the technique for the exact simulation of diffusions and samples from simple sub-densities to obtain perfect samples from the target distribution.

Perfect samples are useful in Bayesian applications either as a complete replacement for MCMC-generated values or as a source of initial values that will guarantee that a conventional MCMC algorithm runs in equilibrium. Perfect samples can also be used as a means of quality control in judging a proposed MCMC implementation when there are questions about the speed of convergence of the MCMC algorithm or whether the chain is capable of exploring all parts of the sample space. Of course, when perfect samples can be obtained speedily, they will be preferred to MCMC values, since they eliminate such doubts. In addition, sophisticated perfect sampling methodology often motivates efficient approximate algorithms and computational techniques. For example, [14] uses regenerated Markov chains to obtain simple standard error estimates for importance sampling under MCMC context. The condition required there will allow us to carry out perfect sampling via multigamma coupling approach [15].

This paper will present a brief review for perfect Monte Carlo sampling and explain the advantages and drawbacks of different types of methods. Section 2 will present rejection sampling techniques, and then CFTP will be covered in Section 3. Recent divide-and-conquer methods will be reviewed in Section 4. The paper ends with a discussion in Section 5.

2. Rejection sampling techniques

Rejection sampling, also known as ‘acceptance-rejection sampling’, generates realizations from a target probability density function $f\left(x\right)$ by using a hat function $\mathit{Mg}\left(x\right)$, where $f\left(x\right)\le \mathit{Mg}\left(x\right)$ and $g\left(x\right)$ are a probability density function from which samples can be readily simulated. The basic rejection sampling algorithm is as follows:

Algorithm 2.1 (Basic rejection sampling)

Sample $x$ from $g\left(x\right)$ and $U$ from Unif$\left[0,1\right]$.	01
If $U\le \frac{f\left(x\right)}{\mathit{Mg}\left(x\right)}$,	02
accept $x$ as a realisation of $f\left(x\right)$ and stop;	03
else	04
reject the value of $x$ and go back to step 01.	05

Many other perfect sampling methods are actually equivalent to rejection sampling. For example, ratio-of-uniform (RoU) method [16] may have to be implemented via a rejection sampling approach.

The efficiency of rejection sampling depends on the acceptance probability, which is $1/M$. To perform rejection sampling efficiently, it is very important to find hat functions which provides higher acceptance probabilities. In other words, we shall choose $M$ as small as possible [17]; however for many complicated problems, there is no easy way to find $M$ small enough to guarantee high acceptance probability. A number of sophisticated rejection sampling methods have been suggested for dealing with complex Bayesian posterior densities, which we discuss below.

2.1 Log-concave densities

A function $h\left(x\right)$ is called log-concave if

$$logh\left(\mathit{\lambda x}+\left(1-\lambda \right)y\right)\ge \lambda logh\left(x\right)+\left(1-\lambda \right)logh\left(y\right),$$

for all $x,y$ and $\lambda \in \left[0,1\right]$. For the special class of log-concave densities, Gilks and wild [3] developed the adaptive rejection sampling (ARS) method. The method constructs an envelope function for the logarithm of the target density, $f\left(x\right)$, by using tangents to $logf\left(x\right)$ at an increasing number, $n$, of points. The envelope function $u_n\left(x\right)$ is the piecewise linear upper hull formed from the tangents. Note that, the envelope can be easily constructed due to the concavity of $logf\left(x\right)$. The method also constructs a squeeze function $l_n\left(x\right)$ which is formed from the chords of the tangent points. The sampling steps of the ARS algorithm are as follows.

Algorithm 2.2 (Adaptive rejection sampling).

Outputs a stream of perfect samples from $f\left(x\right)$.
Initialise $u_n\left(x\right)$ and $l_n\left(x\right)$ by using tangents at several points.	01
Sample $x^{\ast }$ from density $\propto exp\left(u_n\left(x\right)\right)$ and $U\sim \text{Unif}\left(0,1\right)$	02
If $U\le exp\left(l_n\left(x^{\ast }\right)-u_n\left(x^{\ast }\right)\right)$, Output $x^{\ast }$;	03
else if $U\le f\left(x^{\ast }\right)/exp\left(u_n\left(x^{\ast }\right)\right)$,	04
Output $x^{\ast }$; Update $\left(u_n,l_n\right)$ to $\left(u_{n+1},l_{n+1}\right)$ using $x^{\ast }$;	05
Goto 02	06

The ARS algorithm is adaptive and the sampling density is modified whenever $f\left(x^{\ast }\right)$ is evaluated. In this way, the method becomes more efficient as the sampling continues. Leydold [18] extends ARS to log-concave multivariate densities. Leydold’s algorithm is based on decomposing the domain of the density into cones and then computing tangent hyperplanes for these cones. Generic computer code for sampling from a multivariate log-concave density is available on the author’s website; it is only necessary to code a subroutine for the target density. Leydold’s algorithm works well for simple low-dimensional densities. The drawback of ARS algorithm is that it only works for log-concave densities, which is a very small class of posteriors in practice. Also computationally it is very challenging to implement ARS algorithm for high-dimensional densities [19].

2.2 Fill’s rejection sampling algorithm

We consider a discrete Markov chain with transition probability $\mathbf{P}\left(x,y\right)$ and stationary distribution $\pi \left(x\right),x\in S$. Let $\tilde{\mathbf{P}}\left(x,y\right)=\pi \left(y\right)\mathbf{P}\left(y,x\right)/\pi \left(x\right)$ be the transition probability for the time-reversed chain. Suppose that there is a partial ordering on the states $S$, denoted by $x\preccurlyeq y$. Let $\widehat{0}$ and $\widehat{1}$ be unique extremal points of the partial ordering.

To demonstrate the algorithm given by [2], we will assume that there are update functions $\varphi$ and $\tilde{\varphi }$ both mapping $S\times \left[0,1\right]$ to $S$ such that

$$\mathbf{P}\left(x,y\right)=P\left(\varphi \left(x,U\right)=y\right),$$E1

$$\tilde{\mathbf{P}}\left(x,y\right)=P\left(\tilde{\varphi }\left(x,U\right)=y\right),$$E2

where $U\sim$ Unif$\left[0,1\right]$ and

$$x\preccurlyeq y\Rightarrow \tilde{\varphi }\left(x,u\right)\preccurlyeq \tilde{\varphi }\left(y,u\right)a.e.u\in \left[0,1\right].$$

The algorithm is as follows:

Algorithm 2.3 (Fill’s algorithm)

1. Choose an integer $t>0$. Fix $x_0=\widehat{0}$ and $y_0=\widehat{1}$.

2. Generate $x_k=\varphi \left(x_{k-1},U_k\right)$, $k=1,\cdots ,t$, where $\left\{U_k\right\}$ are i.i.d. Unif$\left[0,1\right]$.

3. Generate ${\tilde{U}}_k$ from the conditional distribution of $U$ given $\tilde{\varphi }\left(x_{t-k+1},U\right)=x_{t-k},k=1,\cdots ,t$.

4. Generate $y_k=\tilde{\varphi }\left(y_{k-1},{\tilde{U}}_k\right),k=1,\cdots ,t$.

5. If $y_t=x_0$ then accept $x_t$; else double $t$ and repeat from step 2.

In Algorithm 2.3 (step 2) $z≔x_t$ is sampled from ${\mathbf{P}}^t\left(\widehat{0},\cdot \right)$. If we find an upper bound $M$ for $\pi \left(z\right)/{\mathbf{P}}^t\left(\widehat{0},z\right)$, then we can use rejection sampling. Fill [2] finds a bounding constant given by $M=\pi \left(\widehat{0}\right)/{\tilde{\mathbf{P}}}^t\left(\widehat{1},\widehat{0}\right)$ and proves that steps from 3 to 5 of Algorithm 2.3 are to accept $z$ with probability $\frac{\pi \left(z\right)}{M{\mathbf{P}}^t\left(\widehat{0},z\right)}$. The output of this algorithm is indeed a perfect sample from $\pi$.

From Algorithm 2.3, we can see that rejection sampling can still be possible, even if we do not have a closed form of the hat function. The first limitation of Algorithm 2.3 is that it works only if the time-reversed chain is monotone, but [5] has extended the algorithm theoretically for general chains. The second limitation is that step 3 of Algorithm 2.3 is difficult to perform [2]. For these reasons, Fill’s algorithm is not practical for realistic problems.

3. Coupling from the past

Coupling from the past was introduced in the landmark paper of [4] which showed how to provide perfect samples from the limiting distribution of a Markov chain.

3.1 Basic CFTP algorithms

Let $\left\{X_t\right\}$ be an ergodic Markov chain with state space $\mathcal{X}=\left\{1,\cdots ,n\right\}$, where the probability of going from $i$ to $j$ is $p_{\mathit{ij}}$ and the stationary distribution is $\pi$. Suppose we design an updating function $\varphi \left(\cdot ,U\right)$, which satisfies $P\left[\varphi \left(i,U\right)=j\right]=p_{\mathit{ij}}$, where $\varphi$ is a deterministic function and $U$ is a random variable. To simulate the next state $Y$ of the Markov chain, currently in state $i$, we draw a random variable $U$ and let $Y=\varphi \left(i,U\right)$.

Let $f_t\left(i\right)=\varphi \left(i,U_t\right)$, and define the composition

$$F_{t_1}^{t_2}=f_{t_2-1}\circ f_{t_2-2}\circ \cdots \circ f_{t_1+1}\circ f_{t_1},$$E3

for $t_1<t_2$.

Proposition 3.1 [4] Suppose there exists a time $t=-T$, the backward coupling time, such that chains starting from any state in $\mathcal{X}=\left\{1,\cdots ,n\right\}$, at time $t=-T$, and with the same sequence $\left\{U_t,t=-T,\cdots ,-1\right\}$, arrive at the same state $X_0^{\ast }$. Then it must follow that $X_0^{\ast }$ comes from $\pi$.

This proposition is easy to prove. If we run an ergodic Markov chain from time $t=-\infty$ and with the sequence $\left\{U_t,t=-T,\cdots ,-1\right\}$ after $-T$, the Markov chain will arrive at $X_0^{\ast }$. Then $X_0^{\ast }$ comes exactly from $\pi$ since it is collected at time $0$ and the Markov chain started from $-\infty$. The algorithm is as follows:

Algorithm 3.1 (Basic CFTP)

$t=0$	01
repeat	02
$t=t-1$	03
generate $U_t$	04
$F_t^0=F_{t+1}^0\circ \varphi \left(\cdot ,U_t\right)$	05
until $F_t^0\left(\cdot \right)$ is a constant	06
return $F_t^0\left(\cdot \right)$	07

Propp and Wilson [4] also proved that this algorithm is certain to terminate. The idea of simulating from the past is important. Note that if we collect $F_0^T\left(\cdot \right)$ as the result, where $T$ is the smallest value that makes $F_0^T\left(\cdot \right)$ a constant, we will get a biased sample. This is because $T$ is a stopping time, which is not independent of the Markov chain.

3.2 Read-once CFTP

The basic CFTP algorithm needs to restart the Markov chains at some points in the past if they have not coalesced by time $0$. We must use the same random sequence ${\left\{U_t\right\}}_{-\infty }^{-1}$ when we restart the Markov chains. In Monte Carlo simulations, we usually use pseudorandom number generators, which are deterministic algorithms. So if we give the same random seed, we will get the same random sequence. This means that the same sequence $\left\{U_t\right\}$ can be regenerated in each coupling procedure.

If we can run the Markov chain forward starting at $0$ and collect a perfect sample in the future, we will not have to regenerate $\left\{U_t\right\}$. Wilson [20] developed a read-once CFTP method to implement the forward coupling idea. A simple example is provided by [21]. In fact, the multigamma coupler in [15] can be implemented via the more efficient read-once CFTP algorithm.

3.3 Improvement of CFTP algorithms

Propp and Wilson [4] showed that the computational cost of the algorithm can be reduced if there is a partial order for the state space $\mathcal{X}$ that is preserved by the update function $\varphi$, i.e. if $x\le y$ then $\varphi \left(x,U\right)\le \varphi \left(y,U\right)$. Their procedure is outlined in Algorithm 3.2, whereas before $\widehat{0}$ and $\widehat{1}$ are the unique extremals. Note that their algorithm needs a monotone update function $\varphi$ for the Markov chain, while Algorithm 2.3 requires a monotone update function $\tilde{\varphi }$ for the time-reversed chain.

Algorithm 3.2 (Monotone CFTP)

$T=1$	01
Repeat	02
upper $=\widehat{1}$	03
lower $=\widehat{0}$	04
for $t=-T$ to $t=-1$	05
upper $={\varphi }_t\left(\text{upper},U_t\right)$	06
lower $={\varphi }_t\left(\text{lower},U_t\right)$	07
$T=2T$	08
until upper $=$ lower	09
return upper	10

Algorithm 3.2 is much simpler than Algorithm 3.1, since only two chains have to be run at the same time, but the requirement of monotonicity is very restrictive. Markov chains with transitions given by independent Metropolis-Hastings and perfect slice sampling have been shown to have this property, by [22, 23], respectively. However [23, 24] have also noticed that such independent M-H CFTP is equivalent to simple rejection sampler.

In general it is hard to code perfect slice samplers correctly. For example, Hörmann and Leydold [25] have pointed out that the perfect slice samplers in [26, 27] are incorrect. The challenge of monotone CFTP is usually to construct the detailed updating function with a guarantee of preserving the partial order.

Finding a partial order preserved by the Markov chains is a non-trivial task in many cases. An alternative improvement is to use CFTP with bounding chains, such as that in [28, 29]. If the bounding chains, which bound all the Markov chains, coalesce, then all Markov chains coalesce. Thus if only a few bounding chains are required, the efficiency of the CFTP algorithm can be improved significantly. Sometimes, it may be impossible to define an explicit bounding chain (the maximum of the state space may be infinity, and the upper bound chain cannot start from infinity), but it is possible to use a dominated process to bound all Markov chains [30].

3.4 Applications and challenges

Although CFTP is extremely challenging to be implemented for many practical problems, it did find a few applications in certain discrete state space problems, for example, the Ising model [4]. Also [31] applied CFTP to ancestral selection graph to simulate samples from population genetic models. Refs. [32, 33] applied CFTP to a class of fork-join type queuing system problems. Connor and Kendal [34] applied CFTP for the perfect simulation of M/G/c queues. CFTP also finds its application in signal processing [35].

CFTP for continuous state space Markov chains is very challenging, since a random map from an interval to a finite number of points is required. In recent years, many methods have been developed for unbounded continuous state space Markov chains, such as perfect slice sampler in [23], multigamma coupler and the bounded M-H coupler in [15, 24]. Wilson [36] developed a layered multi-shift coupling, which shifts states in an interval to a finite number of points. However, none of these methods can solve any practical problems.

4. Recent advances in perfect sampling

Recently, a new type of perfect Monte Carlo sampling method based on the decomposition of the target density $f$, as $f\left(\cdot \right)=g_1\left(\cdot \right)g_2\left(\cdot \right)$, was proposed in [12], where $g_1$ and $g_2$ are also (proportional to) proper density functions. Note that here $g_1$ and $g_2$ are continuous density functions which are easy to simulate from. Suppose that $q$-dimensional vector values ${\mathit{x}}_1$ and ${\mathit{x}}_2$ are independently drawn from $g_1$ and $g_2$, respectively. If the two independent samples are equal, i.e. ${\mathit{x}}_1={\mathit{x}}_2=\mathit{y}$ then we have $\mathit{y}$ must be from $f\left(\cdot \right)\propto g_1\left(\cdot \right)g_2\left(\cdot \right)$. Note that such a naive approach may be practical for discrete random variables with low-dimensional state space, but for continuous random variables, it is impossible since $\mathrm{P}\left({\mathit{x}}_1={\mathit{x}}_2\right)=0$. Dai [12] proposed a novel approach to deal with this, which is explained in the following subsection.

4.1 Perfect distributed Monte Carlo without using hat functions

First we introduce the following notations related to the logarithm of the sub-densities:

$$\mathit{\alpha }\left(\mathit{x}\right)={\left({\alpha }^{\left(1\right)},\cdots ,{\alpha }^{\left(q\right)}\right)}^{\mathit{tr}}\left(\mathit{x}\right)=\nabla log{g}_1\left(\mathit{x}\right)$$E4

where $\nabla$ is the partial derivative operator for each component of $\mathit{x}$. Then we consider a $q$-dimensional diffusion process ${\mathit{X}}_t\left(\overrightarrow{\mathit{\omega }}\right),t\in \left[0,T\right]$ ($T<\infty$), defined on the $q$-dimensional continuous function space $\mathbf{\Omega }$, given by:

$$d{\mathit{X}}_t=\mathit{\alpha }\left({\mathit{X}}_t\right)\mathit{dt}+d{\mathit{B}}_t,$$E5

where ${\mathit{B}}_t\left(\overrightarrow{\mathit{\omega }}\right)={\mathit{\omega }}_t$ is a Brownian motion and $\overrightarrow{\mathit{\omega }}=\left\{{\mathit{\omega }}_t,t\in \left[0,T\right]\right\}$ is a typical element of $\mathbf{\Omega }$. Let $\mathbb{W}$ be the probability measure for a Brownian motion with the initial probability distribution ${\mathit{B}}_0={\mathit{w}}_0\sim f_1\left(\cdot \right)=g_1^2\left(\cdot \right)$.

Clearly ${\mathit{X}}_t$ has the invariant distribution $f_1\left(\mathit{x}\right)$ (using the Langevin diffusion results [37]). Let $\mathbb{Q}$ be the probability law induced by ${\mathit{X}}_t,t\in \left[0,T\right]$, with ${\mathit{X}}_0={\mathit{\omega }}_0\sim f_1\left(\cdot \right)$, i.e. under $\mathbb{Q}$ we have ${\mathit{X}}_t\sim f_1\left(\mathit{x}\right)$ for any $t\in \left[0,T\right]$.

The idea in [12] is to use a biased diffusion process $\overline{\mathbf{X}}=\left\{{\overline{\mathit{X}}}_t;0\le t\le T\right\}$ to simulate from the target function $f$. It is defined as follows.

Definition 4.1 (Biased Langevin diffusions) The joint density for the pair $\left({\overline{\mathit{X}}}_0,{\overline{\mathit{X}}}_T\right)$ (the starting and ending points of the biased diffusion process), evaluated at point $\left(\mathit{x},\mathit{y}\right)$, is $f_1\left(\mathit{x}\right){\mathbf{t}}^{\ast }\left(\mathit{y}|\mathit{x}\right)f_2\left(\mathit{y}\right)$, where ${\mathbf{t}}^{\ast }\left(\mathit{y}|\mathit{x}\right)$ is the transition density for the diffusion process $X$ defined in Eq. (6) from ${\mathit{X}}_0=\mathit{x}$ to ${\mathit{X}}_T=\mathit{y}$ and $f_2\left(\mathit{y}\right)=g_2\left(\mathit{y}\right)/g_1\left(\mathit{y}\right)$.

Given $\left({\overline{\mathit{X}}}_0,{\overline{\mathit{X}}}_T\right)$ the process $\left\{{\overline{\mathit{X}}}_t,0<t<T\right\}$ is given by the diffusion bridge driven by Eq. (6).

The marginal distribution for ${\overline{\mathit{X}}}_T$ is $f\left(\mathit{y}\right)$. Therefore, to draw a sample from the target distribution $f\left(\mathit{x}\right)$, we need to simulate a process ${\overline{\mathit{X}}}_t,t\in \left[0,T\right]$ from $\overline{\mathbb{Q}}$ (the law induced by $\overline{\mathit{X}}$) and then ${\overline{X}}_T\sim f\left(\mathit{x}\right)$.

Simulation from $\overline{\mathbb{Q}}$ can be done via rejection sampling. We can use a biased Brownian motion$\left\{{\overline{\mathit{B}}}_t;0\le t\le T\right\}$ as the proposal diffusion:

Definition 4.2 (Biased Brownian motion) The starting and ending points $\left({\overline{\mathit{B}}}_0,{\overline{\mathit{B}}}_T\right)$ follow a distribution with a density $h\left(\mathit{x},\mathit{y}\right)$, and $\left\{{\overline{\mathit{B}}}_t;0<t<T\right\}$ is a Brownian bridge given $\left({\overline{\mathit{B}}}_0,{\overline{\mathit{B}}}_T\right)$.

Under certain mild conditions, Dai [12] proved the following lemma.

Lemma 4.1 Let $\mathbb{Z}$ be the probability law induced by $\left\{{\overline{\mathit{B}}}_t;0\le t\le T\right\}$. By letting

$$h\left({\mathit{\omega }}_0,{\mathit{\omega }}_T\right)=g_2\left({\mathit{\omega }}_T\right)g_1\left({\mathit{\omega }}_0\right)\frac{1}{\sqrt{2\mathit{\pi T}}}{e}^{-\frac{\parallel {\mathit{\omega }}_T-{\mathit{\omega }}_0{\parallel }^2}{2T}}$$E6

we have

$$\frac{d\overline{\mathbb{Q}}}{d\mathbb{Z}}\left(\overrightarrow{\mathit{\omega }}\right)\propto exp\left[-\frac{1}{2}{\int }_0^T\left({‖\mathit{\alpha }‖}^2+div\mathit{\alpha }\right)\left({\mathit{\omega }}_t\right)\mathit{dt}\right]$$E7

where div is the divergence operator.

Condition 4.1 There exists $l>-\infty$ such that

$$\frac{1}{2}\left({‖\mathit{\alpha }‖}^2+\mathbf{div}\mathit{\alpha }\right)\left(\mathit{x}\right)-l\ge 0.$$E8

Under Condition 4.1 the ratio (8) can be rewritten as

$$\frac{d\overline{\mathbb{Q}}}{d\mathbb{Z}}\left(\overrightarrow{\mathit{\omega }}\right)\propto exp\left[-{\int }_0^T\left(\frac{1}{2}\left({‖\mathit{\alpha }‖}^2+\mathbf{div}\mathit{\alpha }\right)\left({\mathit{\omega }}_t\right)-l\right)\mathit{dt}\right],$$E9

which has a value no more than $1$.

Therefore we can use rejection sampling to simulate from $\overline{\mathbb{Q}}$, with proposal measure $\mathbb{Z}$. This acceptance probability (10) can be dealt with using similar methods as that in [9, 11]. The algorithm is presented below; see [12, 38] for more details.

Algorithm 4.1 (Simple distributed sampler)

Simulate $\left({\mathit{\omega }}_0,{\mathit{\omega }}_T\right)$ from density $h$	01
Simulate the biased Brownian bridge $\left({\overline{\mathit{B}}}_t,t\in \left(0,T\right)\right)$	02
Accept ${\mathit{\omega }}_T$ as a sample from $f$, with probability (6); If rejected, go back to step 01.	03

Such a method is a rejection sampling algorithm, but it does not require finding a hat function to bound the target density, which is usually the main challenge of the traditional rejection sampling for complicated target densities. The above algorithm uses $g_2$ as the proposal density function, which does not have to bound the target $f$. However, it requires a bound for the derivatives of the logarithm of the sub-densities (see Condition 4.1). This is usually easier to get in practice, since the logarithm of the posterior is usually easy to deal with. Also [12] noted that we should choose sub-densities $g_1$ and $g_2$ as similar as possible, in order to achieve high acceptance probability.

Dai [12] focused on the simple decomposition of $f=g_1{g}_2$, although it mentioned that for more general decomposition of $f=g_1{g}_2\cdots g_C$, a recursive method can be used. Unfortunately, a naive recursive method is very inefficient. A more sophisticated method is introduced in the following section.

4.2 Monte Carlo fusion for distributed analysis

A more efficient and sophisticated methods were proposed recently in [13], named as Monte Carlo fusion. Suppose that we consider

$$f\left(\mathit{x}\right)\propto g_1\left(\mathit{x}\right)\cdots g_C\left(\mathit{x}\right),$$E10

where each $g_c\left(\mathit{x}\right)$ ($c\in \left\{1,\dots ,C\right\}$) is a density (up to a multiplicative constant). Here $C$ denotes the number of parallel computing cores available in big data problems, and each $g_c\left(\mathit{x}\right)$ means the sub-posterior density based on a subset of the big data. In group decision problems, $C$ means the number of different decisions which should be combined and $g_c\left(\mathit{x}\right)$ stands for the decision from each group member.

Dai et al. [13] considered simulating from the following density on extended space,

$$g\left({\mathit{x}}^{\left(1\right)},\dots ,{\mathit{x}}^{\left(C\right)},\mathit{y}\right)=\prod _{c=1}^C\left[g_c^2\left({\mathit{x}}^{\left(c\right)}\right)\cdot p_c\left(\mathit{y}|{\mathit{x}}^{\left(c\right)}\right)\cdot \frac{1}{g_c\left(\mathit{y}\right)}\right],$$E11

which admits the marginal density $f$ for $\mathit{y}$. Here $p_c\left(\mathit{y}|{\mathit{x}}^{\left(c\right)}\right)$ is the transition density from ${\mathit{x}}^{\left(c\right)}$ to $\mathit{y}$ for the Langevin diffusion defined in Eq. (6) associated with each sub-density $g_c$.

Dai et al. [13] considered a rejection sampling approach with proposal density proportional to the function

$$h\left({\mathit{x}}^{\left(1\right)},\dots ,{\mathit{x}}^{\left(C\right)},\mathit{y}\right)=\prod _{c=1}^C\left[g_c\left({\mathit{x}}^{\left(c\right)}\right)\right]\cdot exp\left(-\frac{C\cdot \parallel \mathit{y}-\overline{\mathit{x}}{\parallel }^2}{2T}\right),$$E12

where $\overline{\mathit{x}}=C^{-1}{\sum }_{c=1}^C{\mathit{x}}^{\left(c\right)}$ and $T$ is an arbitrary positive constant.

Simulation from the proposal $h$ can be achieved directly. In particular, ${\mathit{x}}^{\left(1\right)},\dots {\mathit{x}}^{\left(C\right)}$ are first drawn independently from $g_1,\dots ,g_C$, respectively, and then $\mathit{y}$ is simply a Gaussian random variable centred on $\overline{\mathit{x}}$. This is a distributed analysis or divide-and-conquer approach. Detailed acceptance probabilities and rejection sampling algorithms can be found in [13].

The above fusion approach arises in modern statistical methodologies for ‘big data’. A full dataset will be artificially split into a large number of smaller data sets, and inference is then conducted on each smaller data set and combined (see, for instance, [39, 40, 41, 42, 43, 44, 45, 46]). The advantage for such an approach is that inference on each small data set can be conducted in parallel. Then the heavy computational cost of algorithms such as MCMC will not be a concern. Traditional methods suffer from the weakness that the fusion of the separately conducted inferences is inexact. However, the Monte Carlo fusion in [13] is an exact simulation algorithm and does not have any approximation weakness.

The above fusion approach also arises in a number of other settings, where distributed analysis came naturally. For example, in signal processing, distributed multi-sensor may be used for network fusion systems. Fusion approach arises naturally to combine results from different sensors [47].

5. Conclusion

Although perfect simulation usually refers to correcting the statistical errors for the samples drawn via MCMC, it actually covers a much wider area beyond CFTP. In fact for certain applications, it is often possible to construct other types of perfect sampling methods which are much more efficient than CFTP. For example, for the exact simulation of the posterior of simple mixture models, the geometric-arithmetic mean (GAM) method in [19] is much more efficient than CFTP in [48]. Details of GAM method is provided in Appendix. Also the random walk construction for exact simulation for random spanning trees [7] is much more efficient than the CFTP version.

Bayesian computational algorithms keep evolving, in particular under the current big data era. Although almost all newly developed algorithms are approximate simulation algorithms, perfect sampling is still one of the key wheel-driven forces for new Bayesian computational algorithms, and they usually can quickly motivate new class of ‘mainstream’ algorithms. More focus should be given to methods beyond CFTP, for example, the fusion type of algorithms.

The Monte Carlo fusion method has the potential to be used in many Bayesian big data applications. For example, for large car accident data, the response variable is usually a categorical variable representing the seriousness of the accident, and generalized linear regression model is often used. Under a Bayesian framework, we may estimate the posterior distribution for the regression parameters via such a fusion approach. Then the posterior mean, the posterior median, or other characteristics of the posterior distribution can be estimated using the Monte Carlo samples. Also such an algorithm is perfect sampling algorithm, and no convergence justification is needed, since it always provided realizations exactly from the target distribution.

Observations from a simple mixture model are assumed to be either discrete or continuous. The density function of an individual observation $y$ has the form

$$f\left(y;\mathit{p}\right)=\sum _{k=1}^K{p}_k{f}_k\left(y\right),\text{where}\sum _{k=1}^K{p}_k=1,\text{and}{p}_k>0,k=1,\dots ,K.$$E13

We assume that the component weights $\mathit{p}=\left(p_1,\dots ,p_K\right)$ are unknown parameters and the number of components, $K$, and the component densities, $f_k$, are all known. We focus on the perfect sampling from the posterior distribution of $\mathbf{p}$.

Suppose that we have $N$ observations, $y_1,\dots ,y_N$. Let $L_{\mathit{nk}}=f_k\left(y_n\right)$ and assume that the prior distribution of $\mathbf{p}$ is Dirichlet:

$${\pi }_0\left(\mathbf{p}\right)\propto \prod _{k=1}^K{p}_k^{{\alpha }_k-1},{\alpha }_k>0,k=1,\dots ,K.$$E14

Then the posterior distribution is given by

$$f\left(\mathbf{p}|y\right)\propto {\pi }_0\left(\mathbf{p}\right)\prod _{n=1}^N\left(\sum _{k=1}^K{p}_k{L}_{\mathit{nk}}\right)I_{\mathcal{X}}\left(\mathbf{p}\right),$$E15

where $\mathcal{X}=\left\{\mathbf{p}|{\sum }_{k=1}^K{p}_k=1,p_k>0,k=1,\dots ,K\right\}$.

There are several ways to carry out perfect sampling from Eq. (16). The first method is based on CFTP [48]. An alternative perfect sampling method for simple mixture models is introduced by [49]. The third approach is to use adaptive rejection sampling [3, 18], since the posterior is log-concave. We may also use the ratio-of-uniform method. However, none of these methods are more efficient than the geometric–arithmetic mean method in [19].

A.1 Geometric-arithmetic mean method

Suppose that ${\mathbf{p}}^{\ast }$, the MLE of $\mathbf{p}$ is unique and for simplicity, assume the prior ${\pi }_0\left(\mathbf{p}\right)$ is uniform. Define $a_{\mathit{nk}}=L_{\mathit{nk}}/\sum p_k^{\ast }{L}_{\mathit{nk}}$. The posterior density of $\mathbf{p}$ is then given by

$$f\left(\mathbf{p}|y\right)\propto h\left(\mathbf{p}|y\right)=\prod _{n=1}^N\left(\sum _{k=1}^K{p}_k{a}_{\mathit{nk}}\right)I_{\mathcal{X}}\left(\mathbf{p}\right),$$E16

where $\mathcal{X}$ is defined in (16).

Let $I_n$ be a random element of arg ${max}_k{L}_{\mathit{nk}}$. Define $A_j=\left\{n:I_n=j\right\}$ and let $\mathbf{n}=\left(n_1,\dots ,n_k\right)$ where $n_j$ is the number of elements in $A_j$.

Define $\mathbf{M}=\left\{m_{\mathit{jk}}\right\}$ with $m_{\mathit{jk}}=\left({\sum }_{n\in A_j}{a}_{\mathit{nk}}\right)/n_j$. If $n_j=0$, then set $m_{\mathit{jj}}=1$ and $m_{\mathit{jk}}=0$ for $j\ne k$. We now make two assumptions, which we will return to later on:

A:$\mathbf{M}$ is invertible.

B: The elements of $\mathbf{v}={\left({\mathbf{M}}^T\right)}^{-1}\mathbf{1}$ are all positive.

Under these assumptions, we will show that the following rejection sampler generates simulated values from the posterior distribution of $\mathbf{p}$. First we define $\mathbf{V}$ to be the diagonal matrix with diagonal elements ${\mathbf{v}}^T=\left(v_1\dots ,v_K\right)$.

Algorithm 6.1 (GAM sampler)

Sample $\mathbf{q}$ from the Dirichlet distribution with parameter $\mathbf{n}+\mathbf{1}$.	01
Sample $U$ from Unif$\left[0,1\right]$.	02
Calculate $\mathbf{p}$ with $\mathbf{p}={\mathbf{M}}^{-1}{V}^{-1}\mathbf{q}$.	03
If $U\le h\left(\mathbf{p}|y\right)/{\prod }_{j=1}^K{\left(q_j/v_j\right)}^{n_j}$,	04
Accept $\mathbf{p}$ and stop;	05
else	06
reject $\mathbf{p}$ and go to 01.	07

Proposition 6.1 Under assumptions A and B, Algorithm 6.1 samples $\mathbf{p}$ with probability density (17).

Proof: Since the geometric average is no larger than the arithmetic average, for $\mathbf{p}\in \mathcal{X}$, we have

$$h\left(\mathbf{p}|y\right)=\prod _{n=1}^N\left(\sum _{k=1}^K{p}_k{a}_{\mathit{nk}}\right)=\prod _{j=1}^K\prod _{n\in A_j}\left(\sum _{k=1}^K{p}_k{a}_{\mathit{nk}}\right)$$E17

$$\le \prod _{j=1}^K{\left(\frac{{\sum }_{n\in A_j}\left({\sum }_{k=1}^K{p}_k{a}_{\mathit{nk}}\right)}{n_j}\right)}^{n_j},$$E18

where in the case $n_j=0$, the product term is taken as $1$. So that, for $\mathbf{p}\in \mathcal{X}$, with $m_{\mathit{jk}}$ as previously defined, we have

$$h\left(\mathbf{p}|y\right)\le \prod _{j=1}^K{\left(\sum _{k=1}^K{p}_k{m}_{\mathit{jk}}\right)}^{n_j}$$E19

$$=\left[\prod _{j=1}^K{v}_j^{-n_j}\right]\prod _{j=1}^K{\left(v_j\sum _{k=1}^K{p}_k{m}_{\mathit{jk}}\right)}^{n_j}$$E20

$$=\prod _{j=1}^K{\left(q_j/v_j\right)}^{n_j},$$E21

where $q_j=v_j{\sum }_{k=1}^K{p}_k{m}_{\mathit{jk}},j=1,\dots ,K$ or equivalently $\mathbf{q}\mathbf{=}\mathbf{VMp}$.

Since $v_j>0$ and ${\sum }_{k=1}^K{p}_k{m}_{\mathit{jk}}>0$, it follows that $q_j>0$ for $j=1,\dots ,K$. Furthermore

$$\sum _{j=1}^K{q}_j=\sum _{j=1}^K{v}_j\sum _{k=1}^K{p}_k{m}_{\mathit{jk}}={\mathbf{p}}^T{\mathbf{M}}^T\mathbf{v}\mathbf{=}{\mathbf{p}}^T\mathbf{1}=1,$$

since ${\mathbf{M}}^T\mathbf{v}\mathbf{=}\mathbf{1}$, from the definition of $\mathbf{v}$. It follows that $\mathbf{p}\in \mathcal{X}$ implies $\mathbf{q}\in \mathcal{X}$, so that

$$h\left(\mathbf{p}|y\right)I_{\mathcal{X}}\left(\mathbf{p}\right)\le \prod _{j=1}^K{\left(q_j/v_j\right)}^{n_j}{I}_{\mathcal{X}}\left(\mathbf{q}\right).$$E22

Note that the right-hand side of Eq. (22) is proportional to a Dirichlet distribution with parameters $\left(n_1+1,\dots ,n_K+1\right)$.

Rejection sampling then proceeds as usual:

• A sample $\mathbf{q}$ is drawn from Dirichlet$\left(\mathbf{n}+\mathbf{1}\right)$.

• The value $\mathbf{p}\mathbf{=}{\mathbf{M}}^{-1}{\mathbf{V}}^{-1}\mathbf{q}$ is calculated.

• It is accepted with probability $h\left(\mathbf{p}|y\right)I_{\mathcal{X}}\left(\mathbf{p}\right)/{\prod }_{j=1}^K{\left(q_j/v_j\right)}^{n_j}$.

We now return to assumptions A and B. Suppose that $\mathbf{M}$ is invertible but the elements of $\mathbf{v}={\left({\mathbf{M}}^T\right)}^{-1}\mathbf{1}$ are not all positive. In this case, let

${\alpha }_k=\frac{1}{N}{\sum }_{n=1}^N{a}_{\mathit{nk}},\alpha =\underset{k}{max}\left\{{\alpha }_k\right\},\mathbf{v}={\left(\mathit{\alpha N}\right)}^{-1}\mathbf{n}\text{and}\tilde{\mathbf{M}}=\alpha \mathbf{M}\mathrm{\Delta },$

where $\mathrm{\Delta }$ is a diagonal matrix with diagonal elements $\left(1/{\alpha }_1,\dots ,1/{\alpha }_K\right)$. Note that $\mathbf{v}={\left({\tilde{\mathbf{M}}}^T\right)}^{-1}\mathbf{1}>0$. Algorithm 6.1 and its proof can then be modified by replacing $\mathbf{M}$ by $\tilde{\mathbf{M}}$.

Suppose now that assumption A does not hold, i.e. $\mathbf{M}$ is not invertible. This can be remedied by adding positive quantities to the diagonal elements of $\mathbf{M}$. This also provides an alternative way of ensuring that the elements of $\mathbf{v}$ are positive.

A.2 Dirichlet priors and pseudo data

Suppose that the prior ${\pi }_0\left(\mathbf{p}\right)$ is Dirichlet$\left({\alpha }_1+1,\dots ,{\alpha }_K+1\right)$, where ${\alpha }_i:i=1,\dots ,K$ are positive, integers and let $A={\sum }_{j=1}^K{\alpha }_j$. The prior can be synthesized by introducing pseudo data, ${\tilde{a}}_{\mathit{mk}},m=1,\dots ,A;k=1,\dots ,K$, defined as follows:

$${\tilde{a}}_{\mathit{mk}}=\left\{\begin{array}{ll}1& \text{if}{\sum }_{j=1}^{k-1}{\alpha }_j+1\le m\le {\sum }_{j=1}^k{\alpha }_j\\ 0& \text{otherwise},\end{array}\right.$$E23

since

$${\pi }_0\left(\mathbf{p}\right)\propto \prod _{k=1}^K{p}_k^{{\alpha }_k}=\prod _{m=1}^A\left(\sum _{k=1}^K{\tilde{a}}_{\mathit{mk}}{p}_k\right).$$E24

With the Dirichlet prior, the posterior distribution given by Eq. (16) can be written as

$$f\left(\mathbf{p}|y\right)\propto \prod _{l=1}^{N+A}\left(\sum _{k=1}^K{p}_k{\overline{a}}_{\mathit{lk}}\right)I_{\mathcal{X}}\left(\mathbf{p}\right),$$E25

where $\left\{{\overline{a}}_{\mathit{lk}},l=1,\dots ,N+A\right\}$ contains the real data $\left\{a_{\mathit{nk}},n=1,\dots ,N\right\}$ and the pseudo data $\left\{{\tilde{a}}_{\mathit{mk}},m=1,\dots ,A\right\}$.

The posterior distribution (26) has the same form as Eq. (17). Therefore GAM can be used to sample realizations from the posterior distribution (26).

A.3 Simulation results and discussion

We compare the running time of mixture models with sample sizes ($N$) and different number of components ($K$) in Table 1. The components have specified normal distributions with means $\mu =\left({\mu }_1,\dots ,{\mu }_K\right)$ and variances ${\sigma }^2=\left({\sigma }_1^2,\dots ,{\sigma }_K^2\right)$. The prior on $\mathbf{p}$ is uniform. We sample 10,000 realizations from the posterior of the models.

$K$	$3$	$3$	$4$	$4$	$6$	$6$
$N$	$400$	$1000$	$400$	$1000$	$400$	$1000$
Fearnhead’s	242 s	3610 s	*	*	*	*
Leydold’s	$\le$1 s	3.6 s	*	*	*	*
RoU	16:11 s	28:16 s	31:18 s	68:33 s	88:60 s	152:76 s
GAM	4 s	9 s	11 s	16 s	6 s	11 s
GAM AP	0.7472	0.7509	0.2433	0.3088	0.5325	0.5505

Table 1.

Running times (in s).

Fearnhead’s algorithm, Leydold’s algorithm and ratio-of-uniform. GAM method. GAM acceptance probability. The * indicates that Fearnhead’s method, and Leydold’s method will not run on a standard desktop when $K=4$ and $K=5$.

When $K=3,4$, we simulate $N$ observations from a three-component normal mixture with $\mu =\left(\mathrm{0,0,2}\right)$, ${\sigma }^2=\left(\mathrm{1,4,1}\right)$ and mixture weight ${\mathbf{p}}_0=\left(1/2,1/3,1/6\right)$. We then either sample from the posterior distribution of $\mathbf{p}$ using the same distributional components in the case $K=3$ or sample from the posterior distribution of $\mathbf{p}$ with an additional component having mean ${\mu }_4=4$ and variance ${\sigma }_4^2=4$, in the case $K=4$.

When $K=5$, observations are simulated from the normal mixture distribution with components having means $\mu =\left(-\mathrm{2,0,4,2,3}\right)$, variances ${\sigma }^2=\left(\mathrm{1,1,4,1,4}\right)$ and ${\mathbf{p}}_0=\left(\mathrm{0.35,0.3,0.1,0.2,0.05}\right)$. Samples from the posterior distribution of $\mathbf{p}$ are drawn assuming the same components. Similar calculations are carried out for $K=6$, where $\mu =\left(\mathrm{0,3,2},-2,-4,5\right)$, variances ${\sigma }^2=\left(\mathrm{1,1,1,1,1,4}\right)$ and ${\mathbf{p}}_0=\left(\mathrm{0.05,0.3,0.3,0.1,0.08,0.17}\right)$, again assuming that the component distributions are known.

From Table 1, we can see that the GAM algorithm, while using very little memory, is highly efficient in running time. The last row of the table is the estimated acceptance probability of the GAM algorithm. The algorithm is very efficient when the component densities are known. We can see this not only by simulation but also from theoretical considerations, as follows.

A.3.1 Explanation of efficiency

When $\mathbf{v}={\left({\mathbf{M}}^T\right)}^{-1}\mathbf{1}>0$, we are able to use $\mathbf{M}$ directly without modification to construct the hat function, thereby speeding up the calculations. In the simulations of the previous section, this was always found to be the case. Now we explain why this should be so.

If the maximum likelihood estimator of $\mathbf{p}$ is consistent, then when $N\to \infty$,

$$\frac{1}{n_j}\sum _{n\in A_j}\left|\frac{L_{\mathit{nk}}}{{\sum }_{k=1}^K{p}_k^{\ast }{L}_{\mathit{nk}}}-\frac{L_{\mathit{nk}}}{{\sum }_{k=1}^K{p}_k{L}_{\mathit{nk}}}\right|\stackrel{p}{\to }0.$$E26

Assuming sufficient regularity, we also have

$$m_{\mathit{jk}}=\frac{{\sum }_{n\in A_j}{a}_{\mathit{nk}}}{n_j}$$E27

$$=\frac{1}{n_j}\sum _{n\in A_j}\frac{L_{\mathit{nk}}}{{\sum }_{k=1}^K{p}_k^{\ast }{L}_{\mathit{nk}}}$$E28

$$\stackrel{p}{\to }E\left(\frac{f_k\left(Y\right)}{f\left(Y\right)}|{\mathcal{L}}_j\right)$$E29

$$=\frac{{\int }_{{\mathcal{L}}_j}{f}_k\left(y\right)\mathit{dy}}{{\gamma }_j},\mathrm{as}N\to \infty ,$$E30

where $Y$ has density $f\left(y\right)={\sum }_{k=1}^K{p}_k{f}_k\left(y\right)$, ${\mathcal{L}}_j=\left\{y|f_j\left(y\right)\ge f_k\left(y\right),k=1,\dots ,K\right\}$ and ${\gamma }_j={\int }_{{\mathcal{L}}_j}f\left(y\right)\mathit{dy}$.

Therefore

$${\mathbf{M}}^T\stackrel{p}{\to }{\mathbf{W}}^T,$$E31

where

$$\mathbf{W}=\left\{w_{\mathit{jk}}\right\},w_{\mathit{jk}}=\frac{{\int }_{{\mathcal{L}}_j}{f}_k\left(y\right)\mathit{dy}}{{\gamma }_j}.$$E32

Let $\overline{\mathbf{v}}=\left({\gamma }_1,\dots ,{\gamma }_K\right)$, then

$${\mathbf{W}}^T\overline{\mathbf{v}}=\left[\begin{array}{c}{\sum }_{j=1}^K{\int }_{{\mathcal{L}}_j}{f}_1\left(y\right)\mathit{dy}\\ ⋮\\ {\sum }_{j=1}^K{\int }_{{\mathcal{L}}_j}{f}_K\left(y\right)\mathit{dy}\end{array}\right]=\left[\begin{array}{c}\int f_1\left(y\right)\mathit{dy}\\ ⋮\\ \int f_K\left(y\right)\mathit{dy}\end{array}\right]=\mathbf{1},$$E33

where the second equal sign is because ${\cup }_{j=1}^K{\mathcal{L}}_j=\left(-\infty ,\infty \right)$. So, there exists $\overline{\mathbf{v}}>0$, satisfying ${\mathbf{W}}^T\overline{\mathbf{v}}=\mathbf{1}$. Using Eq. (32), we can conclude that when $N$ is large enough, there also exists $\mathbf{v}\approx \overline{\mathbf{v}}>0$, satisfying ${\mathbf{M}}^T\mathbf{v}=\mathbf{1}$.

Since ${\gamma }_j={\int }_{{\mathcal{L}}_j}f\left(y\right)\mathit{dy}$ and $n_j=\#\left\{A_j\right\}$, we have $n_j/N\stackrel{p}{\to }{\gamma }_j.$ When each $n_j,j=1,\dots ,K$ is large, if a random sample $\mathbf{q}$ is drawn from a Dirichlet distribution with parameter $\mathbf{n}+\mathbf{1}$, then each $q_j$ is approximately equal to $n_j/N\approx {\gamma }_j$. Furthermore, $v_j\approx {\gamma }_j$, so $\mathbf{q}$ satisfies

$${\mathbf{V}}^{-1}\mathbf{q}\approx \mathbf{1},$$E34

and then,

$$\mathbf{p}={\mathbf{M}}^{-1}{\mathbf{V}}^{-1}\mathbf{q}\approx {\mathbf{M}}^{-1}\mathbf{1}={\mathbf{p}}^{\ast }.$$E35

If $\mathbf{p}$ is approximately equal to the mode ${\mathbf{p}}^{\ast }$, the two sides of the inequality,

$$h\left(\mathbf{p}|y\right)=\prod _{n=1}^N\left(\sum _{k=1}^K{p}_k{a}_{\mathit{nk}}\right)\le \left[\prod _{j=1}^K{v}_j^{-n_j}\right]\prod _{j=1}^K{q}_j^{n_j},$$E36

are approximately equal as well. Thus, the closer the sampled realization $\mathbf{p}$ is to ${\mathbf{p}}^{\ast }$, the larger the acceptance probability is. So the algorithm runs very rapidly, since the sampled values of $\mathbf{p}$ are always around the mode ${\mathbf{p}}^{\ast }$.

This algorithm requires calculating the MLE, which can be performed very quickly since the likelihood function is log-concave. In fact an approximate guess for ${\mathbf{p}}^{\ast }$ will suffice. The more accurate the guess is, the more efficient the algorithm will be.

The method performs well when the component densities are correctly specified, as explained in the previous section. For these same reasons, we would expect the algorithm to perform poorly under misspecification. Details of robustness to misspecification can be found in [19].