# Hierarchical Models and Tuning of Random Walk Metropolis Algorithms

Authors: Mylène Bédard
Journal: Journal of Probability and Statistics
 (2019)

---

## Abstract

We obtain weak convergence and optimal scaling results for the random walk Metropolis algorithm with a Gaussian proposal distribution. The sampler is applied to hierarchical target distributions, which form the building block of many Bayesian analyses. The global asymptotically optimal proposal variance derived may be computed as a function of the specific target distribution considered. We also introduce the concept of locally optimal tunings, i.e., tunings that depend on the current position of the Markov chain. The theorems are proved by studying the generator of the first and second components of the algorithm and verifying their convergence to the generator of a modified RWM algorithm and a diffusion process, respectively. The rate at which the algorithm explores its state space is optimized by studying the speed measure of the limiting diffusion process. We illustrate the theory with two examples. Applications of these results on simulated and real data are also presented.


---

## Body

## 1. Introduction

Random walk Metropolis (RWM) algorithms are widely used to sample from complex or multidimensional probability distributions [1, 2]. The simplicity and versatility of these samplers often make them the default option in the MCMC toolbox. Implementing a RWM algorithm involves a tuning step, to ensure that the process explores its state space as fast as possible and that the sample produced be representative of the probability distribution of interest (the target distribution). In this paper, we solve an aspect of the tuning problem for a large class of target distributions with correlated components. This issue has mainly been studied for product target densities, but attention has recently turned towards more complex target models [3, 4]. The specific type of target distribution considered here is formed of components which are related according to a hierarchical structure. These distributions are ubiquitous in several fields of research (finance, biostatistics, and physics, to name a few) and constitute the basis of many Bayesian inferences.Bayesian hierarchical models comprise a likelihood functionfdθ, which is the statistical model for the observed data d. The parameters θ are then modeled using a prior distribution πθρ; since this prior might not be easy to determine, it is common practice to assume that the hyperparameters ρ are themselves distributed according to a noninformative prior distribution πρ. The various models thus represent different levels of hierarchy and give rise to a posterior distribution πθ,ρd, which is often quite complex. Most of the time, this distribution cannot be studied analytically or sampled directly, and thus simulation algorithms such as MCMC methods are required to perform a statistical analysis. Samplers such as the RWM, RWM-within-Gibbs, and Adaptive Metropolis (see [5]) are usually the default algorithms for such targets.The idea behind RWM algorithms is to build a Markov chain having the Bayesian posterior (target) distribution as its stationary distribution. To implement this method, users must select a proposal distribution from which are generated candidates for the Markov chain. This distribution should ideally be similar to the target, while remaining accessible from a sampling viewpoint. A pragmatic choice is to let the proposed moves be normally distributed around the latest value of the sample. Tuning the variance of the normal proposal distribution (σ2) has a significant impact on the speed at which the sampler explores its state space (hereafter referred to as “efficiency”), with extremal variances leading to slow-mixing algorithms. In particular, large variances seldom induce suitable candidates and result in lazy processes; small variances yield hyperactive processes whose tiny steps lead to a time-consuming exploration of the state space. Seeking for an intermediate value that optimizes the efficiency of the RWM algorithm, i.e., a proposal variance σ2 offering sizable steps that are still accepted a reasonable proportion of the time, is called the optimal scaling problem.The optimal scaling issue of the RWM algorithm with a Gaussian proposal has been addressed by many researchers over the last few decades. It has been determined in [6] that target densities formed of independent and identically distributed (i.i.d.) components correspond to an optimal proposal variance σ^2n≈5.66/nElogfX′, where f is the density of one target component and n the number of target components. This optimal proposal variance has also been shown to correspond to an optimal expected acceptance rate of 23.4%, where the acceptance rate is defined as the proportion of candidates that are accepted by the algorithm. Generalizing this conclusion is an intricate task, and further research on the subject has mainly been restricted to the case of target distributions formed of independent components (see [7–12]). In the specific case of multivariate normal target distributions however, the optimal variance and acceptance rate may be easily determined (see [11, 13]). Lately, [3, 4] have also performed scaling analyses of nonproduct target densities. These advances are important, as MCMC methods are mainly used when dealing with complex models, which only rarely satisfy the independence assumption among target components. These results however assume that the correlation structure among target components is known and used in generating candidates for the chain. This is a restrictive assumption that leads, as expected, to an optimal acceptance rate of 23.4% (see [12] for an explanation).In this paper, we focus on solving the optimal scaling problem for a wide class of models that include a dependence relationship, the hierarchical distributions. Weak convergence results are derived without explicitly characterizing the dependency among target components, and thus rely on a Gaussian proposal distribution with diagonal covariance matrix. The optimal proposal variance may then be obtained from these results, i.e., by maximizing the speed measure of the limiting diffusion process. This constitutes significant advances in understanding the theoretical underpinnings of the RWM sampler. More importantly in practice, the results theoretically support the use of RWM-within-Gibbs over RWM samplers and provide a convenient approach for obtaining a new type of proposal variances. These proposal variances are a function of the current state of the Markov chain; they thus evolve with the chain and lead to more appropriate candidates in the RWM-within-Gibbs algorithm.In the next section, we describe the target distribution and introduce some notation related to the RWM sampler. The theoretical optimal scaling results are stated in Section3 and then illustrated with two examples using RWM samplers in Section 4. In Section 5, the potential of RWM-within-Gibbs with local scalings is illustrated in Bayesian contexts through a simulation study and an application on real data. Extensions are briefly discussed in Section 6, while appendices contain proofs.


## 2. Framework

Consider ann-dimensional target distribution consisting of a mixing component X1 and n−1 conditionally i.i.d. components Xi (i=2,…,n) given X1. Suppose that this distribution has a target density π with respect to Lebesgue measure, where(1)πx=f1x1∏i=2nfxix1.To obtain a sample from the target density in (1), we rely on a RWM algorithm with a Gaussian proposal distribution. This sampler builds an n-dimensional Markov chain Xnj;j∈ℕ having πx as its stationary density. Given Xnj=x, the time-j state of the Markov chain, one iteration is performed according to the following steps:(1)
Generate a candidateYnj+1=y from a Nx,Dn, where Dn is a diagonal variance matrix with elements σ12n,σ2n,…,σ2n. In particular, set Dn=ℓ2In/n, where ℓ>0 is a tuning parameter and In the n-dimensional identity matrix(2)
Compute the acceptance probabilityαx,y=min1,πy/πx(3)
GenerateUj+1∼U0,1(4)
IfUj+1≤αx,y, accept the candidate and set Xnj+1=y; otherwise, the Markov chain remains at the same state for another time interval and Xnj+1=xOptimal scaling results widely rely on the use of Gaussian proposal distributions which, due to their symmetry, lead to a simplified form of the acceptance probability. Although generally not emphasized in the literature, we note that the proposal variance could also be a function ofx, which would result in a nonhomogeneous random walk sampler. In that case, there would be no simplification in the Metropolis–Hastings acceptance probability and Step 2 would then replaced by(2)αx,y=min1,πyqnx;yπxqny;x,where qny;x is the density of a Nx,Dnx.In what follows, we work towards finding the optimal value ofℓ, i.e., leading to an optimally mixing chain. The proofs of the theoretical results rely on CLTs and LLNs; as such, the results are obtained by letting n⟶∞. This is a common approach in MCMC theory and does not prevent users from applying the asymptotically optimal value of ℓ in lower dimensional contexts (as small as n=10 or 15). Indeed, a particularity of optimal scaling results is that the asymptotic behaviour kicks in extremely rapidly, as shall be witnessed in the examples of Section 4.The first thought of most MCMC users when facing a target density as in (1) would be to use a RWM-within-Gibbs algorithm, which consecutively updates subgroups of the n components in a given iteration. The tuning of RWM-within-Gibbs algorithms has been addressed in [11], but only for target distributions with i.i.d. components and Gaussian targets with correlation. Focusing on RWM algorithms is thus a good starting point to understand the behaviour of samplers applied to hierarchical target distributions. The results expounded in this paper lead to the concept of local tunings, which is particularly appealing in the context of RWM-within-Gibbs. Incidentally, the proofs in appendices provide a theoretical justification for the use of locally optimal scalings in RWM-within-Gibbs (see [14]). These findings are illustrated in the examples of Section 5.In Sections2.1, 2.2, and 3, we expound how to obtain asymptotically optimal variances Dn and Dnx for RWM and RWM-within-Gibbs, respectively. Section 2.1 describes the regularity conditions imposed on πx, while Section 2.2 explains why the proposal matrix Dn=ℓ2In/n is the optimal choice for obtaining the theoretical results that shall be presented in Section 3.


### 2.1. Assumptions on the Target Density

To characterize the asymptotic behaviour of the conditionally i.i.d. componentsXi (i=2,…,n), we impose some regularity conditions on the densities f1 and f in (1). The density f1 is assumed to be a continuous function on, with X1=x1:f1x1>0 forming an open interval.For all fixedx1∈X1, fxx1 is a positive C2 density on and ∂/∂xlogfxx1 is Lipschitz continuous with constant Kx1 such that EK2X1<∞. Here, C2 denotes the space of real-valued functions with continuous second derivative. For all fixed x∈X=, fxx1 is a C2 function on X1 and ∂/∂xlogfxx1 is Lipschitz continuous with constant Lx such that EL4X<∞. Furthermore,(3)EX∂∂XlogfXx14<∞∀x1∈X1withE∂∂XlogfXX14<∞,and, hereafter, the notation EX⋅ means that the expectation is computed with respect to X conditionally on the other variables in the expression; the first expectation in (3) is thus obtained according to the conditional distribution of X given X1. Where there is no confusion possible, E⋅ shall be used to denote an expectation with respect to all random variables in the expression. The above regularity conditions constitute an extension of those stated in [8] for target distributions with independent components and are weaker than would be a Lipschitz continuity assumption on the bivariate function ∂/∂xlogfxx1. They also imply that the Lipschitz constants Kx1 and Lx themselves satisfy a Lipschitz condition.We now impose further conditions onfxx1 to account for the movements of the coordinate X1 when studying the asymptotic behaviour of a component Xi (i=2,…,n). These movements should not be too abrupt, so for almost all fixed x∈X, ∂/∂x1logfxx1 is Lipschitz continuous with constant Mx such that EM2X<∞ and(4)EX∂∂x1logfXx12<∞∀x1∈X1withE∂∂X1logfXX12<∞.Finally, in order to characterize the asymptotic behaviour of the mixing componentX1, we introduce assumptions that are closely related to the Bernstein von Mises Theorem. Let X2:n=X2,…,Xn, X¯=X2,X3,…, and ⟶p denote convergence in probability. Assume that VX1X2:n⟶p0, and denote μ≡μX¯ such that μn≡μnX2:n=EX1X2:n⟶pμ as n⟶∞, with μ<∞. Hereafter, we make a small abuse of notation by letting μ and μn sometimes denote the random variable or the realisation, depending on the context. Furthermore, define X˜1=nX1−μn; for almost all x2:n∈ℝn−1, the conditional density of X˜1 given X2:n, f1μn+x˜1/nx2:n/n, is assumed to converge almost surely to g1x˜1X¯, a continuous density on ℝ with respect to Lebesgue measure. In fact, the information on X1 increases linearly in n, meaning that the limiting density of X1x2:n is degenerate, but that a standard rescaling leads to a nontrivial density on ℝ (normal distribution).


### 2.2. Form of the Proposal Variance MatrixDn

In Section 3, we focus on deriving weak convergence and optimal scaling results for the RWM algorithm with a Gaussian proposal by lettingn, the dimension of the target density in (1), approach ∞. Traditionally, asymptotically optimal scaling results have been obtained by studying the limiting path of a given component (X2 say) as n⟶∞. In the case of target distributions with i.i.d. components (and some extensions), the components of the RWM algorithm are asymptotically independent of each other and their limiting behaviour is regimented by identical one-dimensional Markovian processes. In the current correlated framework, we expect the presence of an asymptotic dependence relationship among Xi (i∈2,…,n) and X1, in the spirit of (1). In the following section, we thus study the limiting behaviour of components X1 and X2 separately, on their respective conditional space. This approach allows us to quantify the mixing rate of each component Xi conditionally on the others and to propose optimal scalings for the sampler.To obtain nontrivial limiting processes describing the behaviour of the RWM sampler asn⟶∞, we need to fix the form of the proposal scalings σ12n,σ2n. Whilst the proposals are independent, a single accept-reject step is used, which makes the paths of the components dependent. We aim to choose the maximal scalings that avoid a degenerate limit (of either 0 or 1) for this acceptance probability. Since the distribution of X1 conditional on X2:n contracts at a rate of n, if σ1n/n⟶∞, the proposed jumps in X1 will be too large. If σ1n/n⟶0, then the change in X1 makes no contribution to the acceptance probability in the limit; to maximise movements, we, therefore, require σ1n∝1/n. Now, the conditional distribution of X2:n given X1 does not contract with n. Nonetheless, when proposing jumps in X2:n using σ2n=σ2, the odds of rejecting an n-dimensional candidate increase with n and lead to a degenerate (null) acceptance probability. To overcome this problem, we then let the proposal variance be a decreasing function of the dimension. In fact, since Lipschitz conditions control the contribution to the accept-reject ratio coming from the movements of X1, a similar argument to that which leads to σn∝1/n in the case of i.i.d. targets applies again here. We therefore set Dn=ℓ2In/n, where ℓ>0 is a tuning parameter and In the n-dimensional identity matrix.Asn⟶∞, it becomes necessary to speed up time to compensate for the reduced movement along components X2:n. The time interval between each proposed candidate is thus set to 1/n, and we study the continuous-time, sped up version of the initial Markov chain defined as Wnt;t≥0=Xnnt;t≥0, where ⋅ is the floor function. Similarly to the i.i.d. case, a limiting diffusion is obtained for the rescaled one-dimensional process related to Xi (i≥2), but this time its behaviour is conditional on X1.Since the first coordinateX1 converges to a point μ, a transformation X˜1=nX1−μn is required to obtain the limiting behaviour of this component. We thus study the continuous-time process W˜nt;t≥0=X˜1nt,X2:nnt;t≥0; in other words, we are now looking at a magnified, centered version of the path associated to X1. This transformation leads to proposal distributions Y˜1=nY1−μn∼Nx˜1,ℓ2 and Yi∼Nxi,ℓ2/n, i=2,…,n with ℓ>0; it thus cancels the effect of n in σ12n. Without the speed up of time, the limiting process for X˜1 is then a propose-accept-reject on the conditional density for X˜1, given the current values of X2:n; this is made precise in Theorem 1. When considering the diffusion limit for Xi i≥2 with time sped-up, this effectively means that at every instant, X1 is simply a sample from its conditional distribution given the current values of X2:n; this is made precise in Theorem 2.We note that an alternative scaling ofσ1n∝1/n could also be applied. The sped-up limiting process would then be a diffusion for all coordinates and would be easier to describe. However, this would also be a deliberate handicapping of the algorithm since the change in X1 would make no contribution to the acceptance probability in the limit. A suboptimal σ12n, besides altering the movements of X1, would thus also indirectly affect the efficiency according to which X2:n explores its state space.


## 2.1. Assumptions on the Target Density

To characterize the asymptotic behaviour of the conditionally i.i.d. componentsXi (i=2,…,n), we impose some regularity conditions on the densities f1 and f in (1). The density f1 is assumed to be a continuous function on, with X1=x1:f1x1>0 forming an open interval.For all fixedx1∈X1, fxx1 is a positive C2 density on and ∂/∂xlogfxx1 is Lipschitz continuous with constant Kx1 such that EK2X1<∞. Here, C2 denotes the space of real-valued functions with continuous second derivative. For all fixed x∈X=, fxx1 is a C2 function on X1 and ∂/∂xlogfxx1 is Lipschitz continuous with constant Lx such that EL4X<∞. Furthermore,(3)EX∂∂XlogfXx14<∞∀x1∈X1withE∂∂XlogfXX14<∞,and, hereafter, the notation EX⋅ means that the expectation is computed with respect to X conditionally on the other variables in the expression; the first expectation in (3) is thus obtained according to the conditional distribution of X given X1. Where there is no confusion possible, E⋅ shall be used to denote an expectation with respect to all random variables in the expression. The above regularity conditions constitute an extension of those stated in [8] for target distributions with independent components and are weaker than would be a Lipschitz continuity assumption on the bivariate function ∂/∂xlogfxx1. They also imply that the Lipschitz constants Kx1 and Lx themselves satisfy a Lipschitz condition.We now impose further conditions onfxx1 to account for the movements of the coordinate X1 when studying the asymptotic behaviour of a component Xi (i=2,…,n). These movements should not be too abrupt, so for almost all fixed x∈X, ∂/∂x1logfxx1 is Lipschitz continuous with constant Mx such that EM2X<∞ and(4)EX∂∂x1logfXx12<∞∀x1∈X1withE∂∂X1logfXX12<∞.Finally, in order to characterize the asymptotic behaviour of the mixing componentX1, we introduce assumptions that are closely related to the Bernstein von Mises Theorem. Let X2:n=X2,…,Xn, X¯=X2,X3,…, and ⟶p denote convergence in probability. Assume that VX1X2:n⟶p0, and denote μ≡μX¯ such that μn≡μnX2:n=EX1X2:n⟶pμ as n⟶∞, with μ<∞. Hereafter, we make a small abuse of notation by letting μ and μn sometimes denote the random variable or the realisation, depending on the context. Furthermore, define X˜1=nX1−μn; for almost all x2:n∈ℝn−1, the conditional density of X˜1 given X2:n, f1μn+x˜1/nx2:n/n, is assumed to converge almost surely to g1x˜1X¯, a continuous density on ℝ with respect to Lebesgue measure. In fact, the information on X1 increases linearly in n, meaning that the limiting density of X1x2:n is degenerate, but that a standard rescaling leads to a nontrivial density on ℝ (normal distribution).


## 2.2. Form of the Proposal Variance MatrixDn

In Section 3, we focus on deriving weak convergence and optimal scaling results for the RWM algorithm with a Gaussian proposal by lettingn, the dimension of the target density in (1), approach ∞. Traditionally, asymptotically optimal scaling results have been obtained by studying the limiting path of a given component (X2 say) as n⟶∞. In the case of target distributions with i.i.d. components (and some extensions), the components of the RWM algorithm are asymptotically independent of each other and their limiting behaviour is regimented by identical one-dimensional Markovian processes. In the current correlated framework, we expect the presence of an asymptotic dependence relationship among Xi (i∈2,…,n) and X1, in the spirit of (1). In the following section, we thus study the limiting behaviour of components X1 and X2 separately, on their respective conditional space. This approach allows us to quantify the mixing rate of each component Xi conditionally on the others and to propose optimal scalings for the sampler.To obtain nontrivial limiting processes describing the behaviour of the RWM sampler asn⟶∞, we need to fix the form of the proposal scalings σ12n,σ2n. Whilst the proposals are independent, a single accept-reject step is used, which makes the paths of the components dependent. We aim to choose the maximal scalings that avoid a degenerate limit (of either 0 or 1) for this acceptance probability. Since the distribution of X1 conditional on X2:n contracts at a rate of n, if σ1n/n⟶∞, the proposed jumps in X1 will be too large. If σ1n/n⟶0, then the change in X1 makes no contribution to the acceptance probability in the limit; to maximise movements, we, therefore, require σ1n∝1/n. Now, the conditional distribution of X2:n given X1 does not contract with n. Nonetheless, when proposing jumps in X2:n using σ2n=σ2, the odds of rejecting an n-dimensional candidate increase with n and lead to a degenerate (null) acceptance probability. To overcome this problem, we then let the proposal variance be a decreasing function of the dimension. In fact, since Lipschitz conditions control the contribution to the accept-reject ratio coming from the movements of X1, a similar argument to that which leads to σn∝1/n in the case of i.i.d. targets applies again here. We therefore set Dn=ℓ2In/n, where ℓ>0 is a tuning parameter and In the n-dimensional identity matrix.Asn⟶∞, it becomes necessary to speed up time to compensate for the reduced movement along components X2:n. The time interval between each proposed candidate is thus set to 1/n, and we study the continuous-time, sped up version of the initial Markov chain defined as Wnt;t≥0=Xnnt;t≥0, where ⋅ is the floor function. Similarly to the i.i.d. case, a limiting diffusion is obtained for the rescaled one-dimensional process related to Xi (i≥2), but this time its behaviour is conditional on X1.Since the first coordinateX1 converges to a point μ, a transformation X˜1=nX1−μn is required to obtain the limiting behaviour of this component. We thus study the continuous-time process W˜nt;t≥0=X˜1nt,X2:nnt;t≥0; in other words, we are now looking at a magnified, centered version of the path associated to X1. This transformation leads to proposal distributions Y˜1=nY1−μn∼Nx˜1,ℓ2 and Yi∼Nxi,ℓ2/n, i=2,…,n with ℓ>0; it thus cancels the effect of n in σ12n. Without the speed up of time, the limiting process for X˜1 is then a propose-accept-reject on the conditional density for X˜1, given the current values of X2:n; this is made precise in Theorem 1. When considering the diffusion limit for Xi i≥2 with time sped-up, this effectively means that at every instant, X1 is simply a sample from its conditional distribution given the current values of X2:n; this is made precise in Theorem 2.We note that an alternative scaling ofσ1n∝1/n could also be applied. The sped-up limiting process would then be a diffusion for all coordinates and would be easier to describe. However, this would also be a deliberate handicapping of the algorithm since the change in X1 would make no contribution to the acceptance probability in the limit. A suboptimal σ12n, besides altering the movements of X1, would thus also indirectly affect the efficiency according to which X2:n explores its state space.


## 3. Asymptotics of the RWM Algorithm

In this section, we introduce results about the limiting behaviour (asn⟶∞) of the time- and scale-adjusted univariate processes W˜1nt;t≥0 and Wint;t≥0 (i=2,…,n). From these results, we determine the asymptotically optimal scaling (AOS) values and acceptance rate (AOAR) that optimize the mixing of the algorithm.Hereafter, we let⇒ denote weak convergence in the Skorokhod topology and Bt a Brownian motion at time t; the cumulative distribution function of a standard normal random variable is denoted by Φ⋅.Theorem 1.
Consider a RWM algorithm with proposal distributionNx,ℓ2In/n used to sample from a target density π as in (1). Suppose that π satisfies the conditions on f1 and f specified in Section 2.1 and that Xn0 is distributed according to π in (1).If(5)1n∑i=2n∂∂XilogfXiX1=μn+X˜1n2⟶pγ˜μ,with(6)γ˜μ=EX∂∂XlogfXμX¯2=∫ℝ∂∂xlogfxμX¯2fxμX¯dx<∞,then the magnified process W˜1nt;t≥0⟹W˜1t;t≥0. Here, W10 and Wi0 (i=2,3,…) are distributed according to the densities f1 and f, respectively, which implies that W˜10 is distributed according to the density g1 in Section 2.1. Given the time-t state W˜t=x˜1,x¯, the process W˜1t;t>0 evolves as the continuous-time version of a special RWM algorithm applied to the target density g1x˜1x¯; the proposal distribution of this algorithm is a Nx˜1,ℓ2, and the acceptance rule is defined as(7)α∗x˜1,y˜1x¯=Φlogg1y˜1x¯/g1x˜1x−ℓ2/2γ˜μℓγ˜1/2μ+g1y˜1x¯g1x˜1x¯Φ−logg1y˜1x¯/g1x˜1x¯−ℓ2/2γ˜μℓγ˜1/2μ.Proof.
See Appendix A.1.This result describes the limiting path associated to the coordinateX˜1 as n⟶∞, which is Markovian with respect to the history of the multidimensional chain ℱW˜t. We recall that the conditional distribution of X1 given X2:n contracts at a rate of n and that σ1n∝1/n. Conditionally on X¯, the transformed X˜1 thus mixes according to O1 and explores its conditional state space much more efficiently than the other components, as shall be witnessed in Theorem 2. The asymptotic process found can be described as an atypical one-dimensional RWM algorithm, whose acceptance rule α∗x˜1,y˜1x¯ and target density g1x˜1x¯ both vary according to x¯ at every iteration. The acceptance function α∗ in (7) satisfies the reversibility condition with respect to g1x˜1x¯ (see [8] for more details about this acceptance function).Theorem 1 is interesting from a theoretical perspective but cannot be used to optimize the global mixing of the algorithm. Although we could try to determine the value ofℓ leading to the optimal mixing of X1 on its conditional space, it will be wiser to focus instead on optimizing the mixing rate of X2:n on its own conditional space given X1. Since the distribution of X1 contracts about μn, the position of this coordinate heavily depends on the current state of X2:n. We shall also see in Theorem 2 that given X1, the coordinates Xi (i≥2) explore their conditional state space according to On. Since these coordinates take more time exploring their conditional distribution and heavily affect the position of X1, the global performance of the sampler is subjected to the mixing of X2:n conditionally on X1.Theorem 2.
Consider a RWM algorithm with proposal distributionNx,ℓ2In/n used to sample from a target density π as in (1). Suppose that π satisfies the conditions on f1 and f specified in Section 2.1 and that Xn0 is distributed according to π in (1).Fori=2,…,n, we have Wint;t≥0⟹Wit;t≥0, where Wi0 i≥2 is distributed according to f, and W10 according to f1. Conditionally on W1t, the evolution of Wit;t>0 over an infinitesimal interval dt satisfies(8)dWit=υ1/2ℓ,W1tdBt+12υℓ,W1t∂∂WitlogfWitW1tdt,with(9)υℓ,x1=2ℓ2EZ1Φ−ℓ2γ1/2x1,Z1,(10)Z1=nY1−x1ℓ∼N0,1,(11)γx1,z1=z12EX∂∂x1logfXx12+EX∂∂XlogfXx12.Proof.
See Appendix A.2.Equation (8) describes the behaviour of the process at the next instant, (t+dt), given its position at t. This expression should not come as a surprise: each rescaled component Xi (i=2,…,n) asymptotically behaves according to a diffusion process that is Markovian with respect to ℱW1,Wit. Examination of (8) also tells us that fWitW1t is invariant for this diffusion process (see [15], for instance). We finally recall that σn∝1/n and therefore, conditionally on X1, the rescaled Xi mixes according to On. Each coordinate Xi thus requires more iterations than were required by the coordinate X1 to explore its conditional state space.SinceX1 and Xi i≥2 use different time rescaling factors, the asymptotic behaviour of these coordinates cannot be expressed as a bivariate diffusion process. To obtain such a diffusion, we would have to rely on inhomogeneous proposal variances to ensure that X1 also mixes in On iterations; as mentioned at the end of Section 2, this would require setting σ1n=ℓ/n, σn=ℓ/n for ℓ>0. This framework would, of course, be suboptimal as it would restrain the X1 movements. Proposed jumps for X1 would then become insignificant, and so the first term in (11) would be null.Remark 1.
Studying the limiting behaviour ofX1 and Xi (i=2,…,n) separately does not cause information loss. In fact, studying the paths of X1,X2 simultaneously would require letting the test function h of the generator in (A.3) be a function of X1,X2. Such a generator would however be developed as an expression in which cross-derivative terms (e.g., ∂2/∂x1∂x2hx1,x2) are null, which confirms that given the current state of the asymptotic process, one-dimensional moves are performed independently for each coordinate.
The limiting processes in Theorems 1 and 2 indicate that the componentX1 explores its conditional state space at a different (higher) rate than X2:n explores its own. Combined to the specific Markovian forms of the limiting processes obtained (with respect to ℱW˜t and ℱW1,Wit, respectively), this points towards the need for updating X1 and X2:n separately, assessing the superiority of RWM-within-Gibbs samplers for sampling from hierarchical targets. These algorithms update blocks of components successively, a design that allows fully exploiting the characteristics of the target considered. To our knowledge, this is the first time that asymptotic results are used to theoretically validate the superiority of RWM-within-Gibbs over RWM samplers for hierarchical target distributions. This theoretical superiority is obviously tempered in practice by an increased computational effort; the extent of this computational overhead is however difficult to quantify in full generality. To this end, Section 5 presents two examples that illustrate the performance of the RWM-within-Gibbs and compare it to RWM and Adaptive Metropolis samplers.


### 3.1. Optimal Tuning of the RWM Algorithm

To be confident that then-dimensional chain has entirely explored its state space, we must be certain that every one-dimensional path has explored its own space. In the correlated framework considered, the overall mixing rate of the RWM sampler is only as fast as the slowest component. As explained in Section 3, optimal mixing of the algorithm shall be attained by optimizing the mixing of the coordinates Xi, i=2,…,n. In the limit, the only quantity that depends on the proposal variance (i.e., on ℓ) is υℓ,W1t in (9). To optimize mixing, it thus suffices to find the diffusion process that goes the fastest; i.e., the value of ℓ for which the speed measure υℓ,W1t is optimized.The speed measure in (9) is quite intuitive; it is in fact similar to that obtained when studying i.i.d. target densities. The main difference lies in the form of γx,z which, in the i.i.d. case, is given by the constant term γ=E∂/∂XlogfX2. The second term in (11) is thus equivalent to γ and consists in a measure of roughness of the conditional density fxix1 under a variation of xi (i≥2). In the case of hierarchical target distributions, we find an extra term that might be viewed as a measure of roughness of fxix1 under a variation of x1. This term is weighted by z12, the square of the (standardized) candidate increment for the first component; in other words, the further the candidate y1 is from the current x1, the greater is the weight attributed to the associated measure of roughness. Of course, in optimizing the speed measure function, we do not need to know in advance the exact value of the proposed standardized increment z1; the speed measure averages over this quantity.It is interesting to note that optimizing the speed measure leads to local proposal variances of the formℓ^2W1t/n. Such proposal variances would then be used for proposing a candidate at the next instant t+dt, given the position of the mixing coordinate at time t. These local proposal variances thus vary from one iteration to another, by opposition to usual tunings in the literature that are fixed for the duration of the algorithm. Naturally, if both expectations in (11) are constant with respect to x1, then the proposal variance obtained by maximizing the speed measure also is constant.Remark 2.
It turns out that local proposal variances optimizing (9) are bounded above by 2.38/EX1/2∂/∂XlogfXx12, the asymptotically optimal scaling (AOS) values for targets with i.i.d. components given a fixed X1=x1. Indeed, if X1=x1 is fixed across iterations, we find ourselves in an i.i.d. setting and the associated speed measure is expressed as 2ℓ2Φ−ℓEX1/2∂/∂XlogfXx12/2. The mentioned upper bounds then follow from the fact that the function Φ⋅ in (9) decreases faster in ℓ than Φ⋅ in the above expression.
Relying on a local varianceℓ^x1 to propose a candidate for the next time interval is usually time-consuming, as it involves numerically solving for the appropriate local proposal variance at every iteration. Since the process is assumed to start in stationarity and X1 explores its conditional state space faster than the other coordinates, we might determine a value ℓ^ that is fixed across iterations by integrating the speed measure υℓ,⋅ over X1 with respect to the marginal distribution f1. Hence, the global (unconditional) asymptotically optimal scaling value ℓ^ maximizes the function(12)EX1υℓ,X1=2ℓ2EX1,Z1Φ−ℓ2γ1/2X1,Z1=2ℓ2∫X1∫ℝΦ−ℓ2γ1/2x1,z1ϕz1f1x1dz1dx1,where ϕ⋅ is the probability density function of a standard normal random variable.Remark 3.
The asymptotic process introduced in Theorem 2 naturally leads us to the concept of local proposal variances. It is however unclear whether the local tunings obtained by maximizing (9) really optimize the mixing rate of the algorithm. Indeed, the proof of Theorem 2 is carried out with ℓ2 constant; this allows, among other things, relying on the simplified form for the acceptance probability. In order to claim that the local proposal variances obtained are optimal, a weak convergence result would need to be proven using a general proposal variance of the form σ2n,x1=ℓ2x1/n. This extension is not trivial, as the ratio of proposal densities qnx;y/qny;x would then need to be included in the acceptance probability. Since the concept of locally optimal proposal variances is numerically demanding in the current framework, we choose to focus on ℓ2 constant.
In RWM-within-Gibbs, the blocksX1 and X2:n are updated consecutively and the situation is therefore different. In that case, local variances of the form σ2n,x1=ℓ2x1/n obtained by maximizing (9) may be used to update the block X2:n. Since X1 is updated separately, the first term in (11) is null, which makes local variances easier to compute. Furthermore, since local variances only depend on X1 (which is updated separately), the ratio qnx;y/qny;x is equal to 1 and does not need to be included in the acceptance probability. Local variances are thus very appealing in that context and shall be studied in Section 5.
Rather than tuning the sampler using the global AOS value, one may instead monitor the acceptance rate in order to work with an optimally mixing version of the RWM algorithm. To express optimal scaling results in terms of acceptance rates, we introduce the expected acceptance rate of then-dimensional stationary RWM algorithm with a normal proposal:(13)anℓ=∬αx,yℓn−nϕny−xℓ/nπxdydx,where ϕn⋅ denotes the probability density function of an n-dimensional standard normal random variable. Optimal mixing results for the RWM sampler are summarized in the following corollary.Corollary 1.
In the settings of Theorem 2, the global asymptotically optimal scaling valueℓ^ maximizes(14)2ℓ2∫X1∫ℝΦ−ℓ2γ1/2x,zϕzf1xdzdx.Furthermore, we have that(15)limn⟶∞anℓ=aℓ≡2∫X1∫ℝΦ−ℓ2γ1/2x,zϕzf1xdzdx,and the corresponding asymptotically optimal acceptance rate is given by aℓ^.In contrast to the i.i.d. case, the AOAR found is not independent of the densitiesf1 and f. Hence, there is not a huge advantage in choosing to tune the acceptance rate of the algorithm over the proposal variance; in fact, both approaches involve the same effort. Although it would also be possible to compute an overall acceptance rate associated to using local proposal variances, it could not be used to tune the algorithm. Building an optimal Markov chain based on local proposal variances would imply modifying the proposal variance at every iteration, which cannot be achieved by solely monitoring the acceptance rate.For simplicity, the theoretical results expounded in this section attribute the same tuning constantℓ to all n components. In practice, when a RWM algorithm is used to sample from a hierarchical target, users will likely want to use a different proposal variance for the mixing component X1. In fact, the proofs of Theorems 1 and 2 easily generalize to the case of inhomogeneous proposal variances.Corollary 2.
LetY1∼Nx1,ℓ2κ12/n with 0<κ1<∞ and Y2:n∼Nx2:n,ℓ2In−1/n, where Y1,Y2:n are independent. Then, Theorems 1 and 2 hold as stated, except that the limiting proposal distribution in Theorem 1 is Y˜1∼Nx˜1,ℓ2κ12 and the random variable Z1 in Theorem 2 is such that Z1∼N0,κ12.In this paper, we consider the simple, yet useful hierarchical model described in (1) and featuring a single mixing component X1. This is a natural starting point to study weak convergence of RWM algorithms for hierarchical targets, and even for correlated targets in general. There exist many generalizations of (1), just as there are many extensions of the proposal distribution considered. Some extensions of the hierarchical target are considered in the discussion, but we do not aim at presenting a detailed treatment of these cases.


## 3.1. Optimal Tuning of the RWM Algorithm

To be confident that then-dimensional chain has entirely explored its state space, we must be certain that every one-dimensional path has explored its own space. In the correlated framework considered, the overall mixing rate of the RWM sampler is only as fast as the slowest component. As explained in Section 3, optimal mixing of the algorithm shall be attained by optimizing the mixing of the coordinates Xi, i=2,…,n. In the limit, the only quantity that depends on the proposal variance (i.e., on ℓ) is υℓ,W1t in (9). To optimize mixing, it thus suffices to find the diffusion process that goes the fastest; i.e., the value of ℓ for which the speed measure υℓ,W1t is optimized.The speed measure in (9) is quite intuitive; it is in fact similar to that obtained when studying i.i.d. target densities. The main difference lies in the form of γx,z which, in the i.i.d. case, is given by the constant term γ=E∂/∂XlogfX2. The second term in (11) is thus equivalent to γ and consists in a measure of roughness of the conditional density fxix1 under a variation of xi (i≥2). In the case of hierarchical target distributions, we find an extra term that might be viewed as a measure of roughness of fxix1 under a variation of x1. This term is weighted by z12, the square of the (standardized) candidate increment for the first component; in other words, the further the candidate y1 is from the current x1, the greater is the weight attributed to the associated measure of roughness. Of course, in optimizing the speed measure function, we do not need to know in advance the exact value of the proposed standardized increment z1; the speed measure averages over this quantity.It is interesting to note that optimizing the speed measure leads to local proposal variances of the formℓ^2W1t/n. Such proposal variances would then be used for proposing a candidate at the next instant t+dt, given the position of the mixing coordinate at time t. These local proposal variances thus vary from one iteration to another, by opposition to usual tunings in the literature that are fixed for the duration of the algorithm. Naturally, if both expectations in (11) are constant with respect to x1, then the proposal variance obtained by maximizing the speed measure also is constant.Remark 2.
It turns out that local proposal variances optimizing (9) are bounded above by 2.38/EX1/2∂/∂XlogfXx12, the asymptotically optimal scaling (AOS) values for targets with i.i.d. components given a fixed X1=x1. Indeed, if X1=x1 is fixed across iterations, we find ourselves in an i.i.d. setting and the associated speed measure is expressed as 2ℓ2Φ−ℓEX1/2∂/∂XlogfXx12/2. The mentioned upper bounds then follow from the fact that the function Φ⋅ in (9) decreases faster in ℓ than Φ⋅ in the above expression.
Relying on a local varianceℓ^x1 to propose a candidate for the next time interval is usually time-consuming, as it involves numerically solving for the appropriate local proposal variance at every iteration. Since the process is assumed to start in stationarity and X1 explores its conditional state space faster than the other coordinates, we might determine a value ℓ^ that is fixed across iterations by integrating the speed measure υℓ,⋅ over X1 with respect to the marginal distribution f1. Hence, the global (unconditional) asymptotically optimal scaling value ℓ^ maximizes the function(12)EX1υℓ,X1=2ℓ2EX1,Z1Φ−ℓ2γ1/2X1,Z1=2ℓ2∫X1∫ℝΦ−ℓ2γ1/2x1,z1ϕz1f1x1dz1dx1,where ϕ⋅ is the probability density function of a standard normal random variable.Remark 3.
The asymptotic process introduced in Theorem 2 naturally leads us to the concept of local proposal variances. It is however unclear whether the local tunings obtained by maximizing (9) really optimize the mixing rate of the algorithm. Indeed, the proof of Theorem 2 is carried out with ℓ2 constant; this allows, among other things, relying on the simplified form for the acceptance probability. In order to claim that the local proposal variances obtained are optimal, a weak convergence result would need to be proven using a general proposal variance of the form σ2n,x1=ℓ2x1/n. This extension is not trivial, as the ratio of proposal densities qnx;y/qny;x would then need to be included in the acceptance probability. Since the concept of locally optimal proposal variances is numerically demanding in the current framework, we choose to focus on ℓ2 constant.
In RWM-within-Gibbs, the blocksX1 and X2:n are updated consecutively and the situation is therefore different. In that case, local variances of the form σ2n,x1=ℓ2x1/n obtained by maximizing (9) may be used to update the block X2:n. Since X1 is updated separately, the first term in (11) is null, which makes local variances easier to compute. Furthermore, since local variances only depend on X1 (which is updated separately), the ratio qnx;y/qny;x is equal to 1 and does not need to be included in the acceptance probability. Local variances are thus very appealing in that context and shall be studied in Section 5.
Rather than tuning the sampler using the global AOS value, one may instead monitor the acceptance rate in order to work with an optimally mixing version of the RWM algorithm. To express optimal scaling results in terms of acceptance rates, we introduce the expected acceptance rate of then-dimensional stationary RWM algorithm with a normal proposal:(13)anℓ=∬αx,yℓn−nϕny−xℓ/nπxdydx,where ϕn⋅ denotes the probability density function of an n-dimensional standard normal random variable. Optimal mixing results for the RWM sampler are summarized in the following corollary.Corollary 1.
In the settings of Theorem 2, the global asymptotically optimal scaling valueℓ^ maximizes(14)2ℓ2∫X1∫ℝΦ−ℓ2γ1/2x,zϕzf1xdzdx.Furthermore, we have that(15)limn⟶∞anℓ=aℓ≡2∫X1∫ℝΦ−ℓ2γ1/2x,zϕzf1xdzdx,and the corresponding asymptotically optimal acceptance rate is given by aℓ^.In contrast to the i.i.d. case, the AOAR found is not independent of the densitiesf1 and f. Hence, there is not a huge advantage in choosing to tune the acceptance rate of the algorithm over the proposal variance; in fact, both approaches involve the same effort. Although it would also be possible to compute an overall acceptance rate associated to using local proposal variances, it could not be used to tune the algorithm. Building an optimal Markov chain based on local proposal variances would imply modifying the proposal variance at every iteration, which cannot be achieved by solely monitoring the acceptance rate.For simplicity, the theoretical results expounded in this section attribute the same tuning constantℓ to all n components. In practice, when a RWM algorithm is used to sample from a hierarchical target, users will likely want to use a different proposal variance for the mixing component X1. In fact, the proofs of Theorems 1 and 2 easily generalize to the case of inhomogeneous proposal variances.Corollary 2.
LetY1∼Nx1,ℓ2κ12/n with 0<κ1<∞ and Y2:n∼Nx2:n,ℓ2In−1/n, where Y1,Y2:n are independent. Then, Theorems 1 and 2 hold as stated, except that the limiting proposal distribution in Theorem 1 is Y˜1∼Nx˜1,ℓ2κ12 and the random variable Z1 in Theorem 2 is such that Z1∼N0,κ12.In this paper, we consider the simple, yet useful hierarchical model described in (1) and featuring a single mixing component X1. This is a natural starting point to study weak convergence of RWM algorithms for hierarchical targets, and even for correlated targets in general. There exist many generalizations of (1), just as there are many extensions of the proposal distribution considered. Some extensions of the hierarchical target are considered in the discussion, but we do not aim at presenting a detailed treatment of these cases.


## 4. Numerical Studies

To illustrate the theoretical results of Section 3, we consider two toy examples: the first target distribution considered is a normal-normal hierarchical model in which the componentsX2,…,Xn are related through their mean, while the second one is a gamma-normal hierarchical model in which X2,…,Xn are related through their variance. In both cases, we show how to compute the optimal variance ℓ^. We also study the performance of RWM samplers and conclude that even in relatively low-dimensional settings, the samplers behave according to the asymptotic results previously detailed.


### 4.1. Normal-Normal Hierarchical Distribution

Consider ann-dimensional hierarchical target such that X1∼N0,1 and XiX1∼NX1,1 for i=2,…,n. To sample from this distribution, we use a RWM algorithm with a Nx,ℓ2In/n proposal distribution. This simple target shall relate Theorem 2 to the theoretical results derived in [8].Standard calculations lead toX1X2:n∼N∑i=2nXi/n,1/n; as n⟶∞, VX1X2:n⟶0 almost surely. If we let μn=∑i=2nXi/n and X˜1=n1/2X1−μn, then X˜1X2:n∼N0,1. Furthermore, the term ∑i=2nXi−μn−X˜/n2/n is reexpressed as ∑i=2nXi−X12/n=∑i=2nZi2/n and thus converges in probability to EZ2=∫∂/∂xlogfxμ2fxμdx=1, where Z1,…,Zn denote independent standard normal random variables. By Theorem 1, we can thus affirm that the component X˜1 asymptotically behaves according to a one-dimensional RWM algorithm with a standard normal target and acceptance function as in (7); these do not, in the current case, depend on x¯.Evaluating the functionγx1,z1 in (11) is a simple task and leads to γx1,z1=z12+1. The AOS value is then found by maximizing(16)υℓ=2ℓ2EZ1Φ−ℓ2Z12+1,with respect to ℓ, where Z1∼N0,1. This yields an AOS of ℓ^2=4.00 and a corresponding AOAR of υℓ^/ℓ^2=0.205. These values are naturally smaller than those obtained for a target with i.i.d. components (5.66 and 0.234, respectively); indeed, the proposal distribution is formed of i.i.d. components and accordingly better suited for similar targets. Relying on a proposal with correlated components would however require a certain understanding of the target correlation structure, which goes against the general framework we wish to consider.It is worth pointing out that the speed measure of the limiting diffusion process does not depend onX1 in the present case. This holds for arbitrary densities f1 and f satisfying the conditions in Section 2.1, provided that X1 is a location parameter for Xi (i≥2). Since a variation in the location parameter does not perturb the roughness of the distribution, the AOS and AOAR found are valid both locally and globally. This means that ℓ^, which remains fixed across iterations, is the best possible proposal scaling conditionally on the last position of the component X1 (i.e., ℓ^=ℓ^x1).A second peculiarity of this example is that the target distribution is jointly normal with mean0 and n×n covariance matrix Σn given by σ12=1, σj2=2 (j=2,…,n), and σi,j=1 ∀i≠j (i,j=1,…,n). Normal distributions being invariant under orthogonal transformations, we can find a transformation under which the target components become mutually independent. The covariance matrix Σn is thus transformed into a diagonal matrix whose diagonal elements consist in the eigenvalues of Σn. In moderate to large dimensions, the eigenvalues can be approximated by 1/n+1,n+1,1,…,1. It turns out that the optimal scaling problem for target distributions of this sort (i.e., formed of components that are i.i.d. up to a scaling term) has been studied in [13]. Solving for the AOS value and AOAR of the transformed target using Theorem 1 and Corollary 2 in [8] leads to values that are consistent with those obtained using Theorem 2 in Section 3.To illustrate these theoretical results, we consider the 20-dimensional normal-normal target described above and run 50 RWM algorithms that differ by their proposal variance only. For each sampler, we perform 100,000 iterations (sufficient for convergence according to the autocorrelation function) and measure efficiency by recording the average squared jumping distance(17)ASJD=1N∑j=1N∑i=1nxinj−xinj−12,where N is the number of iterations and n is the dimension of the target distribution. We also record the average acceptance rate of each algorithm, expressed as(18)AAR=1N∑j=1N1xnj≠xnj−1.We repeat these steps for 50- and 100-dimensional normal-normal targets and combine all three curves of efficiency versus acceptance rate on a graph along with the theoretical efficiency curve ofυℓ versus the expected acceptance rate υℓ/ℓ2 (Figure 1(b), bottom set of curves). To assess the limiting behaviour of the coordinate X1, we also plot the ASJD of this single component (for the 20-, 50-, and 100-dimensional cases) along with the ASJD for the limiting one-dimensional RWM sampler described in Theorem 1 (Figure 1(a), top set of curves).Figure 1
Efficiency of RWM algorithm against acceptance rate for the normal-normal hierarchical target. (a) Efficiency ofX1 only; the top set of curves corresponds to homogeneous proposal variances. (b) Efficiency of all n components; the top set of curves now corresponds to inhomogeneous proposal variances.
(a)
(b)We now repeat the numerical experiment by taking advantage of the available target variances in the tuning of the proposal distribution. Specifically, we letY1∼Nx1,ℓ2/2n be independent of Y2:n∼Nx2:n,ℓ2/n and run the RWM algorithm in dimensions 20, 50, and 100. The resulting simulated and theoretical efficiency curves are illustrated in Figure 1(a) (bottom set of curves) and Figure 1(b) (top set of curves). Although efficiency curves for X1 are lower when using inhomogeneous proposal variances, this approach still results in a better overall performance (the curves in the right graph are higher than with homogeneous variances). The optimized theoretical efficiency is 0.974, which is related to an AOAR of 0.221. Despite the fact that Theorems 1 and 2 are valid asymptotically, the simulation study yields efficiency curves that are very close together; the theorems thus seem applicable in relatively low-dimensional settings.Each set of curves in Figure1(b) agrees about the optimal acceptance rates 0.205 and 0.221, respectively. These optimal rates have been obtained by running an homogeneous sampler with optimal variance ℓ^2/n=4/n and an inhomogeneous sampler with optimal variance 4.4/n, each optimizing (9). Any other proposal variance leads to a point that is lower on the efficiency curve.According to the shape of these curves, tuning the acceptance rate anywhere between 0.15 and 0.3 would yield a loss of at most 10% in efficiency and would still result in a Markov chain that rapidly explores its state space; in particular, using the usual 0.234 for this target would yield an almost optimal algorithm. Beyond finding the exact AOAR for a specific target distribution, there is thus a need for understanding when and why AOARs significantly differ from 0.234. At the present time, the only way to answer this question is by solving the optimal scaling problem for target distributions of interest.


### 4.2. Gamma-Normal Hierarchical Distribution

As a second example, consider a gamma-normal hierarchical target such thatX1∼Γα,λ and XiX1∼N0,1/X1, i=2,…,n. Although Xi(i≥2) are still normally distributed, the coordinate X1 now acts through the variance of the normal variables. This results in a target that significantly differs from the distribution considered in the previous section, falling slightly outside the framework of Section 2 (∂/∂x1logfxx1 is now only locally Lipschitz continuous). We run the usual RWM algorithm to obtain a sample from this distribution.Standard calculations lead toX1X2:n∼Γα+n−1/2,λ+∑i=2nXi2/2 and as n⟶∞, VX1X2:n⟶p0. The WLLN-type expression in Theorem 1 may be reexpressed as ∑i=2nμn+X˜1/n2Xi2/n  =  μn+X˜1/n∑i=2nZi2/n, where Z1:n are independent standard normal random variables. The condition is thus satisfied as it converges in probability to μX¯=EX∂/∂XlogfXμ2. Using Stirling’s formula, it is not difficult to show that the density of X˜1X2:n converges almost surely to that of a N0,2/μ2X¯. By Theorem 1, the coordinate X˜1 asymptotically behaves according to an atypical one-dimensional RWM algorithm with a normal target; the target variance however varies from one iteration to the next, and so does the acceptance function in (7).To optimize the efficiency of the algorithm, we analyze the speed measure in (9); in the present case, it is expressed as(19)υℓ,x1=2ℓ2EZ1Φ−ℓ212Z12x12+x1,where Z1∼N0,1. Maximizing the function EX1υℓ,X1 in Corollary 1 with respect to ℓ leads to the global AOS value, which is fixed across iterations; when α,λ=3,1 for instance, we find ℓ^2=2.40 and AOAR=0.204.The simulation study described in Section 4.1 has been performed for the gamma-normal target model with variousα and λ. Specifically, for fixed α and λ, we consider a 10-dimensional gamma-normal target distribution and run 50 RWM algorithms possessing their own proposal variance. For each sampler, we perform 1,000,000 iterations (again sufficient for convergence according to the autocorrelation function) and measure efficiency by recording the ASJD of each chain. We then repeat these steps for 20- and 50-dimensional targets. Table 1 presents the optimal efficiency and acceptance rate for various α and λ. Those results are compared to the theoretical optimal values obtained by maximizing EX1υℓ,X1.Table 1
Optimal efficiency and acceptance rate of chains in various dimension (n = 10, 20, and 50) for different parameters α and λ of the gamma distribution for X1. The theoretical optimal efficiency and acceptance rate are also included for comparison.
Parameters
Optimal efficiency
Optimal acceptance rate
Theoretical
n
=
10
n
=
20
n
=
50
Theoretical
n
=
10
n
=
20
n
=
50
α
=
2
,
λ
=
1
0.6381
0.6036
0.6246
0.6456
0.1934
0.1984
0.1968
0.1857
α
=
2
,
λ
=
2
0.8169
0.7430
0.7862
0.8239
0.1815
0.1888
0.1759
0.1886
α
=
2
,
λ
=
3
0.8420
0.7623
0.8170
0.8608
0.1517
0.1682
0.1527
0.1593
α
=
3
,
λ
=
1
0.4889
0.4503
0.4736
0.4926
0.2037
0.2370
0.2158
0.2001
α
=
3
,
λ
=
2
0.7541
0.6739
0.7139
0.7405
0.2038
0.2265
0.2233
0.2040
α
=
3
,
λ
=
3
0.8648
0.7554
0.8075
0.8497
0.1922
0.1931
0.1930
0.1882Although the corresponding graphs are omitted here, they yield curves similar to those obtained in Figure1 for the normal-normal target. We note that even if the gamma-normal departs from a jointly normal distribution assumption and does not yield as nice a target distribution as in the previous example, the AOAR obtained is not too far from the 0.234 found for i.i.d. targets. The AOAR however tends to decrease as λ increases (e.g., 0.152 for α,λ=2,3).In the current example, it also turns out that the agreement between theoretical and simulation results is altered for some valuesα,λ. As mentioned above, one of the Lipschitz conditions is only valid locally and so the change in ∂/∂x1logfxx1 becomes arbitrarily steep as X1⟶0. The amplitude of X1 movements is, therefore, not adequately controlled for some choices of α,λ that yield a density f1 assigning a significant probability close to 0. In cases where regularity assumptions are not all satisfied, the applicability of theoretical results may thus be affected by the choice of hyperparameters.


## 4.1. Normal-Normal Hierarchical Distribution

Consider ann-dimensional hierarchical target such that X1∼N0,1 and XiX1∼NX1,1 for i=2,…,n. To sample from this distribution, we use a RWM algorithm with a Nx,ℓ2In/n proposal distribution. This simple target shall relate Theorem 2 to the theoretical results derived in [8].Standard calculations lead toX1X2:n∼N∑i=2nXi/n,1/n; as n⟶∞, VX1X2:n⟶0 almost surely. If we let μn=∑i=2nXi/n and X˜1=n1/2X1−μn, then X˜1X2:n∼N0,1. Furthermore, the term ∑i=2nXi−μn−X˜/n2/n is reexpressed as ∑i=2nXi−X12/n=∑i=2nZi2/n and thus converges in probability to EZ2=∫∂/∂xlogfxμ2fxμdx=1, where Z1,…,Zn denote independent standard normal random variables. By Theorem 1, we can thus affirm that the component X˜1 asymptotically behaves according to a one-dimensional RWM algorithm with a standard normal target and acceptance function as in (7); these do not, in the current case, depend on x¯.Evaluating the functionγx1,z1 in (11) is a simple task and leads to γx1,z1=z12+1. The AOS value is then found by maximizing(16)υℓ=2ℓ2EZ1Φ−ℓ2Z12+1,with respect to ℓ, where Z1∼N0,1. This yields an AOS of ℓ^2=4.00 and a corresponding AOAR of υℓ^/ℓ^2=0.205. These values are naturally smaller than those obtained for a target with i.i.d. components (5.66 and 0.234, respectively); indeed, the proposal distribution is formed of i.i.d. components and accordingly better suited for similar targets. Relying on a proposal with correlated components would however require a certain understanding of the target correlation structure, which goes against the general framework we wish to consider.It is worth pointing out that the speed measure of the limiting diffusion process does not depend onX1 in the present case. This holds for arbitrary densities f1 and f satisfying the conditions in Section 2.1, provided that X1 is a location parameter for Xi (i≥2). Since a variation in the location parameter does not perturb the roughness of the distribution, the AOS and AOAR found are valid both locally and globally. This means that ℓ^, which remains fixed across iterations, is the best possible proposal scaling conditionally on the last position of the component X1 (i.e., ℓ^=ℓ^x1).A second peculiarity of this example is that the target distribution is jointly normal with mean0 and n×n covariance matrix Σn given by σ12=1, σj2=2 (j=2,…,n), and σi,j=1 ∀i≠j (i,j=1,…,n). Normal distributions being invariant under orthogonal transformations, we can find a transformation under which the target components become mutually independent. The covariance matrix Σn is thus transformed into a diagonal matrix whose diagonal elements consist in the eigenvalues of Σn. In moderate to large dimensions, the eigenvalues can be approximated by 1/n+1,n+1,1,…,1. It turns out that the optimal scaling problem for target distributions of this sort (i.e., formed of components that are i.i.d. up to a scaling term) has been studied in [13]. Solving for the AOS value and AOAR of the transformed target using Theorem 1 and Corollary 2 in [8] leads to values that are consistent with those obtained using Theorem 2 in Section 3.To illustrate these theoretical results, we consider the 20-dimensional normal-normal target described above and run 50 RWM algorithms that differ by their proposal variance only. For each sampler, we perform 100,000 iterations (sufficient for convergence according to the autocorrelation function) and measure efficiency by recording the average squared jumping distance(17)ASJD=1N∑j=1N∑i=1nxinj−xinj−12,where N is the number of iterations and n is the dimension of the target distribution. We also record the average acceptance rate of each algorithm, expressed as(18)AAR=1N∑j=1N1xnj≠xnj−1.We repeat these steps for 50- and 100-dimensional normal-normal targets and combine all three curves of efficiency versus acceptance rate on a graph along with the theoretical efficiency curve ofυℓ versus the expected acceptance rate υℓ/ℓ2 (Figure 1(b), bottom set of curves). To assess the limiting behaviour of the coordinate X1, we also plot the ASJD of this single component (for the 20-, 50-, and 100-dimensional cases) along with the ASJD for the limiting one-dimensional RWM sampler described in Theorem 1 (Figure 1(a), top set of curves).Figure 1
Efficiency of RWM algorithm against acceptance rate for the normal-normal hierarchical target. (a) Efficiency ofX1 only; the top set of curves corresponds to homogeneous proposal variances. (b) Efficiency of all n components; the top set of curves now corresponds to inhomogeneous proposal variances.
(a)
(b)We now repeat the numerical experiment by taking advantage of the available target variances in the tuning of the proposal distribution. Specifically, we letY1∼Nx1,ℓ2/2n be independent of Y2:n∼Nx2:n,ℓ2/n and run the RWM algorithm in dimensions 20, 50, and 100. The resulting simulated and theoretical efficiency curves are illustrated in Figure 1(a) (bottom set of curves) and Figure 1(b) (top set of curves). Although efficiency curves for X1 are lower when using inhomogeneous proposal variances, this approach still results in a better overall performance (the curves in the right graph are higher than with homogeneous variances). The optimized theoretical efficiency is 0.974, which is related to an AOAR of 0.221. Despite the fact that Theorems 1 and 2 are valid asymptotically, the simulation study yields efficiency curves that are very close together; the theorems thus seem applicable in relatively low-dimensional settings.Each set of curves in Figure1(b) agrees about the optimal acceptance rates 0.205 and 0.221, respectively. These optimal rates have been obtained by running an homogeneous sampler with optimal variance ℓ^2/n=4/n and an inhomogeneous sampler with optimal variance 4.4/n, each optimizing (9). Any other proposal variance leads to a point that is lower on the efficiency curve.According to the shape of these curves, tuning the acceptance rate anywhere between 0.15 and 0.3 would yield a loss of at most 10% in efficiency and would still result in a Markov chain that rapidly explores its state space; in particular, using the usual 0.234 for this target would yield an almost optimal algorithm. Beyond finding the exact AOAR for a specific target distribution, there is thus a need for understanding when and why AOARs significantly differ from 0.234. At the present time, the only way to answer this question is by solving the optimal scaling problem for target distributions of interest.


## 4.2. Gamma-Normal Hierarchical Distribution

As a second example, consider a gamma-normal hierarchical target such thatX1∼Γα,λ and XiX1∼N0,1/X1, i=2,…,n. Although Xi(i≥2) are still normally distributed, the coordinate X1 now acts through the variance of the normal variables. This results in a target that significantly differs from the distribution considered in the previous section, falling slightly outside the framework of Section 2 (∂/∂x1logfxx1 is now only locally Lipschitz continuous). We run the usual RWM algorithm to obtain a sample from this distribution.Standard calculations lead toX1X2:n∼Γα+n−1/2,λ+∑i=2nXi2/2 and as n⟶∞, VX1X2:n⟶p0. The WLLN-type expression in Theorem 1 may be reexpressed as ∑i=2nμn+X˜1/n2Xi2/n  =  μn+X˜1/n∑i=2nZi2/n, where Z1:n are independent standard normal random variables. The condition is thus satisfied as it converges in probability to μX¯=EX∂/∂XlogfXμ2. Using Stirling’s formula, it is not difficult to show that the density of X˜1X2:n converges almost surely to that of a N0,2/μ2X¯. By Theorem 1, the coordinate X˜1 asymptotically behaves according to an atypical one-dimensional RWM algorithm with a normal target; the target variance however varies from one iteration to the next, and so does the acceptance function in (7).To optimize the efficiency of the algorithm, we analyze the speed measure in (9); in the present case, it is expressed as(19)υℓ,x1=2ℓ2EZ1Φ−ℓ212Z12x12+x1,where Z1∼N0,1. Maximizing the function EX1υℓ,X1 in Corollary 1 with respect to ℓ leads to the global AOS value, which is fixed across iterations; when α,λ=3,1 for instance, we find ℓ^2=2.40 and AOAR=0.204.The simulation study described in Section 4.1 has been performed for the gamma-normal target model with variousα and λ. Specifically, for fixed α and λ, we consider a 10-dimensional gamma-normal target distribution and run 50 RWM algorithms possessing their own proposal variance. For each sampler, we perform 1,000,000 iterations (again sufficient for convergence according to the autocorrelation function) and measure efficiency by recording the ASJD of each chain. We then repeat these steps for 20- and 50-dimensional targets. Table 1 presents the optimal efficiency and acceptance rate for various α and λ. Those results are compared to the theoretical optimal values obtained by maximizing EX1υℓ,X1.Table 1
Optimal efficiency and acceptance rate of chains in various dimension (n = 10, 20, and 50) for different parameters α and λ of the gamma distribution for X1. The theoretical optimal efficiency and acceptance rate are also included for comparison.
Parameters
Optimal efficiency
Optimal acceptance rate
Theoretical
n
=
10
n
=
20
n
=
50
Theoretical
n
=
10
n
=
20
n
=
50
α
=
2
,
λ
=
1
0.6381
0.6036
0.6246
0.6456
0.1934
0.1984
0.1968
0.1857
α
=
2
,
λ
=
2
0.8169
0.7430
0.7862
0.8239
0.1815
0.1888
0.1759
0.1886
α
=
2
,
λ
=
3
0.8420
0.7623
0.8170
0.8608
0.1517
0.1682
0.1527
0.1593
α
=
3
,
λ
=
1
0.4889
0.4503
0.4736
0.4926
0.2037
0.2370
0.2158
0.2001
α
=
3
,
λ
=
2
0.7541
0.6739
0.7139
0.7405
0.2038
0.2265
0.2233
0.2040
α
=
3
,
λ
=
3
0.8648
0.7554
0.8075
0.8497
0.1922
0.1931
0.1930
0.1882Although the corresponding graphs are omitted here, they yield curves similar to those obtained in Figure1 for the normal-normal target. We note that even if the gamma-normal departs from a jointly normal distribution assumption and does not yield as nice a target distribution as in the previous example, the AOAR obtained is not too far from the 0.234 found for i.i.d. targets. The AOAR however tends to decrease as λ increases (e.g., 0.152 for α,λ=2,3).In the current example, it also turns out that the agreement between theoretical and simulation results is altered for some valuesα,λ. As mentioned above, one of the Lipschitz conditions is only valid locally and so the change in ∂/∂x1logfxx1 becomes arbitrarily steep as X1⟶0. The amplitude of X1 movements is, therefore, not adequately controlled for some choices of α,λ that yield a density f1 assigning a significant probability close to 0. In cases where regularity assumptions are not all satisfied, the applicability of theoretical results may thus be affected by the choice of hyperparameters.


## 5. Applications in Bayesian Contexts

The theoretical results presented in this paper have wide applicability and may be used to improve not only RWM algorithms but other samplers as well (RWM-within-Gibbs, for instance). The examples below study the performance of optimally tuned samplers in the context of hierarchical Bayesian models. They show that the RWM-within-Gibbs sampler with local variances (i.e., variances that are a function of the current state of the chain) is superior to its counterpart with a fixed variance. It is also superior to traditional RWM algorithms and even Adaptive Metropolis (AM) samplers, which use the history of the chain to recursively update the covariance matrix of their proposal distribution (see [5]).


### 5.1. Scottish Secondary School Scores

The dataset ScotsSec in the package mlmRev in R contains the scores attained by 3,435 Scottish secondary school students on a standardized test taken at age 16. The primary schools attended by students are also recorded in this dataset; there aren=148 different primary schools, and the number of students per primary school varies between 1 and 72. We use the following multilevel Bayesian framework to model these data:(20)μ↙↓↘θ1……θnθi∼tνμ,1η,1≤i≤n,↓↓↓Y1,1,…,Y1,r1Yn,1,…,Yn,rnYij∼Nθi,1τ,j=1,…,ri.In this model, the variablesyi,1:ri represent the observed scores obtained by the ri students having attended primary school i, i=1,…148. These observations are modeled according to a normal distribution with mean θi and variance 1/τ. The group sizes range from r148=1 to r61=72. The variables θ1:148, which represent the mean scores of the standardized test for students having attended each of the 148 primary schools, are modeled using a Student distribution with ν=4 degrees of freedom. A translated and scaled Student distribution tνμ,1/η has a density proportional to 1+ηx−μ2/ν−ν+1/2. The mean and precision of the Student distribution, along with the precision of the normally distributed data, are attributed noninformative priors: πμ∝1, πη∝η−1, and πτ∝τ−1.This model leads to then+3-dimensional posterior density:(21)πμ,η,τ,θ1:nYij∝η−1τ−1∏i=1nη1+ηθi−μ2ν−ν+1/2∏i=1n∏j=1riτexp−τ2yij−θi2.The posterior density is too complex for analytic computation, and numerical integration must be ruled out due to the dimensionality of the problem. This distribution is best sampled with MCMC methods, although a classical Gibbs sampler must be ruled out, as the Student distribution destroys conjugacy. In the current setting, we propose to use a RWM-within-Gibbs with four blocks of variables:μ, η, τ, and θ1:n. We are also interested in assessing the performance of full-dimensional RWM and AM algorithms in which μ, η, τ, and θ1:n are updated at once.The RWM-within-Gibbs performs one-dimensional updates ofμ, η, and τ using target densities fμη,τ,θ1:n,Yij, fημ,τ,θ1:n,Yij, and fτμ,η,θ1:n,Yij. It then performs an n-dimensional update of θ1:n with respect to the conditional density fθ1:nμ,η,τ,Yij=∏i=1nfθiμ,η,τ,Yi,1:ri. Since each block of variables is updated individually using a RWM sampler, we may compute local proposal variances for the fourth block using (9) and (11) in Theorem 2. The proposal variances maximizing (9) are adjusted according to the roughness of their corresponding target component’s distribution and should offer a better performance than a fixed proposal variance.The target distribution of the fourth block satisfies(22)fθ1:nμ,η,τYij∝∏i=1n1+ηθi−μ2ν−ν+1/2exp−τ2∑j=1riyij−θi2,and, hence, the partial derivative of the one-dimensional log density with respect to θi is(23)∂∂θilogfθiμ,η,τ,Yi,1:r=τ∑j=1riyij−θi−ν+1νηTi1+Ti2/ν,where Ti=ηθi−μ∼tν0,1, i=1,…,n. Since the variables μ,η,τ are updated separately, the first term in (11) is null, leading to(24)γiμ,η,τ=E∂∂θifθiμ,η,τ,Yi,1:ri2.Optimizing (9) leads to local, inhomogeneous proposal variances of the form 2.382/nγiμ,η,τ.The termsγiμ,η,τ in the proposal variances are not easy to obtain explicitly as the expectation in (24) must be computed with respect to the conditional distribution of θi given μ,η,τ,Yi,1:ri, which is not a Student distribution anymore. However, the terms γiμ,η,τ may be averaged over the random variables Yi,1:ri. Squaring (23) and computing the expectation first with respect to Yi,1:ri and then with respect to θi easily leads to(25)Eγiμ,η,τ=riτ+ην+12νν+2Γν+1/2Γν+4/2Γν/2Γν+5/2.These terms yield local proposal variances that have been averaged over all possible datasets; these are the best local variances for the model under study when no information about the observations is available.The RWM-within-Gibbs is then implemented using Gaussian proposal distributions withσ1=0.95, σ2=0.025, and σ3=0.0005 for μ, η, and τ. This yields acceptance rates in the range 35%–50% for each subalgorithm, as prescribed in the literature for one-dimensional target distributions (see [12]). We update θ1:148 using a Gaussian proposal with local variances 2.382/nEγiμ,η,τ.These steps are then repeated by running a RWM-within-Gibbs in whichθ1:148 is updated using a fixed proposal variance of 52. We also run a 151-dimensional RWM sampler with a Nμ,η,τ,θ1:148,42/151∗diag1,0.01,0.001,1,…,1 proposal distribution, and an AM algorithm in which the tuning factor of the proposal covariance matrix is 8.The ASJD of the chain in (17) offers a reliable way of comparing the four samplers; it is reported in the first column of Table 2. A large value of this measure (relative to other samplers) is indicative of a process that rapidly explores its space, and is equivalent to ordering samplers according to their lag-1 autocorrelations. We also compare the relative efficiency of these samplers by calculating the effective sample size (ESS) of the variables μ, η, τ, and θ2. The effective sample size represents the number of uncorrelated samples that are produced from the output of the sampler. It is also used as a convergence diagnostic: when its value is too small (<100), we may have reasonable doubts that the chain really has converged. It is computed as(26)ESS=N1+2∑k=1∞γk,where N is the number of samples and ∑k=1∞γk is the sum of lag-k sample autocorrelations. An ESS is produced for each variable; since we want to measure the number of samples that are uncorrelated over all variables, we report the minimum ESS (2nd column of Table 2). The ASJD and minimum ESS values are averaged over 10 runs of 100,000 iterations each, with a burn-in period of 1,000. These quantities are then normalized relative to the average running time of samplers (3rd column); this, respectively, yields the average square jumping distance per second (4th column), and the number of uncorrelated samples generated every second (5th column).Table 2
Scottish dataset: efficiency and time-adjusted efficiency measures for the four samplers tested.
Sampler
Efficiency
Time-adjusted efficiency
Mean ASJD
Min ESS
Mean time
a
/
s
e
/
s
(a)
(e)
(s)
×
100
×
100
RWM
2.9712
74.30
145.50
2.0421
51.07
Fixed RWM-w-G
6.4279
157.09
147.77
4.3499
106.31
Local RWM-w-G
8.4108
272.70
148.06
5.6807
184.18
Adapt. Met.
5.2476
473.83
1,081.52
0.4852
43.81According to these results, the RWM-within-Gibbs with local variances is 1.3 times more efficient than the one with a fixed variance; the efficiency gain is even greater (1.7) if we consider the minimum ESS instead of the ASJD. Although the RWM sampler offers a slight improvement in terms of running time, it still results in efficiency measures that are significantly smaller than those of the RWM-within-Gibbs. The Adaptive Metropolis sampler could be an interesting alternative to the RWM-within-Gibbs, if it were not as expensive in terms of computational resources. Indeed, even if its ASJD is smaller than that of the RWM-within-Gibbs, its minimum ESS is greater. This sampler however requires significantly more time than the other samplers to complete its 100,000 iterations. When correcting for computational effort, it thus badly loses ground to its competitors. The results in Table2 thus illustrate that there is an important efficiency gain that is available from preferring a local RWM-within-Gibbs over its constant counterpart. Given that running times for both approaches are equivalent, we should clearly use local proposal variances whenever possible.


### 5.2. Stochastic Volatility Model

As a second example, we wish to study the performance of MCMC samplers in the context of a Bayesian hierarchical model that does not respect the regularity assumptions imposed by the theory of Section 3. We consider a stochastic volatility model in which the latent volatilities form an order-1 autoregressive process. The model, similar to those studied in [16, 17], expresses the mean corrected returns di and log volatilites Xi, for i≥1, as(27)di=εiexpXi2,Xi+1=ϕXi+ηi+1.The variablesεi∼N0,1 and ηi∼N0,τ2 are uncorrelated white noises, and we set X1∼N0,τ2/1−ϕ2. Priors for the parameters τ2 and ϕ are τ2∼IΓδ,λ and ϕ+1/2∼βa,b, where IΓδ,λ is the inverse gamma distribution with density proportional to x−δ+1e−λ/x. This model leads to an n+2-dimensional posterior density πτ2,ϕ,X1,…,Xnd1:n.Before pursuing the analysis, we note thatτ2 and ϕ are constrained to subsets of the real line; since the target density is rather sensitive to changes in these parameters, this will potentially affect the performance of MCMC approaches. To ensure fluidity in the samplers implemented, we apply the transformations τ2=expκ and ϕ=tanhω. The new variables κ and ω take values in and the resulting n+2-dimensional posterior density is given by(28)πκ,ω,x1:nd1:n∝exp−κn2+δe−ω2b+11+e−2ωa+b+1exp−12∑i=1nxi+di2e−xi×exp−e−κ22λ+4e−2ω1+e−2ω2x12+∑i=2nxi−1−e−2ω1+e−2ωxi−12.Using a 100-dimensional datasetd1:100 exhibiting low correlation (obtained from the stochastic volatility model with ϕ=0.1 and τ2=0.75), we sample this posterior density using RWM-within-Gibbs (local and fixed variances), traditional RWM, and AM algorithms. Hyperparameters are set to δ=1, λ=0.75, a=10, and b=6.For the RWM-within-Gibbs, we propose to divide the variables into 3 blocks:κ, ω, and X1:n. The proposal standard deviations associated to κ and ω are set to 0.2 and 0.27, respectively; each subalgorithm thus accepts candidates according to a rate of ≈45%. The n-dimensional update of X1:n is performed according to the conditional target density πx1:nκ,ω,d1:n. In the case of the RWM-within-Gibbs with local variances, the terms(29)γiκ,ω=E∂∂XilogπX1:nκ,ω,d1:n2,i=1,…,n,in (11) are not easy to obtain as the full conditional distribution (given the data) is not normally distributed anymore. As before, we solve this problem by computing the expectation above with respect to d1:n first, and then with respect to X1:n. The resulting proposal variances are thus averaged over all possible datasets; they are the best local proposal variances, independently of the specific dataset considered. Optimizing (9) for i=1,…,n yields the n-dimensional vector:(30)2.382n12+e−κ,12+e−κ1+1−e−2ω1+e−2ω2,…,12+e−κ1+1−e−2ω1+e−2ω2,12+e−κ−1.For the RWM-within-Gibbs with a fixed proposal variance, the proposal standard deviations associated toκ and ω are still 0.2 and 0.27. We then use the theory of Section 3 to obtain an approximately optimal acceptance rate of 0.2 for the block X1:n. We reach a similar conclusion for the traditional RWM sampler. Naturally, we have to keep in mind that regularity assumptions are violated in the current context; the theoretical results might not be robust to a departure from those assumptions. In fact, given that the Xi’s are correlated, we expect the Adaptive Metropolis sampler to better capture this design and to outdo its competitors.The initial covariance matrix of the Adaptive Metropolis algorithm is then+2-dimensional identity matrix. We tune its acceptance rate as close as possible to 0.234, as suggested in the literature. For each sampler, we average the ASJD and minimum ESS over 10 runs of 200,000 iterations each, from which the first 10,000 iterations are discarded as burn-in. Time-adjusted ESJD and minimum ESS are again used a measures of efficiency; their values are reported in Table 3.Table 3
Stochastic volatility: efficiency and time-adjusted efficiency measures for the four samplers tested.
Sampler
Efficiency
Time-adjusted efficiency
Mean ASJD
Minimum ESS
Mean time
a
/
s
e
/
s
(a)
(e)
(s)
(×1,000)
(×1,000)
RWM
0.3994
103.37
367.34
1.0873
281.40
Fixed RWM-w-G
0.6420
116.58
371.93
1.7261
313.45
Local RWM-w-G
0.6740
132.40
371.38
1.8149
356.51
Adapt. Met.
0.6320
347.76
1,149.26
0.5499
302.59In terms of ASJD, the RWM-within-Gibbs with local variances is the best option, although its competitors also offer decent performances. The AM sampler does better, in absolute, for the minimum ESS; when accounting for computational effort however, the AM ends up outdone by the RWM-within-Gibbs (local and fixed). As before, we notice a net efficiency gain when preferring local variances to a fixed one in the RWM-within-Gibbs (net gain between 5% and 13%, depending on the efficiency measure). This modest gain is explained by the fact that, for the specific model studied, variations inκ and ω do not have a huge impact on the value of the local variances in (30). In spite of this, the impact of using local variances remains positive; generally, there does not seem to be a risk associated to using such local variances. Furthermore, the theoretical results seem applicable to contexts where regularity assumptions are violated (to some extent).


## 5.1. Scottish Secondary School Scores

The dataset ScotsSec in the package mlmRev in R contains the scores attained by 3,435 Scottish secondary school students on a standardized test taken at age 16. The primary schools attended by students are also recorded in this dataset; there aren=148 different primary schools, and the number of students per primary school varies between 1 and 72. We use the following multilevel Bayesian framework to model these data:(20)μ↙↓↘θ1……θnθi∼tνμ,1η,1≤i≤n,↓↓↓Y1,1,…,Y1,r1Yn,1,…,Yn,rnYij∼Nθi,1τ,j=1,…,ri.In this model, the variablesyi,1:ri represent the observed scores obtained by the ri students having attended primary school i, i=1,…148. These observations are modeled according to a normal distribution with mean θi and variance 1/τ. The group sizes range from r148=1 to r61=72. The variables θ1:148, which represent the mean scores of the standardized test for students having attended each of the 148 primary schools, are modeled using a Student distribution with ν=4 degrees of freedom. A translated and scaled Student distribution tνμ,1/η has a density proportional to 1+ηx−μ2/ν−ν+1/2. The mean and precision of the Student distribution, along with the precision of the normally distributed data, are attributed noninformative priors: πμ∝1, πη∝η−1, and πτ∝τ−1.This model leads to then+3-dimensional posterior density:(21)πμ,η,τ,θ1:nYij∝η−1τ−1∏i=1nη1+ηθi−μ2ν−ν+1/2∏i=1n∏j=1riτexp−τ2yij−θi2.The posterior density is too complex for analytic computation, and numerical integration must be ruled out due to the dimensionality of the problem. This distribution is best sampled with MCMC methods, although a classical Gibbs sampler must be ruled out, as the Student distribution destroys conjugacy. In the current setting, we propose to use a RWM-within-Gibbs with four blocks of variables:μ, η, τ, and θ1:n. We are also interested in assessing the performance of full-dimensional RWM and AM algorithms in which μ, η, τ, and θ1:n are updated at once.The RWM-within-Gibbs performs one-dimensional updates ofμ, η, and τ using target densities fμη,τ,θ1:n,Yij, fημ,τ,θ1:n,Yij, and fτμ,η,θ1:n,Yij. It then performs an n-dimensional update of θ1:n with respect to the conditional density fθ1:nμ,η,τ,Yij=∏i=1nfθiμ,η,τ,Yi,1:ri. Since each block of variables is updated individually using a RWM sampler, we may compute local proposal variances for the fourth block using (9) and (11) in Theorem 2. The proposal variances maximizing (9) are adjusted according to the roughness of their corresponding target component’s distribution and should offer a better performance than a fixed proposal variance.The target distribution of the fourth block satisfies(22)fθ1:nμ,η,τYij∝∏i=1n1+ηθi−μ2ν−ν+1/2exp−τ2∑j=1riyij−θi2,and, hence, the partial derivative of the one-dimensional log density with respect to θi is(23)∂∂θilogfθiμ,η,τ,Yi,1:r=τ∑j=1riyij−θi−ν+1νηTi1+Ti2/ν,where Ti=ηθi−μ∼tν0,1, i=1,…,n. Since the variables μ,η,τ are updated separately, the first term in (11) is null, leading to(24)γiμ,η,τ=E∂∂θifθiμ,η,τ,Yi,1:ri2.Optimizing (9) leads to local, inhomogeneous proposal variances of the form 2.382/nγiμ,η,τ.The termsγiμ,η,τ in the proposal variances are not easy to obtain explicitly as the expectation in (24) must be computed with respect to the conditional distribution of θi given μ,η,τ,Yi,1:ri, which is not a Student distribution anymore. However, the terms γiμ,η,τ may be averaged over the random variables Yi,1:ri. Squaring (23) and computing the expectation first with respect to Yi,1:ri and then with respect to θi easily leads to(25)Eγiμ,η,τ=riτ+ην+12νν+2Γν+1/2Γν+4/2Γν/2Γν+5/2.These terms yield local proposal variances that have been averaged over all possible datasets; these are the best local variances for the model under study when no information about the observations is available.The RWM-within-Gibbs is then implemented using Gaussian proposal distributions withσ1=0.95, σ2=0.025, and σ3=0.0005 for μ, η, and τ. This yields acceptance rates in the range 35%–50% for each subalgorithm, as prescribed in the literature for one-dimensional target distributions (see [12]). We update θ1:148 using a Gaussian proposal with local variances 2.382/nEγiμ,η,τ.These steps are then repeated by running a RWM-within-Gibbs in whichθ1:148 is updated using a fixed proposal variance of 52. We also run a 151-dimensional RWM sampler with a Nμ,η,τ,θ1:148,42/151∗diag1,0.01,0.001,1,…,1 proposal distribution, and an AM algorithm in which the tuning factor of the proposal covariance matrix is 8.The ASJD of the chain in (17) offers a reliable way of comparing the four samplers; it is reported in the first column of Table 2. A large value of this measure (relative to other samplers) is indicative of a process that rapidly explores its space, and is equivalent to ordering samplers according to their lag-1 autocorrelations. We also compare the relative efficiency of these samplers by calculating the effective sample size (ESS) of the variables μ, η, τ, and θ2. The effective sample size represents the number of uncorrelated samples that are produced from the output of the sampler. It is also used as a convergence diagnostic: when its value is too small (<100), we may have reasonable doubts that the chain really has converged. It is computed as(26)ESS=N1+2∑k=1∞γk,where N is the number of samples and ∑k=1∞γk is the sum of lag-k sample autocorrelations. An ESS is produced for each variable; since we want to measure the number of samples that are uncorrelated over all variables, we report the minimum ESS (2nd column of Table 2). The ASJD and minimum ESS values are averaged over 10 runs of 100,000 iterations each, with a burn-in period of 1,000. These quantities are then normalized relative to the average running time of samplers (3rd column); this, respectively, yields the average square jumping distance per second (4th column), and the number of uncorrelated samples generated every second (5th column).Table 2
Scottish dataset: efficiency and time-adjusted efficiency measures for the four samplers tested.
Sampler
Efficiency
Time-adjusted efficiency
Mean ASJD
Min ESS
Mean time
a
/
s
e
/
s
(a)
(e)
(s)
×
100
×
100
RWM
2.9712
74.30
145.50
2.0421
51.07
Fixed RWM-w-G
6.4279
157.09
147.77
4.3499
106.31
Local RWM-w-G
8.4108
272.70
148.06
5.6807
184.18
Adapt. Met.
5.2476
473.83
1,081.52
0.4852
43.81According to these results, the RWM-within-Gibbs with local variances is 1.3 times more efficient than the one with a fixed variance; the efficiency gain is even greater (1.7) if we consider the minimum ESS instead of the ASJD. Although the RWM sampler offers a slight improvement in terms of running time, it still results in efficiency measures that are significantly smaller than those of the RWM-within-Gibbs. The Adaptive Metropolis sampler could be an interesting alternative to the RWM-within-Gibbs, if it were not as expensive in terms of computational resources. Indeed, even if its ASJD is smaller than that of the RWM-within-Gibbs, its minimum ESS is greater. This sampler however requires significantly more time than the other samplers to complete its 100,000 iterations. When correcting for computational effort, it thus badly loses ground to its competitors. The results in Table2 thus illustrate that there is an important efficiency gain that is available from preferring a local RWM-within-Gibbs over its constant counterpart. Given that running times for both approaches are equivalent, we should clearly use local proposal variances whenever possible.


## 5.2. Stochastic Volatility Model

As a second example, we wish to study the performance of MCMC samplers in the context of a Bayesian hierarchical model that does not respect the regularity assumptions imposed by the theory of Section 3. We consider a stochastic volatility model in which the latent volatilities form an order-1 autoregressive process. The model, similar to those studied in [16, 17], expresses the mean corrected returns di and log volatilites Xi, for i≥1, as(27)di=εiexpXi2,Xi+1=ϕXi+ηi+1.The variablesεi∼N0,1 and ηi∼N0,τ2 are uncorrelated white noises, and we set X1∼N0,τ2/1−ϕ2. Priors for the parameters τ2 and ϕ are τ2∼IΓδ,λ and ϕ+1/2∼βa,b, where IΓδ,λ is the inverse gamma distribution with density proportional to x−δ+1e−λ/x. This model leads to an n+2-dimensional posterior density πτ2,ϕ,X1,…,Xnd1:n.Before pursuing the analysis, we note thatτ2 and ϕ are constrained to subsets of the real line; since the target density is rather sensitive to changes in these parameters, this will potentially affect the performance of MCMC approaches. To ensure fluidity in the samplers implemented, we apply the transformations τ2=expκ and ϕ=tanhω. The new variables κ and ω take values in and the resulting n+2-dimensional posterior density is given by(28)πκ,ω,x1:nd1:n∝exp−κn2+δe−ω2b+11+e−2ωa+b+1exp−12∑i=1nxi+di2e−xi×exp−e−κ22λ+4e−2ω1+e−2ω2x12+∑i=2nxi−1−e−2ω1+e−2ωxi−12.Using a 100-dimensional datasetd1:100 exhibiting low correlation (obtained from the stochastic volatility model with ϕ=0.1 and τ2=0.75), we sample this posterior density using RWM-within-Gibbs (local and fixed variances), traditional RWM, and AM algorithms. Hyperparameters are set to δ=1, λ=0.75, a=10, and b=6.For the RWM-within-Gibbs, we propose to divide the variables into 3 blocks:κ, ω, and X1:n. The proposal standard deviations associated to κ and ω are set to 0.2 and 0.27, respectively; each subalgorithm thus accepts candidates according to a rate of ≈45%. The n-dimensional update of X1:n is performed according to the conditional target density πx1:nκ,ω,d1:n. In the case of the RWM-within-Gibbs with local variances, the terms(29)γiκ,ω=E∂∂XilogπX1:nκ,ω,d1:n2,i=1,…,n,in (11) are not easy to obtain as the full conditional distribution (given the data) is not normally distributed anymore. As before, we solve this problem by computing the expectation above with respect to d1:n first, and then with respect to X1:n. The resulting proposal variances are thus averaged over all possible datasets; they are the best local proposal variances, independently of the specific dataset considered. Optimizing (9) for i=1,…,n yields the n-dimensional vector:(30)2.382n12+e−κ,12+e−κ1+1−e−2ω1+e−2ω2,…,12+e−κ1+1−e−2ω1+e−2ω2,12+e−κ−1.For the RWM-within-Gibbs with a fixed proposal variance, the proposal standard deviations associated toκ and ω are still 0.2 and 0.27. We then use the theory of Section 3 to obtain an approximately optimal acceptance rate of 0.2 for the block X1:n. We reach a similar conclusion for the traditional RWM sampler. Naturally, we have to keep in mind that regularity assumptions are violated in the current context; the theoretical results might not be robust to a departure from those assumptions. In fact, given that the Xi’s are correlated, we expect the Adaptive Metropolis sampler to better capture this design and to outdo its competitors.The initial covariance matrix of the Adaptive Metropolis algorithm is then+2-dimensional identity matrix. We tune its acceptance rate as close as possible to 0.234, as suggested in the literature. For each sampler, we average the ASJD and minimum ESS over 10 runs of 200,000 iterations each, from which the first 10,000 iterations are discarded as burn-in. Time-adjusted ESJD and minimum ESS are again used a measures of efficiency; their values are reported in Table 3.Table 3
Stochastic volatility: efficiency and time-adjusted efficiency measures for the four samplers tested.
Sampler
Efficiency
Time-adjusted efficiency
Mean ASJD
Minimum ESS
Mean time
a
/
s
e
/
s
(a)
(e)
(s)
(×1,000)
(×1,000)
RWM
0.3994
103.37
367.34
1.0873
281.40
Fixed RWM-w-G
0.6420
116.58
371.93
1.7261
313.45
Local RWM-w-G
0.6740
132.40
371.38
1.8149
356.51
Adapt. Met.
0.6320
347.76
1,149.26
0.5499
302.59In terms of ASJD, the RWM-within-Gibbs with local variances is the best option, although its competitors also offer decent performances. The AM sampler does better, in absolute, for the minimum ESS; when accounting for computational effort however, the AM ends up outdone by the RWM-within-Gibbs (local and fixed). As before, we notice a net efficiency gain when preferring local variances to a fixed one in the RWM-within-Gibbs (net gain between 5% and 13%, depending on the efficiency measure). This modest gain is explained by the fact that, for the specific model studied, variations inκ and ω do not have a huge impact on the value of the local variances in (30). In spite of this, the impact of using local variances remains positive; generally, there does not seem to be a risk associated to using such local variances. Furthermore, the theoretical results seem applicable to contexts where regularity assumptions are violated (to some extent).


## 6. Discussion

In this paper, we have studied the tuning of RWM algorithms applied to single-level hierarchical target distributions. The optimal variance of the Gaussian proposal distribution has been found to depend on a measure of roughness of the densityf with respect to x as before, but also with respect to the mixing coordinate x1. This leads to local proposal variances that are a function of the mixing parameter x1. It is however possible to average over the random variable X1 to find a globally optimal proposal variance. In the case where X1 is a location parameter, it does not affect the roughness of the density f, and the optimal proposal scaling is valid both locally and globally.Higher-level hierarchies could be studied using a similar approach. A target featuringp mixing components, expressed as(31)πx=∏j=1pfjxj∏i=p+1nfxix1:p,with x1:p=x1,…,xp would lead to a result similar to Theorem 2, but with the function(32)γx1:p,z1:p=EX∑j=1pzj∂∂xjlogfXx1:p2+EX∂∂XlogfXx1:p2,where z1:p=z1,…,zp come from independent N0,1 random variables. For a target whose mixing component (Xp say) depends itself on higher-level mixing components X1,…,Xp−1, expressed as πx=f1x1:p∏i=p+1nfxixp, the conclusions of Theorem 2 are still valid. These generalizations also hold for Corollary 2, with obvious adjustments (Z1∼N0,ℓ2κ12/n,…,Zp∼N0,ℓ2κp2/n). Similar extensions may be derived for other hierarchical models.In the simulation study of Section 4, we found that the optimal acceptance rate most often lies around 0.2. In the gamma-normal example, there were some values ofα and λ that led to significantly lower optimal acceptance rates (0.15 when α=2,λ=3). The usual 0.234 is thus quite robust and, if preferred, should lead to an efficient version of the sampler. In the case of correlated targets, it would however be wiser to settle for an acceptance rate slightly below 0.234. Since we investigate correlated targets with a proposal distribution featuring a diagonal covariance matrix, it is not surprising to find an AOAR lower than 0.234; the latter is the AOAR for exploring a target distribution with independent components, which is an ideal situation when relying on a proposal distribution with independent components.We conclude by outlining that the concept of locally optimal proposal variances reveals itself to be of interest with other types of samplers, such as RWM-within-Gibbs algorithms. Indeed, the asymptotic results of Section 3 are proof of the theoretical superiority of RWM-within-Gibbs over RWM when sampling from hierarchical targets. The examples of Section 5 illustrate the efficiency gain from using a RWM-within-Gibbs with local variances over some competitors, including an adaptive sampler. Similar ideas may also be applied to different samplers such as Metropolis-adjusted Langevin algorithms (MALA), but this goes beyond the scope of this paper.



---

*Source: 8740426-2019-08-26.xml*