arXiv:1509.06310v2 [math.ST] 11 Aug 2016

Estimating standard errors for importance sampling
estimators with multiple Markov chains
Vivekananda Roy1 , Aixin Tan2 , and James M. Flegal3
1

2

Department of Statistics, Iowa State University
Department of Statistics and Actuarial Science, University of Iowa
3
Department of Statistics, University of California, Riverside
Aug 10, 2016
Abstract

The naive importance sampling estimator, based on samples from a single importance
density, can be numerically unstable. Instead, we consider generalized importance sampling
estimators where samples from more than one probability distribution are combined. We
study this problem in the Markov chain Monte Carlo context, where independent samples
are replaced with Markov chain samples. If the chains converge to their respective target
distributions at a polynomial rate, then under two finite moment conditions, we show a central
limit theorem holds for the generalized estimators. Further, we develop an easy to implement
method to calculate valid asymptotic standard errors based on batch means. We also provide a
batch means estimator for calculating asymptotically valid standard errors of Geyer’s (1994)
reverse logistic estimator. We illustrate the method via three examples. In particular, the
generalized importance sampling estimator is used for Bayesian spatial modeling of binary
data and to perform empirical Bayes variable selection where the batch means estimator
enables standard error calculations in high-dimensional settings.

Key words and phrases: Bayes factors, Markov chain Monte Carlo, polynomial ergodicity,
ratios of normalizing constants, reverse logistic estimator.

1

Introduction

Let π(x) = ν(x)/m be a probability density function (pdf) on X with respect to a measure µ(·).
R
Suppose f : X → R is a π integrable function and we want to estimate Eπ f := X f (x)π(x)µ(dx).
1

Let π1 (x) = ν1 (x)/m1 be another pdf on X such that {x : π1 (x) = 0} ⊂ {x : π(x) = 0}. The
importance sampling (IS) estimator of Eπ f based on independent and identically distributed (iid)
samples X1 , . . . , Xn from the importance density π1 is
Pn
Z
Z
f (Xi )ν(Xi )/ν1 (Xi ) a.s.
f (x)ν(x)/m
ν(x)/m
i=1
Pn
−→
π1 (x) µ(dx)
π1 (x) µ(dx) = Eπ f,
X ν1 (x)/m1
X ν1 (x)/m1
i=1 ν(Xi )/ν1 (Xi )
(1.1)
as n → ∞. This estimator can also be used in the Markov chain Monte Carlo (MCMC) context
when X1 , . . . , Xn are realizations from a suitably irreducible Markov chain with stationary density π1 (Hastings (1970)). Note that (1.1) requires the functions ν and ν1 to be known. On the
other hand, it does not depend on normalizing constants m and m1 , which are generally unknown.
In this article, we consider situations where one wants to estimate Eπ f for all π belonging
to a large collection, say Π. This situation arises in both frequentist and Bayesian statistics.
Although (1.1) provides consistent estimators of Eπ f for all π ∈ Π based on a single Markov
chain {Xn }n≥0 with stationary density π1 , it does not work well when π differs greatly from π1 .
In that case the ratios ν(x)/ν1 (x) can be arbitrarily large for some sample values making the
estimator at (1.1) unstable. In general, there is not a single good importance density π1 which
is close to all π ∈ Π (see e.g. Geyer (1994)). Hence a natural modification is to replace π1 in
(1.1) with a mixture of densities where each density in Π is close to a subset of the k reference
Pk
densities. To this end, denote π ≡
i=1 (ai /|a|)πi , where a = (a1 , . . . , ak ) are k positive
Pk
constants, |a| = i=1 ai , and πi (x) = νi (x)/mi for i = 1, . . . , k are k densities known up to their
normalizing constants. Suppose further that n1 , . . . , nk are positive integers and di := mi /m1 for
i = 2, . . . , k, with d1 ≡ 1. Then define the (k − 1) dimensional vector
d = (m2 /m1 , . . . , mk /m1 ).

(1.2)

(l)

l
be an iid sample from πl or realizations from a positive
Finally for l = 1, . . . , k, let {Xi }ni=1
Harris Markov chain with invariant density πl (for definitions see Meyn and Tweedie (1993)).
Then as nl → ∞, for all l = 1, . . . , k, we have
!, k
!
nl
nl
k
(l)
(l)
(l)
X
X al X
al X
f (Xi )ν(Xi )
ν(Xi )
η̂ ≡
(1.3)
P
P
nl i=1 k as νs (Xi(l) )/ds
nl i=1 k as νs (Xi(l) )/ds
l=1
l=1
s=1
s=1
!, k
!
Z
k
X
X Z
ν(x)
ν(x)
a.s.
−→
al f (x) Pk
πl (x) µ(dx)
al
πl (x) µ(dx)
Pk
a
ν
(x)/d
a
ν
(x)/d
X
X
s
s
s
s
s
s
s=1
s=1
l=1
l=1
Z
Z
ν(x)
ν(x)
π̄(x) µ(dx)
π̄(x) µ(dx) = Eπ f.
= f (x)
π̄(x)
X
X π̄(x)

The generalized IS estimator (1.3) has been discussed widely in the literature, e.g. applications include Monte Carlo maximum likelihood estimation and Bayesian sensitivity analysis. Gill
2

et al. (1988), Kong et al. (2003), Meng and Wong (1996), Tan (2004), and Vardi (1985) consider
estimation using (1.3) based on iid samples. The estimator is applicable to a much larger class
of problems if Markov chain samples are allowed, see e.g. Buta and Doss (2011), Geyer (1994),
and Tan et al. (2015), which is the setting of this paper.
Alternative importance weights have also been proposed. In the case when the normalizing
constants mi ’s are known, the estimator (1.3) resembles the balance heuristic estimator of Veach
and Guibas (1995), which is revisited in Owen and Zhou (2000) as the deterministic mixture.
The standard population Monte Carlo algorithm of Cappe et al. (2004) uses a weighted ratio of
the target π and the proposal πj it was drawn from (evaluated at the sample itself). However,
when iid samples are available from πj , j = 1, 2, . . . , k, Elvira et al. (2015) have shown that
the normalized estimator (mi ’s known) version of (1.3) always has a smaller variance than that
of the population Monte Carlo algorithm. Further, it may be difficult in practice to find fully
known importance densities that approximate the target densities. Indeed, applications such as in
empirical Bayes analysis and Bayesian sensitivity analysis routinely select representatives from
the large number of target posterior densities to serve as proposal densities, and they are known
only up to normalizing constants. See Buta and Doss (2011), Doss (2010), as well as section 5
for examples. Although there is no known proof for the self normalized estimator (Elvira et al.,
2015, p. 16), it is reasonable to assume the superiority of (1.3) over estimators corresponding to
other weighting schemes.
As noted in (1.3), the estimator η̂ converges to Eπ f as the sample sizes increase to infinity,
for iid samples as well as Markov chain samples satisfying the usual regularity conditions. Now
given samples of finite size in practice, it is of fundamental importance to provide some measure
of uncertainty, such as the standard errors (SEs) associated with this consistent estimator. Take
estimators that are sample averages based on iid Monte Carlo samples for example, it is a basic
requirement to report their SEs. But the very same issue is often overlooked in practice when the
estimators have more complicated structure, and when they are based on MCMC samples, largely
due to the difficulty of doing so. See, for e.g. Flegal et al. (2008) on the issue concerning MCMC
experiments and Koehler et al. (2009) for more general simulation studies. For calculating SEs
of η̂ based on MCMC samples, Tan et al. (2015) provide a solution using the method of regenerative simulation (RS). However, this method crucially depends on the construction of a practical
minorization condition, i.e. one where sufficient regenerations are observed in finite simulations
(for definitions and a description of RS see Mykland et al. (1995)). Further, the usual method of
identifying regeneration times by splitting becomes impractical for high-dimensional problems
(Gilks et al., 1998). Hence, successful applications of RS involve significant trial and error and
are usually limited to low-dimensional Gibbs samplers (see e.g. Roy and Hobert (2007); Tan and
3

Hobert (2009)). In this paper we avoid RS and provide SE estimators of η̂ using the batch means
(BM) method, which is straightforward to implement and can be routinely applied in practice.
In obtaining this estimator, we also establish a central limit theorem (CLT) for η̂ that generalizes
some results in Buta and Doss (2011).
The estimator η̂ in (1.3) depends on the ratios of normalizing constants, d, which are unknown
in practical applications. We consider the two-stage scheme studied in Buta and Doss (2011)
where first an estimate dˆ is obtained using Geyer’s (1994) “reverse logistic regression” method
based on samples from πl , and then independently, new samples are used to estimate Eπ f for π ∈
ˆ in (1.3). Buta and Doss (2011) showed that the asymptotic variance
Π using the estimator η̂(d)
ˆ depends on the asymptotic variance of d.
ˆ Thus we study the CLT of dˆ and provide a
of η̂(d)
ˆ Since dˆ involves multiple Markov chain
BM estimator of the asymptotic covariance matrix of d.
samples, we utilize a multivariate BM estimator. Although, the form of the asymptotic covariance
matrix of dˆ is complicated, our consistent BM estimator is straightforward to code.
The problem of estimating d, the ratios of normalizing constants of unnormalized densities
is important in its own right and has many applications in frequentist and Bayesian inference.
For example, when the samples are iid sequences this is the biased sampling problem studied in
Vardi (1985). In addition, the problem arises naturally in the calculations of likelihood ratios in
missing data (or latent variable) models, mixture densities for use in IS, and Bayes factors.
Our work considers the problem of estimating d using Geyer’s (1994) reverse logistic regression method. Specifically, we study the general quasi-likelihood function proposed in Doss
and Tan (2014). Unlike Geyer’s (1994) method, this extended quasi-likelihood function has the
advantage of using user defined weights which are appropriate in situations where the multiple
Markov chains have different mixing rates. We establish the CLT for the resulting estimators of
d and develop the BM estimators of their asymptotic covariance matrix.
Thus we consider two related problems in this paper–firstly, estimating (ratios of) normalizing
constants given samples from k densities, and secondly, estimating expectations with respect to a
large number of (other) target distributions using these samples. In both cases, we establish CLTs
for our estimators and provide easy to calculate SEs using BM methods.
Prior results of Buta and Doss (2011), Doss and Tan (2014), Geyer (1994), and Tan et al.
(2015) all assume that the underlying Markov chains are geometrically ergodic. We weaken this
condition substantially in that we only require the chains to be polynomial ergodic. To this end,
(l)
let Kl (x, ·) be the Markov transition function for the Markov chain Φl = {Xt }t≥1 , so that for

(l)
(l)
any measurable set A, and s, t ∈ {1, 2, . . .} we have P Xs+t ∈ A | Xs = x = Klt (x, A).
Let k · k denote the total variation norm and Πl be the probability measure corresponding to the
density πl . The Markov chain Φl is polynomially ergodic of order m where m > 0 if there exists
4

W : X → R+ with Eπl W < ∞ such that
kKlt (x, ·) − Πl (·)k ≤ W (x)t−m .
There is substantial MCMC literature establishing that Markov chains are at least polynomially
ergodic (see Vats et al. (2015b) and the references therein).
We illustrate the generalized IS method and importance of obtaining SEs through three examples. First, we consider a toy example to demonstrate that BM and RS estimators are consistent
and investigate the benefit of allowing general weights to be used in generalized IS. Second, we
consider a Bayesian spatial model for a root rot disease dataset where we illustrate the importance
of calculating SEs by considering different designs and performing samples size calculations. Finally, we consider a standard linear regression model with a large number of variables and use
the BM estimator developed here for empirical Bayes variable selection.
The rest of the paper is organized as follows. Section 2 is devoted to the important problem
of estimating the ratios of normalizing constants of unnormalized densities, that is estimating
d. Section 3 contains the construction of a CLT for η̂ and describes how valid SEs of η̂ can be
obtained using BM. Sections 4 and Appendix E contain toy examples illustrating the benefits of
different weight functions. Section 5 considers a Bayesian spatial models for binary responses.
The empirical Bayes variable selection example is contained in Appendix F. We conclude with a
discussion in Section 6. All proofs are relegated to the appendices.

2

Estimating ratios of normalizing constants

Consider k densities πl = νl /ml , l = 1, . . . , k with respect to the measure µ, where the νl ’s
are known functions and the ml ’s are unknown constants. For each l we have a positive Harris
(l)
(l)
Markov chain Φl = {X1 , . . . , Xnl } with invariant density πl . Our objective is to estimate all
possible ratios mi /mj , i 6= j or, equivalently, the vector d defined in (1.2).
P
Geyer’s (1994) reverse logistic regression is described as follows. Let n =
nl and set
al = nl /n for now. For l = 1, . . . , k define the vector ζ by
ζl = − log(ml ) + log(al )
and let

νl (x)eζl
pl (x, ζ) = Pk
.
(2.1)
ζs
s=1 νs (x)e
 (l)
Given the value x belongs to the pooled sample Xi , i = 1, . . . , nl , l = 1, . . . , k , pl (x, ζ) is
the probability that x came from the lth distribution. Of course, we know which distribution the
5

sample x came from, but here we pretend that the only thing we know about x is its value and
estimate ζ by maximizing the log quasi-likelihood function
ln (ζ) =

nl
k X
X

(l)

log pl (Xi , ζ)



(2.2)

l=1 i=1

with respect to ζ. Since ζ has a one-to-one correspondence with m = (m1 , . . . , mk ), by estimating ζ we can estimate m.
As Geyer (1994) mentioned, there is a non-identifiability issue regarding ln (ζ): for any constant c ∈ R, ln (ζ) is same as ln (ζ+c1k ) where 1k is the vector of k 1’s. So we can estimate the true
ζ only up to an additive constant. Thus, we can estimate m only up to an overall multiplicative

Pk
constant, that is, we can estimate only d. Let ζ 0 ∈ Rk be defined by [ζ 0 ]l = [ζ]l −
s=1 [ζ]s /k,
the true ζ normalized to add to zero. Geyer (1994) proposed to estimate ζ 0 by ζ̂, the maximizer
of ln subject to the linear constraint ζ > 1k = 0, and thus obtain an estimate of d. The estimator dˆ (written explicitly in Section 2.1), was introduced by Vardi (1985), and studied further by
Gill et al. (1988), who proved that in the iid setting, dˆ is consistent and asymptotically normal,
and established its efficiency. Geyer (1994) proved the consistency and asymptotic normality of
dˆ when Φl , . . . , Φk are k Markov chains satisfying certain mixing conditions. In the iid setting,
Meng and Wong (1996), Kong et al. (2003), and Tan (2004) rederived the estimate under different
computational schemes.
However, none of these articles discuss how to consistently estimate the covariance matrix
ˆ
of d, even in the iid setting. Recently Doss and Tan (2014) address this important issue and
obtain a RS estimator of the covariance matrix of dˆ in the Markov chain setting. Doss and Tan
(2014) also mention the optimality results of Gill et al. (1988) do not hold in the Markov chain
case. In particular, when using Markov chain samples, the choice of the weights aj = nj /n to
the probability density νj /mj in the denominator of (2.1) is no more optimal and should instead
incorporate the effective sample size of different chains as they might have quite different rates
of mixing. They introduce the following more general log quasi-likelihood function
`n (ζ) =

k
X
l=1

wl

nl
X


(l)
log pl (Xi , ζ) ,

(2.3)

i=1

where the vector w ∈ Rk is defined by wl = al n/nl for l = 1, . . . , k for an arbitrary probability
vector a. (Note the change of notation from l to `.) Clearly if al = nl /n, then wl = 1 and (2.3)
becomes (2.2).
When RS can be used, Doss and Tan (2014) proved the consistency (to the true value ζ 0 )
and asymptotic normality of the constrained maximizer ζ̂ (subject to the constraint ζ > 1k = 0) of
6

(2.3) under geometric ergodicity. They also obtain a RS estimator of the asymptotic covariance
matrix and describe an empirical method for choosing the optimal a based on minimizing the
ˆ However, their procedure requires a practical mitrace of the estimated covariance matrix of d.
norization condition for each of the k Markov chains, which can be extremely difficult. Without
a minorization condition, we show dˆ is a consistent estimator of d, show dˆ satisfies a CLT under
significantly weaker mixing conditions, and provide a strongly consistent BM estimator of the
ˆ
covariance matrix of d.

2.1

Central limit theorem and asymptotic covariance estimator

Within each Markov chain l = 1, . . . , k, assume nl → ∞ in such a way that nl /n → sl ∈ (0, 1).
ˆ we first establish a CLT for ζ̂. Note that the function
In order to obtain the CLT result for d,
g : Rk → Rk−1 that maps ζ 0 into d is given by


eζ1 −ζ2 a2 /a1

 ζ1 −ζ3
e
a3 /a1 
,

(2.4)
g(ζ) = 
..

.


eζ1 −ζk ak /a1

and its gradient at ζ 0 (in terms of d) is


d2
d3

−d2 0

D=
 0 −d3
..
 ..
 .
.
0
0

...
...
...
...

dk
0
0
..
.





.




(2.5)

. . . −dk

ˆ
Since d = g(ζ 0 ), and by definition dˆ = g(ζ̂), we can use the CLT result of ζ̂ to get a CLT for d.
First, we introduce the following notations. For r = 1, . . . , k, let

(r,l)
(l)
Yi
= pr (Xi , ζ 0 ) − Eπl pr (X, ζ 0 ) ,
i = 1, . . . , nl .
(2.6)
The asymptotic covariance matrix in the CLT of ζ̂, involves two k × k matrices B and Ω, which
we now define. The matrix B is given by
Brr =

k
X


aj Eπj pr (X, ζ)[1 − pr (X, ζ)] and

j=1

Brs = −

k
X

(2.7)


aj Eπj pr (X, ζ)ps (X, ζ) for r 6= s.

j=1

7

Let Ω be the k × k matrix defined (for r, s = 1, . . . , k) by
∞
∞
k
i
X
X
X
a2l h
(r,l) (s,l)
(r,l) (s,l)
(r,l) (s,l)
Eπl {Y1 Y1 } +
Eπl {Y1 Y1+i } +
Eπl {Y1+i Y1 } .
Ωrs =
sl
i=1
i=1
l=1
(r,l)

(s,l)

(2.8)
(r,l)

(s,l)

Remark 1. The right hand side of (2.8) involves terms of the form Eπl {Y1 Y1+i } and Eπl {Y1+i Y1
(l)
(l)
For any fixed l, r, s and i, the two expectations are the same if X1 and X1+i are exchangeable,
e.g. if the chain Φl is reversible. In general, the two expectations are not equal.
b defined by
The matrix B will be estimated by its natural estimate B
 X

nl
k
X


1
(l)
(l)
brr =
B
al
pr (Xi , ζ̂) 1 − pr (Xi , ζ̂) and
nl i=1
l=1

 X
n
k
X
1 l
(l)
(l)
b
pr (Xi , ζ̂)ps (Xi , ζ̂) for r 6= s.
Brs = −
al
nl i=1
l=1

(2.9)

b suppose we simulate the Markov chain Φl for nl = el bl iterations
To obtain a BM estimate Ω,
(hence el = enl and bl = bnl are functions of nl ) and define for r, l = 1, . . . , k
(r,l)
Z̄m

1
:=
bl

(m+1)bl

X

(l)

pr (Xj , ζ̂)

for m = 0, . . . , el − 1 .

j=mbl +1

>

(l)
(1,l)
(k,l)
Now set Z̄m = Z̄m , . . . , Z̄m
for m = 0, . . . , el − 1. For l = 1, . . . , k, denote Z̄¯ (l) =
>

P l
(l)
¯
¯
(1,l)
(k,l)
pr (Xi , ζ̂)/nl . Let
where Z̄¯ (r,l) = ni=1
Z̄ , . . . , Z̄
b (l) =
Σ
Let

el −1 h
ih
iT
bl X
(l)
(l)
(l)
(l)
¯
¯
Z̄ − Z̄
Z̄m − Z̄
for l = 1, . . . , k.
el − 1 m=0 m


b=
Σ


b (1)
Σ
..

....



0
.....
b (k)
Σ

0
and define the following k × k 2 matrix
 r
r
n
n
An = −
a1 Ik −
a2 Ik
n1
n2

(2.10)

r
...

(2.11)




−

n
ak Ik
nk


,

(2.12)

where Ik denotes the k × k identity matrix. Finally, define
b = An ΣA
b >
Ω
n.
8

(2.13)

}.

We are now ready to describe conditions that ensure strong consistency and asymptotic norˆ The following theorem also provides consistent estimate of the asymptotic covarimality of d.
ance matrix of dˆ using BM method. Consistency of dˆ holds under minimal assumptions, i.e. if
Φ1 , . . . , Φk are positive Harris chains. On the other hand, CLTs and consistency of BM estimator
of asymptotic covariance require some mixing conditions on the Markov chains. For a square
matrix C, let C † denote the Moore-Penrose inverse of C.
(l)

(l)

Theorem 1 Suppose that for each l = 1, . . . , k, the Markov chain {X1 , X2 , . . .} has invariant
distribution πl .
1. If the Markov chains Φ1 , . . . , Φk are positive Harris, the log quasi-likelihood function (2.3)
has a unique maximizer subject to the constraint ζ > 1k = 0. Let ζ̂ denote this maximizer,
a.s.
and let dˆ = g(ζ̂). Then dˆ −→ d as n1 , . . . , nk → ∞.
2. If the Markov chains Φ1 , . . . , Φk are polynomially ergodic of order m > 1, as n1 , . . . , nk →
√
d
∞, n(dˆ − d) → N (0, V ) where V = D> B † ΩB † D.
3. Assume that the Markov chains Φ1 , . . . , Φk are polynomially ergodic of order m > 1 and
b be the matrix D in (2.5)
for all l = 1, . . . , k, bl = bnνl c where 1 > ν > 0. Let D
b and Ω
b be defined by (2.9) and (2.13), respectively. Then,
with dˆ in place of d, and let B
b is a strongly consistent estimator of V .
bB
b†D
b†Ω
b >B
Vb := D

3

IS with multiple Markov chains

This section considers a CLT and SEs for the generalized IS estimator η̂. From (1.3) we see that
η̂ ≡ η̂ [f ] (π; a, d) = v̂ [f ] (π, π1 ; a, d)/û(π, π1 ; a, d), where
nl
k
X
al X
(l)
v̂ ≡ v̂ (π, π1 ; a, d) :=
v [f ] (Xi ; a, d) and
nl i=1
l=1
[f ]

(3.1)

nl
k
X
al X
(l)
u(Xi ; a, d)
û ≡ û(π, π1 ; a, d) :=
n
l i=1
l=1

with
v [f ] (x; a, d) := f (x)u(x; a, d)

and u(x; a, d) := Pk

s=1

ν(x)
as νs (x)/ds

.

Note that û converges almost surely to
Pk
Z
k
X
al νl (x)/ml
m
al Eπl u(X; a, d) =
ν(x) µ(dx) =
,
Pk l=1
m1
X
s=1 as νs (x)/(ms /m1 )
l=1
9

(3.2)

(3.3)

as n1 , . . . , nk → ∞. Thus û itself is a useful quantity as it consistently estimates the ratios of
normalizing constants {u(π, π1 ) ≡ m/m1 |π ∈ Π}. Unlike the estimator dˆ in Section 2, û does
not require a sample from each density π ∈ Π. Thus û is well suited for situations where one
wants to estimate the ratios u(π, π1 ) for a very large number of π’s based on samples from a small
number of skeleton densities, say k. Thus this method is particularly efficient when obtaining
samples from the target distributions is computationally demanding and the distributions within
Π are similar.
In the context of Bayesian analysis, let π(x) = lik(x)p(x)/m be the posterior density corresponding to the likelihood function lik(x) and prior p(x) with normalizing constant m. In this
case, u(π, π1 ) is the so-called Bayes factor between the two models, which is commonly used in
model selection.
The estimators û and v̂ in (3.1) depend on d, which is generally unknown in practice. Here
we consider a two-stage procedure for evaluating û. In the 1st stage, d is estimated by its reverse
l
logistic regression estimator dˆ described in Section 2 using Markov chains Φ̃l ≡ {X˜il }N
i=1 with
stationary densities πl , for l = 1, . . . , k. Note the change of notation from Section 2 where we
used nl ’s to denote the length of the Markov chains. In order to avoid more notations, we use
Φ̃l ’s and Nl ’s to denote the stage 1 chains and their length respectively. Once dˆ is formed, new
l
MCMC samples Φl ≡ {Xil }ni=1
, l = 1 . . . , k are obtained and u(π, π1 )(Eπ f ) is estimated using
[f ]
ˆ
ˆ
û(π, π1 ; a, d) (η̂ (π; a, d)) based on these 2nd stage samples. Buta and Doss (2011) propose
this two-stage method and quantify its benefits over the method where the same MCMC samples
are used to estimate both d and u(π, π1 ).

3.1

Estimating ratios of normalizing constants

ˆ we require some notation. Let
Before we state a CLT for û(π, π1 ; a, d),
τl2 (π; a, d)

=

(l)
Varπl (u(X1 ; a, d))

+2

∞
X

(l)

(l)

Covπl (u(X1 ; a, d), u(X1+g ; a, d))

(3.4)

g=1

P
and τ 2 (π; a, d) = kl=1 (a2l /sl )τl2 (π; a, d). Further, define c(π; a, d) as a vector of length k − 1
with (j − 1)th coordinate as
Z
aj νj (x)
u(π, π1 )
π(x)dx for j = 2, . . . , k
(3.5)
[c(π; a, d)]j−1 =
Pk
2
dj
a
ν
(x)/d
X
s
s
s
s=1
and ĉ(π; a, d) as a vector of length k − 1 with (j − 1)th coordinate as
[ĉ(π; a, d)]j−1

nl
k
(l)
(l)
X
1 X
aj al ν(Xi )νj (Xi )
=
for j = 2, . . . , k.
P
nl i=1 ( k as νs (Xi(l) )/ds )2 d2j
l=1
s=1

10

(3.6)

Assuming nl = el bl , then let
τ̂l2 (π; a, d) =

el −1
bl X
[ūm (a, d) − ū¯(a, d)]2 ,
el − 1 m=0
(l)

(3.7)
(l)

where ūm (a, d) is the average of the (m + 1)st block {u(Xmbl +1 ; a, d), · · · , u(X(m+1)bl ; a, d)},
(l)
(l)
and ū¯(a, d) is the overall average of {u(X1 ; a, d), · · · , u(Xnl ; a, d)}. Here, bl and el are the
P
block sizes and the number of blocks respectively. Finally let τ̂ 2 (π; a, d) = kl=1 (a2l /sl )τ̂l2 (π; a, d).
ˆ
Theorem 2 Suppose that for the stage 1 chains, conditions of Theorem 1 holds such that N 1/2 (d−
P
d
d) → N (0, V ) as N ≡ kl=1 Nl → ∞. Suppose there exists q ∈ [0, ∞) such that n/N → q
Pk
where n =
l=1 nl is the total sample size for stage 2. In addition, let nl /n → sl for l =
1, · · · , k.
1. Assume that the stage 2 Markov chains Φ1 , . . . , Φk are polynomially ergodic of order m,
and for some δ > 0 Eπl |u(X; a, d)|2+δ < ∞ for each l = 1, · · · , k where m > 1 + 2/δ.
Then as n1 , . . . , nk → ∞,
√
d
ˆ − u(π, π1 )) →
n(û(π, π1 ; a, d)
N (0, qc(π; a, d)> V c(π; a, d) + τ 2 (π; a, d)). (3.8)
2. Let Vb be the consistent estimator of V given in Theorem 1 (3). Assume that the Markov
chains Φ1 , . . . , Φk are polynomially ergodic of order m ≥ (1 + )(1 + 2/δ) for some
, δ > 0 such that Eπl |u(X; a, d)|4+δ < ∞, and for all l = 1, . . . , k, bl = bnνl c where
ˆ + τ̂ 2 (π; a, d))
ˆ is a strongly consistent estimator
ˆ > Vb ĉ(π; a, d)
1 > ν > 0. Then qĉ(π; a, d)
of the asymptotic variance in (3.8).
Note that the asymptotic variance in (3.8) has two components. The second term is the variance of û when d is known. The first term is the increase in the variance of û resulting from
using dˆ instead of d. Since we are interested in estimating u(π, π1 ) for a large number of π’s and
for every π the computational time needed to calculate û in (3.1) is linear in the total sample size
n, this can not be very large. If generating MCMC samples is not computationally demanding,
then long chains can be used in the 1st stage (that is, large Nl ’s can be used) to obtain a precise
estimate of d, and thus greatly reducing the first term in the variance expression (3.8).

3.2

Estimation of expectations using generalized IS

This section discusses estimating SEs of the generalized IS estimator η̂ given in (1.3). In order to
state a CLT for η̂ we define the following notations:
∞
X
(l)
(l)
(l)
γl11 ≡ γl11 (π; a, d) = Varπl (v [f ] (X1 ; a, d)) + 2
Covπl (v [f ] (X1 ; a, d), v [f ] (X1+g ; a, d)),
g=1

11

(l)

(l)

γl12 ≡ γl12 (π; a, d) = γl21 ≡ γl21 (π; a, d) = Covπl (v [f ] (X1 ; a, d), u(X1 ; a, d))
∞
X
(l)
(l)
(l)
(l)
+
[Covπl (v [f ] (X1 ; a, d), u(X1+g ; a, d)) + Covπl (v [f ] (X1+g ; a, d), u(X1 ; a, d))],
g=1

γl22

≡

γl22 (π; a, d)

=

(l)
Varπl (u(X1 ; a, d))

+2

∞
X

(l)

(l)

Covπl (u(X1 ; a, d), u(X1+g ; a, d)),

g=1

(note γl22 is the same as τl2 (π; a, d) defined in (3.4)) and
!
k
X
γ 11 γ 12
a2l
Γl (π; a, d) =
;
Γ(π;
a,
d)
=
Γl (π; a, d).
s
γ 21 γ 22
l
l=1

(3.9)

Since η̂ has the form of a ratio, to establish a CLT for it, we apply the Delta method on the
function h(x, y) = x/y, with ∇h(x, y) = (1/y, −x/y 2 )0 . Let
ρ(π; a, d) = ∇h(Eπ f u(π, π1 ), u(π, π1 ))0 Γ(π; a, d)∇h(Eπ f u(π, π1 ), u(π, π1 )),
e(π; a, d) is a vector of length k − 1 with (j − 1)th coordinate as
Z
aj
[f (x) − Eπ f ]νj (x)
π(x)dx, j = 2, . . . , k,
[e(π; a, d)]j−1 = 2
Pk
dj X
s=1 as νs (x)/ds

(3.10)

(3.11)

and ê(π; a, d) is a vector of length k − 1 with (j − 1)th coordinate as
Pk
[ê(π; a, d)]j−1 ≡

al
l=1 nl

(l)

(l)

(l)

aj f (Xi )ν(Xi )νj (Xi )
i=1 d2 (Pk as νs (X (l) )/ds )2
s=1
j
i

Pnl

û(π, π1 ; a, d)

−

[c(π; a, d)]j−1 η̂ [f ] (π; a, d)
,
û(π, π1 ; a, d)

(3.12)

where [c(π; a, d)]j−1 is defined in (3.6). Assuming nl = el bl , let
"
!
! #"
!
! #>
eX
[f ]
l −1
¯[f ]
¯[f ]
v̄
v̄
v̄
v̄
b
m
m
l
bl (π; a, d) =
−
−
Γ
el − 1 m=0
ūm
ū¯
ūm
ū¯


h
i
h
i
2
Pel −1 [f ]
Pel −1 [f ]
[f ]
[f ]
¯
¯
¯
v̄
−
v̄
v̄
−
v̄
[ū
−
ū
]
m
m
bl 
m
m=0
m=0

h
i
=
P
P
[f
]
e
−1
e
−1
2
l
l
[f ]
el − 1
¯
¯
¯
[ūm − ū]
m=0 v̄m − v̄
m=0 [ūm − ū]
!
γ̂ 11 (π; a, d) γ̂ 12 (π; a, d)
=
,
γ̂ 21 (π; a, d) γ̂ 22 (π; a, d)
[f ]

(l)

(l)

where v̄m is the average of the (m + 1)st block {v [f ] (Xmbl +1 ; a, d), · · · , v [f ] (X(m+1)bl ; a, d)},
(l)
(l)
v̄¯[f ] is the overall average of {v [f ] (X1 ; a, d), · · · , v [f ] (Xnl ; a, d)} and ūm ≡ ūm (π, a, d), ū¯ ≡
b a, d) = Pk (a2 /sl )Γ
bl (π; a, d), and
ū¯(π, a, d) defined in Section 3.1. Finally let Γ(π;
l=1 l
[f ] ˆ
b a, d)∇h(v̂
ˆ = ∇h(v̂ [f ] (d),
ˆ û(d))
ˆ 0 Γ(π;
ˆ
ˆ
ρ̂(π; a, d)
(d), û(d)).

12

ˆ
Theorem 3 Suppose that for the stage 1 chains, conditions of Theorem 1 hold such that N 1/2 (d−
P
d
d) → N (0, V ) as N ≡ kl=1 Nl → ∞. Suppose there exists q ∈ [0, ∞) such that n/N → q
Pk
where n =
l=1 nl is the total sample size for stage 2. In addition, let nl /n → sl for l =
1, · · · , k.
1. Assume that the stage 2 Markov chains Φ1 , . . . , Φk are polynomially ergodic of order m,
and for some δ > 0 Eπl |u(X; a, d)|2+δ < ∞ and Eπl |v [f ] (X; a, d)|2+δ < ∞, for each
l = 1, · · · , k where m > 1 + 2/δ. Then as n1 , . . . , nk → ∞,
√

d
ˆ − Eπ f ) →
n(η̂ [f ] (π; a, d)
N (0, qe(π; a, d)> V e(π; a, d) + ρ(π; a, d)).

(3.13)

2. Let Vb be the consistent estimator of V given in Theorem 1 (3). Assume that the Markov
chains Φ1 , . . . , Φk are polynomially ergodic of order m ≥ (1+)(1+2/δ) for some , δ > 0
such that Eπl |u(X; a, d)|4+δ < ∞, Eπl |v [f ] (X; a, d)|4+δ < ∞, and for each l = 1, · · · , k,
ˆ + ρ̂(π; a, d))
ˆ is a strongly
ˆ > Vb ê(π; a, d)
bl = bnνl c where 1 > ν > 0. Then qê(π; a, d)
consistent estimator of the asymptotic variance in (3.13).
Remark 2. Part (1) of Theorems 2 and 3 extend Buta and Doss’s (2011) Theorems 1 and 3,
respectively. Specifically, they require al = nl /n which is a non-optimal choice for a (Tan et al.,
2015). Our results also substantially weaken the Markov chain mixing conditions.
Remark 3. Theorems 2 and 3 prove consistency of the BM estimators of the variances of û and
η̂ for a general a. This extends results in Tan et al. (2015), which provides RS based estimators
ˆ With this particular
of the asymptotic variance of û and η̂ in the special case when a = (1, d).
ˆ and v [f ] (x; a, d)
ˆ in (3.2) become free of dˆ leading to independence among
choice, u(x; a, d)
ˆ for any user specified fixed vector w,
certain quantities. However, one can set a = w ∗ (1, d)
which allows the expression in (2.18) of Tan et al. (2015) to be free of dˆ and thus the necessary
independence. Hence, their RS estimator can also be applied to an arbitrary vector a (details are
given in Appendix D).
Remark 4. A sufficient
for the moment
assumptions for u in Theorems 2 and 3 is that,
n condition
o
 Pk
for any π ∈ Π, supx π(x)
s=1 as πs (x) < ∞. That is, in any given direction, the tail of at
least one of {πs , s = 1, . . . , k} is heavier than that of π. This is not hard to achieve in practice by
properly choosing {πs } with regard to Π (see e.g. Roy , 2014). Further, if Eπ |f |4+δ < ∞ then
the moment assumptions for v [f ] are satisfied.

13

4

Toy example

In this section, we employee a toy example to confirm that both the BM and the RS estimators
are consistent, as well as demonstrate the benefit of allowing general weights to be used in the
generalized IS estimator. Let tr,µ denote the t-distribution with degree of freedom r and central
parameter µ. We set π1 (·) = ν1 (·) and π2 (·) = ν2 (·), which are the density functions for a
t5,µ1 =1 and t5,µ2 =0 , respectively. Pretending that we do not know the value of the ratio between
the two normalizing constants, d = m2 /m1 = 1/1, we estimate it by the stage 1 estimator dˆ from
section 2, and compare the BM and the RS method in estimating the asymptotic variance. As for
the stage 2 estimators from section 3, the choice of weight and performance of the BM and the
RS methods in assessing estimators’ uncertainty are studied in Appendix E.
We draw iid samples from π1 and Markov chain samples from π2 using the independent
(x)
>
Metropolis Hastings algorithm with proposal density t5,1 . It is simple to show inf x tt5,µ
5,0 (x)
0, which implies the algorithm is uniformly ergodic (Mengersen and Tweedie (1996) Theorem
2.1) and hence polynomially ergodic and geometrically ergodic. For RS, our carefully tuned
minorization condition enables the Markov chain for π2 to regenerate about every 3 iterations. In
contrast, the BM method proposed here requires no such theoretical development.
We evaluated the variance estimators at various sample sizes with different choices of weight.
ˆ in
Figure 1 displays traces of the BM and the RS estimates of the asymptotic variance of d,
dashed and solid lines, respectively. Overall, the BM and the RS estimates approach the empirical asymptotic variance as the sample size increases, suggesting their consistency. Due to the
frequency of regenerations, BM estimates are generally more variable than RS estimates. Further, the left panel of Figure 1 is for estimators based on the naive weight, a = (0.5, 0.5), that is
proportional to the sample sizes; and the right panel is for estimators based on a = (0.82, 0.18),
that emphasizes the iid sample more than the Markov chain sample. Indeed, the latter weight is
a close-to-optimal weight obtained with a small pilot study (see Appendix E for details). Using
such a method to choose weight can lead to big improvement in the efficiency of dˆ if the mixing
rate of the multiple samples differ a lot.

5

Bayesian spatial models for binary responses

In this section, we analyze a root rot disease dataset collected from a 90-acre farm in the state
of Washington (Zhang, 2002). All computations are done in R, using the package geoBayes
(Evangelou and Roy, 2015). Recorded at M = 100 chosen sites are the longitude and the latitude
si , the root counts `(si ), and the number of infected roots y(si ), i = 1, . . . , M . Of interest is
14

0.12

0.16

BM
RS
empirical

0.04

0.08

estimates of asy var

0.12
0.08
0.04

estimates of asy var

0.16

BM
RS
empirical

3.0

3.5

4.0

4.5

5.0

3.0

log10 of chain 1 length in stage 1

3.5

4.0

4.5

5.0

log10 of chain 1 length in stage 1

Figure 1: Plots of BM and RS estimates of the asymptotic variance of dˆ in stage 1 for 100
randomly chosen replications. The left panel is based on the naive weight, a[1] = (0.5, 0.5)
and the right panel is based on a close-to-optimal weight, a[1] = (0.82, 0.18). Horizontal lines
represent the empirical asymptotic variance of dˆ obtained over all replications.
a map of the disease rate over the entire area for precision farming. We consider the following
spatial generalized linear mixed model (SGLMM), similar to that used by Zhang (2002) and Roy
et al. (2016). Taking `(si ) and si as fixed, let
ind

y(si )|z(si ) ∼ Binomial (`(si ), Φ (z(si ))) , i = 1, . . . , M.
Here z = (z(s1 ), . . . , z(sM )) is a vector of latent variables, which is assumed to be a subvector
of a Gaussian random field (GRF) {zs , s ∈ S}, that has a constant mean µ, and a covariance
function
Cov (z(s), z(s0 )) = σ 2 ρφ (ks − s0 k) + ω σ 2 Is (s0 ) .
Here, σ 2 is the partial sill, k · k denotes the Euclidean distance, and ρφ is a correlation function
 3
from the spherical family with range parameter φ. That is, ρφ (u) = 1 − 23 φu + 21 φu for
u ∈ (0, φ). Next, Is (s0 ) is an indicator that equals 1 if s = s0 , and equals 0 otherwise. Finally, ωσ 2
is the nugget effect, accounting for any remaining variability at site s such as measurement error,
while ω ∈ R+ is the relative size of the nugget to the partial sill. Following Roy et al. (2016)
we assign a non-informative Normal-inverse-Gamma prior to (µ, σ 2 ) which is (conditionally)
conjugate for the model. Specifically,


 1
1
2
2
2
2 − 2 −1
µ|σ ∼ N(0, 100σ ), and f (σ ) ∝ σ
exp − 2 .
2σ
15

Assigning priors for h = (φ, ω) in the correlation function of the Gaussian random field is usually
difficult, and the choice of prior may influence the inference (Christensen, 2004). Hence we
perform a sensitivity analysis focusing on obtaining the Bayes factor (BF) of the model at h
relative to a baseline h0 for a range of values h ∈ H. Note that for a fixed h = (φ, ω), the
Bayesian model described above has parameters ψ = (µ, σ 2 ). Conditioning on the observed data
y = (y(s1 ), . . . , y(sM )), inference is based on the posterior density
πh (ψ|y) =

Lh (ψ|y)π(ψ)
,
mh (y)

(5.1)

R
where Lh (ψ|y) = RM f (y|z)fh (z|ψ)dz is the likelihood, π(ψ) is the prior on ψ, and mh (y) =
R
L (ψ|y)π(ψ)dψ is the normalizing constant, also called the marginal likelihood. The
R×R+ h
Bayes factor between any two models indexed by h and h0 is mh (y)/mh0 (y), and the empirical
Bayes choice of h is arg max mh (y) = arg max [mh (y)/mh0 (y)]. Our plan is to get MCMC
h∈H

h∈H

samples for a small reference set of h, to estimate the BF among them using the reverse logistic
estimator, and then get new samples to estimate {mh (y)/mh0 (y), h ∈ H} using the generalized
IS method. Below, we describe the MCMC algorithms and the practical concern of how long to
run them, which illustrates the importance of calculating a SE.
The two high-dimensional integrals lend the posterior density in (5.1) intractable. But there
exist MCMC algorithms to sample from the augmented posterior distribution, that is
πh (ψ, z|y) =

f (y|z)fh (z|ψ)π(ψ)
.
mh (y)

(5.2)

R
Note that R πh (ψ, z|y)dz = πh (ψ|y). Hence, a two-component Gibbs sampler that updates ψ
and z in turn from their respective conditional distributions based on (5.2) yields a Markov chain
{ψ (i) , z (i) }i≥1 with stationary distribution πh (ψ, z|y). As a result, the marginal {ψ (i) }i≥1 is also
a Markov chain with stationary distribution πh (ψ|y) (Tanner and Wong, 1987).
As a starting point, we use a small pilot study to identify a range for h = (φ, ω) that corresponds to reasonably large BF values. This step is carried out by obtaining the reverse logistic
estimator of BF at a coarse grid of h values over a wide area, based on short runs of Markov
chains. Specifically, (φ, ω) ∈ [80, 200] × [0.2, 2] and within this range the minimum BF is about
1% the size of the maximum. To more carefully estimate BF over this range, we examine a fine
grid H that consists of 130 different h values, with increments of size 10 for the φ component,
and that of .2 for the w component.
A natural choice for the set of skeleton points is {80, 140, 200}×{0.5, 1, 2}, with an arbitrarily
chosen baseline at (200, 2). We first experiment with samples of sizes n1 = · · · = n9 = 500 at the
16

skeleton points (after a burn-in period of 500 iterations and a thinning procedure that keeps one
sample every 10 iterations), of which the first 80% are used in stage 1, and the remaining in stage 2
of the generalized IS procedure. BF estimates at all h ∈ H are obtained, though not shown.
Given the current Monte Carlo sample sizes, it is natural to consider how trustworthy these BF
estimates are. To this end, Figure 2 shows the point-wise SEs over the entire plot obtained via the
BM method. In this setting for some h, the magnitude of the SE is over 63% of the corresponding
BF estimate. Hence, inference based on these BF estimates could be questioned because of the
high computational uncertainty.
Suppose we wish to reduce the relative SE to 5% or less for all h ∈ H, then roughly
(60%/5%)2 = 144 times as many samples would be required under the current design. Instead, we consider an alternative design using a bigger set of skeleton points, {80, 140, 200} ×
{0.5, 1, 1.5, 2} keeping the baseline unchanged at (200, 0.2). In this alternative design, with sample sizes n1 = · · · = n12 = 500, the largest relative SE is 17%. This is almost one fourth that
of the previous design, but the computing time only increased by 40%. Accordingly, running the
alternative design would achieve the computing goal much faster. In this case, we increase the
sample sizes to n1 = · · · = n12 = 6000, which is approximately (17%/5%)2 times 500. Overall,
the new process takes 2 minutes to run on a 3.4GHz Intel Core i7 running linux. The resulting
BF estimates are shown in Figure 3, with maximum relative SE reduced to 5.0%. For the sake of
comparison, running the original design for the same amount of time allows a common sample
size of ni = 8000 resulting in a maximum relative SE of 14.4%. In short, easily obtainable SE
estimates allows us to experiment, choose among different designs, and perform samples size
calculations.
The simplicity of the method matters when it comes to estimating SEs in practice. Using
the BM method to obtain SE requires no extra input beyond what is needed for obtaining the
generalized IS estimates. Indeed, as long as one can run existing software to obtain the Markov
chain samples, there is no need to know the Markov transition kernels utilized in the background.
Unlike the BM method, the RS method depends on identifying regeneration times, typically
through constructing minorization conditions for the Markov transition kernels (see Mykland
et al. (1995) for details). Despite the fact that minorization conditions can be established for any
Markov transition kernel, we demonstrate for the current example the amount of effort needed
to obtain a regeneration can be prohibitively high. Recall the MCMC scheme involves sampling
from πh (ψ|z, y) and πh (z|ψ, y) in turn. The former is a standard distribution hence easy to sample from. The latter is not, and we followed Diggle et al. (1998) that updates zj , j = 1, · · · , M in
turn, each using a one-dimensional Metropolis-Hastings step that keeps invariant the conditional
posterior distribution of zj given all other components. Denote the transition density of these MH
17

SE of BF estimates
2.0

−4

−6

−4.5

1.5

−3

.5

1.5

−5

−4.5
−5.5

−5

−5
.5

−5.5

−4.5

−4

−5

−5.5

.5

.5

−5

−4

−3

2.0

SE of BF estimates

−4.5
−2.5

−3.5

−3.5

1.0

−3

w

−3.5

−2.5

−3

−3.5

1.0

w

−4

−4

−4.5

−3

80

0.5
100

120

−4

−3

140

160

180

200

80

φ

−4

−4

−4.5

−3.5

100

−5

−4

.5

−4
−3

0.5

−4

−4

−3.5

120

140

160

180

200

φ

Figure 2: Contour plots of the SEs (in log scale) evaluated for the BF estimates. The two plots
are based on the original and alternative designs. Skeleton points used in each design are marked
by crosses, with the baselines marked by boxed crosses.
steps as f1 , · · · , fM , and suppressing the notations of their dependence on y, the transition kernel
of the Markov chain can be represented as
0
0
0 0
p(z 0 , ψ 0 |z, ψ) = f1 (z10 |z2 , · · · , zM , ψ)f2 (z20 |z10 , z3 , · · · , zM , ψ) · · · fn (zM
|z10 , · · · , zM
−1 , ψ)πh (ψ |z ) .

According to a common method described in Jones and Hobert (2004), one can build a minorization condition by finding D ⊂ RM × R × R+ ,  > 0, and k(·) such that,
p(z 0 , ψ 0 |z, ψ) ≥ ID (z, ψ) k(z 0 , ψ 0 ) for all (z 0 , ψ 0 ) ∈ RM × R × R+ .
Further, the above condition can be established if
0
0 0
0
|z10 , · · · , zM
f1 (z10 |z2 , · · · , zM , ψ)f2 (z20 |z10 , z3 , · · · , zM , ψ) · · · fM (zM
−1 , ψ)πh (ψ |z )
0
)πh (ψ 0 |z 0 ) for all (z 0 , ψ 0 ) ∈ RM × R × R+ ,
≥ID (z, ψ) 1 k1 (z10 ) 2 k2 (z10 , z20 ) · · · M kM (z10 , · · · , zM

where the common term πh on both sides of the inequality will cancel, and hence the work
is in finding 1 , · · · , M , and k1 (·), · · · , kM (·). It’s easy to see that the smaller the set D, the
larger  = ΠM
i=1 j can possibly be, where  can be interpreted as the conditional regeneration
rate given D is visited. Suppose we take D to be small enough such that j takes on a very
18

2.0

SE/BF

2.0

BF

−5.5

−6

−7

−7

5

−5.

1
1.5

−7

−6

5

0

.
−6

1.5

0

−6

0.5

w

1.0

0.5

160
phi

0.5
180

0

−0.5

−2

−1

.5

−6

−5.5

140

−7

−6

0.5

120

−5.5

1

0.5

100

−5.5

1.0

2.0
1.5

w

−3
80

1.5

1.0

w
−2

−8

BF

−1

−5

−4

−4.5

−3.5

200
80

100

120

140
φ

160

180

200

80

100

120

140

160

180

200

φ

Figure 3: Based on samples of size 6000 of each of the 16 Markov chains, the left and the middle
panels display a surface plot and a contour plot for the BF estimates (in log scale), respectively.
The right panel shows the ratio of the SE to the BF estimates (in log scale), where the SEs are
evaluated using the BM method.
large value of 0.8 for each j, then the probability of getting a regeneration given a visit to D is
 = (0.8)100 ≈ 2 × 10−10 . Being overoptimistic that the Markov chain visits D with probability
close to 1, it would still take 100 billon iterations for the chain to regenerate about twenty times,
barely enough for the RS method to be effective.
Using the EB estimate ĥ of h, estimation of the remaining parameters ψ and prediction of the
spatial random field can be done in the standard method using MCMC samples from πĥ (ψ|y)
(see e.g. Roy et al., 2016, section 3.2). Thus we can produce the root rot disease prediction map
similar to that in Roy et al. (2016, Web Fig. 10).

6

Discussion

In this paper we consider two separate but related problems. The first problem is estimating the
ratios of unknown normalizing constants given Markov chain samples from each of the k > 1
probability densities. The second problem is estimating expectations of a function with respect
to a large number of probability distributions. These problems are related in the sense that generalized IS estimators used for the latter utilize estimates derived when solving the first problem.
The first situation also arises in a variety of contexts other than the generalized IS estimators.
For both problems, we derive estimators with flexible weights and thus these estimators are
19

appropriate for Markov chains with different mixing behaviors. We establish CLTs for these
estimators and develop BM methods for consistently estimating their SEs. These easy to calculate
SEs are important for at least three reasons. First, SEs are needed to assess the quality of the
estimates. Second, our ability to calculate SEs allows us to search for optimal weights a for both
stage 1 and 2. And last but not least, SEs form the basis for comparison of generalized IS with
other available methods for estimating large number of (ratios of) normalizing constants.
Although we compare BM and RS in this paper, spectral estimators can also be derived for
variance estimation using the results in Vats et al. (2015b). However, estimation by spectral
methods is generally more expensive computationally. Further, Flegal and Jones (2010) compare
the performance of confidence intervals produced by BM, RS and spectral methods for the time
average estimator, and they conclude that if tuning parameters are chosen appropriately, all these
three methods perform equally well. Control variates can be used to further improve the accuracy
of our generalized IS estimators (Doss, 2010; Owen and Zhou , 2000). A direction of future
research would be to establish a BM estimator of the SEs for control variate based methods.

A

Proof of Theorem 1

In the appendices, we denote Doss and Tan (2014) by D&T. The proof of the consistency of
dˆ follows from D&T section A.1 and is omitted. Establishing a CLT for dˆ is analogous to
section A.2 of D&T, but there are significant differences. Below we establish the CLT for dˆ and
finally show that Vb is a consistent estimator of V .
We begin by considering n1/2 (ζ̂ − ζ 0 ). As before, let ∇ represents the gradient operator. We
expand ∇`n at ζ̂ around ζ 0 , and using the appropriate scaling factor, we get

− n−1/2 ∇`n (ζ̂) − ∇`n (ζ 0 ) = −n−1 ∇2 `n (ζ ∗ )n1/2 (ζ̂ − ζ 0 ),

(A.1)

where ζ ∗ is between ζ̂ and ζ 0 . Consider the left side of (A.1), which is just n1/2 n−1 ∇`n (ζ 0 ),
since ∇`n (ζ̂) = 0. There are several nontrivial components to the proof, so we first give an
outline.
1. Following D&T we show that each element of the vector n−1 ∇`n (ζ 0 ) can be represented
as a linear combination of mean 0 averages of functions of the k chains.
2. Based on Step 1, applying CLT for each of the k Markov chain averages, we obtain a CLT
d
for the scaled score vector. In particular, we show that n1/2 n−1 ∇`n (ζ 0 ) → N (0, Ω), where
Ω defined in (2.8) involves infinite sums of auto-covariances of each chain.
20

† a.s.
a.s.
3. Following Geyer (1994) it can be shown that −n−1 ∇2 `n (ζ ∗ ) −→ B and that −n−1 ∇2 `n (ζ ∗ ) −→
B † , where B is defined in (2.7).
d

4. We conclude that n1/2 (ζ̂ − ζ 0 ) → N (0, B † ΩB † ).
5. Since d = g(ζ 0 ) and dˆ = g(ζ̂), where g is defined in (2.4), by the Delta method it follows
d
that n1/2 (dˆ − d) → N (0, V ) where V = D> B † ΩB † D.
We now provide the details.
1. Start by considering n−1 ∇`n (ζ 0 ). For r = 1, . . . , k, from D&T we have
nl
nr
k
X
X
 X
∂`n (ζ 0 )
(l)
(r)
pr (Xi , ζ 0 )
= wr
wl
1 − pr (Xi , ζ 0 ) −
∂ζr
i=1
i=1
l=1
l6=r

(can be shown to) = wr

nr 
X
i=1
k
X

−

l=1
l6=r



(r)
1 − pr (Xi , ζ 0 ) − 1 − Eπr pr (X, ζ 0 )
wl

nl
X


(A.2)


(l)
pr (Xi , ζ 0 ) − Eπl pr (X, ζ 0 ) .

i=1

That is, (A.2) can be used to view n−1 ∂`n (ζ 0 )/∂ζr as a linear combination of mean 0
averages of functions of the k chains.
2. Next, we need a CLT for the vector ∇`n (ζ 0 ) = (∂`n (ζ 0 )/∂ζ1 , · · · , ∂`n (ζ 0 )/∂ζk )T , that is,
d
to show that n−1/2 ∇`n (ζ 0 ) → N (0, Ω) as n → ∞. Note that,
nl
k


1 X X
1 ∂`n (ζ 0 )
(l)
√
= −√
wl
pr (Xi , ζ 0 ) − Eπl pr (X, ζ 0 )
n ∂ζr
n l=1
i=1
r
nl
k
X n


1 X
(l)
=−
al √
pr (Xi , ζ 0 ) − Eπl pr (X, ζ 0 )
nl
nl i=1
l=1

=−

k
X
√
nal Ȳ¯ (r,l) ,
l=1

P l
(r,l)
(r,l)
where Ȳ¯ (r,l) := n1l ni=1
Yi
and Yi
is as defined in (2.6). Since pr (x, ζ) ∈ (0, 1) for

all x, r and ζ, we have Eπl |pr (X, ζ 0 ) − Eπl (pr (X, ζ 0 ))|2+δ < ∞ for any δ > 0. Then
since Φl is polynomially ergodic of order m > 1, we have asymptotic normality for the

21

√
univariate quantities nl Ȳ¯ (r,l) (see e.g. Corollary 2 of Jones (2004)). Since nl /n → sl for
l = 1, . . . , k and al ’s are known, by independence of the k chains, we conclude that
1 ∂`n (ζ 0 ) d
√
→ N (0, Ωrr ) as n → ∞,
n ∂ζr
where Ω is defined in (2.8). Next, we extend the component-wise CLT to a joint CLT.
Consider any t ∈ (t1 , · · · , tk ) ∈ Rk , we have
1 ∂`n (ζ 0 )
1 ∂`n (ζ 0 )
+ · · · + tk √
t1 √
n ∂ζ1
n ∂ζk
!
P
Pnl
k
(1,l)
(k,l)
nl
X
√
√
Y
Y
=−
t1 nal i=1 i
+ · · · + tk nal i=1 i
n
nl
l
l=1


Pnl
(1,l)
(k,l)
k r
+ · · · + tk Yi
X
i=1 t1 Yi
n
d
=−
al
→ N (0, tT Ωt) as n → ∞.
√
nl
nl
l=1
Hence, the Cramér-Wold device implies the joint CLT,
d

n−1/2 ∇`n (ζ 0 ) → N (0, Ω)

as n → ∞.

(A.3)

Steps 3-5 are omitted since the derivations are basically the same as in D&T.
Next we provide a proof of the consistency of the estimate of the asymptotic covariance
a.s.
b −→
bB
b†D
b†Ω
b >B
V ≡ D> B † ΩB † D as n → ∞. Since
matrix V , that is, we show that Vb ≡ D
a.s.
a.s.
a.s.
a.s.
b −→
b −→
ζ̂ −→ ζ 0 and dˆ −→ d, it implies that D
D. From D&T, we know that B
B and using
†
† a.s.
b
b
the spectral representation of B and of B, it follows that B −→ B .
a.s.
b −→
b is defined in
To complete the proof, we now show that Ω
Ω where the BM estimator Ω
(2.13). This will be proved in couple of steps. First, we consider a single chain Φl used to calculate k quantities and establish a multivariate CLT. We use the results in Vats et al. (2015a) who
obtain conditions for the nonoverlapping BM estimator to be strongly consistent in multivariate
b is
settings. Second, we combine results from the k independent chains. Finally, we show that Ω
a strongly consistent estimator of Ω.

>
Denote Ȳ¯ (l) = Ȳ¯ (1,l) , Ȳ¯ (2,l) , . . . , Ȳ¯ (k,l) . Similar to deriving (A.3) via the Cramér-Wold
device, we have the following joint CLT for Wl :,
is a k × k covariance matrix with
Σ(l)
rs

=

(r,l) (s,l)
Eπl {Y1 Y1 }

+

∞
X

√ ¯ (l) d
nl Ȳ → N (0, Σ(l) ) as nl → ∞, where Σ(l)

(r,l) (s,l)
Eπl {Y1 Y1+i }

i=1

+

∞
X
i=1

22

(r,l)

(s,l)

Eπl {Y1+i Y1

}.

(A.4)

The nonoverlapping BM estimator of Σ(l) is given in (2.10). We now prove the strong consis(r,l)
b (l) . Note that Σ
b (l) is defined using the terms Z̄m
tency of Σ
’s which involve the random quantity
(l)
(l)
b (ζ 0 ) to be Σ
b with ζ 0 substituted for ζ̂, that is,
ζ̂. We define Σ
b (l) (ζ 0 ) =
Σ

el −1 h
ih
i>
bl X
Ȳm(l) − Ȳ¯ (l) Ȳm(l) − Ȳ¯ (l)
for l = 1, . . . , k,
el − 1 m=0

>
P(m+1)bl (r,l)
a.s.
(1,l)
(k,l)
(r,l)
b (l) −→
Yj /bl . We prove Σ
Ȳm , . . . , Ȳm
with Ȳm := j=mbl +1
Σ(l)
a.s.
a.s.
b (l) (ζ 0 ) −→
b (l) − Σ
b (l) (ζ 0 ) −→
in two steps: (1) Σ
Σ(l) and (2) Σ
0. Strong consistency of
(l)
b (ζ 0 ) requires both el → ∞ and bl → ∞. Since for all r,
the multivariate BM estimator Σ

Eπl |pr (X, ζ 0 ) − Eπl (pr (X, ζ 0 ))|4+δ < ∞ for any δ > 0, Φl is polynomially ergodic of order
a.s.
b (l) (ζ 0 ) −→
m > 1, and bl = bnνl c where 1 > ν > 0, it follows from Vats et al. (2015a) that Σ
Σ(l)
a.s.
b (l)
b (l)
b (l)
b (l)
as nl → ∞. We show Σ
rs − Σrs (ζ 0 ) −→ 0 where Σrs and Σrs (ζ 0 ) are the (r, s)th elements
b (l)
b (l)
of the k × k matrices Σ
rs and Σrs (ζ 0 ) respectively. By the mean value theorem (in multiple
variables), there exists ζ ∗ = tζ̂ + (1 − t)ζ 0 for some t ∈ (0, 1), such that
(l)

where Ȳm =



b (l) − Σ
b (l) (ζ 0 ) = ∇Σ
b (l) (ζ ∗ ) · (ζ̂ − ζ 0 ),
Σ
rs
rs
rs

(A.5)

where · represents the dot product. Note that
b (l)
Σ
rs (ζ) =

el −1
bl X
(s,l)
[Z̄ (r,l) (ζ) − Z̄¯ (r,l) (ζ)][Z̄m
(ζ) − Z̄¯ (s,l) (ζ)],
el − 1 m=0 m

P(m+1)bl
P l
(l)
(r,l)
(l)
pr (Xj , ζ)/nl . Some calcuwhere Z̄m (ζ) := j=mbl +1
pr (Xj , ζ)/bl and Z̄¯ (r,l) (ζ) := nj=1
lations show that for t 6= r
(r,l)

∂ Z̄m (ζ)
1
=−
∂ζ t
bl

(m+1)bl

X

(l)

(l)

pr (Xj , ζ)pt (Xj , ζ)

j=mbl +1

and
(r,l)

∂ Z̄m (ζ)
1
=
∂ζ r
bl

(m+1)bl

X

(l)

(l)

pr (Xj , ζ)(1 − pr (Xj , ζ)).

j=mbl +1

(r,l)
r
We denote Ūm
:= Z̄m (ζ) − Eπl [pr (X, ζ)], Ū¯ r := Z̄¯ (r,l) (ζ) − Eπl [pr (X, ζ)], and similarly the
(r,l)
(r,t)
centered versions of ∂ Z̄m (ζ)/∂ζ t and ∂ Z̄¯ (r,l) (ζ)/∂ζ t by V̄m and V̄¯ (r,t) respectively. Since
pr (X, ζ) is uniformly bounded by 1 and Φl is polynomially ergodic of order m > 1, there exist
√ r d
√ (r,t) d
√
d
2
2
σr2 , τr,t
< ∞ such that bl Ūm
→ N (0, σr2 ), nl Ū¯ r → N (0, σr2 ), bl V̄m → N (0, τr,t
), and

23

√ ¯ (r,t) d
2
nl V̄
). We have
→ N (0, τr,t
el −1
b (l)
p
p
p
p
1 X
∂Σ
rs (ζ)
r
s
=
[ bl (Ūm
− Ū¯ r ) bl (V̄m(s,t) − V̄¯ (s,t) ) + bl (V̄m(r,t) − V̄¯ (r,t) ) bl (Ūm
− Ū¯ s )]
∂ζ t
el − 1 m=0

=

el −1 h
i
p
p
p
p
1 h√ ¯ r √ ¯ (s,t) √ ¯ (r,t) √ ¯ s i
1 X
r
s
bl Ūm
bl V̄m(s,t) + bl V̄m(r,t) bl Ū¯m
−
nl Ū nl V̄
+ nl V̄
nl Ū .
el − 1 m=0
el − 1

It is easy to see that the negative term in the above expression goes to zero as el → ∞. Further,
since
p
p
 1

1
r 2
r
) + bl (V̄m(s,t) )2 ,
bl (Ūm
bl Ūm
bl V̄m(s,t) ≤
2
2
we have
el −1
el −1
l −1
p
p

 1 1 eX

 a.s. 1
1 X
1 1 X
1 2
r
(s,t)
r 2
≤
bl Ūm bl V̄m
bl (Ūm ) +
bl (V̄m(s,t) )2 −→ σr2 + τs,t
,
el − 1 m=0
2 el − 1 m=0
2 el − 1 m=0
2
2

where the last step above is due to strong consistency of the BM estimators for the asymptotic
(l)
(l)
variances of the sequences {pr (Xj , ζ), j = 1, · · · , nl } and {∂ps (Xj , ζ)/∂ζ t , j = 1, · · · , nl }
respectively. Similarly, we have
"
#
el −1
eX
eX
l −1
l −1
p
p


1 X
1
1
1
1
1
a.s. 1 2
s
s 2
≤
−→ τr,t
bl V̄m(r,t) bl Ūm
bl (V̄m(r,t) )2 +
bl (Ūm
)
+ σs2 .
el − 1 m=0
2 el − 1 m=0
2 el − 1 m=0
2
2
(r,t)

2
r
, etc, above actually depends on ζ, and we are indeed conVm , σr2 , τr,t
Note that the terms Um
a.s.
cerned with the case where ζ takes on the value ζ ∗ , lying between ζ̂ and ζ 0 . Since, ζ̂ −→ ζ 0 ,
a.s.
ζ ∗ −→ ζ 0 as nl → ∞. Let kuk denotes the L1 norm of a vector u ∈ Rk . So from (A.5), and the
b (l)
fact that ∂ Σ
rs (ζ)/∂ζ t is bounded with probability one, we have
(
)
∗
b (l)
∂
Σ
(ζ
)
rs
a.s.
b (l)
b (l)
|Σ
kζ̂ − ζ 0 k −→ 0 as n → ∞.
rs − Σrs (ζ 0 )| ≤ max
1≤t≤k
∂ζ t
a.s.
a.s.
b (l) −→
b −→
b is defined in (2.11)
Since Σ
Σ(l) , for l = 1, . . . , k, it follows that Σ
Σ where Σ
b
and Σ is the corresponding k 2 × k 2 covariance matrix, that is, Σ is a block diagonal matrix as Σ
b (l) , l = 1, . . . , k. Since nl /n → sl for l = 1, . . . , k, we have An → As
with Σ(l) substituted for Σ
as n → ∞ where An is defined in (2.12) and
r
r

 r
1
1
1
a1 Ik −
a2 Ik . . . −
ak Ik .
As = −
s1
s2
sk
a.s.
b ≡ An ΣA
b Tn −→
Finally from (2.8) and (A.4) we see that Ω = As ΣATs . So from (2.13) we have Ω
As ΣATs = Ω as n → ∞.

24

B

Proof of Theorem 2

As in Buta and Doss (2011) we write
√

√
ˆ
n(û(π, π1 ; a, d)−û(π,
π1 ; a, d))+ n(û(π, π1 ; a, d)−u(π, π1 )).
(B.1)
First, consider the 2nd term, which involves randomness only from the 2nd stage. From (3.3)
P
note that kl=1 al Eπl u(X; a, d) = u(π, π1 ). Then from (3.1) we have
ˆ
n(û(π, π1 ; a, d)−u(π,
π1 )) =

√

√

r Pnl
k
(l)
X
n i=1 (u(Xi ; a, d) − Eπl u(X; a, d))
.
n(û(π, π1 ; a, d) − u(π, π1 )) =
al
√
n
n
l
l
l=1

Since Φl is polynomially ergodic of order m and Eπl |u(X; a, d)|2+δ is finite where m > 1 + 2/δ,
P l
√ d
(l)
(u(Xi ; a, d) − Eπl u(X; a, d))/ nl → N (0, τl2 (π; a, d)) where τl2 (π; a, d)
it follows that ni=1
is defined in (3.4). As nl /n → sl and the Markov chains Φl ’s are independent, it follows that
√
d
n(û(π, π1 ; a, d) − u(π, π1 )) → N (0, τ 2 (π; a, d)).
Now we consider the 1st term in the right hand side of (B.1). Letting F (z) = û(π, π1 ; a, z),
by Taylor series expansion of F about d we have
√
√
√
>
ˆ − F (d)) = n∇F (d) (dˆ − d) + n (dˆ − d)> ∇2 F (d∗ )(dˆ − d),
(B.2)
n(F (d)
2
ˆ Simple calculations show that
where d∗ is between d and d.
[∇F (d)]j−1

nl
k
(l)
(l)
X
al X
aj νj (Xi )ν(Xi )
a.s.
=
−→ [c(π; a, d)]j−1
Pk
(l)
2 2
nl i=1 (
l=1
s=1 as νs (Xi )/ds ) dj

(B.3)

where [c(π; a, d)]j−1 is defined in (3.5). We know that n/N → q. Using similar arguments as
in Buta and Doss (2011), it follows that ∇2 F (d∗ ) is bounded in probability. Thus from (B.2) we
have
r
r
√
√
√
1
n
n √ ˆ
>
ˆ
ˆ
n(F (d) − F (d)) =
∇F (d) N (d − d) + √
[ N (d − d)]> ∇2 F (d∗ )[ N (dˆ − d)]
N
2 N N
√
√
= qc(π; a, d)> N (dˆ − d) + op (1).
Then Theorem 2 (1) follow from (B.1) and the independence of the two stages of Markov chain
sampling.
Next to prove Theorem 2 (2), note that, we already have a consistent BM estimator Vb of
a.s.
V . From (B.3), we have [ĉ(π; a, d)]j−1 = [∇F (d)]j−1 −→ [c(π; a, d)]j−1 . Applying mean
25

value theorem on [∇F (d)]j−1 and the fact that ∇2 F (d∗ ) is bounded in probability, it follows that
P
Pk−1
a.s.
ˆ j−1 − [ĉ(π; a, d)]j−1 −→
[ĉ(π; a, d)]
0. Writing c(π; a, d)> V c(π; a, d) as k−1
i=1
j=1 ci Vij cj , it
a.s.
>
>
ˆ Vb ĉ(π; a, d)
ˆ −→ c(π; a, d) V c(π; a, d).
then follows that ĉ(π; a, d)
ˆ is a consistent estimator of τ 2 (π; a, d) where τ 2 and τ̂ 2 are defined
We now show τ̂l2 (π; a, d)
l
l
l
(l) l
in (3.4) and (3.7), respectively. Since the Markov chains {Xi }ni=1
are independent, it then folˆ completing the proof of Theorem 2
lows that τ 2 (π; a, d) is consistently estimated by τ̂ 2 (π; a, d)
(2).
If d is known from the assumptions of Theorem 2 (2) and the results in Vats et al. (2015a),
we know that τl2 (π; a, d) is consistently estimated by its BM estimator τ̂l2 (π; a, d). Note that,
(l)
ˆ −
τ̂l2 (π; a, d) is defined in terms of the quantities u(Xi ; a, d)’s. We now show that τ̂l2 (π; a, d)
a.s.
τ̂l2 (π; a, d) −→ 0.
Denoting τ̂l2 (π; a, z) by G(z), by the mean value theorem (in multiple variables), there exists
ˆ − G(d) = ∇G(d∗ ) · (dˆ − d). For any
d∗ = tdˆ + (1 − t)d for some t ∈ (0, 1), such that G(d)
k−1
j ∈ {2, · · · , k}, and z ∈ R+ ,
"e −1

#
l
X
∂G(z)
∂ ūm (a, z) ∂ ū¯(a, z)
bl
2(ūm (a, z) − ū¯(a, z))
=
−
(B.4)
∂zj
el − 1 m=0
∂zj
∂zj
¯ := ū¯(a, z) − E (u(X; a, z)). Note that, there
Let W̄m := ūm (a, z) − Eπl (u(X; a, z)) and W̄
πl
√
√ ¯ d
d
2
2
exists, σ < ∞ such that bl W̄m → N(0, σ ), and nl W̄ → N(0, σ 2 ). Simple calculations
show that


(m+1)bl
(l)
(l)
aj 1 X  ν(Xi )νj (Xi )
∂ ūm (a, z)

= 2
 P
 2 .
∂zj
zj bl i=mb +1
(l)
l
s as νs (Xi )/zs
P
Hence, letting αj = Eπl [ν(X)νj (X)/ ( s as νs (X)/zs )2 ], we write
∂ ūm (a, z) ∂ ū¯(a, z)
aj 
aj n o
−
≡ 2 Z̄m,j − 2 Z̄¯j ,
∂zj
∂zj
zj
zj
Pl
P
(l)
(l)
(l)
where Z̄1,j = (1/bl ) bi=1
[ν(Xi )νj (Xi )/{ s as νs (Xi )/zs }2 ] − αj and Z̄¯j is similarly de√
√
d
d
fined. Note that, there exists τj2 < ∞, such that bl Z̄m,j → N(0, τj2 ), and nl Z̄¯j → N(0, τj2 ).

26

From (B.4) we have
el −1 h
i
p
p 
∂G(z) aj 2 X
¯ ) b Z̄ − Z̄¯
= 2
bl (W̄m − W̄
l
m,j
j
∂zj
zj el − 1 m=0
el −1 h
i
p
p
aj 2 X
bl W̄m bl Z̄m,j
= 2
zj el − 1 m=0
"
#
eX
eX
l −1
l −1
aj
1
e
1
l
¯ Z̄¯
¯
− 2 2bl Z̄¯j
W̄m + W̄
Z̄m,j −
W̄
j
zj
el − 1 m=0
el − 1 m=0
el − 1

=

el −1 h
i a
p
p
2 h√ ¯ √ ¯ i
aj 2 X
j
b
W̄
b
Z̄
−
nl W̄ nl Z̄j .
l m
l m,j
zj2 el − 1 m=0
zj2 el − 1

Then using similar arguments as in the proof of Theorem 1, it can be shown that ∂G(z)/∂zj is
bounded with probability one. Then it follows that


∂G(d∗ )
a.s.
ˆ
|G(d) − G(d)| ≤ max
kdˆ − dk −→ 0.
1≤j≤k−1
∂zj

C

Proof of Theorem 3

As in the proof of Theorem 2 we write
√ [f ]
√
√
ˆ − Eπ f ) = n(η̂ [f ] (π; a, d)
ˆ − η̂ [f ] (π; a, d)) + n(η̂ [f ] (π; a, d) − Eπ f ). (C.1)
n(η̂ (π; a, d)
First, consider the 2nd term, which involves randomness only from the 2nd stage. Since
a.s.

v̂ −→

k
X
l=1

[f ]

Z

al Eπl v (X; a, d) =
X

f (x)
Pk
s=1

Pk

l=1

al νl (x)/ml

as νs (x)/(ms /m1 )

Pk

ν(x) µ(dx) =

m
Eπ f,
m1

al Eπl v [f ] (X; a, d) = Eπ f u(π, π1 ). Then from (3.1) we have
!
!
r
nl
k
(l)
X
√
v̂ [f ] (π; a, d) − Eπ f u(π, π1 )
v [f ] (Xi ; a, d) − Eπl v [f ] (X; a, d)
n 1 X
n
=
al
.
√
(l)
n
n
û(π, π1 ; a, d) − u(π, π1 )
u(X
;
a,
d)
−
E
u(X;
a,
d)
l
l i=1
π
i
l
l=1
(C.2)
From the conditions of Theorem 3 and the fact that the Markov chains Φl , l = 1, . . . , k are independent, it follows that the above vector (C.2) converges in distribution to the bivariate normal
distribution with mean 0 and covariance matrix Γ(π; a, d) defined in (3.9). Then applying the
Delta method to the function g(x, y) = x/y we have a CLT for the ratio estimator η̂ [f ] (π; a, d),
we have

l=1

27

√
d
that is, we have n(η̂ [f ] (π; a, d) − Eπ f ) → N (0, ρ(π; a, d)) where ρ(π; a, d)) is defined in
(3.10).
Next letting L(z) = η̂ [f ] (π; a, z), by Taylor series expansion of L about d we have
√
√
√
>
ˆ − L(d)) = n∇L(d) (dˆ − d) + n (dˆ − d)> ∇2 L(d∗ )(dˆ − d),
n(L(d)
(C.3)
2
ˆ Simple calculations show that
where d∗ is between d and d.
a.s.

[∇L(d)]j−1 = [ê(π; a, d)]j−1 −→ [e(π; a, d)]j−1

(C.4)

where [e(π; a, d)]j−1 and [ê(π; a, d)]j−1 are defined in (3.11) and (3.12) respectively. It can be
√
ˆ − L(d)) =
shown that ∇2 L(d∗ ) is bounded in probability. Thus from (C.3) we have n(L(d)
√
√
qe(π; a, d)> N (dˆ− d) + op (1). Then Theorem 3 (1) follow from (C.1) and the independence
of the two stages of Markov chain sampling.
Next to prove Theorem 3 (2), note that, we already know that Vb is a consistent BM estimator
a.s.
of V . From (C.4), we have [ê(π; a, d)]j−1 −→ [e(π; a, d)]j−1 . Applying mean value theorem on
ˆ j−1 −
[∇L(d)]j−1 and the fact that ∇2 L(d∗ ) is bounded in probability, it follows that [ê(π; a, d)]
a.s.
[ê(π; a, d)]j−1 −→ 0.
a.s.
a.s.
ˆ −→
ˆ −→
From (3.8) we know that û(π, π1 ; a, d)
u(π, π1 ). From (3.13) we know η̂ [f ] (π; a, d)
a.s.
ˆ −→
Eπ f . Since v̂ [f ] (π, π1 ; a, d) = η̂ [f ] (π; a, d)û(π, π1 ; a, d), it follows that v̂ [f ] (π, π1 ; a, d)
a.s.
ˆ û(π, π1 ; a, d))
ˆ −→
Eπ f u(π, π1 ). Thus ∇h(v̂ [f ] (π, π1 ; a, d),
∇h(Eπ f u(π, π1 ), u(π, π1 )). Thus
a.s.
bl (π; a, d)
ˆ −→
to prove Theorem 3 (2), we only need to show that Γ
Γl (π; a, d).
If d is known from the assumptions of Theorem 3 (2) and the results in Vats et al. (2015a),
bl (π; a, d). We now show
we know that Γl (π; a, d) is consistently estimated by its BM estimator Γ
a.s.
bl (π; a, d)
bl (π; a, d) −→ 0.
ˆ −Γ
that Γ
a.s.
ˆ − γ̂ 22 (π; a, d) −→
From Theorem 2 (2), we know that γ̂l22 (π; a, d)
0. We now show
l
a.s.
11
11
11
ˆ − γ̂ (π; a, d) −→ 0. Letting γ̂ (π; a, z) by H(z), by the mean value theorem,
γ̂l (π; a, d)
l
l
∗
ˆ
ˆ − H(d) = ∇H(d∗ ) · (dˆ− d).
there exists d = td + (1 − t)d for some t ∈ (0, 1), such that H(d)
k−1
For any j ∈ {2, · · · , k}, and z ∈ R+ ,
"e −1
!#
[f ]
l
[f ]
X
¯
∂H(z)
∂v̄
(a,
z)
bl
∂
v̄
(a,
z)
m
[f ]
.
=
2(v̄m
(a, z) − v̄¯[f ] (a, z))
−
∂zj
el − 1 m=0
∂zj
∂zj
¯ [f ] := v̄¯[f ] (a, z) − E (v [f ] (X; a, z)). Note
Let W̄m := v̄m (a, z) − Eπl (v [f ] (X; a, z)) and W̄
πl
√
√ ¯ [f ] d
[f ] d
2
2
that, there exists, σf < ∞ such that bl W̄m → N(0, σf ), and nl W̄ → N(0, σf2 ). Simple
[f ]

[f ]

28

calculations show that
[f ]
∂v̄m (a, z)

∂zj

(m+1)bl

=





(l)
(l)
(l)
aj 1 X  f (Xi )ν(Xi )νj (Xi ) 
 P
 2 .
zj2 bl i=mb +1
(l)
a
ν
(X
)/z
l
s
i
s s s

[f ]

Hence, letting αj = Eπl [f (X)ν(X)νj (X)/ (

P

s

as νs (X)/zs )2 ], we write

[f ]
aj n [f ] o aj n ¯ [f ] o
∂v̄m (a, z) ∂ ū¯[f ] (a, z)
−
≡ 2 Z̄m,j − 2 Z̄j ,
∂zj
∂zj
zj
zj

as νs (Xi )/zs }2 ] − αj and Z̄¯j is sim√
√
d
d
2
2
), and nl Z̄¯j →
< ∞, such that bl Z̄m,j → N(0, τj,f
ilarly defined. Note that, there exists τj,f
2
N(0, τj,f
). The rest of the proof is analogous to Theorem 2, in that we have
[f ]

where Z̄1,j = (1/bl )

Pbl

(l)

(l)

(l)

i=1 [f (Xi )ν(Xi )νj (Xi )/{

(l)

P

[f ]

[f ]

s

el −1 h
i
p
p [f ] i aj 2 h√
∂H(z)
aj 2 X
¯ [f ] √n Z̄¯ [f ] .
= 2
bl W̄m[f ] bl Z̄m,j − 2
nl W̄
l j
∂zj
zj el − 1 m=0
zj el − 1
a.s.

a.s.

11
12
12
ˆ
ˆ
Then it can be shown γ̂l11 (π; a, d)−γ̂
l (π; a, d) −→ 0 and finally γ̂l (π; a, d)−γ̂l (π; a, d) −→
0.

D

Regeneration with general weights

Tan et al. (2015) provide a regeneration based central limit theorem (CLT) for the estimators η̂
and û defined in 1.3 and 3.1 respectively in the main text. In the case when d is unknown, they
ˆ for their results to hold,
allow only a special choice for the weight vector, namely a = (1, d)
where dˆ is the estimator of d based on the Stage 1 chains discussed in Section 2 of the main text.
In this section, we establish a regeneration based CLT for η̂ and û with any choice of the weight
vector a.
We will refer to the following conditions.
(l)

(l)

A1 For each l = 1, . . . , k, the Markov chain Φl = {X0 , X1 , . . .} is geometrically ergodic
and has πl as its invariant density.
A2 Let kl : X×X → [0, ∞) be the Markov transition density for Φl , so that for any measurable

R
(l)
(l)
set A we have P Xn+1 ∈ A | Xn = x = A kl (y|x)µ(dy). Suppose that for each l =
1, . . . , k, kl satisfies the following minorization condition:
kl (y|x) ≥ sl (x) ql (y)
29

for all x, y ∈ X,

(D.1)

where the function sl : X → [0, 1) with Eπl sl > 0, and ql is a probability density function
on X.
A3 Recall the functions u(X; a, d) and v [f ] (X; a, d) defined in (3.2) of our paper. There exists
 > 0 such that Eπl |v [f ] (X; a, d)|2+ and Eπl |u(X; a, d)|2+ are finite.
A4 Suppose Φl is simulated for Rl regenerative tours for l = 1, . . . , k. Assume Rl /R1 → bl ∈
(0, ∞) as R1 → ∞.
(l)

Following Tan et al. (2015), let the regeneration times for the lth Markov chain be τ0 =


(l)
(l)
0, τ1 , τ2 , . . .. Accordingly, the chain Φl is broken up into “tours” Xτ (l) , . . . , Xτ (l) −1 , t =
t

t−1

1, 2, . . . that are independent stochastic replicas of each other. Suppose we simulate Rl tours of
(l)
the lth Markov chain for l = 1, . . . , k, so the length of the lth chain is nl = τRl . Also as in Tan
et al. (2015), for t = 1, 2, . . . , Rl define
(l)

(l)

Vt

=

(l)

(l)

τt −1

τt −1

τt −1

X

X

X

(l)

i=τt−1

(l)

v [f ] (Xi ; a, d),

(l)

Ut =

(l)

u(Xi ; a, d),

(l)

and

(l)

Tt =

(l)

(l)

1 = τt − τt−1 ,

(l)

i=τt−1

i=τt−1

(D.2)
where the sums range over the values of i that constitute the tth tour.
ˆ
Recall from Remark 4 in Section 3 of our paper, when d is unknown, we set a = w ∗ (1, d)
where ∗ denotes component-wise multiplication. That is, (a1 , . . . , ak ) = (w1 , w2 , . . . , wk ) ∗
(1, dˆ2 , . . . , dˆk ) for any pre-determined weight w. With this choice, the expressions for u and v [f ]
in (3.2) become

ˆ dˆ = P ν(x)
u x; w ∗ (1, d),
k
l=1 wl νl (x)


ˆ dˆ = Pf (x)ν(x) .
and v [f ] x; w ∗ (1, d),
k
l=1 wl νl (x)

(D.3)


ˆ and consequently for each l, the triples Vt(l) , Ut(l) , Tt(l) , t =
The above quantities do not involve d,
0, 1, 2, . . . defined in (D.2) are iid, and we have independence across l’s. The estimator for η re-

30

duces to
η̂ = η̂N,n

nl
nl
k
k
(l)
(l)  X
(l)
ˆl X
ˆl X
 X
f
(X
)ν(X
)
w
d
ν(Xi )
w
d
l
l
i
i
ˆ dˆ =
w ∗ (1, d),
P
P
nl i=1 k wl νs (Xi(l) ) l=1 nl i=1 k ws νs (Xi(l) )
l=1
s=1
s=1

(D.4)
=

Rl
k
X
wl dˆl X

nl

l=1

=

k
X
l=1

(l)
Vt

X
k

t=1

V̄ (l)
wl d̂l (l)
T̄

l=1

X
k

wl d̂l

l=1

Rl
wl dˆl X
(l)
Ut
nl t=1

Ū (l)
,
T̄ (l)

where
(l)

τt −1
(l)
Ut

=

(l)
i=τt−1

(l)

τt −1

(l)

ν(Xi )

X

(l)

Pk

s=1 ws νs (Xi )

(l)
Vt

and

=

X
(l)

i=τt−1

(l)

(l)

f (Xi )ν(Xi )
,
Pk
(l)
w
ν
(X
)
s
s
i
s=1

P l (l)
P l (l)
T̄ (l) = Rl−1 R
be the average tour length and, analogously, V̄ (l) = Rl−1 R
and
t=1 Tt
t=1 Vt
PRl (l)
−1
(l)
Ū = Rl
t=1 Ut . Similarly, the estimator for m/m1 reduces to
nl
k
(l)
 X
wl dˆl X
ν(Xi )
ˆ dˆ =
û = ûN,n w ∗ (1, d),
P
nl i=1 k ws νs (Xi(l) )
l=1
s=1

=

Rl
k
X
wl dˆl X
l=1

nl

(l)

Ut =

t=1

k
X

(D.5)

(l)

wl d̂l

l=1

Ū
.
T̄ (l)

Theorem 4 below gives the asymptotic distributions of η̂ and û. It extends Tan et al.’s (2015)
Theorem 2 to the general choice of weight vector a. To state the theorem, we first need to define
some notation. Let M̃ and L̃ be the vectors of length k − 1 for which the (j − 1)th coordinates
are, for j = 2, . . . , k,
M̃j−1 = wj Eπj u and
wj Eπj v [f ]

L̃j−1 = Pk

l=1

wl dl Eπl u

Pk

l=1

−



wl dl Eπl v [f ] wj Eπj u
.
2
Pk
w
d
E
u
l
l
π
l
l=1

(D.6)

As in Tan et al. (2015), assume that in Stage 1, for l = 1, . . . , k, the lth chain has been run for ρl
(l)
(l)
regenerations. So the length of the lth chain, Nl = T1 + . . . + Tρl , is random. We assume that
ρ1 , . . . , ρk → ∞ in such a way that ρl /ρ1 → cl ∈ (0, ∞), for l = 1, . . . , k.
31

Theorem 4 Suppose that for the Stage 1 chains, conditions A1 and A2 hold, and that for the
Stage 2 chains, conditions A1–A4 hold. If ρ1 → ∞ and R1 → ∞ in such a way that R1 /ρ1 →
q ∈ [0, ∞), then
 d

1/2
R1 û − m/m1 → N 0, q M̃ > W M̃ + κ2
and
1/2

R1

 d

η̂ − η → N 0, q L̃> W L̃ + τ 2 ,

where M̃ , L̃ are given in equations (D.6), W , κ2 and τ 2 are given in equations (2.15), (2.8), and
(2.10) of Tan et al. (2015), respectively. In their (2.8) and (2.10), a is taken to be a = w ∗ (1, d).
c,
Furthermore, we can form strongly consistent estimates of the asymptotic variances if we use W
κ̂2 , and τ̂ 2 defined in (2.16) and (2.11) of Tan et al. (2015), respectively, and use the standard
empirical estimates of M̃ and L̃.

D.1

Proof of Theorem 4

We first prove the CLT for η̂. Note that



 
 

1/2 
ˆ dˆ −η = R11/2 η̂ w∗(1, d),
ˆ dˆ −η̂ w∗(1, d), d +R11/2 η̂ w∗(1, d), d −η .
R1 η̂ w∗(1, d),
(D.7)
The second term on the right side of (D.7) involves randomness coming only from Stage 2
sampling, and its distribution is given by Theorem 1 of Tan et al. (2015): it is asymptotically
normal with mean 0 and variance τ 2 . The first term involves randomness from both Stage 1 and
Stage 2 sampling. However, as in the proofs of Theorem 2 and 3, we can show that for this term,
the randomness from Stage 2 is asymptotically negligible, so that only Stage 1 sampling contributes to its asymptotic distribution. Finally, the asymptotic normality of the left side of (D.7)
follows since the two stages of sampling are independent. We now provide the details of the
proof.
Consider the first term on the right side of (D.7). Recall that if a = w ∗ (1, d), then
f (x)ν(x)
v [f ] (x) := v [f ] (x; a, d) = Pk
l=1 wl νl (x)

and u(x) := u(x; a, d) = Pk

ν(x)

l=1

wl νl (x)

.

With (D.4) and (D.5) in mind, define the function


A(z) = η̂ w ∗ (1, z), z =

nl
k
X
wl zl X
l=1

nl

i=1

v

[f ]

(l)
(Xi )

X
k
l=1

nl
wl zl X
(l)
u(Xi )
nl i=1

for z = (z2 , . . . , zk )> , with zl > 0 for l = 2, . . . , k, and z1 = 1. Note that setting z = d gives
ˆ = η̂(w ∗ (1, d),
ˆ d).
ˆ
A(d) = η̂(w ∗ (1, d), d), and setting z = dˆ gives A(d)
32

By a Taylor series expansion of A about d we get
1/2


R1 ˆ
ˆ
ˆ dˆ −η̂ w∗(1, d), d = R11/2 ∇A(d)> (d−d)+
ˆ
(d−d)> ∇2 A(d∗ )(d−d)
η̂ w∗(1, d),
2
1/2
>

R
1/2
1/2
1/2
= R1 ∇A(d)> (dˆ − d) + 1 ρ1 (dˆ − d) ∇2 A(d∗ ) ρ1 (dˆ − d) ,
2ρ1
ˆ As R1 → ∞, nl → ∞ for each l. We first show that the
where d∗ is between d and d.
gradient ∇A(d) converges almost surely to a finite constant vector by proving that each one of
its components, [A(d)]j−1 , j = 2, . . . , k, converges almost surely as R1 → ∞. As nl → ∞ for
l = 1, . . . , k, for j = 2, . . . , k, we have
Pnj [f ] (j)
(wj /nj ) i=1
v (Xi )
[∇A(d)]j−1 = Pk
Pnl
(l)
i=1 u(Xi )
l=1 (wl dl /nl )
Pnj
Pnl [f ] (l) 
Pk
(j) 
i=1 u(Xi )
i=1 v (Xi ) (wj /nj )
l=1 (wl dl /nl )
−
Pnl
Pk
(l) 2
i=1 u(Xi )
l=1 (wl dl /nl )


Pk
[f ]
wj Eπj u
wj Eπj v [f ]
a.s.
l=1 wl dl Eπl v
−→ Pk
.
−
2
Pk
u
w
d
E
u
w
d
E
l
l
π
l
l
π
l
l=1
l
l=1
1/2 

R1

The expression above corresponds to L̃j−1 , which is defined in (D.6), and it is finite by assumption A3. Next, we show that the random Hessian matrix ∇2 A(d∗ ) is bounded in probability, i.e., each element of this matrix is Op (1). As nl → ∞ for l = 1, . . . , k, for any
j, t ∈ {2, . . . , k}, j 6= t, we have
P t
(j) 
(t) 
wj Pnj
v [f ] (Xi ) wntt ni=1
u(Xi )
i=1
n
j
[∇2 A(d∗ )]t−1,j−1 = −
Pk wl d∗l Pnl
(l) 2
l=1 nl
i=1 u(Xi )
#
Pk wl d∗l Pnl [f ] (l)  wt Pnt
Pnt [f ] (t)
 X
"
(t) 
nj
wt
v
(X
)
u(X
)
)
v
(X
wj
i
i
i
l=1
i=1
i=1
i=1
(j)
n
n
nt
t
l
−
u(Xi )
−2
Pk wl d∗l Pnl
Pk wl d∗l Pnl
(l) 2
(l) 3
nj i=1
u(X
u(X
)
)
i
i
l=1 nl
i=1
l=1 nl
i=1
"
#

Pk
[f ]
(wt Eπt u)
(wj Eπj v [f ] )(wt Eπt u)
wt Eπt v [f ]
a.s.
l=1 wl dl Eπl v
−→ − Pk
,
2 − (wj Eπj u) Pk
2 − 2
3
Pk
l=1 wl dl Eπl u
l=1 wl dl Eπl u
l=1 wl dl Eπl u
where the limits are also finite.
Now, we can rewrite (D.7) as





1/2 
ˆ dˆ − η = (R1 /ρ1 )1/2 ∇A(d)> ρ11/2 (dˆ − d) + R11/2 η̂ w ∗ (1, d), d − η
R1 η̂ w ∗ (1, d),

 1/2
>
 1/2
1
+ 1/2 (R1 /ρ1 )1/2 ρ1 (dˆ − d) ∇2 A(d∗ ) ρ1 (dˆ − d)
2ρ1


1/2
1/2 
= q 1/2 [∇A(d)]> ρ1 (dˆ − d) + R1 η̂ w ∗ (1, d), d − η + op (1).
33

d
1/2
Since from Tan et al. (2015) we have ρ1 (dˆ − d) → N (0, W ) and the two sampling stages are
assumed to be independent, we conclude that

 d

1/2 
ˆ dˆ − η →
R1 η̂ w ∗ (1, d),
N 0, q L̃> W L̃ + τ 2 .

The proof of the CLT for û is similar. As in (D.7), we have





1/2 
ˆ dˆ − m/m1 = R11/2 û w ∗ (1, d),
ˆ dˆ − û w ∗ (1, d), d
R1 û w ∗ (1, d),


1/2 
+ R1 û w ∗ (1, d), d − m/m1 .

(D.8)

The asymptotic distribution of the second term in (D.8) is given in Tan et al.’s (2015) Theorem
1. The first term is linear in dˆ − d:

 X
nj
k

 X
1
(j)
ˆ
ˆ
(D.9)
u(Xi ) (dˆj − dj ).
û w ∗ (1, d), d − û w ∗ (1, d), d =
wj
n
j
i=1
j=2
For j = 2, . . . , k, the coefficient of (dˆj − dj ) in (D.9) converges almost surely to wj Eπj u, which
is the term M̃j−1 defined in (D.6).
Finally, from the independence of the two terms in (D.8) we conclude that

 d

1/2 
ˆ dˆ − m/m1 →
R1 û (1, d),
N 0, q M̃ > W M̃ + κ2 .

E

Toy example

In this section, we follow up on the simulation studies that involve t distributions from Section 4
of the main text to verify Theorems 1-3. We also discuss different weights in forming generalized
IS estimators and their effects on estimates of expectations and ratios of normalizing constants.
Let tr,µ denote the t-distribution with degree of freedom r and central parameter µ. We consider π1 (·) and π2 (·) as the density functions for a t5,µ1 =1 and t5,µ2 =0 , respectively. For simplicity,
let νi (·) = πi (·) for i = 1, 2. Our plan is to first estimate the ratio between the two normalizing
constants, d = m2 /m1 . Then we will study a sea of t-distributions Π = {t5,µ : µ ∈ M } where
M is a fine grid over [0, 1], say M = {0, .01, · · · , .99, 1}. For each µ ∈ M , we assume that
νµ (·) = πµ (·) and we estimate the ratio between its normalizing constant and m1 , denoted by
µ
. We also estimate the expectation of each distribution in Π, denoted Et5,µ X or Eµ X
dµ := m
m1
for short. Clearly, the exact answers are d = dµ = 1 and Eµ X = µ for any µ ∈ M . Nevertheless,
we follow the two-stage procedure from Sections 2 and 3 to generate Markov chains from π1
and π2 and build MCMC estimators from Theorems 1-3. The primary goal is to compare the
performance of BM and RS estimators.
34

We draw iid samples from π1 and Markov chain samples from π2 using the so called independent Metropolis Hastings algorithm with proposal density t5,1 . For RS, we follow the idea of
Mykland et al. (1995, Section 4.1) on constructing minorization conditions to identify regeneration times. Based on a carefully tuned minoration condition, the Markov chain for π2 regenerates
about every 3 iterations on average. In contrast, for users of the BM method proposed in this paper, no such theoretical development is needed. For i = 1, 2 we draw Ni observations from πi in
stage 1 and ni observations from πi in stage 2. We set N1 = N2 and n1 = n2 = N1 /10 = N2 /10.
Recall the reason for smaller stage 2 sample sizes is due to computing cost. For completeness,
note generating Markov chain samples using RS results in a random chain length so these chains
were run in such a way that N1 ∼ N2 and n1 ∼ n2 .
For estimators based on stage 1 samples, Theorem 1 allows any choice of weight, a[1] . For
estimators based on stage 2 samples, Theorem 2 and 3 allow any choice of weight, a[2] , in constructing consistent BM estimators of the asymptotic variances. RS based estimators in stage 2
are calculated using Theorems stated in Doss and Tan (2014) and Tan et al. (2015) with a general
weight choice noted in Remark 4. This is an important generalization in that now any nonnegative numerical weight vector can be used. We discuss the choice of weights and their impact
on the estimators later in this section.
The following details the simulation study presented in the main text. We consider increasing
sample sizes from N1 = 103 to 105 in order to examine trace plots for BM and RS estimates.
The two stage procedure is repeated 1000 times independently. The unknown true value of the
asymptotic variance of dˆ is estimated by its empirical asymptotic variance over the 1000 replications at N1 = 105 . We consider the naive weight, a[1] = (0.5, 0.5), that is proportional to
the sample sizes, and an alternative a[1] = (0.82, 0.18) that weighs the iid sample more than
the Markov chain sample. As illustrated in Figure 1 of the main text, both the BM and the RS
estimates approach the empirical asymptotic variance as the sample size increases suggesting
consistency. Similarly for stage 2, Figure 4 shows convergence of the BM and the RS estimates
to the corresponding empirical asymptotic variances of dˆµ and Êµ (X). Plots for other µ ∈ M
show similar results, but are not included here.
Overall, the simulation study suggests BM and RS methods provide consistent estimators for
the true asymptotic variance. RS estimators enjoy smaller mean squared error in most cases.
Nevertheless, when the number of regenerations is not great, BM estimators could be the more
stable estimator. For example, in the top left panel of Figure 4, at stage 2 sample size n2 ≈ n1 =
100, or about 35 regenerations for chain 2, the RS method substantially over-estimated the target
in about 5% of the replications. Further, in the cases where regeneration is unavailable or the
number of regenerations is extremely small, then BM would be the more viable estimator.
35

µ= 0.5

µ= 2
BM
RS
empirical

2.5

0.06

BM
RS
empirical

2.0

2.5

3.0

3.5

4.0

0.5

0.02

5

0.04

1.5

10

15

BM
RS
empirical

0

^
estimates of asy var for d

µ= −1

2.0

2.5

3.0

3.5

4.0

2.0

2.5

3.0

3.5

4.0

log10 of chain 1 length in stage 2

7

BM
RS
empirical

6

BM
RS
empirical

1

20

2

2

3

40

4

4

5

8

µ= 2

60

BM
RS
empirical

6

80 100

µ= 0.5

0

0

^
estimates of asy var for Eµ(X)

µ= −1

2.0

2.5

3.0

3.5

4.0

2.0

2.5

3.0

3.5

4.0

2.0

2.5

3.0

3.5

4.0

log10 of chain 1 length in stage 2

Figure 4: Estimates of the asymptotic variance of dˆµ (upper panels) and Êµ (X) (lower panels) in
stage 2, with naive weight a[2] = (.5, .5).

E.1

Choice of stage 1 weights

For stage 1, we recommend obtaining a close-to-optimal weight â[1,opt] using a pilot study described in Doss and Tan (2014). In short, one can generate samples of small size from π1 and
π2 , estimate dˆ and its asymptotic variance based on Theorem 1 for a grid of weights, and then
identify the weight that minimizes the estimated variance. With a small pilot study based on samples of size 1000 from both distributions, we obtained â[1,opt] = (0.82, 0.18). As depicted by the
horizontal lines accross the pictures in Figure 1 of the main text, the asymptotic variance of the
estimator dˆ based on â[1,opt] is approximately 0.07, which is more than 30% smaller than of the
estimator based on the naive choice a[1] = (.5, .5). Note that the naive weight is proportional to
the sample sizes from π1 and π2 , which is asymptotically optimal if both samples were independent. However, since sample 2 is from a Markov chain sample, using a weight that appropriately
favors the independent sample has lead to smaller error in the estimator. The gain in efficiency
using a close-to-optimal weight will be more pronounced if the difference in the mixing rates of
36

the two samples is larger.

E.2

Choice of stage 2 weights

In stage 2, for each µ ∈ M , the asymptotic variance of dˆµ and Êµ (X) are minimized at different
weights. Instead of searching for each of the 2|M | optimal weights in a pilot study, it is more
practical to set sub-optimal weights using less costly strategies. Below, we perform a simulation
study to examine three simple weighting strategies:
1. naive: a[2] ∝ (n1 , n2 ),
2. inverse distance (inv-dist): a[2] (µ) ∝



n1
, n2
|µ−µ1 | |µ−µ2 |



,
[2]

3. effective sample size (ess) by inverse distance (inv-dist): a (µ) ∝



ess1
, ess2
|µ−µ1 | |µ−µ2 |



.

Using each of the three strategies, we construct generalized IS estimators for dµ and Eµ (X)
for a grid of µ values between −1.5 and 4. Note that samples are drawn from two reference
distributions indexed by µ = 1 and µ = 0 respectively. Hence our simulation study concerns
both interpolation and extrapolation. A summary of their performance is provided in Figure 5,
and detailed results for selected simulation setups are shown in Figures 4, 6, and 7 for strategies 1,
2, and 3, respectively. Figure 5 suggests that none of the three strategies is uniformly better than
the others. In particular, we observe the following.
1. For estimating dµ
(a) For µ ∈ (0, 1), strategy 2 works the best.
(b) For µ = 0, strategies 2 and 3 work better than strategy 1. Indeed, both of them simply
ˆ This would be a better
set their stage 2 estimates dˆ0 to be the stage 1 estimate, d.
choice than strategy 1 because in a two-step procedure, stage 1 chains are often much
longer than stage 2 chains, and hence dˆ is already a very accurate estimate for d0 = d.
(c) For µ ∈
/ [0, 1], strategies 2 and 3 generally lead to more stable estimates of dµ . However, all strategies lead to very large asymptotic variances for µ < 0. Hence, one
needs to be mindful when doing extrapolation with IS estimators — always obtain an
estimate of the standard error, or reconsider the placement of the reference points.
2. For estimating Eµ (X)
(a) For µ ∈ (0, 1), strategy 2 works the best in general, while strategy 3 is very unstable.
37

^
relative asy var for dµ

naive
inv−dist
ess by inv−dist

1.0

0.5

1.5

1.0

2.0

1.5

2.5

^
relative asy var for Eµ(X)

0.0

naive
inv−dist
ess by inv−dist
−1

0

1

2

3

4

−1

2

3

weight for chain 1 (the iid sample)

●

●

0.8

0.8

●

●

●

4

1.0

weight for chain 1 (the iid sample)
1.0

µ

●

●
●

●
●

●

●

●

●
●

●
●

●

●

●

●

●

●

0.4

●

0.4

●

●
●

0.6

0.6

1

µ

●

●

0

0.0

●

−1

0

0.2

naive
inv−dist
ess by inv−dist
1

2

3

naive
inv−dist
ess by inv−dist

0.0

0.2

●
●

4

−1

µ

0

1

2

3

4

µ

Figure 5: Comparisons of three weight strategies in terms of the asymptotic variance of the corresponding estimators dˆµ and Êµ (X). The solid dots show which strategy achieves the smallest
asymptotic variance among the three at any given µ (ties awarded to the more basic strategy).

(b) For either µ = 0 or 1, strategy 2 and 3 are the same, and they only utilize the reference
chain from µ. This was a wise choice for estimating dµ as explained before, but not
so for other quantities of interest.
(c) For µ ∈
/ [0, 1], all strategies lead to fairly large asymptotic variances, especially for
µ < 0.
Overall in stage 2, strategy 2 has an advantage when the estimands are ratios between normalizing constants. However, when estimating Eµ (X), the situation is more complicated. Our
impression is that assigning any extreme weight will lead to high variability in the estimator. So
it is reasonable to simply use the naive weight, or other strategies that bound the weights away
38

µ= 2
BM
RS
empirical

2.0

BM
RS
empirical

0.02

0.5

5

1.0

0.04

10

1.5

15

BM
RS
empirical

0.06

20

µ= 0.5

0

^
estimates of asy var for d

µ= −1

2.0

2.5

3.0

3.5

4.0

2.0

2.5

3.0

3.5

4.0

2.0

2.5

3.0

3.5

4.0

log10 of chain 1 length in stage 2

BM
RS
empirical

4

0

2

2

60

8

BM
RS
empirical

6

8

µ= 2

4

100

BM
RS
empirical

6

140

µ= 0.5

0 20

^
estimates of asy var for Eµ(X)

µ= −1

2.0

2.5

3.0

3.5

4.0

2.0

2.5

3.0

3.5

4.0

2.0

2.5

3.0

3.5

4.0

log10 of chain 1 length in stage 2

Figure 6: Estimates of the asymptotic variance of dˆµ (upper panels) and Êµ (X) (lower panels) in
stage 2, with weight a[2] (µ) chosen by strategy 2.
from 0 and 1.

F

Bayesian variable selection models

Here, we consider a class of Bayesian variable selection (BVS) models for linear regression with
independent normal priors on the regression coefficients. This model involves a 2-dimensional
prior hyperparameter that influences inference, yet no default choice guarantees good performance in practice. Hence, displaying the effect of different hyperparameter values on the posterior distribution would greatly benefit users of the model. When the number of predictors, q, is
large, the computing is challenging. Our solution is to obtain MCMC samples for a small number of models with different hyperparameter values, based on which generalized IS estimates can
be obtained for BFs and other posterior expectations for a large number of models. Again, an
important problem in practice is how long the Markov chains need to be run? In this context, the
39

2.0
0.5

2

0.05

4

1.0

6

BM
RS
empirical

1.5

0.15

BM
RS
empirical

0.10

8

10

BM
RS
empirical

µ= 2
2.5

µ= 0.5

0

^
estimates of asy var for d

µ= −1

2.0

2.5

3.0

3.5

4.0

2.0

2.5

3.0

3.5

4.0

2.0

2.5

3.0

3.5

4.0

log10 of chain 1 length in stage 2

10

µ= 2
BM
RS
empirical

2

1

4

2

6

8

BM
RS
empirical

4

BM
RS
empirical

3

10 20 30 40 50 60

µ= 0.5

0

0

^
estimates of asy var for Eµ(X)

µ= −1

2.0

2.5

3.0

3.5

4.0

2.0

2.5

3.0

3.5

4.0

2.0

2.5

3.0

3.5

4.0

log10 of chain 1 length in stage 2

Figure 7: Estimates of the asymptotic variance of dˆµ (upper panels) and Êµ (X) (lower panels) in
stage 2, with naive weight a[2] (µ) chosen by strategy 3.
only affordable method that we are aware of is to estimate the SE of these IS estimators using the
proposed BM method.
As introduced by Mitchell and Beauchamp (1988), let Y = (Y1 , . . . , Ym )> denote the vector
of responses and X1 , . . . , Xq denote q potential predictors, each a vector of length m. The predictors are standardized, so that for j = 1, . . . , q, 1Tm Xj = 0 and XjT Xj = m, where 1m is the
vector of m 1’s. The BVS model is given by:
given γ, σ 2 , β0 , βγ ,
given γ, σ 2 , β0 ,
given γ,

Y ∼ Nm (1m β0 + Xγ βγ , σ 2 I),
γj 
ind
βj ∼ N 0, σ 2 for j = 1, . . . , q,
λ
(σ 2 , β0 ) ∼ p(β0 , σ 2 ) ∝ 1/σ 2 ,
γ ∼ p(γ) = wqγ (1 − w)q−qγ .

(F.1a)
(F.1b)
(F.1c)
(F.1d)

The binary vector γ = (γ1 , . . . , γq )> ∈ {0, 1}q identifies a subset of predictors, such that Xj is
40

P
included in the model if and only if γj = 1, and |γ| = qj=1 γj denotes the number of predictors
included. So (F.1a) says that each γ corresponds to a model given by Y = 1m β0 + Xγ βγ +  ,
where Xγ is an n × |γ| sub-matrix of X that consists of predictors included by γ, βγ is the vector
that contains corresponding coefficients, and  ∼ Nm (0, σ 2 I). It is sometimes more convenient
to use the notation, Y = X0γ β0γ + , where X0γ has one more column of 1’s than Xγ and
T
β0γ
= (β0 , βγT ). Unknown parameters are θ = (γ, σ, β0 , βγ ) for which we set a hierarchical prior
in (F.1b) to (F.1d). In (F.1d), an independent Bernoulli prior is set for γ, where w ∈ (0, 1) is a
hyperparameter that reflects the prior inclusion probability of each predictor. In (F.1c), a noninformative prior is set for (σ 2 , β0 ). In (F.1b), an independent normal prior is assigned to βγ ,
where λ > 0 is a second hyperparameter, that controls the precision of the prior. Overall, θ is
given an improper prior due to (F.1c) but the posterior of θ is indeed proper.
One can actually integrate out (βγ , β0 , σ 2 ) and arrive at the following model with parameter
γ only:
Z Z Z
f (Y |γ, σ 2 , β0 , βγ )f (βγ |γ, σ 2 , β0 )f (σ 2 , β0 )dβγ dβ0 dσ 2
Y |γ ∼ `h (γ; Y ) =
R+

= cm λ
γ ∼ ph (γ)

=

R
|γ|
2

R|γ|

A0γ

− 21 

(Y − Y )T (Y − Y ) − βeγT A0γ βeγ

−(m−1)/2

,

(F.2)

wqγ (1 − w)q−qγ .

T
Here, cm is a constant depending only on the sample size m. Further, A0γ = X0γ
X0γ + Λ0γ ,
where Λ0γ is a diagonal matrix, the main diagonal of which is the (1 + |γ|)-dimensional vector
T
(0, λ, · · · , λ). Finally, βeγ = A−1
0γ X0γ Y .
Using the model at (F.2) requires specification of the hyperparameter h = (w, λ). Smaller w
values assign high prior probabilities to models with fewer predictors, and priors with smaller λ
values allow selected predictors to have large coefficients. It is common to set w = 0.5 (a uniform
prior on the model space) and λ = 1 (a unit information prior for uncorrelated predictors, see e.g.
Kass and Raftery (1995)). One can also choose h adaptively, say according to the marginal likeliP
hood mh = γ `h (γ; Y )ph (γ). A small value of mh indicates that the prior ph is not compatible
with the observed data, while hEB = arg max mh is defined to be the empirical Bayesian choice
of h. The empirical Bayes idea has been successfully applied to various models with variable
selection components (see e.g. George and Foster (2000); Yuan and Lin (2005)). However, we
have not seen this idea being carried out for the model in (F.1), except where n = p and the design matrices are orthogonal (Clyde and George (2000); Johnstone and Silverman (2005)). Due
to the improper prior in (F.1d), mh is not uniquely defined. Nevertheless, the Bayes factor among
any two models, say mh /mh0 , is well-defined because the same improper prior is assigned to the
shared parameters of the two models (see e.g. Kass and Raftery (1995, sec.5) and Liang et al.

41

(2008, sec.2)).
Here, we concentrate on two goals. The first is to evaluate {mh /mh1 , h ∈ H}, the marginal
likelihood of model h relative to a reference model h1 , which allows us to identify the empirical
Bayesian choice of h. The second is to evaluate the posterior mean of the vector of coefficients
β for each h ∈ H, which we denote by bh . Predictions can then be made for new observations
using Y (new) = (x(new) )T bh .
For model (F.2) with a fixed h, a Metropolis Hastings random-swap algorithm (Clyde et al.
(2011)) can be used to generate Markov chains of γ from its posterior distribution. In each
iteration, with probability ρ(γ), we propose flipping a random pair of 0 and 1 in γ, and with
probability 1 − ρ(γ), we propose changing γj to 1 − γj for a random j while leaving other
coordinates untouched. We set ρ(γ) = 0 when γ corresponds to the null model or the full model,
and ρ(γ) = .5 otherwise. Finally, the proposal is accepted with an appropriate probability. Since
this Markov chain lies on a finite state space, it is uniformly ergodic and hence polynomially
ergodic as well. Further, moment conditions in Theorems 2 and 3 are satisfied because they
reduce to summations of 2q terms, a large but finite number. To achieve the goal, we generate
Markov chains of γ with respect to model (F.2) at several h values that scatters in H, from which
we build generalized IS estimators, and estimate their standard errors.

F.1

Cookie dough data

We demonstrate the aforementioned sensitivity analysis using the biscuit dough dataset (Brown
et al. (2001); Osborne et al. (1984)). The dataset, available in the R package ppls (Kraemer and
Boulesteix (2012)), contains a training set of 39 observations and a test set of 31 observations.
These data were obtained from a near-infrared spectroscopy experiment that study the composition of biscuit dough pieces. For each biscuit, the reflectance spectrum is measured at 700 evenly
spaced wavelengths. We use these measurements as covariates to predict the response variable,
the percentage of water in each dough. We follow previous studies (Hans (2011)) and thin the
spectral to q = 50 evenly spaced wavelengths.
Figure 8 provides a general picture of the sensitivity analysis. The left plots provide two ways
to visualize estimates of the BFs. To form the plot, we took the 12 reference values of h = (w, λ)
to be such that (w, − log(λ)) ∈ {0.1, 0.2, 0.3, 0.4} × {1, 3, 5}. In stage 1 we ran each of the
ˆ In stage 2, we ran the
12 Markov chains at the above values of h for 105 iterations to obtain d.
same 12 Markov chains for 50, 000 iterations each, to form the estimates ûn over a fine grid that
consists of 475 different h values, with the w component ranging from 0.05 to 0.5 in increments
of 0.025 and the − log(λ) component ranging from 0 to 6 in increments of 0.25.
42

standard error relative to estimates of BF
6

BF estimates

0.03
4
0.0

0.02

5

8

3
0.2

2
1

0.4

0

0.3
w

− lo
g(
λ)

0.1

0.0

1

6
5
4

0.05

2

0.01

3

4

2

− log (λ)

3

4

6

0.1

0.5 0

0.2

0.3

0.4

0.5

w

6

contour of prediction mse
0.4

1

4

5

5

2
0.1
25

0.5

0.35

4
0.5

2

2

5

0.25

0.125

0.4

0.3

3

3

− log (λ)

4

0.45

3
0.

− log (λ)

0.55

5

8

0.4

6

contour of BF estimates

0.4

0

0

0.6

0.45

1

0.5

1

0.4

0.1

0.2

0.3

0.4

0.5

0.1

w

0.2

0.3

0.4

0.5

w

Figure 8: Left panels provide a surface plot and a contour plot for BF estimates. The upper-right
panel displays standard errors with respect to the BF estimates. The lower-right panel shows
pmse over the test set.
How trustworthy are these BF estimates? Their estimated standard errors are obtained using
the BM method, based on Theorem 2. We choose to display the relative SE with respect to the
BF estimates, as shown in the upper right panel of Figure 8. The relative SEs are smaller than or
equal to 5%, and we believe the BF estimates are accurate enough. Finally, the lower-right panel
of Figure 8 shows the prediction mean squared error (pmse) over the test set for all h.
Based on our estimation, the BF attains the maximum value 9.75 at hEB = (0.075, e−5 ).
Recall when comparing any two models indexed by h and h0 respectively, the BF between them
m /m
is given by BFh,h0 = m h0 /mhh1 . Also, according to Jeffreys (1998) and Kass and Raftery (1995),
h

1

43

5

Standard method, (w, − log (λ) ) = (0.5, 0)

5

EB method, (w, − log (λ) ) = (0.075, 5)

−5

−5

β

0

pmse=0.431

0

pmse=0.411

1200

1600

2000

2400

1200

wavelength

1600

2000

2400

wavelength

Figure 9: Estimated posterior mean of regression coefficients (with 95% point-wise confidence
intervals) for the empirical Bayesian and the standard choice of h, respectively.
the evidence for h over h0 is considered to be strong only if BFh,h0 is greater than 10 or 20. Hence,
all h with BF over 1/10 or 1/20 times the maximum BF can be considered as reasonably well
supported by the data as that of the empirical Bayesian choice. Comparing the lower two plots
of Figure 8, we see that the set Ac := {h ∈ H : BF > c} for c = 1 and 0.5 do overlap with an
area that corresponds to relatively small pmse. Outside A0.5 , a region that consists of larger w
and smaller − log(λ) also enjoys small pmse values, at around 0.3 to 0.4. This region includes
the common choice of h0 = (0.5, e0 ). These suggest that hEB and its vicinity might not be the
only area of h that has good prediction performances.
To better compare the effect of hEB = (0.075, e−5 ) and the commonly used h = (0.5, e0 ),
Figure 9 displays the estimated posterior mean of regression coefficients at both choices of h,
together with the point-wise 95% confidence intervals for the posterior means. Due to the small
size of wEB and λEB , the empirical Bayesian method yields a model with a few covariates that
have big coefficients. In comparison, the common choice has larger w and λ values, leading to
a regression model that combines more covariates each having smaller effects. It turns out these
two opposite strategy of modeling both predict the test dataset well, with pmse being 0.411 and
0.431 respectively. For comparison, pmses were calculated for several frequentist penalized linear
regression methods with their respective penalty parameters chosen by 10-fold cross validation.
The resulting pmses for the ridge, the lasso and the elastic net method are 4.675, 0.633 and 0.536,
respectively.
The BM method for estimating SE is carried out above without the need of further user
44

input. Theoretically, its competitor RS can be developed too, if enough regeneration times can
be identified for each Markov chain. Recall that with the random-swap algorithm, each Markov
chain lives on the discrete state space Γ of size 2q . A naive way to introduce regeneration is to
specify a single point γ1 , then each visit of the Markov chain to γ1 marks a regeneration time.
Note that the chance of visiting γ1 converges to π(γ1 ), the posterior probability of γ1 . In our
BVS model with 250 ≈ 1.1 × 1015 states of γ, even maxγ∈Γ π(γ) could be very small. Take for
example the Markov chain for the BVS model with h = (w, log(g)) = (.4, 1), the point with
the highest frequency appeared only 8 out of a run of 104 iterations. And that the waiting times
between consecutive regenerations are highly variable, which ranges from less than ten iterations
to a few thousand iterations. To obtain alternative schemes of identifying regeneration times, one
can take the general minorization condition approach. It could potentially increase the chance of
regeneration and reduce variability of the waiting times. Specifically, for any α ∈ {1, 2, · · · , 2q },
one could define Dα to contain the α points with the highest posterior probabilities, and find
α ∈ (0, 1] and a probability mass function kα (·) such that p(γ 0 |γ) ≥ α IDα (γ)kα (γ 0 ) for all
γ 0 ∈ Γ. Note that as α increases, the chance of visiting Dα improves, but α , the conditional rate
of regeneration given the current state γ is in Dα , would decrease sharply. Finding a good α to
maximize the overall chance of regeneration requires tuning that is specific for both the model
specification h and the dataset. Even if we can find the optimal α for each Markov chain used in
the example, it is unlikely that all of them would regenerate often enough for the RS estimator to
be stable.

References
Brown, P. J., Fearn, T. and Vannucci, M. (2001). Bayesian wavelet regression on curves with
application to a spectroscopic calibration problem. J. Amer. Statist. Assoc. 96 398–408.
Buta, E. and Doss, H. (2011). Computational approaches for empirical Bayes methods and
Bayesian sensitivity analysis. Ann. Statist. 39 2658–2685.
Clyde, M. and George, E. I. (2000). Flexible empirical Bayes estimation for wavelets. J. R. Stat.
Soc. Ser. B. Stat. Methodol. 62 681–698.
Clyde, M. A., Ghosh, J. and Littman, M. L. (2011). Bayesian adaptive sampling for variable
selection and model averaging. J. Comput. Graph. Statist. 20.
Cappé, O. and Guillin, A. and Marin, J. M. and Robert, C. P.(2004). Population Monte Carlo. J.
Comput. Graph. Statist. 13 907–929.
45

Christensen, O. F. (2004). Monte Carlo maximum likelihood in model based geostatistics. J.
Comput. Graph. Statist. 13 702–718.
Diggle, P. J., Tawn, J. A. and Moyeed, R. A. (1998). Model-based geostatistics. Applied Statistics 47, 299–350.
Doss, H. (2010). Estimation of large families of Bayes factors from Markov chain output. Statist.
Sinica, 20 537–560.
Doss, H. and Tan, A. (2014). Estimates and standard errors for ratios of normalizing constants
from multiple Markov chains via regeneration. J. R. Stat. Soc. Ser. B. Stat. Methodol. 76 683–
712.
Elvira, V. and Martino, L. and Luengo, D. and Bugallo, M. F. (2015). Generalized multiple importance sampling. ArXiv e-prints.
Evangelou, E. and Roy, V. (2015). geoBayes. http://cran.r-project.org/web/packages/geoBayes.
R package version 0.3-3.
Flegal, J. M. and Haran, M. and Jones, G. L. (2008). Markov chain Monte Carlo: Can we trust
the third significant figure? Statist. Sci. 23 250–260.
Flegal, J. M. and Jones, G. L. (2010). Batch means and spectral variance estimators in Markov
chain Monte Carlo. Ann. Statist., 38:1034–1070.
George, E. and Foster, D. (2000). Calibration and empirical Bayes variable selection. Biometrika
87 731.
Geyer, C. J. (1994). Estimating normalizing constants and reweighting mixtures in Markov chain
Monte Carlo. Tech. Rep. 568r, Department of Statistics, University of Minnesota.
Gilks, W. R., Roberts, G. O. and Sahu, S. K. (1998) Adaptive Markov chain Monte Carlo through
regeneration. J. Amer. Statist. Assoc. 93 1045–1054.
Gill, R. D., Vardi, Y. and Wellner, J. A. (1988). Large sample theory of empirical distributions in
biased sampling models. Ann. Statist. 16 1069–1112.
Hans, C. (2011). Elastic net regression modeling with the orthant normal prior. J. Amer. Statist.
Assoc. 106 1383–1393.

46

Hastings, W. K. (1970). Monte Carlo sampling methods using Markov chains and their applications. Biometrika 57 97–109.
Jeffreys, H. (1998). The theory of probability. 3rd ed. Oxford.
Johnstone, I. M. and Silverman, B. W. (2005). Empirical Bayes selection of wavelet thresholds.
Ann. Statist. 33 1700–1752.
Jones, G. L. (2004). On the Markov chain central limit theorem. Probab. Surv., 1:299–320.
Johns, G. and Hobert, J. (2004). Sufficient Burn-in for Gibbs Samplers for a Hierarchical Random
Effects Model. Ann. Statist. 32 784–817.
Kass, R. E. and Raftery, A. E. (1995). Bayes factors. J. Amer. Statist. Assoc. 90 773–795.
Koehler, E. and Brown, E. and Haneuse, S. (2009). On the assessment of Monte Carlo error in
simulation-based statistical analyses. Amer. Statist. 63 155–162.
Kong, A., McCullagh, P., Meng, X.-L., Nicolae, D. and Tan, Z. (2003). A theory of statistical
models for Monte Carlo integration (with discussion). J. R. Stat. Soc. Ser. B. Stat. Methodol.
65 585–618.
Kraemer, N. and Boulesteix, A. (2012). ppls: Penalized partial least squares, r package version
1.05.
Liang, F., Paulo, R., Molina, G., Clyde, M. A. and Berger, J. O. (2008). Mixtures of g-priors for
Bayesian variable selection. J. Amer. Statist. Assoc. 103 410–423.
Meng, X.-L. and Wong, W. H. (1996). Simulating ratios of normalizing constants via a simple
identity: A theoretical exploration. Statist. Sinica 6 831–860.
Mengersen, K. L. and Tweedie, R. L. (1996). Rates of convergence of the Hastings and Metropolis algorithms. Ann. Statist. 24 101–121.
Mitchell, T. and Beauchamp, J. (1988). Bayesian variable selection in linear regression. J. Amer.
Statist. Assoc. 83 1023–1036.
Meyn, S. P. and Tweedie, R. L. (1993). Markov Chains and Stochastic Stability. Springer-Verlag,
New York, London.

47

Mykland, P., Tierney, L. and Yu, B. (1995). Regeneration in Markov chain samplers. J. Amer.
Statist. Assoc. 90 233–41.
Osborne, B. G., Fearn, T., Miller, A. R. and Douglas, S. (1984). Application of near infrared
reflectance spectroscopy to the compositional analysis of biscuits and biscuit doughs. Journal
of the Science of Food and Agriculture 35 99–105.
Owen A. and Zhou, Y. (2000). Safe and effective importance sampling. J. Amer. Statist. Assoc.
95 135-143.
Roy, V. (2014) Efficient estimation of the link function parameter in a robust Bayesian binary
regression model. Comput. Stat. Data. Anal., 73 87–102.
Roy, V. and Evangelou, E. and Zhu, Z. (2016). Efficient estimation and prediction for the
Bayesian binary spatial model with flexible link functions. Biometrics 72 289–298.
Roy, V. and Hobert, J. P. (2007) Convergence rates and asymptotic standard errors for MCMC
algorithms for Bayesian probit regression. J. R. Stat. Soc. Ser. B. Stat. Methodol. 69 607–623.
Tan, Z. (2004). On a likelihood approach for Monte Carlo integration. J. Amer. Statist. Assoc. 99
1027–1036.
Tan, A. and Doss, H. and Hobert, J. P. (2015). Honest importance sampling with multiple Markov
chains. J. Comput. Graph. Statist., 24 792–826.
Tan, A. and Hobert, J. P. (2009). Block Gibbs sampling for Bayesian random effects models with
improper priors: convergence and regeneration. J. Comput. Graph. Statist. 18 861-878.
Tanner, M. A. and Wong, W. H. (1987) The calculation of posterior distributions by data augmentation(with discussion). J. Amer. Statist. Assoc. 82 528–550.
Vardi, Y. (1985). Empirical distributions in selection bias models. Ann. Statist. 13 178–203.
Vats, D., Flegal, J. M., and Jones, G. L.. (2015a). Multivariate output analysis for Markov chain
Monte Carlo. ArXiv e-prints.
Vats, D., Flegal, J. M., and Jones, G. L.. (2015b). Strong consistency of the multivariate spectral
variance estimator in Markov chain Monte Carlo. ArXiv e-prints.
Veach, E. and Guibas, L. (1995). Optimally combining sampling techniques for Monte Carlo
rendering. SIGGRAPH 95 Conference Proceedings, Reading MA. Addison-Wesley, 419-428.
48

Yuan, M. and Lin, Y. (2005). Efficient empirical Bayes variable selection and estimation in linear
models. J. Amer. Statist. Assoc. 100 1215–1225.
Zhang, H. (2002). On estimation and prediction for spatial generalized linear mixed models.
Biometrics 58 129–136.
Department of Statistics, Iowa State University
E-mail: vroy@iastate.edu
Department of Statistics and Actuarial Science, University of Iowa
E-mail: aixin-tan@uiowa.edu
Department of Statistics, University of California, Riverside
E-mail: jflegal@ucr.edu

49

