Variance reduction of estimators arising from Metropolis–Hastings algorithms

The Metropolis–Hastings algorithm is one of the most basic and well-studied Markov chain Monte Carlo methods. It generates a Markov chain which has as limit distribution the target distribution by simulating observations from a different proposal distribution. A proposed value is accepted with some particular probability otherwise the previous value is repeated. As a consequence, the accepted values are repeated a positive number of times and thus any resulting ergodic mean is, in fact, a weighted average. It turns out that this weighted average is an importance sampling-type estimator with random weights. By the standard theory of importance sampling, replacement of these random weights by their (conditional) expectations leads to more efficient estimators. In this paper we study the estimator arising by replacing the random weights with certain estimators of their conditional expectations. We illustrate by simulations that it is often more efficient than the original estimator while in the case of the independence Metropolis–Hastings and for distributions with finite support we formally prove that it is even better than the “optimal” importance sampling estimator.

The authors wish to thank the two anonymous referees for their helpful comments and suggestions which considerably improved the paper.

Authors and Affiliations

1. Department of Statistics and Insurance Science, University of Piraeus, 80 Karaoli & Dimitriou str., 18534, Piraeus, Greece

   George Iliopoulos

2. Department of Engineering Sciences, University of Patras, 26500, Rio, Patras, Greece

   Sonia Malefaki

Corresponding author

Correspondence to George Iliopoulos.

Appendix: Proof of Theorem 1

(a) Assume without loss of generality that the state space is \(\mathcal {X}=\{1, \ldots, m\}\) for some m≥2 and set for convenience h _k =h(k), π _k =π(k), g _k =g(k), g _kl =g(l|k) and w _k =w(k). Then,

(5)

By the ergodic theorem it holds \(n^{-1}\sum_{i=1}^{n}I(X_{i}=k)\stackrel {\mathrm{a.s.}}{\longrightarrow}g_{k}\) and \(\sum_{j=1}^{n}\frac{\xi_{j}}{\sum_{l=1}^{n}\xi_{l}}I(X_{j}=\ell )\stackrel{\mathrm{a.s.}}{\longrightarrow}\pi _{\ell}\). Since (5) contains only finite sums we conclude that it converges almost surely to

because \(g_{k} = \kappa\sum_{i=1}^{m} \min\{\pi_{k} g_{k\ell}, \pi _{\ell }g_{\ell k}\}\). By taking h≡1 we get that \(n^{-1}\sum _{i=1}^{n}\hat{w}(X_{i})\stackrel{\mathrm{a.s.}}{\longrightarrow}\kappa \) and thus

(b) Let us define

and

Billingsley (1961) proved that

where \(\bar{\boldsymbol{U}}= (\bar{U}_{1}, \ldots, \bar{U}_{m})^{T}\), g=(g ₁,…,g _m )^T and the ij entry of Σ ₁₁ is

Here \(g_{ij}^{(n)}\) denotes the n-step transition probability from state i to state j and δ _ij is Kronecker’s delta. Using similar arguments it can be proven that

where the ij entry of the submatrix Σ ₁₂ is σ ₁₂(ij)=κw _j σ ₁₁(ij) while the ij entry of the submatrix Σ ₂₂ is σ ₂₂(ij)=δ _ij g _i κw _i (κw _i −1)+κ ² w _i w _j σ ₁₁(ij).

It is clear that the asymptotic distribution of \(n^{1/2}\{\hat {h}_{\hat {w}}-E_{\pi}(h)\}\) can be found via the standard delta method. Observe that \(\hat{h}_{\hat{w}}\) can be expressed as

(6)

say. By differentiation we get

and

The variance of the asymptotic normal distribution of \(n^{1/2}\{\hat {h}_{\hat{w}}-E_{\pi}(h)\}\) is

where ∇ _u f ₁, ∇ _v f ₁ are the vectors containing the derivatives with respect to u, v, respectively, evaluated at (u,v)=(g,κ π). We will see below that the above expression is in fact as stated in the theorem. To this end, let us also consider the IS “estimator” \(\hat {h}_{IS}\) and evaluate its asymptotic variance too. After some algebra we get that

say. By (6) and (7) we see that the only difference between the two “estimators” is the replacement of π _l by its unbiased estimate \(\kappa^{-1}\bar{V}_{l}\) in \(\hat {h}_{\hat{w}}\). By differentiation we get

Since the corresponding variance of the importance sampling estimator is

it is clear that

Set now \(A_{1} = - 2\nabla_{\boldsymbol{u}}f_{1}^{T} \boldsymbol{\varSigma } _{12} \nabla_{\boldsymbol{v}}f_{1}\), \(A_{2} = -\nabla_{\boldsymbol{u}}f_{1}^{T} \boldsymbol{\varSigma }_{22} \nabla_{\boldsymbol{u}}f_{1}\), B _n =w(X _n )h(X _n ) for n=0,1,…, and b _i =w _i h _i , μ _i =E(B ₁|X ₀=i) for i=1,…,m. Consider further without loss of generality that E _π (h)=0 and note that under stationarity, it holds that E _π (h)=E _g (B _n ) for all n. Then,

But

and since X is reversible, i.e., it holds g _i g _ij =g _j g _ji and more generally \(g_{i} g_{ij}^{(n)} = g_{j} g_{ji}^{(n)}\) for all n, we conclude that

Moreover,

and for all n, we have that

Thus,

On the other hand,

But

and for all n,

Hence,

From (8) and (9) we get

By the standard asymptotic theory for Markov chains \(\sigma ^{2}_{IS}(h)=E_{g}(B_{0}^{2})+2\sum_{n=1}^{\infty}E_{g}(B_{0}B_{n})\). Since \(E_{g}(B_{0}^{2})=\operatorname{Var}_{g}\{w(X_{0})h(X_{0})\}\) part (b) of the theorem follows.

(c) In order to compare the two variances note first that the first term of the sum in (10) is clearly positive. Moreover, we know that \(\sum_{n=2}^{\infty}E_{g}(B_{0}B_{n}) \equiv\sum _{n=2}^{\infty}\mathrm{Cov}_{g}(B_{0},B_{n}) \geq 0\) since each residual sum of autocovariances starting from an even integer is nonnegative (see for example Geyer 1992). Thus a sufficient condition for the variances difference to be positive is E _g (B ₀ B ₁)≥0. This expectation can be expressed as a quadratic form, namely,

where \(\tilde{g}_{ij}=g_{i}g_{ij}\). We will show that in the case of independence MH, the matrix \(\tilde{\boldsymbol{G}}\) is nonnegative definite. Indeed, in this case,

Consider without loss of generality that

and set γ _ij =κ ⁻¹ q _i π _j . Then clearly \(\tilde {g}_{ij}=\gamma_{i\wedge j,i\vee j}\) where i∧j=min{i,j} and i∨j=max{i,j} so \(\tilde{\boldsymbol{G}}\) is nonnegative definite. To see that, proceed by induction to show that all its principal minors are nonnegative. Let \(\tilde{\boldsymbol{G}}_{kk}\) denote the matrix consisting of the first k rows and columns of \(\tilde{\boldsymbol{G}}\). Then, \(|\tilde{\boldsymbol{G}}_{11}|\equiv\gamma_{1\wedge1}=\kappa^{-1}q_{1}\pi_{1}>0\). Suppose now that \(|\tilde{\boldsymbol{G}}_{kk}|\geq 0\) for some k≥1. In order to prove that \(|\tilde{\boldsymbol{G}}_{k+1,k+1}|\geq 0\), multiply the k-th row of \(\tilde{\boldsymbol{G}}_{k+1,k+1}\) by π _k+1/π _k and subtract it from the (k+1)-th. Then the resulting matrix has all last row’s elements equal to zero except from the last one which is \(\kappa^{-1}q_{k+1}\pi_{k+1}-\kappa^{-1}q_{k}\pi _{k+1}^{2}/\pi _{k}=\kappa^{-1}\pi_{k+1}^{2}(q_{k+1}/\pi_{k+1}-q_{k}/\pi_{k})\). By (11) this quantity is nonnegative thus, \(|\tilde {\boldsymbol{G}}_{k+1,k+1}|= \kappa^{-1}\pi_{k+1}^{2}(q_{k+1}/\pi_{k+1}- q_{k}/\pi _{k})|\tilde{\boldsymbol{G}}_{kk}|\geq 0\).

Reprints and permissions