Plateau proposal distributions for adaptive component-wise multiple-try metropolis

Lau, F. Din-Houn; Krumscheid, Sebastian Article — Published Version Plateau proposal distributions for adaptive componentwise multiple-try metropolis METRON Provided in Cooperation with: Springer Nature Suggested Citation: Lau, F. Din-Houn; Krumscheid, Sebastian (2022) : Plateau proposal distributions for adaptive component-wise multiple-try metropolis, METRON, ISSN 2281-695X, Springer Milan, Milano, Vol. 80, Iss. 3, pp. 343-370, https://doi.org/10.1007/s40300-022-00235-y This Version is available at: https://hdl.handle.net/10419/310010 Standard-Nutzungsbedingungen: Die Dokumente auf EconStor dürfen zu eigenen wissenschaftlichen Zwecken und zum Privatgebrauch gespeichert und kopiert werden. Sie dürfen die Dokumente nicht für öffentliche oder kommerzielle Zwecke vervielfältigen, öffentlich ausstellen, öffentlich zugänglich machen, vertreiben oder anderweitig nutzen. Sofern die Verfasser die Dokumente unter Open-Content-Lizenzen (insbesondere CC-Lizenzen) zur Verfügung gestellt haben sollten, gelten abweichend von diesen Nutzungsbedingungen die in der dort genannten Lizenz gewährten Nutzungsrechte. Terms of use: Documents in EconStor may be saved and copied for your personal and scholarly purposes. You are not to copy documents for public or commercial purposes, to exhibit the documents publicly, to make them publicly available on the internet, or to distribute or otherwise use the documents in public. If the documents have been made available under an Open Content Licence (especially Creative Commons Licences), you may exercise further usage rights as specified in the indicated licence. https://creativecommons.org/licenses/by/4.0/ METRON (2022) 80:343–370 https://doi.org/10.1007/s40300-022-00235-y Plateau proposal distributions for adaptive component-wise multiple-try metropolis F. Din-Houn Lau1 ·Sebastian Krumscheid2 Received: 15 June 2021 / Accepted: 18 June 2022 / Published online: 15 July 2022 © The Author(s) 2022 Abstract Markov chain Monte Carlo (MCMC) methods are sampling methods that have become a commonly used tool in statistics, for example to perform Monte Carlo integration. As a consequence of the increase in computational power, many variations of MCMC methods exist for generating samples from arbitrary, possibly complex, target distributions. The performance of an MCMC method, in particular that of a Metropolis–Hastings MCMC method, is predominately governed by the choice of the so-called proposal distribution used. In this paper, we introduce a new type of proposal distribution for the use in Metropolis–Hastings MCMC methods that operates component-wise and with multiple trials per iteration. Specifically, the novel class of proposal distributions, called Plateau distributions, does not overlap, thus ensuring that the multiple trials are drawn from different regions of the state space. Furthermore, the Plateau proposal distributions allow for a bespoke adaptation procedure that lends itself to a Markov chain with efficient problem dependent state space exploration and favourable burn-in properties. Simulation studies show that our novel MCMC algorithm outperforms competitors when sampling from distributions with a complex shape, highly correlated components or multiple modes. Keywords Component-wise Metropolis–Hastings ·Multiple-try Metropolis ·Adaptive Markov chain Monte Carlo ·Plateau proposal distribution 1 Introduction Markov chain Monte Carlo (MCMC) methods are essentially used to perform Monte Carlo integration, which has become a standard statistical tool. Specifically, MCMC methods produce samples from a target distribution πby constructing an ergodic Markov chain whose stationary distribution is π. Typically, MCMC methods are used when it is difficult to sample BSebastian Krumscheid [email protected] F. Din-Houn Lau [email protected] 1Cheshunt, Hertfordshire, UK 2Department of Mathematics, RWTH Aachen University, 52062 Aachen, Germany 123 344 F. D-H. Lau, S. Krumscheid from the target distribution directly, e.g. when the normalisation constant is unknown. There are many ways to construct this Markov chain, which have led to many variations of MCMC methods; see, e.g., [2]. The classic MCMC method is the Metropolis-Hastings algorithm [21]. At each iteration, the Metropolis-Hastings algorithm is designed to update the entire current state (i.e., all components of the random vector generated at the previous iteration) at once. However, updating individual components, or subsets of components, is possible. Indeed, this type of component-wise updating was initially proposed in [21], but did not receive much attention at first. Recently, component-wise updates in MCMC methods have been more established in the literature; see, e.g., [14] where the convergence rates of such methods are studied. In this paper, we focus on updating individual components. Using individual component updates is a way of sampling from the (lower dimensional) conditional distributions of the target, provided the conditional distributions are known. However, the conditional distributions are unknown in practise. Modelling the conditional distributions using parametric models leads to an inflexible approximation, whereas non-parametric models do not scale well with the number of dimensions. For these reasons, we focus our attention on independent componentwise updates. Another variant of MCMC sampling, separate from component-wise updating, is the multiple-try method [15], where several proposals or trials are suggested at each iteration. The motivation behind the multiple-try method is that more of the space is explored at the expense of an increased computational cost (i.e., proposal generation and evaluation of acceptance criterion). Several variations of multiple-try MCMC algorithms exist in the literature; see [20]. Recently, in [27] the authors introduce a component-wise, multiple-try MCMC method with Gaussian proposals for each trial, where each univariate proposal has the same mean but different variances. In this work, we introduce a new class of proposal distributions for use in a component-wise, multiple-try MCMC method. These proposals, called Plateau distributions, do not overlap to exploit the multiple-try nature of the method. Indeed, by using proposals that do not overlap for each trial, the Markov chain is forced to explore different parts of the state space. Conversely, using proposals that overlap, such as Gaussians with the same mean and different variances, can lead to an inefficient algorithm, as the trials tend to be from a similar region of the state space. The idea of Plateau proposals is intuitive, easy to implement and leads to good results, in the sense of exploring the state space. Moreover, using the Plateau proposals leads to a reversible Markov chain with the target distribution as its invariant distribution, see, e.g., [27]. As is common for most proposals used in MCMC methods, the Plateau distributions depend on parameters that need to be selected with respect to the target distribution to obtain an effective algorithm. Adaptation of MCMC methods typically entails tuning the parameter(s) of a class of proposal distributions, e.g., the variance parameter in a Normal distribution, in order to improve the convergence properties of the Markov chain. For instance, in [11], a Gaussian proposal distribution is used whose variance (or covariance) is adapted using the previously generated states of the Markov chain. In this work, we propose an adaptation procedure that is designed for use with the nonoverlapping Plateau proposals. The Plateau proposals together with the associated adaptation procedure are illustrated in Fig. 1. Suppose the MCMC algorithm is initiated with multiple-try proposals whose distributions are presented by the different coloured lines in Fig. 1a. These proposals operate independently on each (univariate) component in the target space. As the Markov chain evolves, the number of times each Plateau’s proposed candidate state is selected 123 Plateau proposal distributions... 345 −20 −10 0 10 20 0.0 0.2 0.4 0.6 0.8 1.0 x y (a) −20 −10 0 10 20 0.0 0.2 0.4 0.6 0.8 1.0 x y (b) −20 −10 0 10 20 0.0 0.2 0.4 0.6 0.8 1.0 x y (c) Fig. 1 Illustration of 4 Plateau proposal distributions (probability density functions) represented by the different colours and the adaptation procedure. The initial Plateau distributions are presented in (a); If the outermost (blue) proposals are selected frequently, the proposals are adapted to (b); If the innermost (black) proposals are selected frequently, the proposals are adapted to (c) (colour figure online) is recorded. If the innermost or outermost Plateaus are overly selected, then the set of Plateau proposals are re-organised as depicted in Figs. 1band1c. More precisely, if the innermost (black) proposal is overly selected, then one halves the width of each plateau and shifts the set of proposals to remove the gaps. A similar procedure is performed if the outermost (blue) proposals is overly selected by doubling the width of the plateaus. This procedure appropriately scales the set of Plateau proposals to the components of the target distribution. This intuitive adaptation procedure makes explicit use of the non-overlapping feature of the Plateau proposals. Specifically, if a sample is selected from a particular Plateau distribution, then that sample could not have been obtained by sampling from any other Plateau distribution since their supports do not overlap. The advantages of using non-overlapping proposals, as opposed to overlapping proposals, is explored in Sects. 5and 7. The mathematical definition of the Plateau proposals and the pseudocode for the adaptation procedure are presented in Sects. 3and 4respectively. The works in [9]and[5] share a similar objective with our proposed work, namely: draw samples that are far apart from each other to facilitate an effective exploration of the state space. For example, in [9] the authors introduce an independent (i.e., componentwise), yet single-try, Metropolis–Hastings method using a Normal mixture distribution as proposal. The mixture distribution is adapted using a k-means clustering approach in order to explore the space more efficiently. In [5] an improved state space exploration is achieved by combining variance reduction techniques with multiple trials. Specifically, it involves drawing multiple trails from a single, well-chosen proposal distribution that is based on, e.g., Latin Hypercube sampling. The performance of the approaches in [9]and[5] depends on 123 346 F. D-H. Lau, S. Krumscheid a well-chosen transformation. Conversely, the approach introduced in this work does not require such a transformation, as the Plateau proposals reside in the state space. Further, our adaptation procedure tunes the Plateau proposals parameters as the MCMC runs and does not require additional clustering (i.e., optimisation) algorithms to be performed. To summarise, the setting of this work is to propose a new MCMC algorithm along with an adaptation procedure, which is simple to implement, intuitive and can be used without using explicit derivative information (e.g., Hamiltonian Monte Carlo [22]) or conditional distributions (e.g., Gibbs sampling [4]) of the target. Therefore, a general-purpose MCMC algorithm is required. To this end, we introduce Plateau proposals as a general class of proposal distributions for use in a component-wise multiple-try MCMC methods. The combination of the non-overlapping characteristic of Plateau proposals with the multiple-try approach and a bespoke adaptation procedure leads to good results for a variety of target distributions. Examples of target distributions that will benefit from an adaptive multiple-try MCMC based on Plateau proposals include, in particular, multimodal distributions, such as mixture distributions, and, more generally, target distributions with complex shapes, such as those with highly correlated components or Bayesian posteriors. The remainder of this paper is organised as follows. In Sect. 2.1 a generic component-wise multiple-try algorithm is presented. The novel class of Plateau proposals is introduced and discussed in Sect. 3. In Sect. 4we discuss how to adaptively select the parameter of the Plateau proposals and offer a detailed algorithmic description of the complete method. The performance of our new method is then compared with other MCMC methods in Sect. 6. Finally, a commentary on improvements and a summary are provided in Sect. 8. 2 Component-wise update with multiple trials Let πbe a probability density function, π:X→R+,whereX⊆Rd. Our main interest is to sample from π; this is the target distribution. We assume that sampling directly from π is difficult or impossible, for example because πmay only be known up to a multiplicative constant. In order to sample from πwe use MCMC methods. In the Metropolis-Hastings algorithm, one of the simplest MCMC methods, a candidate Yfor the chain’s next state Xn+1 is drawn from the probability density function (PDF) T(x,·):X→R+, based on the current state Xn=x∈Xof the Markov chain. The proposal T(x,·)is typically easier to sample from than the target distribution. The density T(x,·)is the conditional density given the current value x. This density is commonly known as the proposal distribution. For example, the random-walk proposal [12,21] uses a multivariate Normal distribution, which we will write as T(x,·)=N(x,σ2Id)with σ>0 given and Idbeing the d×didentity matrix. Another example is Nx+τ∇lnπ(x),2τId)with τ>0 fixed, which is the proposal used in the Metropolis-adjusted Langevin algorithm [25]. The realisation of the Markov chain’s next state is then selected according to an accept-reject procedure designed such that the resulting Markov chain’s stationary distribution is π. Constructing an MCMC method that produces good results typically depends on the choice of proposal distribution. In particular, the choice of the proposal may significantly affect the properties of the MCMC method, including the speed of convergence to equilibrium and mixing properties [26]. Typically, the proposal distribution is selected from some family of distributions, e.g., from the family of Normal distributions. The performance of MCMC methods depends on appropriately scaling of the proposal distribution, especially for highdimensional target distributions [23]. 123 Plateau proposal distributions... 347 Instead of proposing a single multivariate candidate from Y∼T(x,·)by updating all components of the current state Xn=xsimultaneously via the (global) proposal distribution T(x,·), it is also possible to split the state space Xinto its individual components (or small groups of components) and propose candidates for each component (or group of components) independently. This local (or projected) approach is intuitive, computationally efficient, and reduces the problem of selecting a multidimensional proposal into lower dimensional proposals that are easier to handle. MCMC methods using these types of updates are called component-wise MCMC methods, which may use a potentially different proposal distribution per component of the state [8, Ch. 1]. In this work, we focus on one-dimensional, component-wise proposals. Considering independent one-dimensional, component-wise proposals is equivalent to the proposal distribution T(x,·)given x=(x1,...,xd)∈Rdbeing separable, in the sense that T(x,y)≡T1(x1,y1)×···×Td(xd,yd)= d  k=1 Tk(xk,yk), (1) for any y=(y1,...,yd)∈Rd. Here, Tk(x,·):R→R+,k=1,...,dwith x∈R, denotes the one-dimensional proposal density given xused to draw the candidate for component k. That is, the proposed (global) candidate Y=(Y1,...,Yd)∈Rdis obtained by sampling each Yk∼Tk(xk,·)mutually independent of any other Y. By construction proposing candidates component-wise does not account for correlations between the components. Consequently, the component-wise proposed candidates Ymay not be good representatives of the target distribution π(i.e., most candidates Ywill be rejected, resulting in a very low acceptance rate), if πhas highly correlated components. Therefore, these candidates with independently sampled components may lead to a poor state space exploration and thus to a poor performance of the MCMC method. To mitigate this issue, we will combine component-wise proposals with multiple trials. The multiple-try technique proposes many candidates from a proposal distribution, amongst which the “best” one is selected. Each trial may be proposed from a different proposal. Thus, in combination with component-wise proposals, we generate Mindependent candidates (or Mtrials) for each separate component k=1,...,d.LetTj,k(x,·):R→R+,j=1,...,M, denote the proposal PDF of the jth trial for the kth component given xk=x. 2.1 Generic procedure of component-wise multiple-try metropolis A generic component-wise multiple-try algorithm is now described – the full pseudocode for generating Nsamples from the target distribution πis presented in Algorithm 1. Each iteration within the MCMC algorithm draws multiple trials and then performs an acceptancerejection step for each component sequentially. We now describe the intuition behind the steps of the procedure. Suppose that the state of the chain at the beginning of the nth iteration is x=(x1,...,xd). For the first component (i.e., k=1), Mindependent trials, z1,...,zM∈R, are drawn from Tj,1(x1,·),j=1,...,M. These trials are then weighted according to wj,1(zj,x)=π(zj;x[−1])Tj,1(x1,zj)λj,1(x1,zj), where (z;x[−i])∈Rddenotes the vector that is identical to xexcept for its ith component which is replaced by z∈R,thatis(z;x[−i])=(x1,...,xi−1,z,xi+1,...,xd). 123 348 F. D-H. Lau, S. Krumscheid Algorithm 1: Generic Component-wise Multiple-Try Metropolis Input: number of trials M; number of MCMC realisations N; starting position x0∈Rd; target distribution π(possibly un-normalised); proposal distributions Tj,k 1: Let X0=x0=x. 2: for n=1,...,Ndo 3: for k=1,...,ddo 4: Propose Mtrials: zj∼Tj,k(xk,·)for j=1,...,M. 5: Compute the trial weights wj,k(zj,x)=π(zj;x[−k])Tj,k(xk,zj)λ j,k(xk,zj), j=1,...,M. 6: Draw y∈{z1,...,zM}randomly with probability proportional to w1,k,...,wM,k. 7: Draw x∗ j∼Tj,k(y,·)for j=1,...,M−1 and let x∗ M=xk. 8: Let y=(y;x[−k])and compute α=min 1,w1,k(z1,x)+···+wM,k(zM,x) w1,k(x∗ 1,y)+···+wM,k(x∗ M,y). Draw r∼Uniform(0,1). 9: if r<αthen 10: Accept Xn=y=(y;x[−k])and set x=y. 11: else 12: Xn=x. 13: end if 14: end for 15: end for 16: return X1,...,XN The functions λj,k(x,y)with x,y∈Rfor any k=1,...,d, are non-negative, symmetric functions in xand ywhich are selected by the user. It is required that λj,k(x,y)>0 whenever Tj,k(x,y)>0. Each trial zj,j=1,...,M, has an associated weight wj,1(zj,x).A candidate for the first component of the chain’s next state is then randomly selected amongst all trials zj(j=1,...,M)according to these weights. The selected candidate is then accepted or rejected. The remaining components k=2,3,...,dof the chain’s state are updated in order in a similar fashion; Algorithm 1details the full pseudocode of the MCMC method. In step 5 in Algorithm 1, the accepted candidate for the kth component is selected randomly from Mtrails z1,...,zMwith weights proportional to wj,k(zj,x)=π(zj;x[−k])Tj,k(xk,zj)λj,k(xk,zj), j=1,...,M. The weight depends on the function λj,k, for which there are a number of choices in the multiple-try literature [15], such as Tj,k(x,y)+Tj,k(y,x) 2−1 ,Tj,k(x,y)Tj,k(y,x)−βand 1. In work [27] the authors suggest using λj,k(x,y)=Tj,k(x,y)y−xα,(2) where α=2.9 based upon a simulation study focusing on large moves in the state spaces. For our proposal algorithm we use (2) with λj,k. In simulations, not reported here, we found that 123 Plateau proposal distributions... 349 α=2.5 performed best in terms of the mean squared error for a variety of target distributions. Consequently, we use λj,kin (2) with α=2.5 for the remainder of the paper unless stated otherwise. We note that it is beyond the scope of this paper to propose a particular form for the class of functions λj,kdue to the problem dependent nature of this choice. Instead, we advocate that users perform a trial run with a variety of λj,kto determine which is best suit for their application and performance metric. 3 Non-overlapping proposal distributions In principle, any family of proposals Tj,kcan be used in Algorithm 1. However, a careful choice of the proposals can lead to a more efficient algorithm. Recall that Tj,kis the jth trial proposal distribution for the kth component. The motivation of using multiple trials is to explore a larger region of the state space than is achieved by using a single proposal. Therefore, it would not be beneficial to sample trials (from T1,k,...,TM,k) that are similar. To illustrate this point, suppose that for a fixed component kwe use the proposal distributions Tj,k(x,·)=N(x,(jσ)2)for j=1,...,M=5 with known σ>0 and take x=0 without loss of generality. The probability density functions of these proposals with σ=1 are presented in Fig. 2. As illustrated in Fig. 2, these proposals are very similar. Indeed, 99% of T1,k’s density mass lies within the interval J=−σ−1(0.995), σ−1(0.995), where −1is the inverse of the cumulative distribution function of a standard Normal so that −1(0.995)≈2.6. For a candidate from the second proposal distribution Y2,k∼T2,k(0,·), we then have P(Y2,k∈J)≈0.8foranyσ>0. That is, draws from T1,kand T2,kwill be located in the same region with high probability. Similar arguments hold for the wider Gaussian proposals. Thus, draws from these Gaussian proposals will tend to be similar, thus leading to an inefficient use of the multiple-try technique. To avoid proposing candidates from similar regions, we seek densities which do not overlap (or overlap to a small degree). Specifically, we advocate using proposals of the type illustrated in Fig. 3b. That is, each trial for each component-wise proposal distribution combines uniform distributions with exponentially decaying tails. The amount of overlap between different proposals is controlled −10 −5 0 5 10 0.0 0.2 0.4 0.6 y Tj, k(0,y) j=1 j=2 j=3 j=4 j=5 Fig. 2 Probability density function of 5 Normal distributions with zero mean and standard deviations jσfor σ=1, j=1,...,5 123 350 F. D-H. Lau, S. Krumscheid y f(y;μ, δ, σ1, σ2) 1 C 2δ μ−δ μ μ+δ (a) −10 −5 0 5 10 0.0 0.2 0.4 0.6 y Tj, k(0, y) j=1 j=2 j=3 j=4 j=5 (b) Fig. 3 PDF of Plateau proposal distributions with M=5. The five different proposals each have a different colour through how fast the tails decay. Further, note that there is no gap between two contiguous distributions, i.e., no interval where neither distribution is likely to be sampled from. The family of Plateau proposals are now defined. We begin by introducing the PDF (see Fig. 3a) f(y;μ, δ, σ1,σ 2)=1 C ⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ exp −1 2σ2 1 [y−(μ −δ)]2for y<μ−δ 1forμ−δ≤y≤μ+δ exp −1 2σ2 2 [y−(μ +δ)]2for y>μ+δ (3) where C=2πσ2 1 2+2πσ2 2 2+2δ denotes the normalisation constant. For each component kwith given value xk=x,wethen set the PDF of each trial proposal as Tj,k(x,y)=⎧ ⎪ ⎨ ⎪ ⎩ f(y;x,δ 1,σ,σ) j=1 1 2f(y;x−2(j−1)δ1−δ,δ,σ,σ) +1 2f(y;x+2(j−1)δ1+δ,δ,σ,σ) j=2,...,M−1, 1 2f(y;x−2(M−1)δ1−δ,δ,σ0,σ)+1 2f(y;x+2(M−1)δ1+δ,δ,σ,σ1)j=M for some values of δ1,δ,σ,σ 0,σ 1>0. The M=5 trial proposals shown in Fig. 3bcor- respond to δ1=δ=1, σ=0.05 and σ0=σ1=0.5. We shall refer to the proposals of this type as Plateau proposals, given the shape of their PDFs. The δ1parameter controls the width of the central Plateau centred at the current state x.Theδparameter is the width of the other Plateaus. The σvalue controls the decay of the tails either side of the inner Plateaus. The outer tails for the Mth proposal are described by σ0and σ1. To compare with the earlier calculations for coverage probabilities for the Gaussian proposals; 99% of the density of T1,kwith x=0, δ=1, and σ=0.5 lies in the interval J=(−2.11,2.11). Suppose that Y2,k∼T2,k(0,·),thenP(Y2,k∈J)≈0.43, which is reduced by almost a factor of two compared to the overlapping Gaussian proposal. Further, if σ=0.25 then P(Y2,k∈J)≈0.31 and if σ=0.05 then P(Y2,k∈J)≈0.06. Thus, the Plateau proposals overlap less than the Gaussian proposals and further, the extent of the overlapping of the proposals is controlled by the values of σ,σ1,andσ2. Note that each Plateau proposal distribution has a support on R. This is to ensure that the support of the target distribution is included within the support of the proposals. In theory, this allows the Markov chain to explore the entire support of the target distribution. In a practical setting, however, by selecting the value of σappropriately, the tails of the distribution decay 123 Plateau proposal distributions... 357 200 400 600 800 1000 AP AG2 First Iteration within 95% ellipse Fig. 7 Correlated 2-dimensional Gaussian Target, 2: Violin plot of first hitting time of the Markov chain into 95% ellipse. AP results represented in red and AG2 results in green (colour figure online) some iterations, the AP method always (5,000 out of 5,000 iterations) entered the highdensity region in less than 381 iterations. In comparison, in 517 out of 5,000 iterations, the AG2 took more than 381 iterations to enter the ellipse. In summary, the AP proposals adapt well to the scale of the target components during the initial steps of the algorithm. This was demonstrated in Sect. 5.1 using independent and correlated target components. Moreover, in Sect. 5.2 we showed that the AP method is able to effectively move into high-probability regions when initiated from low probability regions. 6 Results In this section, we compare the long-term performance of the proposed adaptive componentwise, multiple-try MCMC method with other MCMC methods in simulations. The adaptive Plateau MCMC (AP) is compared to a Metropolis-Hastings algorithm with Gaussian proposals (MH) and two versions of an adaptive Gaussian MCMC (AG1 and AG2) as introduced in [27]. The difference between AG1 and AG2 is that AG1 uses (2) with α=2.5 (i.e., as in AP), while AG2 uses α=2.9 as is suggested in [27]. The proposal distributions in both AG1 and AG2 are adapted in the fashion outlined in Sect. 5. For all simulations and methods with multiple trials, we fix the number of trials M=5. Investigation of the method’s performance for differing values of Mis beyond the scope of this paper, which we leave for future work. However, interesting discussions in that direction already exist, see [19] for example. The proposal standard deviation in the AG1 and AG2 methods are initialised at 2j−2for j=1,...,M. The Plateau parameters values are initialised at δ=δ1=1, σ=0.05 and ς=3 with η1=η2=0.4 for the AP method. The proposal distributions used in the MH method depend on the particular target distribution and the choices used are summarised in Table 2. The various target distributions are now introduced. 123 358 F. D-H. Lau, S. Krumscheid 6.1 Target distributions In order to compare the aforementioned methods, we investigate their performances by applying them to sample from a variety of target distributions. Mixture of Gaussians Consider a mixture of two 4-dimensional Gaussians 1 2N(μ1, 1)+1 2N(μ2, 2), where μ1=(5,5,0,0)T,μ2=(15,15,0,0)T and 1=diag(6.25,6.25,6.25,0.01), 2=diag(6.25,6.25,0.25,0.01). We refer to this target distribution as π1. Banana distribution Consider the 8-dimensional “banana-shaped” distribution [10], which is defined as follows. Let fbe the density of the 8D Normal distribution N(0, 3)with covariance given by 3=diag(100,1,...,1). The density function of the banana distribution with non-linearity parameter b>0isgivenby fb=f◦φbwhere the function φbis φb(x)=(x1,x2+bx2 1−100b,x3,...,x8)for x∈R8. The value of bdetermines the amount of non-linearity of φb. Here, we consider the target distribution π2=f0.03, which leads to the unusual banana-shape of the first two components as shown in Fig. 8a. Distributions perturbed by oscillations Another target we consider is the perturbed 2-dimensional Gaussian, whose probability density function is given by π3(x)∝exp −xTAx−cos x1 0.1−0.5cosx2 0.1=: π3(x)for x∈R2 where A=11 13/2. Figure 8c displays the un-normalised function π3(x). Lastly, we also consider the following perturbed version of the 1D bi-stable distribution x→ Z−1e−x2+5x2, whose PDF is given by π4(x)∝exp −x4+5x2−cos x 0.02 for x∈R. Figure 8d displays the PDF of π4(x)where the normalising constant is approximated by numerical integration. 123 Plateau proposal distributions... 359 Bayesian posterior distribution A more complex target is given by the Bayesian posterior distribution originating from a localisation problem using multiple noisy sensor measurements in a wireless sensor network [18]. The problem involves using multiple observations from NSsensor locations to estimate a target location z∈R2. Following [18], the true target location is z∗=(2.5,2.5)Tand the NS=6 sensor locations are h1=(3,−8)T,h2=(8,10)T,h3=(−4,−6)T,h4= (−8,1)T,h5=(10,0)T,andh6=(0,10)T.TheNO=10 observations at each sensor location are modelled as Yk,l∼N20 log (z−hl2),ζ2 l,k=1,...,NO;l=1,...,NS.(4) The observations Yk,lare assumed to be all independent of each other, and are drawn using the true location z=z∗and standard deviations ζ∗=(1,2,1,0.5,3,0.2)T∈R6. The objective is to estimate the target location zand the unknown sensor noise standard deviations ζ=(ζ1,...,ζNS)from the given observations Y=(Yk,l)k=1,...NO l=1,...NS .Thisis achieved using a Bayesian approach to compute the posterior distribution of zand ζgiven the observations Y. The prior for the target location, z∈R2, is uniform on Rz=(−30,30)2 and similarly the prior for the standard deviations ζis uniform on Rζ=(0,20)NS.The MCMC target distribution is the posterior distribution of x:= (z,ζ)T∈R2+NSgiven the observations Y, π5(x):= π(z,ζ|Y=y)∝ NO  k=1 NS  l=1 1 2πζ2 l exp −1 2ζ2 lyk,l−20 log (z−hl2)2 ×I(z∈Rz)I(ζ∈Rζ), where I(·)denotes the indicator function. (a) (c) (b) (d) Fig. 8 Selected marginal density plots of simulation target distributions 123 360 F. D-H. Lau, S. Krumscheid 6.2 Run parameters Each simulation run of the MCMC methods was independently repeated R=200 times and for each run a burn-in period of 50% of the MCMC iterations was used. During the burn-in period, the AP, AG1 and AG2 were allowed to adapt their proposals. For each repetition, all methods started at the same random initial position x0. The number of MCMC iterations, N, used for each method is presented in Table 1, which was determined by a trial run. In order to make fair comparisons, the number of MCMC iterations performed by Metropolis-Hastings algorithm is d×Mtimes larger than the multiple-try versions. This is because multiple-try methods cycle over all dcomponents and evaluate the target for Mtrials each iteration. This will ensure that the number of times the target distribution is evaluated by each MCMC method is the same and that the computational effort is approximately the same. The proposal distribution in the MH method for each target is presented in Table 2.Note that these proposals are based on the target distribution, which would be typically be unknown in practice. The MH method should have an unfair advantage as its proposal is tuned to the target distribution. The particular scaling of 2.4/√dfollows from [6]. The exception is the proposal distribution used for π5, where we use the same proposal as used in [18]. 6.3 Simulation results For each target distribution, we compare the performance of the MCMC methods using measures based on the autocorrelations and jumping distances of the chain. We now define these measures mathematically. Denote the Markov chain produced by one of the MCMC methods for the rth independent repetition as X(r) 0,...,X(r) Nwhere X(r) i=(X(r) i,1,...,X(r) i,d)Tfor r=1,...,R. Denote the component-wise variances of the target as σ2 1,...,σ2 d. Table 1 Number of MCMC iterations used in simulations for each target distribution Target Dimensions (d) Adaptive MH π14 4,000 80,000 π28 10,000 400,000 π32 3,000 30,000 π41 3,000 15,000 π58 10,000 400,000 Table 2 Proposal distributions used in the Metropolis-Hastings algorithm Target Proposal Distribution π1 2.4 √4[0.5N(0, 1)+0.5N(0, 2)] π2 2.4 √8N(0, 3) π3 2.4 √2N(0,A−1) π42.4N(0,1) π5N(0,I8) 123 Plateau proposal distributions... 361 We will use the integrated autocorrelation time (ACT) of the MCMC methods as a measure of performance. The ACT for the chain’s kth component is given by ACTk=1+2 σ2 k N  i=1 cov(X0,k,Xi,k), provided that the chain is stationary so that X0∼π. For every repetition r=1,...,Rof the MCMC method, the integrated autocorrelation times are estimated based on the observed Markov chain X(r) 0,...,X(r) Ncomponent-wise using the initial sequence estimator introduced in [7]. With slight abuse of notation, we will denote the resulting chain-based autocorrelation times by ACT(r) k. Smaller autocorrelation times indicate that consecutive samples have lower correlation. Autocorrelation times are inversely proportional to the effective sample size [16, 17], which is another commonly used measure of performance. In fact, the effective sample size is often interpreted as the number of samples that would need to be (directly) drawn from the target in order to achieve the same variance as that from an estimator of interest using independent samples. Higher effective sample sizes and therefore lower autocorrelation times are desirable. Another way of interpreting the ACTs is through the accuracy of a chain-based Monte Carlo integration. Moreover, the mean squared error of a Monte Carlo estimator can be expressed as a sum of the component-wise ACTs weighted by the component-wise variance. Consequently, a method that produces a chain with low ACTs is preferable. Further, the ACTs of an MCMC characterise the asymptotic variance (so-called time-averaged variance constant in this context) of a Monte Carlo estimator in the central limit theorem for Markov chains [1]. In fact, lower ACTs will lead to smaller time-averaged variance constants. In practice, the target distribution is intractable and therefore the variances of the components are unknown. However, these variances are estimated by the initial sequence estimator method. Another measure of performance is the chain’s average squared jump distance (ASJD), which, for the kth component and repetition r,wedefineas ASJD(r) k=1 N N  i=1|X(r) i,k−X(r) i−1,k|2. The average squared jump distance measures the movement of the chain and also is linked with the acceptance rate of the MCMC method. Higher values of average squared jump distances are desired as it indicates larger moves and therefore more exploration of the space. We use the ASJD as a measure of the ability of an MCMC method to move around the state space. In summary, in the following results, we present the ACTs and the ASJD per component. The distribution of the ACTs and ASJDs over the repetitions, i.e., {ACT(r) k:r=1,...,R} and {ASJD(r) k:r=1,...,R}, will be presented using violin plots. 6.3.1 Mixture of Gaussians The results for the 4-dimensional mixture of two Gaussians target, π1, are presented in Fig. 9. Fig. 9a indicates that the AP method achieves lower ACTs than the other methods for all components (including the MH method, which is not included in this figure due to very high ACTs). Further, the range of ACT values suggest that the AP method consistently produces MCMC chains with lower ACTs. 123 362 F. D-H. Lau, S. Krumscheid Component 1 Component 2 0 30 70 110 160 ACT Component 3 Component 4 02468 (a) Autocorrelation times Component 1 Component 2 Component 3 0246811 14 ASJD 0.001 0.006 0.011 0.016 Component 4 (b) Average square jump distance Fig. 9 Mixture of Gaussians, π1: Distribution of the ACTs and ASJDs of MCMC methods for target π1.Red is AP, blue is AG1, green is AG2 and pink is MH (colour figure online) In terms of the movement of the MCMC chains, the AP method outperforms the other methods for this target distribution. In fact, the ASJDs presented in Fig. 9b show that the AP method moving around the state space in larger jumps than the other methods. Since the AP and AG1 methods use the same weight function as discussed in Sect. 2.1 this advantage is due to using the Plateau proposals in contrast to Gaussian proposals. These results suggest that the AP method is able to move between the two Gaussians in the target density efficiently. 6.3.2 Banana distribution The target, π2, is a difficult distribution to sample from due to the wide-ranging variances across the components and the unusual banana-shape of the first two components. The ACTs for the methods, presented in Fig. 10a, show similar results across the AP, AG1 and AG2 for the first two components. However, for the remaining components, the AP is achieving notably smaller ACTs. The ACT results for the MH are substantially larger for all components, and thus are not included in this figure. As an indication, the median ACTs for components 1 to 8 respectively are: 1131.74, 2066.35, 54.24, 54.37, 54.34, 54.03, 54.74, 54.47 for the MH method to 2 decimal places. The ASJDs for the methods presented in Fig. 10b. The AP method again outperforms the other methods by achieving higher ASJDs for all components. Note that for the first component, the wide range of jumping distance produced when using the Plateau proposals. This suggests that the AP method is able to navigate the banana-shape in the first component easily. 6.3.3 2D perturbed distribution For the perturbed 2-dimensional distribution, π3, the perturbations represent local modes where MCMC methods may potentially get stuck. Again, the AP method’s ability to move 123 Plateau proposal distributions... 363 Component 1 Component 2 0100 250 400 550 ACT Component 3 Component 5 Component 7 00.511.522.53 (a) Autocorrelation times Component 1 0816 26 36 46 56 ASJD Component 2 Component 4 Component 6 Component 8 00.5 11.522.5 (b) Average square jump distance Fig. 10 Banana Distribution, π2: Distribution of the ACTs and ASJDs of MCMC methods for target π2.Red is AP, blue is AG1, green is AG2 and pink is MH (colour figure online) slightly larger distances, as shown in Fig. 11b, gives it a slight advantage over the other methods. This ability to jump further may explain the lower ACTs for the AP method as depicted in Fig. 11a. 6.3.4 1D perturbed distribution Similar to π3, the oscillations in π4are areas where an MCMC may get stuck. The ACTs for the AP, AG1 and AG2 method are presented in Fig. 12a. The AP achieves the lowest ACTs; however, there are a few outliers which may indicate a few runs where the sampler got stuck in the local modes. This may also be the case for the AG1 method. For the MH method, the ACTs (not presented in the figure) are extremely large in comparison to the other methods – with a median of 178.54 and a range of 111.62 to 457.56 to 2 decimal places. The ASJDs for the AP method, is on average jumping also twice the distance of the AG1 and AG2 methods – see Fig. 12b. Again, there are some outlying ASJDs for the AP method, which may indicate some repetitions where the MCMC got stuck in local modes. 6.3.5 Posterior distribution The posterior distribution, π5, poses a difficult distribution to sample from due to the nonlinearity of the parameters. The ACTs for the MCMC methods (except MH) are presented in Fig. 13a. The AP method consistently achieves a lower ACT than the AG1 and AG2 method, especially for the variance parameters (components 3 to 8). In terms of ASJDs, all methods give similar performance – see Fig. 13b, with a slight advantage to using the AP method which achieves higher ASJD albeit infrequently. The ACT and ASJD results for the MH method are not presented in Fig. 13a. The median ACT for components 1 to 8 range from 19730 to 37430. Further, the median ASJDs range from 2.9×10−8to 4.4×10−6over 123 364 F. D-H. Lau, S. Krumscheid Component 1 Component 2 0 5 10 15 20 ACT (a) Autocorrelation times Component 1 Component 2 0 0.15 0.3 0.45 0.6 0.75 0.9 ASJD (b) Average square jump distance Fig. 11 2D Perturbed Distribution, π3: Distribution of the ACTs and ASJDs of MCMC methods for target π3. Red is AP, blue is AG1, green is AG2 and pink is MH (colour figure online) Component 1 0 5 10 15 20 25 ACT (a) Autocorrelation times Component 1 00.3 0.6 0.9 1.2 1.5 1.8 2.1 2.4 ASJD (b) Average square jump distance Fig. 12 1D Perturbed Distribution, π4: Distribution of the ACTs and ASJDs of MCMC methods for target π4. Red is AP, blue is AG1, green is AG2 and pink is MH (colour figure online) 123 Plateau proposal distributions... 365 Component 1 Component 2 Component 3 Component 4 Component 5 Component 6 Component 7 Component 8 0123456 ACT (a) Autocorrelation times Component 1 Component 2 Component 3 Component 4 Component 5 Component 6 Component 7 Component 8 0 0.5 1 1.5 2 2.5 3 3.5 ASJD (b) Average square jump distance Fig. 13 Posterior Distribution, π5: Distribution of the ACTs and ASJDs of MCMC methods for target π5. Red is AP, blue is AG1, green is AG2 all 8 components. Both the high ACTs and very low ASJDs indicate that the MH method performed poorly. 6.4 Comparing computational speed In this section, we compare the computational effort of the AP and AG methods. For this comparison, we use the 8-dimensional banana target distribution, π2and investigate the time to run the MCMC methods with a different number of trials Mand MCMC steps N. Specially, we run a single simulation with a unique combination of M∈{1,4,5,10}and N∈{5000,10000,50000,100000}. Each simulation is repeated 1,000 times. We compare the AP and AG methods with α=2.9 both without and with the adaptation. The methods are written in C++. The results presented in Tables 3and 4reports the average effort to run each MCMC methods over 1,000 iterations. The reported effort is relative to the effort to generate the MCMC samples from the AG method. Table 3shows that the AP and AG methods, with and without adapting, are all comparable in terms of effort for all values of N.Table4suggests that keeping Nfixed and increasing the number of trial Mdoes not lead to any one method being more efficient than the others. We Table 3 Computational effort of MCMC algorithms using M=5 trials divided by the effort using the AG method N=5,000 10,000 50,000 100,000 AP 1.0092 1.0208 1.0118 1.0064 AP (no adapt) 1.0079 1.0168 1.0069 1.0051 AG 1.0000 1.0000 1.0000 1.0000 AG (no adapt) 0.9975 1.0043 1.0038 1.0048 123 366 F. D-H. Lau, S. Krumscheid Table 4 Computational effort of MCMC algorithms using N=1000 MCMC steps divided by the effort using the AG method M=13 5 10 AP 1.1197 1.0320 1.0208 0.9939 AP (no adapt) 1.1228 1.0343 1.0168 0.9909 AG 1.0000 1.0000 1.0000 1.0000 AG (no adapt) 1.0421 0.9990 1.0043 1.0005 conclude that both the AP and AG methods are comparable in terms of computation time, as there are negligible differences between the two methods. 7 Generalisation of plateau proposals In this section we explain why the family of Plateau proposals are suitable within multiple-try MCMC methods. As mentioned in Sect. 3, the Plateau proposals are designed to be nonoverlapping and simultaneously have no (or small) gaps between the individual proposals. The non-overlapping feature ensures that trial samples are drawn from separate areas meaning that the chain explores the area efficiently. The feature of having no gaps ensures that there are no regions between proposals that will not be explored – meaning that the entire combined support of the proposals is being explored with a non-negligible probability. We now introduce measures of overlap and gap between two trial distributions. In this section, we simplify the setting to two contiguous trial distributions rather than using M distributions, to make clear the overlap and gap measures. The intuition carries over to the full case with Mtrial distributions. The overlap and gap between distributions is illustrated in (i) the Gaussian case and (ii) the Plateau case which are now defined. Let Aand Btwo independent random variables. In the Gaussian case let A∼N(0,s2)and B∼N(, s2) for some >0. In the Plateau case let A∼f(·,μ =0,δ =/2,σ 1=σ, σ2=σ) and B∼f(·,μ=,δ =/2,σ 1=σ, σ2=σ)where the PDF fis defined by Eq. (3). Recall that in the PDF f, the parameter μis the location of the centre of a single plateau and σ1(σ2) determines the rate of exponential decay to the left (right) side of the plateau. As suggested earlier, and used in all numerical simulation conducted in this paper, we set σ1=σ2=σ. The Gaussian and Plateau cases are illustrated in Fig.14. We measure the overlap between the distributions of Aand Bas poverlap := P(B<mp)+P(A>mp), 0mpΔ ( a ) 0mpΔ ( b ) Fig. 14 PDF illustrations for random variables A(black) and B(red) for athe Gaussian case and bthe Plateau case. Shaded grey area represents the overlap measure, poverlap (colour figure online) 123