# Selection of proposal distributions for generalized importance sampling
estimators

Vivekananda Roy
Department of Statistics, Iowa State University, USA
and
Evangelos Evangelou
Department of Mathematical Sciences, University of Bath, UK

###### Abstract

The standard importance sampling (IS) estimator, generally does not work well
in examples involving simultaneous inference on several targets as the
importance weights can take arbitrarily large values making the estimator
highly unstable. In such situations, alternative generalized IS estimators
involving samples from multiple proposal distributions are preferred. Just
like the standard IS, the success of these multiple IS estimators crucially
depends on the choice of the proposal distributions. The selection of these
proposal distributions is the focus of this article. We propose three methods
based on (i) a geometric space filling coverage criterion, (ii) a minimax
variance approach, and (iii) a maximum entropy approach. The first two methods
are applicable to any multi-proposal IS estimator, whereas the third approach
is described in the context of Doss,’s (2010) two-stage IS estimator. For the
first method we propose a suitable measure of coverage based on the symmetric
Kullback-Leibler divergence, while the second and third approaches use
estimates of asymptotic variances of Doss,’s (2010) IS estimator and Geyer,’s
(1994) reverse logistic estimator, respectively. Thus, we provide consistent
spectral variance estimators for these asymptotic variances. The proposed
methods for selecting proposal densities are illustrated using various
detailed examples.

Key words: Bayes factor, Markov chain CLT, marginal likelihood, multiple
importance sampling, polynomial ergodicity, reverse logistic regression.

##  1 Introduction

Importance sampling (IS) is a popular Monte Carlo procedure where samples from
one distribution are weighted to estimate features of other distributions. To
set the framework, let $\Pi$ be the family of target densities on X with
respect to a measure $\mu$ where if $\pi(x)=\nu(x)/\theta\in\Pi$ then $\nu(x)$
is known but $\theta=\int_{\textsf{X}}\nu(x)\mu(dx)$ is unknown. Suppose the
goal is to estimate the normalizing constants $\theta$ for all $\pi\in\Pi$.
Let $q_{1}(x)=\varphi_{1}(x)/c_{1}$ be another density on X with respect to
$\mu$ such that $\nu(x)=0$ whenever $\varphi_{1}(x)=0$. If
$\\{X_{i}\\}_{i=1}^{n}$ is either iid samples from $q_{1}$ or a positive
Harris Markov chain with invariant density $q_{1}$, then the IS estimator

$\frac{1}{n}\sum_{i=1}^{n}\frac{\nu(X_{i})}{\varphi_{1}(X_{i})}\buildrel\text{a.s.}\over{\longrightarrow}\int_{\textsf{X}}\frac{\nu(x)}{\varphi_{1}(x)}q_{1}(x)\,\mu(dx)=\frac{\theta}{c_{1}}\int_{\textsf{X}}\frac{\nu(x)/\theta}{\varphi_{1}(x)/c_{1}}q_{1}(x)\,\mu(dx)=\frac{\theta}{c_{1}}.$
(1.1)

To illustrate the method, consider the popular autologistic models (Besag,,
1974), which are particular Markov random field models for binary
observations. Let $s_{i}$ denote the $i$th spatial location, and let
$\mbox{nb}_{i}\equiv\\{s_{j}:s_{j}\mbox{ is a neighbor of}\;s_{i}\\}$ denote
the neighborhood set of $s_{i},i=1,\dots,m$. Markov random field models for
$\boldsymbol{x}=\\{x(s_{i}),i=1,\dots,m\\}$ are formulated by specifying the
conditional probabilities $p_{i}=P(x(s_{i})=1|\\{x(s_{j}):j\neq
i\\})=P(x(s_{i})=1|\\{x(s_{j}):s_{j}\in\mbox{nb}_{i}\\}),i=1,\dots,m$. For
simplicity, we assume that all neighborhoods have the same size
$w=|\mbox{nb}_{i}|,i=1,\dots,m$. Here, we consider a centered parameterization
(Kaiser et al.,, 2012) given by
$\mbox{logit}(p_{i})=\mbox{logit}(\kappa)+(\gamma/w)\sum_{s_{j}\in\mbox{nb}_{i}}(x(s_{j})-\kappa)$,
where $\mbox{logit}(z)=\log(z/[1-z])$, $\gamma$ is the dependence parameter,
and $\kappa$ represents the probability of observing one in the absence of
statistical dependence. The autologistic model is defined via conditional
distributions, and the joint pmf $\pi(\boldsymbol{x}|\gamma,\kappa)$ of
$\boldsymbol{x}$ is given by (derived in Section J of the supplementary
materials)

$\pi(\boldsymbol{x}|\gamma,\kappa)\propto\exp\big{\\{}(\mbox{logit}(\kappa)-\gamma\kappa)\sum\nolimits_{i=1}^{m}x(s_{i})+\frac{\gamma}{2w}\sum\nolimits_{i=1}^{m}\sum\nolimits_{s_{j}\in\mbox{nb}_{i}}x(s_{i})x(s_{j})\big{\\}}.$
(1.2)

The normalizing constant $\theta\equiv\theta(\gamma,\kappa)$ in
$\pi(\boldsymbol{x}|\gamma,\kappa)$ is intractable when $\gamma\neq 0$.
Sherman et al., (2006) mention that ‘there is no known simple way to
approximate this normalizing constant’.

We consider a $10\times 10$ square lattice on torus, with four-nearest (east-
west, north-south) neighborhood structure. We consider the family of
autologistic pmfs
$\Pi=\\{\pi(\boldsymbol{x}|\gamma,\kappa=0.5),\gamma=-4,-3.9,\dots,3.9,4\\}$.
We estimate the ratios of normalizing constants
$\theta(\gamma,0.5)/\theta(0,0.5)$, by the naive IS (NIS) estimator (1.1) with
$n=60,000$ (Gibbs sampler) samples from
$q_{1}\equiv\pi(\boldsymbol{x}|\gamma=0,\kappa=0.5)$. Figure 1 provides the
relative standard error (SE) of (1.1). We use the spectral variance (SV)
method provided in section 3.2 for estimating the SE. The NIS estimator (1.1)
suffers from high relative SE when the target pmf is not ‘close’ to the
proposal pmf (Geyer,, 2011). Indeed, when $\gamma$ is away from $0$, the ratio
$\nu(X_{i})/\varphi_{1}(X_{i})$ takes large values for some samples $X_{i}$’s,
and it results in a highly variable IS estimator.

Figure 1: Relative standard errors of naive and multiple proposal IS
estimators corresponding to different proposal distributions for the
autologistic model. The $\gamma$ values corresponding to the proposal
distributions are denoted by $+$.

There are several modifications of the simple IS estimator (1.1) available in
the literature, e.g., multiple importance sampling (Veach and Guibas,, 1995;
Owen and Zhou,, 2000; Elvira et al.,, 2015), umbrella sampling (Geyer,, 2011;
Doss,, 2010), parallel, serial or simulated tempering (George and Doss,, 2018;
Geyer and Thompson,, 1995; Marinari and Parisi,, 1992) which utilize samples
from more than one importance distribution such that any pdf in the given set
of target densities is close to at least one of the importance densities. In
IS estimation based on multiple proposal densities, the single density $q_{1}$
is generally replaced with a linear combination of $k$ densities (Elvira et
al.,, 2015; Geyer,, 2011). In particular, let $q_{i}(x)=\varphi_{i}(x)/c_{i}$
for $i=1,\dots,k$ be $k$ proposal densities, where $\varphi_{i}$’s are known
but $c_{i}$’s may be unknown, $\boldsymbol{a}=(a_{1},\dots,a_{k})$ be $k$
positive constants such that $\sum_{i=1}^{k}a_{i}=1$, and
$\overline{q}\equiv\sum_{i=1}^{k}a_{i}q_{i}$. Let $d_{i}=c_{i}/c_{1}$ for
$i=1,2,\dots,k$ with $d_{1}=1$, and
$\boldsymbol{d}\equiv(c_{2}/c_{1},\ldots,c_{k}/c_{1})$. For $l=1,\dots,k$, let
$\\{X_{i}^{(l)}\\}_{i=1}^{n_{l}}$ be either iid samples from $q_{l}$ or a
positive Harris Markov chain with invariant density $q_{l}$. Then as
$n_{l}\rightarrow\infty,$ for all $l=1,\dots,k$, we have

$\displaystyle\hat{u}\equiv\sum_{l=1}^{k}\frac{a_{l}}{n_{l}}\sum_{i=1}^{n_{l}}\frac{\nu(X_{i}^{(l)})}{\sum_{j=1}^{k}a_{j}\varphi_{j}(X_{i}^{(l)})/d_{j}}$
$\displaystyle\buildrel\text{a.s.}\over{\longrightarrow}\sum_{l=1}^{k}a_{l}\int_{\textsf{X}}\frac{\nu(x)}{\sum_{j=1}^{k}a_{j}\varphi_{j}(x)/d_{j}}q_{l}(x)\,\mu(dx)$
(1.3)
$\displaystyle=\frac{1}{c_{1}}\int_{\textsf{X}}\frac{\nu(x)}{\overline{q}(x)}\overline{q}(x)\,\mu(dx)=\frac{\theta}{c_{1}}.$

Gill et al., (1988), Kong et al., (2003), Meng and Wong, (1996), Tan, (2004),
and Vardi, (1985) considered 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; Tan et al.,,
2015), which is the setting of the current paper. There are alternative
weighting schemes proposed in the literature, although none is as widely
applicable as (1.3). For example, if the normalizing constants $c_{i}$’s are
known, the estimator (1.3) resembles the balance heuristic estimator of Veach
and Guibas, (1995), which is discussed in Owen and Zhou, (2000) as the
deterministic mixture. Alternatively, the standard population Monte Carlo
algorithm (Cappé et al.,, 2004) uses a weighted ratio of the target $\pi$ and
the proposal $q_{j}$ it was drawn from (evaluated at the sample itself).
However, if iid samples are available from $q_{j},j=1,2,\dots,k$, Elvira et
al., (2015) have shown that the normalized estimator ($c_{i}$’s known) version
of (1.3) always has a smaller variance than that of the population Monte Carlo
algorithm.

The effectiveness of (1.3) depends on the choice of $k$, $\boldsymbol{a}$,
$n_{l}$, and the importance densities $\\{q_{j},j=1,\dots,k\\}$. This article
focuses on the choice of the importance densities because it is the most
crucial, and the multiple IS estimator (1.3), just like the NIS estimator
(1.1), is useless if the importance densities are ‘off targets’. Although
increasing $k$ or $n_{l}$, may lead to estimators with less variance, it
results in higher computational cost, therefore these are often determined
based on the available computational resources. On the other hand, for fixed
$k$, $\boldsymbol{a}$, and $n_{l}$, efficiency and stability of the estimator
(1.3) can be highly improved by appropriately choosing the $k$ importance
densities $\boldsymbol{q}=\\{q_{1},\dots,q_{k}\\}$ from the set of potential
densities $Q\equiv\\{q(x)=\varphi(x)/c\\}$ on X with respect to $\mu$.
Unfortunately, in the literature, there is not much discussion on the choice
of the importance densities in multiple IS methods. One exception is Buta and
Doss, (2011) who described an ad-hoc method in the important special case when
$Q=\Pi$, that is, when samples from a subset of densities of $\Pi$ are used to
estimate the normalizing constants for the entire family via (1.3). Routine
applications of IS estimation with $Q=\Pi$ can be found in Monte Carlo maximum
likelihood estimation, Bayesian sensitivity analysis and model selection
(Geyer and Thompson,, 1992; Buta and Doss,, 2011; Doss,, 2010). We denote the
asymptotic variance of (1.3) by $\sigma^{2}_{u}(\pi;\boldsymbol{q})$, which
depends on the choice of the proposal distributions $\boldsymbol{q}$. An ideal
choice of $\boldsymbol{q}$ would be the one that minimizes the maximum of some
measure of discrepancy, $\phi(\boldsymbol{q})$, e.g., the asymptotic variance
or relative variance. Buta and Doss, (2011) stated that solving the minimax
variance design problem, that is, the one that minimizes
${\mbox{max}}_{\pi\in\Pi}\sigma^{2}_{u}(\pi;\boldsymbol{q})$ exactly is
‘hopeless’, and instead described a procedure where starting from some ‘trial’
proposal pdfs, $\hat{\sigma}^{2}_{u}(\pi;\boldsymbol{q})$ is computed for all
$\pi\in\Pi$. Here $\hat{\sigma}^{2}_{u}$ is a consistent estimator of the
asymptotic variance $\sigma^{2}_{u}$. Then proposal densities are either moved
to regions of $\Pi$ where $\hat{\sigma}^{2}_{u}$ is large, or new proposal
densities from these high variance regions are added increasing $k$. We
implemented Buta and Doss,’s (2011) method with relative SE criterion for the
autologistic example mentioned before. Since the relative SE for the IS
estimator with $k=1$ takes large values when $|\gamma|=4$, following Buta and
Doss, (2011), we added two new proposal pmfs
$\\{\pi(\boldsymbol{x}|\gamma,\kappa=0.5),\gamma=\pm 4\\}$ increasing $k$ to
3\. In this case, 20,000 samples were used from each of the 3 proposal pmfs to
keep the total sample size the same as in the NIS method. From Figure 1 we see
that Buta and Doss,’s (2011) method substantially decreases the maximum
relative SE. A contribution of this paper is the derivation of a consistent
estimator of $\sigma^{2}_{u}(\pi,\boldsymbol{q})$. This allows us to directly
compute a minimax optimal design based on this estimator, i.e., the one that
minimizes
${\mbox{max}}_{\pi\in\Pi}[\hat{\sigma}_{u}(\pi;\boldsymbol{q})/\hat{u}(\pi;\boldsymbol{d})]$
where $\hat{u}(\pi;\boldsymbol{d})\equiv\hat{u}$ defined in (1.3). The
proposed design with $k=3$, shown in Figure 1, reduces the relative SE
(uniformly over $\gamma$) compared to Buta and Doss,’s (2011) method for this
example.

This paper is the first where systematic methods of selection of proposal
distributions for IS are developed and tested. We propose three approaches.
Our first approach is based on a geometric spatial design method. In
particular, among all subsets $\boldsymbol{q}\subset Q$ with
$|\boldsymbol{q}|=k$, the one with minimum distance from $\Pi$ is chosen. The
choice of the distance is crucial, and here we propose to use the symmetric
Kullback-Leibler divergence. The second approach is the minimax approach
mentioned before. The third approach is applicable when $\boldsymbol{d}$ in
(1.3) is unknown, which is the case in the autologistic example as the
proposal pmfs are chosen from the family of autologistic (target) models. If
$\boldsymbol{d}$ is unknown, Doss, (2010) proposed a two-stage method, where
in the first step, based on (iid or Markov chain) samples from
$q_{i},i=1,\dots,k$, $\boldsymbol{d}$ is estimated using Geyer,’s (1994)
reverse logistic estimator $\boldsymbol{\hat{d}}$ or Meng and Wong,’s (1996)
bridge sampling method. Then, independent of step one, new samples are used to
calculate (1.3) with $\boldsymbol{d}$ replaced by $\boldsymbol{\hat{d}}$.
Doss, (2010) suggested that $\boldsymbol{a}$ in (1.3) need not be proportional
to sample sizes, which was termed the generalized IS estimator by Roy et al.,
(2018). In our third approach, following the maximum entropy criterion of
experimental design, proposal distributions are chosen by maximizing
determinant of the asymptotic covariance matrix of $\boldsymbol{\hat{d}}$.
Note that, in Doss,’s (2010) two-stage IS estimation, the estimator
$\boldsymbol{\hat{d}}$ based on the first stage samples from $k$ proposal
densities is used to explore large number of $\pi$’s in the second stage. We
describe and compare these three methods in details in Section 3. These
methods are illustrated using several detailed examples involving autologistic
models, Bayesian robit models and spatial generalized linear mixed models.

As mentioned above, two of the three approaches proposed here utilize
asymptotic standard errors of $\boldsymbol{\hat{d}}$ and $\hat{u}$. A
contribution of this paper is the development of SV estimators for
$\boldsymbol{\hat{d}}$ and $\hat{u}$. Recently, Roy et al., (2018) provided
standard errors estimators of $\boldsymbol{\hat{d}}$ and $\hat{u}$ using the
batch means (BM) method. For the simulation study considered in the supplement
to the paper, we observe that the proposed SV estimators are generally less
variable than the BM estimators. This observation is in line with Flegal and
Jones, (2010) who showed that, for estimating means of scalar valued
functions, certain SV estimators are less variable than BM estimators by a
factor of 1.5.

The rest of the paper is organized as follows. In Section 2, we describe both
the generalized IS estimation as well as the reverse logistic regression
estimation. The proposed methods of selecting proposal densities for IS
estimators are described in Section 3. An illustrative example is given in
Section 4. Section 5 contains conclusions of the paper. Proofs of theorems and
several examples are relegated to the supplementary materials.

##  2 Generalized IS estimation of normalizing constants and expectations

Recall that $\Pi$ is a family of target densities on X. Let
$f:\textsf{X}\rightarrow\mathbb{R}$ be the function of interest. Assume that
$f$ is $\pi$-integrable for all $\pi\in\Pi$. Given samples
$\Phi_{l}\equiv\\{X_{i}^{(l)}\\}_{i=1}^{n_{l}},l=1\ldots,k$ from a small
number of proposal densities $\\{q_{l}=\varphi_{l}(x)/c_{l},l=1,\dots,k\\}$,
one wants to estimate $\theta$ (or, rather $\theta/c_{1}$) and
$E_{\pi}f:=\int_{\textsf{X}}f(x)\pi(x)\mu(dx)$ for all $\pi\in\Pi$. We
estimate $u(\pi,q_{1})\equiv\theta/c_{1}$ and $E_{\pi}f$ by
$\hat{u}(\boldsymbol{d})\equiv\hat{u}(\pi;\boldsymbol{d})$ defined in (1.3)
and
$\hat{\eta}^{[f]}(\pi;\boldsymbol{d})\equiv\hat{v}^{[f]}(\pi;\boldsymbol{d})/\hat{u}(\boldsymbol{d})$
respectively, where

$\hat{v}^{[f]}(\pi;\boldsymbol{d}):=\sum_{l=1}^{k}\frac{a_{l}}{n_{l}}\sum_{i=1}^{n_{l}}\frac{f(X_{i}^{(l)})\nu(X_{i}^{(l)})}{\sum_{j=1}^{k}a_{j}\varphi_{j}(X_{i}^{(l)})/d_{j}}.$

We also consider the more general setting when $\boldsymbol{d}$ is unknown,
which is the case if $Q=\Pi$. In such situations, we use the two-stage IS
procedure of Doss, (2010), where first, $\boldsymbol{d}$ is estimated using
Geyer,’s (1994) reverse logistic (RL) regression method (described in Section
2.1) based on Markov chain samples
$\tilde{\Phi}_{l}\equiv\\{X_{i}^{(l)}\\}_{i=1}^{N_{l}}$ with stationary
density $q_{l}$, for $l=1,\ldots,k$. Once $\boldsymbol{\hat{d}}$ is formed,
independent of stage 1, new samples $\Phi_{l},l=1\ldots,k$ are obtained to
estimate $u(\pi,q_{1})$ and $E_{\pi}f$ by $\hat{u}(\boldsymbol{\hat{d}})$ and
$\hat{\eta}^{[f]}(\pi;\boldsymbol{\hat{d}})$, respectively. Buta and Doss,
(2011) quantify benefits of the two-stage scheme as opposed to using the same
MCMC samples to estimate both $\boldsymbol{d}$ and $u(\pi,q_{1})$.

###  2.1 Reverse logistic regression estimator of $\boldsymbol{d}$

Let $N=\sum_{l=1}^{k}N_{l}$ and $a_{l}\in[0,1]$ for $l=1,\ldots,k$ such that
$\sum_{l=1}^{k}a_{l}=1$. Define

$\zeta_{l}=-\log(c_{l})+\log(a_{l}),\;l=1,\ldots,k,$ (2.1)

and

$p_{l}(x,\boldsymbol{\zeta})=\frac{\varphi_{l}(x)e^{\zeta_{l}}}{\sum_{s=1}^{k}\varphi_{s}(x)e^{\zeta_{s}}},\;l=1,\ldots,k,$
(2.2)

where $\boldsymbol{\zeta}=(\zeta_{1},\dots,\zeta_{k})$. (Note that, if
$a_{l}=N_{l}/N$, given that $x$ belongs to the pooled sample
$\big{\\{}X_{i}^{(l)},\,i=1,\ldots,N_{l},\,l=1,\ldots,k\big{\\}}$,
$p_{l}(x,\boldsymbol{\zeta})$ is the probability that $x$ comes from the
$l^{\text{th}}$ distribution.) Following Doss and Tan, (2014), consider the
log quasi-likelihood function

$\ell_{N}(\boldsymbol{\zeta})=\sum_{l=1}^{k}a_{l}\frac{N}{N_{l}}\sum_{i=1}^{N_{l}}\log\bigl{(}p_{l}(X_{i}^{(l)},\boldsymbol{\zeta})\bigr{)}.$
(2.3)

Note that adding the same constant to all $\zeta_{l}$’s leaves (2.3)
invariant. Let $\boldsymbol{\zeta}^{0}\in\mathbb{R}^{k}$ denote the true
$\boldsymbol{\zeta}$ normalized to add to zero, that is,
$\boldsymbol{\zeta}^{0}_{l}=\boldsymbol{\zeta}_{l}-\bigl{(}\sum_{j=1}^{k}\boldsymbol{\zeta}_{j}\bigr{)}/k$.
Here, $\boldsymbol{\zeta}_{l}$ denotes the $l$th element of
$\boldsymbol{\zeta}$. We estimate $\boldsymbol{\zeta}^{0}$ by
$\boldsymbol{\hat{\zeta}}$, the maximizer of $\ell_{N}$ defined in (2.3)
subject to the linear constraint $\sum_{j=1}^{k}\zeta_{j}=0$, and thus obtain
$\boldsymbol{\hat{d}}$, an estimate of $\boldsymbol{d}$ from (2.1).

##  3 Selection of proposal distributions

In this section we propose some methods for selecting the reference
distributions for efficient use of the multiple IS estimators.

###  3.1 Space filling approach

In this method, among all subsets $\boldsymbol{q}=\\{q_{1},\dots,q_{k}\\}$ of
$Q$, the one with the minimum distance from $\Pi$ is chosen. For $\pi\in\Pi$,
$q\in Q$, let $\Upsilon(\pi,q)$ be a suitably chosen metric. Define

$\psi_{p}(\boldsymbol{q},\pi)=\Big{(}\sum_{q\in\boldsymbol{q}}\Upsilon(\pi,q)^{p}\Big{)}^{1/p},$

as a measure of ‘coverage’ of $\pi$ by $\boldsymbol{q}$. The design criterion
is to choose $\boldsymbol{q}$ to minimize

$\Psi_{p,\tilde{p}}(\boldsymbol{q})=\Big{(}\sum_{\pi\in\Pi}\psi_{p}(\boldsymbol{q},\pi)^{\tilde{p}}\Big{)}^{1/\tilde{p}}$

over all subsets $\boldsymbol{q}$ with $|\boldsymbol{q}|=k$. In the limit
($p\rightarrow-\infty,\tilde{p}\rightarrow\infty$), $\Psi_{p,\tilde{p}}$ is
related to the minimax design. We use $p=-30,\tilde{p}=30$ in the examples
considered here. We use an efficient point swapping algorithm of Royle and
Nychka, (1998) to implement this space filling (SF) method. This point
swapping algorithm is described in Section H.1 of the supplement of this
paper, and is available in the R package fields (Nychka et al.,, 2017). We can
also use the simulated annealing algorithm described in the supplement for
performing the SF method.

The choice of the metric $\Upsilon(\pi,q)$ is crucial. For instance, in the
binomial robit model with degrees of freedom (df) parameter, $\xi$ example in
Section I of the supplemental materials, the family of target densities
$\Pi\equiv\\{\pi_{\xi}(x)=\nu_{\xi}(x)/\theta_{\xi}:\xi\in\Xi\\}$ is indexed
by $\xi$. In this example the IS estimator (1.3) is used with samples from
some selected target densities. Thus, $Q=\Pi$, and the choice of proposal
distributions boils down to choosing $k$ skeleton points
${\boldsymbol{\xi}}=\\{\xi_{1},\dots,\xi_{k}\\}$ inside $\Xi$. Here the
relevant geometry (with respect to $\xi$) in $\mathbb{R}$ is not Euclidean.
Indeed, df $\xi=10^{2}$ and $10^{3}$ are close, but $\xi=0.5$ and $\xi=1$ are
not. Thus the SF based on the Euclidean distance metric (SFE) may not be
appropriate unless the indexing variable is a location parameter. Another
choice is the information metric (Kass,, 1989; Rao,, 1982) which measures the
distance between two parametric distributions using asymptotic standard
deviation units of the best estimator. The Kullback-Leibler divergence
generates the information number through the information metric (Ghosh et
al.,, 2007). The information metric although seems to be appropriate for the
context, it may be difficult to implement in practice. Here, we use the
symmetric Kullback-Leibler divergence (SKLD), although it is not a metric and
denote the corresponding method by SFS. Thus,

$\Upsilon(\pi,q)=\int_{\textsf{X}}\pi(x)\log\frac{\nu(x)}{\varphi(x)}\mu(dx)-\int_{\textsf{X}}q(x)\log\frac{\varphi(x)}{\nu(x)}\mu(dx).$
(3.1)

In the special case when
$\Pi\equiv\\{\pi_{\xi}(x)=\nu_{\xi}(x)/c_{\xi}:\xi\in\Xi\\}$, that is, the
target family is indexed by some variable $\xi$, and $Q=\Pi$, the SKLD between
$\pi_{\xi_{1}}(x)$ and $\pi_{\xi_{2}}(x)$, is

$\displaystyle\Upsilon(\xi_{1},\xi_{2})$
$\displaystyle=\int_{\textsf{X}}\pi_{\xi_{1}}(x)\log\frac{\nu_{\xi_{1}}(x)}{\nu_{\xi_{2}}(x)}\mu(dx)-\int_{\textsf{X}}\pi_{\xi_{2}}(x)\log\frac{\nu_{\xi_{1}}(x)}{\nu_{\xi_{2}}(x)}\mu(dx)$
(3.2a)
$\displaystyle=\frac{\int_{\textsf{X}}\nu_{\xi_{1}}(x)\log\frac{\nu_{\xi_{1}}(x)}{\nu_{\xi_{2}}(x)}\mu(dx)}{\int_{\textsf{X}}\nu_{\xi_{1}}(x)\mu(dx)}-\frac{\int_{\textsf{X}}\nu_{\xi_{2}}(x)\log\frac{\nu_{\xi_{1}}(x)}{\nu_{\xi_{2}}(x)}\mu(dx)}{\int_{\textsf{X}}\nu_{\xi_{2}}(x)\mu(dx)}.$
(3.2b)

The SKLD (3.1) is generally not available in closed form. We use a modified
Laplace method (Evangelou et al.,, 2011) to approximate (3.2b), and we
describe the method in the supplement to this article. The second order
approximation described in the supplement is exact for two Gaussian densities.
If X is discrete, or the target distributions are far from Gaussian, a Monte
Carlo estimate of (3.2a) can be used with samples from $\pi_{\xi_{1}}$ and
$\pi_{\xi_{2}}$. Indeed, for the autologistic example, we use the Monte Carlo
estimate of (3.2a) to implement the SF approach.

The SF method does not involve the form of the IS estimator (1.3), and can be
used for selection of proposal distributions for any multiple IS estimators.
Although the method may not always be ideal as it is independent of the form
of the estimator and the particular estimand of interest, in our experience,
the SF design with a properly chosen metric consistently provides desirable
results. When $Q=\Pi$, the uniform (with respect to the chosen metric)
selection of the proposal distributions attempts to guarantee that each target
density is close to at least one reference distribution. The SF method is
attractive as generally, the IS estimator is used to simultaneously estimate
several quantities of interest resulting in different optimal design criteria.
In every step of Royle and Nychka,’s (1998) point swapping algorithm, the
number of computation of the pairwise metric between design and candidate
points is linear in $|\Xi|$, thus efficient computation of the metric is the
key facet of the algorithm. The variable $\xi$ can be multidimensional and the
range of $\xi$, in every direction, can be infinite. Thus, for computational
purpose, it may be required to narrow down the potential region of search (for
skeleton points) inside $\Xi$. This region may depend on the application.
Evangelou and Roy, (2019) considered the problem of maximizing (1.3) with
respect to $\xi$ so they used iterative, nested Laplace approximations to
identify the region where the maximizer may lie. Thus using Laplace
approximations, as in Evangelou and Roy, (2019), we can narrow $\Xi$ down to a
search set $\tilde{\Xi}$. In Section 4, we calculate the quantity of interest
(the multi-chain IS estimators $\hat{u}$ or $\hat{\eta}^{[f]}$) based on
samples from reference distributions corresponding to a preliminary SF set
chosen over a wide range of parameters values. Then using a cut-off criterion,
we decide a search set $\tilde{\Xi}$ (see Section 4 for details).

###  3.2 Minimax approach

Our second method is the minimax (MNX) design mentioned in the Introduction.
Thus the design is based on minimizing the maximum (estimated) asymptotic
variance or relative SE of $\hat{u}(\pi,\boldsymbol{\hat{d}})$ or
$\hat{\eta}^{[f]}(\pi;\boldsymbol{d})$ over $\pi\in\Pi$. Consistency and
asymptotic normality of $\boldsymbol{\hat{d}}$,
$\hat{u}(\pi;\boldsymbol{\hat{d}})$ and
$\hat{\eta}^{[f]}(\pi;\boldsymbol{\hat{d}})$ are described in Theorems 1, 2
and 3, respectively of Roy et al., (2018). In the following, we discuss
estimation of the asymptotic variances of these estimators.

Spectral variance estimation in reverse logistic regression and generalized IS
methods: First, we provide an SV estimator of the asymptotic covariance matrix
of $\boldsymbol{\hat{d}}$, as it is needed for the asymptotic variances of
$\hat{u}(\pi;\boldsymbol{\hat{d}})$ and
$\hat{\eta}^{[f]}(\pi;\boldsymbol{\hat{d}})$. Also, SV estimator of
Var$(\boldsymbol{\hat{d}})$ is important in its own right, and is used in
section 3.3 in our third approach to selection of proposal distributions.

As in Roy et al., (2018), we assume that the Markov chains $\Phi_{l}$,
$\tilde{\Phi}_{l}$ are polynomially ergodic for $l=1,\dots,k$. (The definition
of polynomial ergodicity of Markov chains can be found in Roy et al., (2018).)
They showed that if the Markov chain $\tilde{\Phi}_{l}$ is polynomially
ergodic of order $t>1$ for $l=1,\dots,k$, then $\boldsymbol{\hat{\zeta}}$ and
$\boldsymbol{\hat{d}}$ defined in section 2.1 are consistent and
asymptotically normal as $N_{1},\ldots,N_{k}\rightarrow\infty$, that is, there
exist matrices $B,\Omega\in\mathbb{R}^{k,k}$ and $D\in\mathbb{R}^{k,k-1}$ such
that

$\sqrt{N}(\boldsymbol{\hat{\zeta}}-\boldsymbol{\zeta})\buildrel
d\over{\rightarrow}{\cal
N}(0,U)\quad\mbox{and}\quad\sqrt{N}(\boldsymbol{\hat{d}}-\boldsymbol{d})\buildrel
d\over{\rightarrow}{\cal N}(0,V),$

where $U=B^{\dagger}\Omega B^{\dagger}$ and $V=D^{\top}UD$. Here, for a square
matrix $C$, $C^{\dagger}$ denotes its Moore-Penrose inverse. The matrices $B$,
$\Omega$ and $D$ are as defined in (2.7), (2.8), and (2.5) respectively in Roy
et al., (2018). Theorem 1 below provides consistent SV estimators of the
asymptotic variances of $\boldsymbol{\hat{\zeta}}$ and $\boldsymbol{\hat{d}}$.

We now introduce some notations. Assume $N_{l}\rightarrow\infty$ such that
$\lim N_{l}/N\in(0,1)$ for $l=1,\dots,k$. Note that the function
$g\colon{\mathbb{R}}^{k}\rightarrow{\mathbb{R}}^{k-1}$ that maps
$\boldsymbol{\zeta}^{0}$ into $\boldsymbol{d}$ is given by

$g(\boldsymbol{\zeta})=(e^{\zeta_{1}-\zeta_{2}}a_{2}/a_{1},e^{\zeta_{1}-\zeta_{3}}a_{3}/a_{1},\dots,e^{\zeta_{1}-\zeta_{k}}a_{k}/a_{1})^{\top}.$
(3.3)

Since $\boldsymbol{d}=g(\boldsymbol{\zeta}^{0})$, by definition
$\boldsymbol{\hat{d}}=g(\boldsymbol{\hat{\zeta}})$, and its gradient at
$\boldsymbol{\hat{\zeta}}$ (in terms of $\boldsymbol{\hat{d}}$) is

$\widehat{D}=\begin{pmatrix}\hat{d}_{2}&\hat{d}_{3}&\ldots&\hat{d}_{k}\\\
-\hat{d}_{2}&0&\ldots&0\\\ 0&-\hat{d}_{3}&\ldots&0\\\
\vdots&\vdots&\ddots&\vdots\\\ 0&0&\ldots&-\hat{d}_{k}\end{pmatrix}.$ (3.4)

As in Roy et al., (2018), the $k\times k$ matrix $\widehat{B}$ is defined by

$\begin{split}\widehat{B}_{rr}&=\sum_{l=1}^{k}a_{l}\biggl{(}\frac{1}{N_{l}}\sum_{i=1}^{N_{l}}p_{r}(X_{i}^{(l)},\boldsymbol{\hat{\zeta}})\bigl{[}1-p_{r}(X_{i}^{(l)},\boldsymbol{\hat{\zeta}})\bigr{]}\biggr{)}\text{
and}\\\
\widehat{B}_{rs}&=-\sum_{l=1}^{k}a_{l}\biggl{(}\frac{1}{N_{l}}\sum_{i=1}^{N_{l}}p_{r}(X_{i}^{(l)},\boldsymbol{\hat{\zeta}})p_{s}(X_{i}^{(l)},\boldsymbol{\hat{\zeta}})\biggr{)}\text{
for }r\neq s,\end{split}$ (3.5)

that is, $\widehat{B}$ denotes the matrix of second derivatives of
$-\ell_{N}(\boldsymbol{\zeta})/N$ evaluated at $\boldsymbol{\hat{\zeta}}$,
where $\ell_{N}(\boldsymbol{\zeta})$ is defined in (2.3). Set
$Z^{(l)}_{i}=\left(p_{1}(X_{i}^{(l)},\boldsymbol{\hat{\zeta}}),\ldots,p_{k}(X_{i}^{(l)},\boldsymbol{\hat{\zeta}})\right)^{\top}$
for $i=1,\dots,N_{l}$ and $\bar{Z}^{(l)}=\sum_{i=1}^{N_{l}}Z^{(l)}_{i}/N_{l}$.
Define the lag $j$ sample autocovariance as

$\gamma_{N}^{(l)}(j)=\frac{1}{N_{l}}\sum_{i\in
S_{j,N}}\left[Z^{(l)}_{i}-\bar{Z}^{(l)}\right]\left[Z^{(l)}_{i+j}-\bar{Z}^{(l)}\right]^{\top}\;\mbox{
for}\;\;l=1,\ldots,k,$ (3.6)

where $S_{j,N}=\\{1,\dots,N-j\\}$ for $j\geq 0$ and
$S_{j,N}=\\{(1-j),\dots,N\\}$ for $j<0$. Let

$\widehat{\Sigma}^{(l)}=\sum_{j=-(b_{N_{l}}-1)}^{b_{N_{l}}-1}w_{N_{l}}(j)\gamma_{N}^{(l)}(j),$
(3.7)

where $w_{N_{l}}(\cdot)$ is the lag window, $b_{N_{l}}$’s are the truncation
points for $l=1,\ldots,k$. Finally, define

$\widehat{\Omega}=\sum_{l=1}^{k}\frac{N}{N_{l}}a^{2}_{l}\widehat{\Sigma}^{(l)}.$
(3.8)

######  Theorem 1.

Assume that the Markov chains $\tilde{\Phi}_{1},\ldots,\tilde{\Phi}_{k}$ are
polynomially ergodic of order $t>1$, and for all $l=1,\ldots,k$, $w_{N_{l}}$
and $b_{N_{l}}$ satisfy conditions 1-4 in Vats et al., (2018, Theorem 2). Let
$\widehat{D}$, $\widehat{B}$ and $\widehat{\Omega}$ be the matrices defined by
(3.4), (3.5) and (3.8), respectively. Then,
$\widehat{U}:=\widehat{B}^{\dagger}\widehat{\Omega}\widehat{B}^{\dagger}$ and
$\widehat{V}:=\widehat{D}^{\top}\widehat{U}\widehat{D}$ are strongly
consistent estimators of $U$ and $V$ respectively.

Next, we consider estimation of the asymptotic variances of
$\hat{u}(\pi;\boldsymbol{\hat{d}})$ and
$\hat{\eta}^{[f]}(\pi;\boldsymbol{\hat{d}})$. Roy et al., (2018) showed that,
under certain conditions, there exist $\sigma^{2}_{u},\sigma^{2}_{\eta}>0$
such that, as $n_{1},\ldots,n_{k}\rightarrow\infty$,

$\sqrt{n}(\hat{u}(\pi;\boldsymbol{\hat{d}})-u(\pi,q_{1}))\buildrel
d\over{\rightarrow}N(0,\sigma^{2}_{u})\quad\mbox{and}\quad\sqrt{n}(\hat{\eta}^{[f]}(\pi;\boldsymbol{\hat{d}})-E_{\pi}f)\buildrel
d\over{\rightarrow}N(0,\sigma^{2}_{\eta}),$ (3.9)

where $n=\sum_{l=1}^{k}n_{l}$. In Theorem 2 we provide consistent SV
estimators of $\sigma^{2}_{u}$ and $\sigma^{2}_{\eta}$. We first introduce
some notations. Let

$u^{\pi}(x;\boldsymbol{d}):=\frac{\nu(x)}{\sum_{s=1}^{k}a_{s}\varphi_{s}(x)/d_{s}}\quad\text{and}\quad
v^{[f],\pi}(x;\boldsymbol{d}):=f(x)u^{\pi}(x;\boldsymbol{d}).$ (3.10)

Define the vectors $c(\pi;\boldsymbol{d})$ and $e(\pi;\boldsymbol{d})$ of
length $k-1$ with $(j-1)$th coordinate as

$\displaystyle[c(\pi;\boldsymbol{d})]_{j-1}=\frac{u(\pi,q_{1})}{d_{j}^{2}}\int_{\textsf{X}}\frac{a_{j}\varphi_{j}(x)}{\sum_{s=1}^{k}a_{s}\varphi_{s}(x)/d_{s}}\pi(x)\mu(dx)$
(3.11)
$\displaystyle[e(\pi;\boldsymbol{d})]_{j-1}=\frac{a_{j}}{d_{j}^{2}}\int_{\textsf{X}}\frac{[f(x)-E_{\pi}f]\varphi_{j}(x)}{\sum_{s=1}^{k}a_{s}\varphi_{s}(x)/d_{s}}\pi(x)\mu(dx),$
(3.12)

for $j=2,\dots,k$, and their estimators $\hat{c}(\pi;\boldsymbol{d})$ and
$\hat{e}(\pi;\boldsymbol{d})$ as

$\displaystyle[\hat{c}(\pi;\boldsymbol{d})]_{j-1}=\sum_{l=1}^{k}\frac{1}{n_{l}}\sum_{i=1}^{n_{l}}\frac{a_{j}a_{l}\nu(X_{i}^{(l)})\varphi_{j}(X_{i}^{(l)})}{(\sum_{s=1}^{k}a_{s}\varphi_{s}(X_{i}^{(l)})/d_{s})^{2}d_{j}^{2}},$
(3.13)
$\displaystyle[\hat{e}(\pi;\boldsymbol{d})]_{j-1}=\frac{\sum_{l=1}^{k}\frac{a_{l}}{n_{l}}\sum_{i=1}^{n_{l}}\frac{a_{j}f(X_{i}^{(l)})\nu(X_{i}^{(l)})\varphi_{j}(X_{i}^{(l)})}{d_{j}^{2}(\sum_{s=1}^{k}a_{s}\varphi_{s}(X_{i}^{(l)})/d_{s})^{2}}}{\hat{u}(\pi;\boldsymbol{d})}-\frac{[\hat{c}(\pi;\boldsymbol{d})]_{j-1}\hat{\eta}^{[f]}(\pi;\boldsymbol{d})}{\hat{u}(\pi;\boldsymbol{d})}.$
(3.14)

Suppose $b_{n_{l}}$’s are the truncation points, $w_{n_{l}}(j)$’s are lag
window, $u_{i}\equiv u_{i}(\boldsymbol{d})\equiv
u^{\pi}(X_{i}^{(l)};\boldsymbol{d})$, $v^{[f]}_{i}\equiv
v^{[f]}_{i}(\boldsymbol{d})\equiv v^{[f],\pi}(X_{i}^{(l)};\boldsymbol{d})$,
and $\bar{u}\equiv\bar{u}(\boldsymbol{d})$,
$\bar{v}^{[f]}\equiv\bar{v}^{[f]}(\boldsymbol{d})$ are the averages of
$\\{u^{\pi}(X_{1}^{(l)};\boldsymbol{d}),\cdots,u^{\pi}(X_{n_{l}}^{(l)};\boldsymbol{d})\\}$
and
$\\{v^{[f],\pi}(X_{1}^{(l)};\boldsymbol{d}),\cdots,v^{[f],\pi}(X_{n_{l}}^{(l)};\boldsymbol{d})\\}$,
respectively. (Note that, abusing notations, the dependence on $l$ is ignored
in $u_{i},v^{[f]}_{i},\bar{u}$ and $\bar{v}^{[f]}$.) Let

$\hat{\tau}^{2}_{l}(\pi;\boldsymbol{d})=\frac{1}{n_{l}}\sum_{j=-(b_{n_{l}}-1)}^{b_{n_{l}}-1}w_{n_{l}}(j)\sum_{i\in
S_{j,n}}\left[u_{i}-\bar{u}\right]\left[u_{i+j}-\bar{u}\right],\;\mbox{and}$
(3.15)
$\widehat{\Gamma}_{l}(\pi;\boldsymbol{d})=\frac{1}{n_{l}}\sum_{j=-(b_{n_{l}}-1)}^{b_{n_{l}}-1}w_{n_{l}}(j)\sum_{i\in
S_{j,n}}\Bigg{[}\left(\begin{array}[]{c}v_{i}^{[f]}\\\ u_{i}\\\
\end{array}\right)-\left(\begin{array}[]{c}\bar{v}^{[f]}\\\ \bar{u}\\\
\end{array}\right)\Bigg{]}\Bigg{[}\left(\begin{array}[]{c}v_{i+j}^{[f]}\\\
u_{i+j}\\\ \end{array}\right)-\left(\begin{array}[]{c}\bar{v}^{[f]}\\\
\bar{u}\\\ \end{array}\right)\Bigg{]}^{\top}.$

Finally, let
$\hat{\tau}^{2}(\pi;\boldsymbol{d})=\sum_{l=1}^{k}(a_{l}^{2}n/n_{l})\hat{\tau}^{2}_{l}(\pi;\boldsymbol{d})$,
$\widehat{\Gamma}(\pi;\boldsymbol{d})=\sum_{l=1}^{k}(a_{l}^{2}n/n_{l})\widehat{\Gamma}_{l}(\pi;\boldsymbol{d})$,
and

$\hat{\rho}(\pi;\boldsymbol{\hat{d}})=\nabla
h(\hat{v}^{[f]}(\pi;\boldsymbol{\hat{d}}),\hat{u}(\boldsymbol{\hat{d}}))^{\top}\widehat{\Gamma}(\pi;\boldsymbol{\hat{d}})\nabla
h(\hat{v}^{[f]}(\pi;\boldsymbol{\hat{d}}),\hat{u}(\boldsymbol{\hat{d}})),$

where $\nabla h(x,y)=(1/y,-x/y^{2})^{\top}$.

######  Theorem 2.

Suppose that for $\tilde{\Phi}_{l},l=1,\dots,k$, conditions of Theorem 1 hold
and $\widehat{V}$ is the consistent SV estimator of $V$. Suppose there exists
$r\in[0,\infty)$ such that $n/N\rightarrow r$ where $N=\sum_{l=1}^{k}N_{l}$
and $n=\sum_{l=1}^{k}n_{l}$ are the total sample sizes for stage 1 and stage
2, respectively. In addition, let $n_{l}/n\rightarrow s_{l}\in(0,1)$ for
$l=1,\cdots,k$. Assume that the Markov chains $\Phi_{1},\ldots,\Phi_{k}$ are
polynomially ergodic of order $t\geq(1+\epsilon)(1+2/\delta)$ for some
$\epsilon,\delta>0$ such that
$E_{q_{l}}|u^{\pi}(X;\boldsymbol{d})|^{4+\delta}<\infty$, and for each
$l=1,\cdots,k$, $w_{n_{l}}$ and $b_{n_{l}}$ satisfy conditions 1-4 in Vats et
al., (2018, Theorem 2).

  1. (a)

Then
$\hat{\sigma}^{2}_{u}=(n/N)\hat{c}(\pi;\boldsymbol{\hat{d}})^{\top}\widehat{V}\hat{c}(\pi;\boldsymbol{\hat{d}})+\hat{\tau}^{2}(\pi;\boldsymbol{\hat{d}})$
is a strongly consistent estimator of $\sigma^{2}_{u}$.

  2. (b)

In addition, suppose that
$E_{q_{l}}|v^{[f],\pi}(X;\boldsymbol{d})|^{4+\delta}<\infty$. Then
$\hat{\sigma}^{2}_{\eta}=(n/N)\hat{e}(\pi;\boldsymbol{\hat{d}})^{\top}\widehat{V}\hat{e}(\pi;\boldsymbol{\hat{d}})+\hat{\rho}(\pi;\boldsymbol{\hat{d}})$
is a strongly consistent estimator of the asymptotic variance
$\sigma^{2}_{\eta}$.

Since samples are obtained by running the Markov chains with the stationary
densities in $\boldsymbol{q}$ we denote the corresponding RL estimator of
$\boldsymbol{d}\equiv\boldsymbol{d}_{\boldsymbol{q}}$ by
$\boldsymbol{\hat{d}}_{\boldsymbol{q}}$ and its asymptotic variance as
${V}_{\boldsymbol{q}}$. Similarly, in this case, we denote the SV estimators
of the asymptotic variances (3.9) of
$\hat{u}(\pi;\boldsymbol{\hat{d}}_{\boldsymbol{q}})$ and
$\hat{\eta}^{[f]}(\pi;\boldsymbol{\hat{d}}_{\boldsymbol{q}})$ as
$\hat{\sigma}^{2}_{u}(\pi;\boldsymbol{q})$ and
$\hat{\sigma}^{2}_{\eta}(\pi;\boldsymbol{q})$, respectively. Thus, if the goal
is to estimate the normalizing constants in $\Pi$ via (1.3), then the MNX set
$\boldsymbol{q}$ minimizes
$\phi(\boldsymbol{q})\equiv\mbox{max}_{\pi\in\Pi}\hat{\sigma}^{2}_{u}(\pi;\boldsymbol{q})$.
In general, solving this design problem is difficult. Here, we use simulated
annealing as described in Section H.2 of the supplementary materials, to
estimate the MNX set $\boldsymbol{q}$. If the goal is to estimate $E_{\pi}f$
for $\pi\in\Pi$, then simulated annealing is used to find the MNX set
$\boldsymbol{q}$ that minimizes
$\mbox{max}_{\pi\in\Pi}\hat{\sigma}^{2}_{\eta}(\pi;\boldsymbol{q})$.

When $Q=\Pi$, a less computationally demanding approach is a method in which
densities are chosen sequentially from $\Pi$ where
$\hat{\sigma}^{2}_{u}(\pi;\boldsymbol{q})$ is the largest. Specifically,
starting with an initial density $\boldsymbol{q}_{0}=\\{\tilde{q}\\}$, suppose
that we have completed the $i$th step with the set $\boldsymbol{q}_{i}$ chosen
along with (Markov chain) samples from each density in $\boldsymbol{q}_{i}$.
If $\boldsymbol{d}$ is unknown, part of this sample (stage 1) is used for
calculating the RL estimator $\boldsymbol{\hat{d}}$, and the remaining sample
is used to compute $\hat{\sigma}^{2}_{u}(\pi;\boldsymbol{q}_{i})$ for the
remaining densities $\pi\in\Pi\setminus\boldsymbol{q}_{i}$. Then
$\boldsymbol{q}_{i+1}=\boldsymbol{q}_{i}\cup\\{\pi_{j}\\}$ where
$\pi_{j}=\operatorname*{argmax}_{\pi\in\Pi\setminus\boldsymbol{q}_{i}}\hat{\sigma}^{2}_{u}(\pi;\boldsymbol{q}_{i})$,
and the existing (Markov chain) sample is augmented with samples from
$\pi_{j}$. Thus, at each step, the density corresponding to the largest
(estimated) asymptotic variance is chosen. The process is repeated until $k$
densities have been selected. This sequential method is a principled approach
formalizing the ad-hoc procedure of Buta and Doss, (2011). The initial
$\tilde{q}$ can be the density where the generalized IS estimator (1.3) or any
other interesting quantity based on samples from a preliminary SF set is
maximized (see Section 4 for details). Note that, both the MNX approach as
well as the sequential (SEQ) approach can be used to select proposal
distributions for any multi-density IS estimator for which we have valid
estimates of its asymptotic variance.

###  3.3 Maximum entropy approach

The third method uses maximum entropy sampling (Shewry and Wynn,, 1987) for
selecting $\boldsymbol{q}$. This method is applicable when $\boldsymbol{d}$ is
unknown and is developed in the context of Doss,’s (2010) two-stage IS
estimation scheme described before. We use the notation $\mbox{Ent}(\cdot)$ to
denote the Boltzmann-Shannon entropy of the expression inside the brackets.
Under a suitable decomposition of the joint distribution of the estimator, let
that be $\hat{u}$ or $\hat{\eta}^{[f]}$, over the whole set $\Pi$, and
$\boldsymbol{\hat{d}}_{\boldsymbol{q}}$, the maximum entropy (ENT) approach
chooses $\boldsymbol{q}$ that maximizes
$\mbox{Ent}(\boldsymbol{\hat{d}}_{\boldsymbol{q}})$. Note that since
$\boldsymbol{\hat{d}}_{\boldsymbol{q}}$ depends on the reference density
$q_{1}$, it is assumed that $q_{1}$ remains fixed, which can be the density
$\tilde{q}$ discussed in Section 3.2. In the following we assume that the
objective is to estimate ratios of normalizing constants. In the supplementary
materials we derive similar results under the objective of estimating means
$E_{\pi}f$.

Let $\hat{\bf u}(\boldsymbol{\pi};\boldsymbol{\hat{d}}_{\boldsymbol{q}})$ be
the vector of length $|\Pi|$ consisting of
$\hat{u}(\pi;\boldsymbol{\hat{d}}_{\boldsymbol{q}})$’s, $\pi\in\Pi$ in a (any)
fixed order. Indeed, we refer to this fixed ordering whenever we write $\Pi$
in this section. Similarly define the vector of true (ratios of) normalizing
constants ${\bf u}(\boldsymbol{\pi},q_{1})$. Let
$C(\boldsymbol{\pi};\boldsymbol{d}_{\boldsymbol{q}})$ be the
$|\Pi|\times(k-1)$ matrix with rows $c(\pi;\boldsymbol{d}_{\boldsymbol{q}})$
(defined in (3.11)), $\pi\in\Pi$. Similarly, define
$\widehat{C}(\boldsymbol{\pi};\boldsymbol{d}_{\boldsymbol{q}})$ with rows
$\hat{c}(\pi;\boldsymbol{d}_{\boldsymbol{q}})$ (defined in (3.13)),
$\pi\in\Pi$. Let ${\bf
u}^{\boldsymbol{\pi}}(x;\boldsymbol{d}_{\boldsymbol{q}})$ be the $|\Pi|$
dimensional vector consisting of
$u^{\pi}(x;\boldsymbol{d}_{\boldsymbol{q}})$’s defined in (3.10). Let
$T_{l}(\boldsymbol{d}_{\boldsymbol{q}})$ be the $|\Pi|\times|\Pi|$ matrix with
elements

$\displaystyle\tau^{2}_{l}(\pi,\pi^{\prime};\boldsymbol{d}_{\boldsymbol{q}})$
$\displaystyle=\operatorname{Cov}_{q_{l}}(u^{\pi}(X_{1}^{(l)};\boldsymbol{d}_{\boldsymbol{q}}),u^{\pi^{\prime}}(X_{1}^{(l)};\boldsymbol{d}_{\boldsymbol{q}}))$
(3.16)
$\displaystyle+\sum_{g=1}^{\infty}\operatorname{Cov}_{q_{l}}(u^{\pi}(X_{1}^{(l)};\boldsymbol{d}_{\boldsymbol{q}}),u^{\pi^{\prime}}(X_{1+g}^{(l)};\boldsymbol{d}_{\boldsymbol{q}}))+\sum_{g=1}^{\infty}\operatorname{Cov}_{q_{l}}(u^{\pi}(X_{1+g}^{(l)};\boldsymbol{d}),u^{\pi^{\prime}}(X_{1}^{(l)};\boldsymbol{d}_{\boldsymbol{q}})).$

Finally, let

$\widehat{T}_{l}(\boldsymbol{d}_{\boldsymbol{q}})=\frac{1}{n_{l}}\sum_{j=-(b_{n_{l}}-1)}^{b_{n_{l}}-1}w_{n_{l}}(j)\sum_{i\in
S_{j,n}}\left[{\bf
u}^{\boldsymbol{\pi}}(X_{i}^{(l)};\boldsymbol{d}_{\boldsymbol{q}})-\bar{{\bf
u}}(\boldsymbol{d}_{\boldsymbol{q}})\right]\left[{\bf
u}^{\boldsymbol{\pi}}(X_{i+j}^{(l)};\boldsymbol{d}_{\boldsymbol{q}})-\bar{{\bf
u}}(\boldsymbol{d}_{\boldsymbol{q}})\right]^{\top},$ (3.17)

where $b_{n_{l}}$’s are the truncation points, $w_{n_{l}}(j)$’s are the lag
windows, and $\bar{\bf
u}(\boldsymbol{d}_{\boldsymbol{q}})=\sum_{i=1}^{n_{l}}{\bf
u}^{\boldsymbol{\pi}}(X_{i}^{(l)};\boldsymbol{d}_{\boldsymbol{q}})/n_{l}$.

######  Theorem 3.

Suppose there exists $r\in[0,\infty)$ such that $n/N\rightarrow r$ where
$N\equiv\sum_{l=1}^{k}N_{l}$ and $n=\sum_{l=1}^{k}n_{l}$ are the total sample
sizes for stages 1 and 2, respectively. In addition, let $n_{l}/n\rightarrow
s_{l}\in(0,1)$ for $l=1,\cdots,k$.

  1. (a)

Assume that the stage 1 Markov chains $\tilde{\Phi}_{l},l=1,\dots,k$ are
polynomially ergodic of order $t>1$. Further, assume that the stage 2 Markov
chains $\Phi_{l},l=1,\dots,k$ are polynomially ergodic of order $t$, and for
some $\delta>0$
$E_{q_{l}}|u^{\pi}(X;\boldsymbol{d}_{\boldsymbol{q}})|^{2+\delta}<\infty$ for
each $\pi\in\Pi$ and $l=1,\cdots,k$ where $t>1+2/\delta$. Then as
$n_{1},\ldots,n_{k}\rightarrow\infty$,

$\sqrt{n}\left(\begin{array}[]{c}\boldsymbol{\hat{d}}_{\boldsymbol{q}}-\boldsymbol{d}_{\boldsymbol{q}}\\\
\hat{\bf u}(\boldsymbol{\pi};\boldsymbol{\hat{d}}_{\boldsymbol{q}})-{\bf
u}(\boldsymbol{\pi},q_{1})\\\ \end{array}\right)\buildrel
d\over{\rightarrow}N\Bigg{(}0,\left(\begin{array}[]{cc}rV_{\boldsymbol{q}}&\Sigma_{12}\\\
\Sigma_{21}&\Sigma_{22}\\\ \end{array}\right)\Bigg{)},$ (3.18)

where
$\Sigma_{21}=rC(\boldsymbol{\pi};\boldsymbol{d}_{\boldsymbol{q}})V_{\boldsymbol{q}}$,
$\Sigma_{12}=\Sigma_{21}^{\top}$, and
$\Sigma_{22}=rC(\boldsymbol{\pi};\boldsymbol{d}_{\boldsymbol{q}})V_{\boldsymbol{q}}C(\boldsymbol{\pi};\boldsymbol{d}_{\boldsymbol{q}})^{\top}+\sum_{l=1}^{k}(a_{l}^{2}/s_{l})T_{l}(\boldsymbol{d}_{\boldsymbol{q}})$.

  2. (b)

Suppose that the conditions of Theorem 1 hold for the stage 1 Markov chains.
Let $\widehat{V}_{\boldsymbol{q}}$ be the consistent estimator of
$V_{\boldsymbol{q}}$ given in Theorem 1. Assume that the Markov chains
$\Phi_{l},l=1,\dots,k$ are polynomially ergodic of order
$t\geq(1+\epsilon)(1+2/\delta)$ for some $\epsilon,\delta>0$ such that
$E_{q_{l}}\|{\bf
u}^{\boldsymbol{\pi}}(X;\boldsymbol{d}_{\boldsymbol{q}})\|^{4+\delta}<\infty$,
($\|\cdot\|$ denotes the Euclidean norm) for all $l=1,\ldots,k$, and
$w_{n_{l}}$ and $b_{n_{l}}$ satisfy conditions 1-4 in Vats et al., (2018,
Theorem 2). Then
$(n/N)\widehat{C}(\boldsymbol{\pi};\boldsymbol{\hat{d}}_{\boldsymbol{q}})\widehat{V}_{\boldsymbol{q}}\widehat{C}(\boldsymbol{\pi};\boldsymbol{\hat{d}}_{\boldsymbol{q}})^{\top}+\sum_{l=1}^{k}(a_{l}^{2}/s_{l})\widehat{T}_{l}(\boldsymbol{\hat{d}}_{\boldsymbol{q}})$
is a strongly consistent estimator of $\Sigma_{22}$ and
$(n/N)\widehat{C}(\boldsymbol{\pi};\boldsymbol{\hat{d}}_{\boldsymbol{q}})\widehat{V}_{\boldsymbol{q}}$
is a consistent estimator of $\Sigma_{21}$.

Let $Y\equiv(Y^{T}_{\boldsymbol{q}},Y^{T}_{\Pi})^{T}$ be a random vector
having the normal distribution in (3.18). The Boltzmann-Shannon entropy of $Y$
is $\mbox{Ent}(Y)=\mbox{constant}+\frac{1}{2}\log\mbox{det}(\Sigma)$, where
$\Sigma$ is the covariance matrix of $Y$. Note that

$\displaystyle\log\mbox{det}(\Sigma)=\log\mbox{det}(rV_{\boldsymbol{q}})+\log\mbox{det}(\Sigma_{22}-rC(\boldsymbol{\pi};\boldsymbol{d}_{\boldsymbol{q}})V_{\boldsymbol{q}}C(\boldsymbol{\pi};\boldsymbol{d}_{\boldsymbol{q}})^{\top}),$

where the second matrix on the right side is the covariance matrix of the
conditional distribution of $Y_{\Pi}|Y_{\boldsymbol{q}}$. Since Theorem 3 (b)
provides a consistent estimator of this conditional covariance matrix, we can
minimize the determinant of this estimator matrix to choose $\boldsymbol{q}$.

As mentioned in Shewry and Wynn, (1987), great computational benefit can be
achieved by converting this conditional problem to an unconditional problem.
In particular, as noted in Shewry and Wynn, (1987), minimization of the second
term is equivalent to maximization of $\log\mbox{det}(V_{\boldsymbol{q}})$
Caselton and Zidek, (1984). In practice, we would replace $V_{\boldsymbol{q}}$
by its estimator given in Theorem 1, i.e. $\widehat{V}_{\boldsymbol{q}}$,
using Markov chain samples from densities in $\boldsymbol{q}$.

It is known that the entropy sampling problem is NP-hard (Ko et al.,, 1995).
Here, we use simulated annealing with
$\phi(\boldsymbol{q})=-\log\det(\widehat{V}_{\boldsymbol{q}})$ to choose the
set $\boldsymbol{q}$. Unlike the SF, MNX and SEQ methods, the ENT approach is
applicable only in the context of Doss,’s (2010) two-stage IS estimation
scheme. In contrast, if the generalized IS estimator (1.3) is used, since ENT
avoids the second stage IS estimation, it needs fewer samples than the MNX and
SEQ methods which require enough samples to be used for both stages. However,
one advantage of the MNX and SEQ methods is that at the end of the procedure,
we already have available samples from densities in $\boldsymbol{q}$ which can
be used in the two-stage IS estimation scheme. Another benefit of ENT is that
it avoids computing the target density $\nu$. The ENT set of proposal
densities works no matter whether one is interested in means or normalizing
constants, whereas the MNX approach uses different criteria for these two
problems. Finally, among the proposed methods, the SF geometric approach is
the most straightforward to implement as it uses only distance between the
densities, and estimation of any asymptotic variance is not required for the
SF method.

We compare the different methods in the context of the autologistic model
described in the Introduction. The family of importance densities $Q=\Pi$ in
this case, therefore, choosing the importance densities amounts to choosing
the parameters $(\gamma,\kappa)$. We compare six methods, namely NIS (naive
IS), SFE, SFS, SEQ, MNX, and ENT. For NIS, samples are taken from
$\pi(\boldsymbol{x}|\gamma=0,\kappa=0.5)$ and for the multiple IS methods the
samples are taken from $k=5$ different densities, one of which must be the one
used for NIS. We consider two cases, one where $\kappa=0.5$ and only $\gamma$
can vary, and one where both $(\gamma,\kappa)$ vary. Further details on the
computational methods used can be found in Section J of the supplementary
material. The relative SE computed based on the SV estimator corresponding to
the optimal skeleton set of each method is obtained. This estimator is based
on a total of 100,000 samples using the full Gibbs sampler. From Figure 2, we
see that NIS estimators, as expected, result in large relative SE when the
target pmf is farther away from the proposal pmf corresponding to $\gamma=0$
and $\kappa=0.5$. In particular, the relative SE for NIS is about 10 times
larger than that for MNX at $\gamma=-4$ ($\kappa=0.5$ fixed case) meaning that
NIS requires 100 times more samples to achieve the same accuracy. Similarly,
the relative SE for NIS is about five times larger than that for MNX at
$\gamma=-4$, $\kappa=0.9$ ($\kappa$ varying case). MNX, SFE, and SFS turn out
to be the best performers in both cases ($\kappa$ known and unknown) overall,
whereas SEQ and ENT have the worst performance in selecting proposal
distributions for the generalized IS estimator. Indeed, as can be seen from
Figure 2, the relative SEs for SEQ and ENT are about two times larger than
that for MNX uniformly over all values of $\gamma$ when $\kappa=0.9$.

Figure 2: Relative SE of the estimator (1.3) for the naive and multiple
proposal IS estimators corresponding to different proposal distributions for
the autologistic model, plotted against $\gamma$. The left plot is with known
$\kappa=0.5$. The middle and right plots are with both $\gamma,\kappa$ unknown
and the plots show the relative SE at $\kappa=0.4$ and $\kappa=0.9$. The
$\gamma$ values corresponding to the proposal distributions are denoted by
$+$.

##  4 Analysis of radionuclide concentrations using spatial GLMM

The dataset consists of spatial measurements of $\gamma$-ray counts $y_{i}$
observed during $\ell_{i}$ seconds at the $i$th coordinate on the Rongelap
island, $i=1,\ldots,157$. These data were analyzed by Diggle et al., (1998)
and Evangelou and Roy, (2019) among others using a spatial generalized linear
mixed model (SGLMM). We consider a Poisson SGLMM using a parametric link
function for the $\gamma$-ray counts, that is, we assume
$y_{i}|\mu_{i}\buildrel\text{ind}\over{\sim}\mathrm{Po}(\ell_{i}\mu_{i})$ with
$g_{\lambda}(\mu_{i})=z_{i}$ for $i=1,\dots,157,$ where $g_{\lambda}(\cdot)$
is a modified Box-Cox link given in Evangelou and Roy, (2019) with parameter
$\lambda$, and $z_{i}$’s are the latent variables. Let $\boldsymbol{y}$ and
$\boldsymbol{\mu}$ denote the vectors of $y_{i}$’s and $\mu_{i}$’s,
respectively. Then $\boldsymbol{z}=(z_{1},\dots,z_{157})$ is modeled by a
multivariate Gaussian distribution corresponding to a Gaussian random field
(GRF) $\mathcal{Z}$ at the sampled locations. In particular, we assume
$\mathcal{Z}|\beta,\sigma^{2}\sim\mathrm{GRF}(\beta,\sigma^{2},\phi,\omega,\kappa)$,
the GRF with constant mean $\beta$, Matérn correlation, variance $\sigma^{2}$,
range $\phi$, relative nugget $\omega$ and smoothness $\kappa$. The partial
sill parameter $\sigma^{2}$ is assigned a scaled-inverse-chi-square prior
($\mathrm{ScInv}\mathcal{X}^{2}(1,1)$), and conditioned on $\sigma^{2}$, the
mean parameter $\beta\sim N(0,100\sigma^{2})$. Let
$\xi=(\lambda,\phi,\omega,\kappa)$. We consider estimating the marginal
likelihood for $\xi$ (relative to an arbitrary reference point
$\xi_{1}=\tilde{\xi}$ to be defined later) by (1.3). Note that, the empirical
Bayes estimate of $\xi$ is the point where the marginal likelihood function is
maximized (see e.g. Roy et al.,, 2016). Since $\beta,\sigma^{2}$ can be
analytically integrated out, one can work with the posterior density of
$\boldsymbol{z}$, $\pi_{\xi}(\boldsymbol{z}|\boldsymbol{y})$. Here, we
consider generalized IS estimator (1.3) based on samples from the
(transformed) density $\pi_{\xi}(\boldsymbol{\mu}|\boldsymbol{y})$ (see
Evangelou and Roy,, 2019, for the reasons for considering the transformed
samples). Thus, here
$Q=\Pi=\\{\pi_{\xi}(\boldsymbol{\mu}|\boldsymbol{y}),\xi\in\Xi\\}$ for some
$\Xi$ defined later.

In order to narrow down the potential region of search, we initially choose a
wide range of values for each component of $\xi$, and form a large grid,
denoted by $\Xi$, by combining discrete points within these ranges. This gives
us the set consisting of the following $9^{4}$ points:

$\Xi=\\{0,0.5,\ldots,4\\}\times\\{100,425,\ldots,2700\\}\times\\{0,0.75,\ldots,6\\}\times\\{0.1,0.35,\ldots,2.1\\}.$

The SFE method, after each range is scaled in $[0,1]$, was applied to choose
$k=5$ points from $\Xi$. Markov chain samples from the $k$ densities
$\pi_{\xi_{i}}(\boldsymbol{\mu}|\boldsymbol{y}),i=1,\dots,k$ corresponding to
this preliminary skeleton set are generated. We evaluate (1.3) with
$N_{l}=1000,n_{l}=1000,l=1,\dots,5$ for all points in $\Xi$, and retain only
those points for which the value of (1.3) is not less than 60% of the maximum
value. The maximum value is attained at $\tilde{\xi}=(1,425,2.25,0.6)$ and
there are 33 points satisfying this criterion. These points form the search
set $\tilde{\Xi}$ which is a subset of

$\tilde{\Xi}\subset\\{1\\}\times\\{100,425,\ldots,1400\\}\times\\{1.50,2.25,\ldots,4.50\\}\times\\{0.35,0.60,\ldots,2.10\\}.$

Our aim is to choose $k=5$ elements from $\tilde{\Xi}$, one of which must be
$\tilde{\xi}$, to form the skeleton sets, using the methods discussed in
Section 3, with the objective of estimating the ratios of marginal densities
in $\Xi$ relative to $\tilde{\xi}$. The NIS method with samples from the
posterior density $\pi_{\tilde{\xi}}(\boldsymbol{\mu}|\boldsymbol{y})$ is
considered for comparison.

The SFE optimal set is computed on $\tilde{\Xi}$ after each dimension is
scaled in $[0,1]$. For SFS, we write (3.2b) as an integral over
$(\boldsymbol{z},\log\sigma^{2})$, because the prior for $\boldsymbol{z}$ is
multivariate normal and $\sigma^{2}>0$, and use the approximation given in
Section A of the supplementary materials. The SEQ, MNX, and ENT optimal sets
are computed iteratively. At each iteration, an estimate of the asymptotic
relative SE is computed using Theorem 2(a) based on 1000 samples for the first
stage and 1000 samples for the second stage. Further computational details
about our implementation are provided in Section K.1.1 and plots of the
optimal sets are given in Figures 4-9 of the supplementary materials.

Finally, for each obtained skeleton set, we generate new samples which we use
to estimate the ratio of the marginal likelihoods and its relative SE for all
$\xi\in\Xi$. We generate a total of 50,000 samples which are equally divided
between each proposal density. Further computational details are given in
Section K.1.2 of the supplementary materials. The maximum relative SE
estimates corresponding to one component of $\xi$ fixed across the other
components are shown in Figure 3. It can be seen that across all parameters,
NIS and SEQ have the highest variance, and that SFS, MNX, ENT have the lowest
maximum variance.

Figure 3: Profile relative SE estimates for the Rongelap data. One parameter
is fixed and the maximum relative SE across the other parameters is plotted
against the fixed parameter. The crosses indicate points included in the
skeleton set.

##  5 Discussions

We consider situations where one is simultaneously interested in large number
of target distributions as in model selection and sensitivity analysis
examples. The multiple IS estimators with more than one proposal distributions
are particularly useful in this context. We address a crucial issue in
multiple IS estimation, namely the selection of proposal distributions in a
systematic way — which has not received much attention in the literature.
Three methods are described for choosing representative proposal
distributions. The first approach based on a geometric ‘coverage’ criterion
and the second method based on the maximum asymptotic variance can be used for
any multiple IS estimators. The maximum entropy method is designed for the
two-stage generalized IS estimators of Doss, (2010). We compare the
performance of these three methods in terms of Monte Carlo standard errors for
the estimators at all target pdfs. In particular, we undertake these
comparisons in the context of three detailed examples involving autologistic
model, Bayesian robit model, and Bayesian Poisson spatial GLMM. Our results
show that multiple IS estimators with proposals selected by the described
methods outperform IS based on a single proposal density.

The proposed methods except the one based on the geometric coverage criterion
use standard errors for the generalized IS method and the RL regression
method. We construct consistent SV estimators for these standard errors.
Although these standard errors are used here for selection of proposal
distributions, these are important in their own rights as they are required
for assessing the quality of the generalized IS estimators and the reverse
logistic estimator. The SV estimation of standard errors is illustrated using
a negative binomial GLM example.

Supplementary materials

##  Appendix A A modified Laplace approximation for Kullback-Leibler
divergence

In this section, we describe a modified Laplace approximation for SKLD defined
in the paper. Let $\textsf{X}=\mathbb{R}^{r}$, for some $r\geq 1$, and $\mu$
be the Lebesgue measure. Consider the SKLD between two densities
$\pi_{\xi_{1}}(x)=\nu_{\xi_{1}}(x)/c_{\xi_{1}}$ and
$\pi_{\xi_{2}}(x)=\nu_{\xi_{2}}(x)/c_{\xi_{2}}$, with the assumption
$\log\nu_{\xi_{i}}(x)=O(M)$ for some $M,i=1,2$. Note that

$\Upsilon(\xi_{1},\xi_{2})=\frac{M\int_{\textsf{X}}J(x)\exp(G(x))\mu(dx)}{\int_{\textsf{X}}\exp(G(x))\mu(dx)}-\frac{M\int_{\textsf{X}}J(x)\exp(H(x))\mu(dx)}{\int_{\textsf{X}}\exp(H(x))\mu(dx)},$
(A.1)

where $G(x)=\log\nu_{\xi_{1}}(x)$, $H(x)=\log\nu_{\xi_{2}}(x)$, and
$J(x)=(G(x)-H(x))/M$. We apply Laplace approximation on each integral in (A.1)
separately. Specifically, we expand the integrals in the first term around
$\hat{x}=\operatorname*{argmax}_{x\in\textsf{X}}G(x)$ and the integrals in the
second term around $\tilde{x}=\operatorname*{argmax}_{x\in\textsf{X}}H(x)$.
Let $\hat{J}$ and $\tilde{J}$ denote $J$ evaluated at $\hat{x}$ and
$\tilde{x}$ respectively. We denote $\hat{G}_{i}=\frac{\partial}{\partial
x_{i}}G(x)|_{x=\hat{x}}$ and similarly $\hat{G}_{ij}$ for second order partial
derivatives and so on. We also denote $\hat{G}^{-1}_{ij}$ to be the $(i,j)$th
element of the inverse of the matrix with elements $\hat{G}_{ij}$’s. Then by
an application of (17) from Evangelou et al., (2011), we have

$\frac{\int_{\textsf{X}}J(x)\exp(G(x))\mu(dx)}{\int_{\textsf{X}}\exp(G(x))\mu(dx)}\approx\hat{J}+\frac{1}{2}\hat{J}_{i_{1}}\hat{G}_{i_{2}i_{3}i_{4}}\hat{G}^{-1}_{i_{1}i_{2}}\hat{G}^{-1}_{i_{3}i_{4}}-\frac{1}{2}\hat{J}_{i_{1}i_{2}}\hat{G}^{-1}_{i_{1}i_{2}}$

with an implicit summation $i_{1},\dots,i_{4}\in\\{1,\dots,r\\}$. A similar
approximation is derived for the second term:

$\frac{\int_{\textsf{X}}J(x)\exp(H(x))\mu(dx)}{\int_{\textsf{X}}\exp(H(x))\mu(dx)}\approx\tilde{J}+\frac{1}{2}\tilde{J}_{i_{1}}\tilde{H}_{i_{2}i_{3}i_{4}}\tilde{H}^{-1}_{i_{1}i_{2}}\tilde{H}^{-1}_{i_{3}i_{4}}-\frac{1}{2}\tilde{J}_{i_{1}i_{2}}\tilde{H}^{-1}_{i_{1}i_{2}}.$

The first order approximation to $\mbox{SKLD}(\xi_{1},\xi_{2})$ is
$M(\hat{J}-\tilde{J})$, which may be sufficient, but not if
$\hat{x}=\tilde{x}$. Note that, the second order approximation is exact for
two Gaussian densities.

##  Appendix B Proof of Theorem 1

From Roy et al., (2018), we only need to show
$\widehat{\Omega}\buildrel\text{a.s.}\over{\longrightarrow}\Omega$ where the
SV estimator $\widehat{\Omega}$ is defined in (3.8) and the $k\times k$ matrix
$\Omega$, following Roy et al., (2018), is defined through

$\Omega_{rs}=\sum_{l=1}^{k}\frac{a_{l}^{2}}{\tilde{s}_{l}}\Big{[}E_{q_{l}}\\{Y_{1}^{(r,l)}Y_{1}^{(s,l)}\\}+\sum_{i=1}^{\infty}E_{q_{l}}\\{Y_{1}^{(r,l)}Y_{1+i}^{(s,l)}\\}+\sum_{i=1}^{\infty}E_{q_{l}}\\{Y_{1+i}^{(r,l)}Y_{1}^{(s,l)}\\}\Big{]},$

for $r,s=1,\ldots,k$, where, $N_{l}/N\rightarrow\tilde{s}_{l}$ and for
$r,l=1,\dots,k$,

$Y_{i}^{(r,l)}\equiv
p_{r}(X_{i}^{(l)},\boldsymbol{\zeta}^{0})-E_{q_{l}}\bigl{(}p_{r}(X,\boldsymbol{\zeta}^{0})\bigr{)},\qquad
i=1,\ldots,n_{l}.$

As in Roy et al., (2018), this will be proved in couple of steps. First, we
consider a single chain $\tilde{\Phi}_{l}$ used to calculate $k$ quantities.
We use the results in Vats et al., (2018) who obtain conditions for the
multivariate SV estimator to be strongly consistent. Second, we combine
results from the $k$ independent chains. Finally, we show that
$\widehat{\Omega}$ is a strongly consistent estimator of $\Omega$.

Denote
$\bar{Y}^{(l)}=\left(\bar{Y}^{(1,l)},\bar{Y}^{(2,l)},\dots,\bar{Y}^{(k,l)}\right)^{\top}$
where $\bar{Y}^{(r,l)}=\sum_{i=1}^{N_{l}}Y_{i}^{(r,l)}/N_{l}$. From Roy et
al., (2018) we have $\sqrt{N_{l}}\bar{Y}^{(l)}\buildrel
d\over{\rightarrow}{\cal N}(0,\Sigma^{(l)})$ as $N_{l}\rightarrow\infty$,
where $\Sigma^{(l)}$ is a $k\times k$ covariance matrix with

$\Sigma^{(l)}_{rs}=E_{q_{l}}\\{Y_{1}^{(r,l)}Y_{1}^{(s,l)}\\}+\sum_{i=1}^{\infty}E_{q_{l}}\\{Y_{1}^{(r,l)}Y_{1+i}^{(s,l)}\\}+\sum_{i=1}^{\infty}E_{q_{l}}\\{Y_{1+i}^{(r,l)}Y_{1}^{(s,l)}\\}.$
(B.1)

The SV estimator of $\Sigma^{(l)}$ is given in (3.7). We now prove the strong
consistency of $\widehat{\Sigma}^{(l)}$. Note that $\widehat{\Sigma}^{(l)}$ is
defined using the terms $\bar{Z}^{(l)}_{i}$’s which involve the random
quantity $\boldsymbol{\hat{\zeta}}$. We define
$\widehat{\Sigma}^{(l)}(\boldsymbol{\zeta}^{0})$ to be
$\widehat{\Sigma}^{(l)}$ with $\boldsymbol{\zeta}^{0}$ substituted for
$\boldsymbol{\hat{\zeta}}$, that is,

$\widehat{\Sigma}^{(l)}(\boldsymbol{\zeta}^{0})=\frac{1}{N_{l}}\sum_{j=-(b_{N_{l}}-1)}^{b_{N_{l}}-1}w_{N_{l}}(j)\sum_{i\in
S_{j,N}}\left[Y^{(l)}_{i}-\bar{Y}^{(l)}\right]\left[Y^{(l)}_{i+j}-\bar{Y}^{(l)}\right]^{\top}\mbox{
for}\;\;l=1,\ldots,k,$

where $Y^{(l)}_{i}=\left(Y^{(1,l)}_{i},\ldots,Y^{(k,l)}_{i}\right)^{\top}$. We
prove
$\widehat{\Sigma}^{(l)}\buildrel\text{a.s.}\over{\longrightarrow}\Sigma^{(l)}$
in two steps: (1)
$\widehat{\Sigma}^{(l)}(\boldsymbol{\zeta}^{0})\buildrel\text{a.s.}\over{\longrightarrow}\Sigma^{(l)}$
and (2)
$\widehat{\Sigma}^{(l)}-\widehat{\Sigma}^{(l)}(\boldsymbol{\zeta}^{0})\buildrel\text{a.s.}\over{\longrightarrow}0$.
Under the conditions of Theorem 1, it follows from Vats et al., (2018) that
$\widehat{\Sigma}^{(l)}(\boldsymbol{\zeta}^{0})\buildrel\text{a.s.}\over{\longrightarrow}\Sigma^{(l)}$
as $N_{l}\rightarrow\infty$. We show
$\widehat{\Sigma}_{rs}^{(l)}-\widehat{\Sigma}_{rs}^{(l)}(\boldsymbol{\zeta}^{0})\buildrel\text{a.s.}\over{\longrightarrow}0$
where $\widehat{\Sigma}_{rs}^{(l)}$ and
$\widehat{\Sigma}_{rs}^{(l)}(\boldsymbol{\zeta}^{0})$ are the $(r,s)$th
elements of the $k\times k$ matrices $\widehat{\Sigma}_{rs}^{(l)}$ and
$\widehat{\Sigma}_{rs}^{(l)}(\boldsymbol{\zeta}^{0})$ respectively. By the
mean value theorem (in multiple variables), there exists
$\boldsymbol{\zeta}^{*}=t\boldsymbol{\hat{\zeta}}+(1-t)\boldsymbol{\zeta}^{0}$
for some $t\in(0,1)$, such that

$\widehat{\Sigma}_{rs}^{(l)}-\widehat{\Sigma}_{rs}^{(l)}(\boldsymbol{\zeta}^{0})=\nabla\widehat{\Sigma}_{rs}^{(l)}(\boldsymbol{\zeta}^{*})\cdot(\boldsymbol{\hat{\zeta}}-\boldsymbol{\zeta}^{0}),$
(B.2)

where $\cdot$ represents the dot product. Note that

$\widehat{\Sigma}_{rs}^{(l)}(\boldsymbol{\zeta})=\frac{1}{N_{l}}\sum_{j=-(b_{N_{l}}-1)}^{b_{N_{l}}-1}w_{N_{l}}(j)\sum_{i}\left[Z^{(r,l)}_{i}(\boldsymbol{\zeta})-\bar{Z}^{(r,l)}(\boldsymbol{\zeta})\right]\left[Z^{(s,l)}_{i+j}(\boldsymbol{\zeta})-\bar{Z}^{(s,l)}(\boldsymbol{\zeta})\right],$

where
$Z^{(r,l)}_{i}(\boldsymbol{\zeta}):=p_{r}(X_{i}^{(l)},\boldsymbol{\zeta})$ and
$\bar{Z}^{(r,l)}(\boldsymbol{\zeta}):=\sum_{j=1}^{N_{l}}p_{r}(X_{j}^{(l)},\boldsymbol{\zeta})/N_{l}$.
Some calculations show that for $t\neq r$

$\frac{\partial
Z^{(r,l)}_{j}(\boldsymbol{\zeta})}{\partial\boldsymbol{\zeta}_{t}}=-p_{r}(X_{j}^{(l)},\boldsymbol{\zeta})p_{t}(X_{j}^{(l)},\boldsymbol{\zeta})$

and

$\frac{\partial
Z^{(r,l)}_{j}(\boldsymbol{\zeta})}{\partial\boldsymbol{\zeta}_{r}}=p_{r}(X_{j}^{(l)},\boldsymbol{\zeta})(1-p_{r}(X_{j}^{(l)},\boldsymbol{\zeta})).$

Simplifying the notations, we denote $U^{(r,t)}_{j}:=\partial
Z^{(r,l)}_{j}(\boldsymbol{\zeta})/\partial\boldsymbol{\zeta}_{t}$,
$\bar{U}^{r}:=\partial\bar{Z}^{(r,l)}(\boldsymbol{\zeta})/\partial\boldsymbol{\zeta}_{t}$
and simply write $Z^{(r,l)}_{j}$ and $\bar{Z}^{(r,l)}$ for
$Z^{(r,l)}_{j}(\boldsymbol{\zeta})$ and $\bar{Z}^{(r,l)}(\boldsymbol{\zeta})$
respectively. Thus we have

$\displaystyle\frac{\partial\widehat{\Sigma}_{rs}^{(l)}(\boldsymbol{\zeta})}{\partial\boldsymbol{\zeta}_{t}}=\frac{1}{N_{l}}\sum_{j=-(b_{N_{l}}-1)}^{b_{N_{l}}-1}w_{N_{l}}(j)\sum_{i}\left[(Z^{(r,l)}_{i}-\bar{Z}^{(r,l)})(U^{(s,t)}_{i+j}-\bar{U}^{(s,t)})+(U^{(r,t)}_{i}-\bar{U}^{(r,t)})(Z^{(s,l)}_{i+j}-\bar{Z}^{(s,l)})\right]$
$\displaystyle=$
$\displaystyle\frac{1}{N_{l}}\sum_{j=-(b_{N_{l}}-1)}^{b_{N_{l}}-1}w_{N_{l}}(j)\sum_{i}\left[(Z^{(r,l)}_{i}-\bar{Z}^{(r,l)})(U^{(s,t)}_{i+j}-\bar{U}^{(s,t)})\right]$
(B.3) $\displaystyle\hskip
72.26999pt+\frac{1}{N_{l}}\sum_{j=-(b_{N_{l}}-1)}^{b_{N_{l}}-1}w_{N_{l}}(j)\sum_{i}\left[(U^{(r,t)}_{i}-\bar{U}^{(r,t)})(Z^{(s,l)}_{i+j}-\bar{Z}^{(s,l)})\right],$
(B.4)

Let $V_{i}^{(l)}:=(Z^{(r,l)}_{i},U^{(s,t)}_{i})^{T}$ and

$\widehat{\Sigma}^{(l)}_{V}(\boldsymbol{\zeta})=\frac{1}{N_{l}}\sum_{j=-(b_{N_{l}}-1)}^{b_{N_{l}}-1}w_{N_{l}}(j)\sum_{i}\left[V^{(l)}_{i}-\bar{V}^{(l)}\right]\left[V^{(l)}_{i+j}-\bar{V}^{(l)}\right]^{\top}.$

Since $p_{r}(X,\boldsymbol{\zeta})$ is uniformly bounded by 1 and
$\tilde{\Phi}_{l}$ is polynomially ergodic of order $m>1$, from Vats et al.,
(2018) we know that
$\widehat{\Sigma}^{(l)}_{V}(\boldsymbol{\zeta})\buildrel\text{a.s.}\over{\longrightarrow}\Sigma^{(l)}_{V}(\boldsymbol{\zeta})$
where $\Sigma^{(l)}_{V}(\boldsymbol{\zeta})$ is the covariance matrix of the
asymptotic distribution of $\sqrt{N_{l}}(\bar{V}^{(l)}-E_{q_{l}}V)$. Since the
expression in (B.3) is the off-diagonal elements of
$\widehat{\Sigma}^{(l)}_{V}(\boldsymbol{\zeta})$, it is bounded with
probability one. We can similarly see that the expression in (B.4) is bounded
with probability one. Note that, the proof to show that
$\partial\widehat{\Sigma}_{rs}^{(l)}(\boldsymbol{\zeta})/\partial\boldsymbol{\zeta}_{t}$
is bounded with probability one is quite different from the proof in Roy et
al., (2018).

Note that the terms $Z^{(r,l)}_{i},U_{i}^{(r,t)}$, etc, above actually depends
on $\boldsymbol{\zeta}$, and we are indeed concerned with the case where
$\boldsymbol{\zeta}$ takes on the value $\boldsymbol{\zeta}^{*}$, lying
between $\boldsymbol{\hat{\zeta}}$ and $\boldsymbol{\zeta}^{0}$. Since,
$\boldsymbol{\hat{\zeta}}\buildrel\text{a.s.}\over{\longrightarrow}\boldsymbol{\zeta}^{0}$,
we have
$\boldsymbol{\zeta}^{*}\buildrel\text{a.s.}\over{\longrightarrow}\boldsymbol{\zeta}^{0}$
as $N_{l}\rightarrow\infty$. Let $\|u\|_{L_{1}}$ denotes the $L_{1}$ norm of a
vector $u\in\mathbb{R}^{k}$. So from (B.2), and the fact that
$\partial\widehat{\Sigma}_{rs}^{(l)}(\boldsymbol{\zeta})/\partial\boldsymbol{\zeta}_{t}$
is bounded with probability one, we have

$|\widehat{\Sigma}_{rs}^{(l)}-\widehat{\Sigma}_{rs}^{(l)}(\boldsymbol{\zeta}^{0})|\leq\underset{1\leq
t\leq
k}{\max}\left\\{\left|\frac{\partial\widehat{\Sigma}_{rs}^{(l)}(\boldsymbol{\zeta}^{*})}{\partial\boldsymbol{\zeta}_{t}}\right|\right\\}\|\boldsymbol{\hat{\zeta}}-\boldsymbol{\zeta}^{0}\|_{L_{1}}\buildrel\text{a.s.}\over{\longrightarrow}0\;\;\text{
as}\;n\rightarrow\infty.$

Let

$\widehat{\Sigma}=\begin{pmatrix}\widehat{\Sigma}^{(1)}\vtop{\hbox
to0.0pt{\smash{\kern 6.0pt\raise-1.5pt\hbox to0.0pt{\hss.}\kern
6.0pt\raise-3.0pt\hbox to0.0pt{\hss.}\kern 6.0pt\raise-4.5pt\hbox
to0.0pt{\hss.}\kern 6.0pt\raise-6.0pt\hbox to0.0pt{\hss.}\kern
6.0pt\raise-7.5pt\hbox to0.0pt{\hss.}\kern 6.0pt\raise-9.0pt\hbox
to0.0pt{\hss.}\kern 6.0pt\raise-10.5pt\hbox to0.0pt{\hss.}\kern
6.0pt\raise-12.0pt\hbox to0.0pt{\hss.}\kern 6.0pt\raise-13.5pt\hbox
to0.0pt{\hss.}\kern 6.0pt\raise-15.0pt\hbox to0.0pt{\hss.}\kern
6.0pt\raise-16.5pt\hbox to0.0pt{\hss.}\kern 6.0pt}\hss}}&\quad&\quad&\quad
0&\\\ &&&&\\\ &0&&&\widehat{\Sigma}^{(k)}\\\ \end{pmatrix}.$

Since
$\widehat{\Sigma}^{(l)}\buildrel\text{a.s.}\over{\longrightarrow}\Sigma^{(l)},$
for $l=1,\dots,k$, it follows that
$\widehat{\Sigma}\buildrel\text{a.s.}\over{\longrightarrow}\Sigma$ where
$\Sigma$ is the corresponding $k^{2}\times k^{2}$ covariance matrix, that is,
$\Sigma$ is a block diagonal matrix as $\widehat{\Sigma}$ with $\Sigma^{(l)}$
substituted for $\widehat{\Sigma}^{(l)},l=1,\dots,k$. Define the following
$k\times k^{2}$ matrix

$A_{N}=\left(-\sqrt{\frac{N}{N_{1}}}a_{1}I_{k}\quad-\sqrt{\frac{N}{N_{2}}}a_{2}I_{k}\quad\dots\quad-\sqrt{\frac{N}{N_{k}}}a_{k}I_{k}\right)\;,$

where $I_{k}$ denotes the $k\times k$ identity matrix. Then we have
$\widehat{\Omega}\equiv
A_{N}\widehat{\Sigma}A_{N}^{T}\buildrel\text{a.s.}\over{\longrightarrow}\Omega$
as $N\rightarrow\infty$.

##  Appendix C Proof of Theorem 2 (a)

From Roy et al., (2018) we know that
$\sigma^{2}_{u}=rc(\pi;\boldsymbol{d})^{\top}Vc(\pi;\boldsymbol{d})+\tau^{2}(\pi;\boldsymbol{d})$,
where
$\tau^{2}(\pi;\boldsymbol{d})=\sum_{l=1}^{k}(a_{l}^{2}/s_{l})\tau^{2}_{l}(\pi;\boldsymbol{d})$,
and

$\tau^{2}_{l}(\pi;\boldsymbol{d})=\mbox{Var}_{q_{l}}(u^{\pi}(X_{1}^{(l)};\boldsymbol{d}))+2\sum_{g=1}^{\infty}\mbox{Cov}_{q_{l}}(u^{\pi}(X_{1}^{(l)};\boldsymbol{d}),u^{\pi}(X_{1+g}^{(l)};\boldsymbol{d})).$
(C.1)

To prove Theorem 2 (a), note that, we already have a consistent SV estimator
$\widehat{V}$ of $V$. From Roy et al., (2018) it follows that
$\hat{c}(\pi;\boldsymbol{\hat{d}})^{\top}\widehat{V}\hat{c}(\pi;\boldsymbol{\hat{d}})\buildrel\text{a.s.}\over{\longrightarrow}c(\pi;\boldsymbol{d})^{\top}Vc(\pi;\boldsymbol{d})$.

We now show $\hat{\tau}_{l}^{2}(\pi;\boldsymbol{\hat{d}})$ is a consistent
estimator of $\tau_{l}^{2}(\pi;\boldsymbol{d})$ where $\hat{\tau}^{2}_{l}$ is
defined in (3.15). Since the Markov chains
$\\{X_{i}^{(l)}\\}_{i=1}^{n_{l}},l=1,\dots,k$ are independent, it then follows
that $\tau^{2}(\pi;\boldsymbol{d})$ is consistently estimated by
$\hat{\tau}^{2}(\pi;\boldsymbol{\hat{d}})$ completing the proof of Theorem 2
(a).

If $\boldsymbol{d}$ is known from the assumptions of Theorem 2 and the results
in Vats et al., (2018), we know that $\tau_{l}^{2}(\pi;\boldsymbol{d})$ is
consistently estimated by its SV estimator
$\hat{\tau}_{l}^{2}(\pi;\boldsymbol{d})$. Note that,
$\hat{\tau}_{l}^{2}(\pi;\boldsymbol{d})$ is defined in terms of the quantities
$u^{\pi}(X_{i}^{(l)};\boldsymbol{d})$’s. We now show that
$\hat{\tau}_{l}^{2}(\pi;\boldsymbol{\hat{d}})-\hat{\tau}_{l}^{2}(\pi;\boldsymbol{d})\buildrel\text{a.s.}\over{\longrightarrow}0.$
Let

${\partial U}_{i}^{m}(\boldsymbol{z}):=\frac{\partial
u_{i}(\boldsymbol{z})}{\partial
z_{m}}=\frac{a_{m}}{z_{m}^{2}}\frac{\nu(X_{i}^{(l)})\varphi_{m}(X_{i}^{(l)})}{\left(\sum_{s}a_{s}\varphi_{s}(X_{i}^{(l)})/z_{s}\right)^{2}},$

and $\bar{\partial U}^{m}(\boldsymbol{z})$ be the averages of $\\{{\partial
U}_{i}^{m}(\boldsymbol{z}),i=1,\dots,n_{l}\\}$. Denoting
$\hat{\tau}_{l}^{2}(\pi;\boldsymbol{z})$ by $G(\boldsymbol{z})$, by the mean
value theorem (in multiple variables), there exists
$\boldsymbol{d}^{*}=t\boldsymbol{\hat{d}}+(1-t)\boldsymbol{d}$ for some
$t\in(0,1)$, such that $G(\boldsymbol{\hat{d}})-G(\boldsymbol{d})=\nabla
G(\boldsymbol{d}^{*})\cdot(\boldsymbol{\hat{d}}-\boldsymbol{d})$. For any
$m\in\\{2,\cdots,k\\}$, and $\boldsymbol{z}\in{R^{+}}^{k-1}$,

$\displaystyle\frac{\partial G(\boldsymbol{z})}{\partial
z_{m}}=\sum_{j=-(b_{n_{l}}-1)}^{b_{n_{l}}-1}w_{n_{l}}(j)\sum_{i}\left[u_{i}(\boldsymbol{z})-\bar{u}(\boldsymbol{z})\right]\left[{\partial
U}^{m}_{i+j}(\boldsymbol{z})-\bar{\partial U}^{m}(\boldsymbol{z})\right]$
$\displaystyle\hskip
43.36243pt+\sum_{j=-(b_{n_{l}}-1)}^{b_{n_{l}}-1}w_{n_{l}}(j)\sum_{i}\left[{\partial
U}^{m}_{i}(\boldsymbol{z})-\bar{\partial
U}^{m}(\boldsymbol{z})\right]\left[u_{i+j}(\boldsymbol{z})-\bar{u}(\boldsymbol{z})\right]$

Then using similar arguments as in the proof of Theorem 1, it can be shown
that $\partial G(\boldsymbol{z})/\partial z_{m}$ is bounded with probability
one. Then it follows that

$|G(\boldsymbol{\hat{d}})-G(\boldsymbol{d})|\leq\underset{1\leq m\leq
k-1}{\max}\left\\{\left|\frac{\partial G(\boldsymbol{d}^{*})}{\partial
z_{m}}\right|\right\\}\|\boldsymbol{\hat{d}}-\boldsymbol{d}\|_{L_{1}}\buildrel\text{a.s.}\over{\longrightarrow}0.$

##  Appendix D Proof of Theorem 2 (b)

From Roy et al., (2018) we know that
$\sigma^{2}_{\eta}=re(\pi;\boldsymbol{d})^{\top}Ve(\pi;\boldsymbol{d})+\rho(\pi;\boldsymbol{d})$,
where

$\rho(\pi;\boldsymbol{d})=\nabla
h(E_{\pi}fu(\pi,\pi_{1}),u(\pi,\pi_{1}))^{\top}\Gamma(\pi;\boldsymbol{d})\nabla
h(E_{\pi}fu(\pi,\pi_{1}),u(\pi,\pi_{1})),$
$\Gamma(\pi;\boldsymbol{d})=\sum_{l=1}^{k}\frac{a_{l}^{2}}{s_{l}}\Gamma_{l}(\pi;\boldsymbol{d});\Gamma_{l}(\pi;\boldsymbol{d})=\left(\begin{array}[]{lr}\gamma^{11}&\gamma^{12}\\\
\gamma^{21}&\gamma^{22}\\\ \end{array}\right),$

with

$\displaystyle\gamma^{11}\equiv\gamma^{11}_{l}(\pi;\boldsymbol{d})$
$\displaystyle=\mbox{Var}_{q_{l}}(v^{[f],\pi}(X_{1}^{(l)};\boldsymbol{d}))+2\sum_{g=1}^{\infty}\mbox{Cov}_{q_{l}}(v^{[f],\pi}(X_{1}^{(l)};\boldsymbol{d}),v^{[f],\pi}(X_{1+g}^{(l)};\boldsymbol{d})),$
$\displaystyle\gamma^{12}\equiv\gamma^{12}_{l}(\pi;\boldsymbol{d})$
$\displaystyle=\gamma^{21}\equiv\gamma^{21}_{l}(\pi;\boldsymbol{d})$
$\displaystyle=\mbox{Cov}_{q_{l}}(v^{[f],\pi}(X_{1}^{(l)};\boldsymbol{d}),u^{\pi}(X_{1}^{(l)};\boldsymbol{d}))$
$\displaystyle\quad{}+\sum_{g=1}^{\infty}[\mbox{Cov}_{q_{l}}(v^{[f],\pi}(X_{1}^{(l)};\boldsymbol{d}),u^{\pi}(X_{1+g}^{(l)};\boldsymbol{d}))+\mbox{Cov}_{q_{l}}(v^{[f],\pi}(X_{1+g}^{(l)};\boldsymbol{d}),u^{\pi}(X_{1}^{(l)};\boldsymbol{d}))]$
$\displaystyle\gamma^{22}_{l}\equiv\gamma^{22}_{l}(\pi;\boldsymbol{d})$
$\displaystyle=\mbox{Var}_{q_{l}}(u^{\pi}(X_{1}^{(l)};\boldsymbol{d}))+2\sum_{g=1}^{\infty}\mbox{Cov}_{q_{l}}(u^{\pi}(X_{1}^{(l)};\boldsymbol{d}),u^{\pi}(X_{1+g}^{(l)};\boldsymbol{d})).$

From Roy et al., (2018) we know that
$\hat{e}(\pi;\boldsymbol{\hat{d}})^{\top}\widehat{V}\hat{e}(\pi;\boldsymbol{\hat{d}})\buildrel\text{a.s.}\over{\longrightarrow}e(\pi;\boldsymbol{d})^{\top}Ve(\pi;\boldsymbol{d})$.
Thus, to prove Theorem 2 (b), we only need to show that
$\widehat{\Gamma}_{l}(\pi;\boldsymbol{\hat{d}})\buildrel\text{a.s.}\over{\longrightarrow}\Gamma_{l}(\pi;\boldsymbol{d})$.
Note that,

$\displaystyle\widehat{\Gamma}_{l}(\pi;\boldsymbol{d})$
$\displaystyle=\frac{1}{n_{l}}\sum_{j=-(b_{n_{l}}-1)}^{b_{n_{l}}-1}w_{n_{l}}(j)\sum_{i\in
S_{j,n}}\left(\begin{array}[]{cc}\left[v^{[f]}_{i}-\bar{v}^{[f]}\right]\left[v^{[f]}_{i+j}-\bar{v}^{[f]}\right]&\left[v^{[f]}_{i}-\bar{v}^{[f]}\right]\left[u_{i+j}-\bar{u}\right]\\\
\left[v^{[f]}_{i+j}-\bar{v}^{[f]}\right]\left[u_{i}-\bar{u}\right]&\left[u_{i}-\bar{u}\right]\left[u_{i+j}-\bar{u}\right]\\\
\end{array}\right).$
$\displaystyle=\left(\begin{array}[]{cc}\hat{\gamma}^{11}(\pi;\boldsymbol{d})&\hat{\gamma}^{12}(\pi;\boldsymbol{d})\\\
\hat{\gamma}^{21}(\pi;\boldsymbol{d})&\hat{\gamma}^{22}(\pi;\boldsymbol{d})\\\
\end{array}\right).$

If $\boldsymbol{d}$ is known, from the assumptions of Theorem 2 (b) and the
results in Vats et al., (2018), we know that $\Gamma_{l}(\pi;\boldsymbol{d})$
is consistently estimated by its SV estimator
$\widehat{\Gamma}_{l}(\pi;\boldsymbol{d})$. We now show that
$\widehat{\Gamma}_{l}(\pi;\boldsymbol{\hat{d}})-\widehat{\Gamma}_{l}(\pi;\boldsymbol{d})\buildrel\text{a.s.}\over{\longrightarrow}0.$

From Theorem 2 (a), we know that
$\hat{\gamma}_{l}^{22}(\pi;\boldsymbol{\hat{d}})-\hat{\gamma}_{l}^{22}(\pi;\boldsymbol{d})\buildrel\text{a.s.}\over{\longrightarrow}0$
as $\gamma^{22}_{l}$ is the same as $\tau^{2}_{l}(\pi;\boldsymbol{d})$ defined
in (C.1). We now show
$\hat{\gamma}_{l}^{11}(\pi;\boldsymbol{\hat{d}})-\hat{\gamma}_{l}^{11}(\pi;\boldsymbol{d})\buildrel\text{a.s.}\over{\longrightarrow}0$.

Let

${\partial V}_{i}^{[f],m}(\boldsymbol{z}):=\frac{\partial
v^{[f]}_{i}(\boldsymbol{z})}{\partial
z_{m}}=\frac{a_{m}}{z_{m}^{2}}\frac{f(X_{i}^{(l)})\nu(X_{i}^{(l)})\varphi_{m}(X_{i}^{(l)})}{\left(\sum_{s}a_{s}\varphi_{s}(X_{i}^{(l)})/z_{s}\right)^{2}},$

and $\bar{\partial V}^{[f],m}(\boldsymbol{z})$ be the averages of
$\\{{\partial V}_{i}^{[f],m}(\boldsymbol{z}),i=1,\dots,n_{l}\\}$.

Letting $\hat{\gamma}_{l}^{11}(\pi;\boldsymbol{z})$ by $H(\boldsymbol{z})$, by
the mean value theorem, there exists
$\boldsymbol{d}^{*}=t\boldsymbol{\hat{d}}+(1-t)\boldsymbol{d}$ for some
$t\in(0,1)$, such that $H(\boldsymbol{\hat{d}})-H(\boldsymbol{d})=\nabla
H(\boldsymbol{d}^{*})\cdot(\boldsymbol{\hat{d}}-\boldsymbol{d})$. For any
$m\in\\{2,\cdots,k\\}$, and $\boldsymbol{z}\in{R^{+}}^{k-1}$,

$\displaystyle\frac{\partial H(\boldsymbol{z})}{\partial
z_{m}}=\sum_{j=-(b_{n_{l}}-1)}^{b_{n_{l}}-1}w_{n_{l}}(j)\sum_{i}\left[v^{[f]}_{i}(\boldsymbol{z})-\bar{v}^{[f]}(\boldsymbol{z})\right]\left[{\partial
V}^{[f],m}_{i+j}(\boldsymbol{z})-\bar{\partial
V}^{[f],m}(\boldsymbol{z})\right]$ $\displaystyle\hskip
43.36243pt+\sum_{j=-(b_{n_{l}}-1)}^{b_{n_{l}}-1}w_{n_{l}}(j)\sum_{i}\left[{\partial
V}^{[f],m}_{i}(\boldsymbol{z})-\bar{\partial
V}^{[f],m}(\boldsymbol{z})\right]\left[v^{[f]}_{i+j}(\boldsymbol{z})-\bar{v}^{[f]}(\boldsymbol{z})\right]$

The rest of the proof is analogous to Theorem 2 (a) and thus we have
$\hat{\gamma}_{l}^{11}(\pi;\boldsymbol{\hat{d}})-\hat{\gamma}_{l}^{11}(\pi;\boldsymbol{d})\buildrel\text{a.s.}\over{\longrightarrow}0$.
Finally, using similar arguments as before we can show
$\hat{\gamma}_{l}^{12}(\pi;\boldsymbol{\hat{d}})-\hat{\gamma}_{l}^{12}(\pi;\boldsymbol{d})\buildrel\text{a.s.}\over{\longrightarrow}0$.

##  Appendix E Proof of Theorem 3

Since the Markov chains used in stage 1 are polynomially ergodic of order
$m>1$, from Roy et al., (2018, Theorem 1), we have
$N^{1/2}(\boldsymbol{\hat{d}}_{\boldsymbol{q}}-\boldsymbol{d}_{\boldsymbol{q}})\buildrel
d\over{\rightarrow}{\cal N}(0,V_{\boldsymbol{q}})$. Since $n/N\rightarrow r$,
it follows that
$\sqrt{n}(\boldsymbol{\hat{d}}_{\boldsymbol{q}}-\boldsymbol{d}_{\boldsymbol{q}})\buildrel
d\over{\rightarrow}{\cal N}(0,rV_{\boldsymbol{q}})$. Following Roy et al.,
(2018, Proof of Theorem 2) we write

$\begin{split}\sqrt{n}(\hat{\bf
u}(\boldsymbol{\pi};\boldsymbol{\hat{d}}_{\boldsymbol{q}})-{\bf
u}(\boldsymbol{\pi},q_{1}))=\sqrt{n}(\hat{\bf
u}(\boldsymbol{\pi};\boldsymbol{\hat{d}}_{\boldsymbol{q}})-\hat{\bf
u}(\boldsymbol{\pi};\boldsymbol{d}_{\boldsymbol{q}}))+\sqrt{n}(\hat{\bf
u}(\boldsymbol{\pi};\boldsymbol{d}_{\boldsymbol{q}})-{\bf
u}(\boldsymbol{\pi},q_{1})).\end{split}$ (E.1)

Note that the 2nd term involves randomness only from the 2nd stage Markov
chains. Since $\sum_{l=1}^{k}a_{l}E_{q_{l}}{\bf
u}^{\boldsymbol{\pi}}(X;\boldsymbol{d}_{\boldsymbol{q}})={\bf
u}(\boldsymbol{\pi},q_{1})$, we have

$\sqrt{n}(\hat{\bf u}(\boldsymbol{\pi};\boldsymbol{d}_{\boldsymbol{q}})-{\bf
u}(\boldsymbol{\pi},q_{1}))=\sum_{l=1}^{k}a_{l}\sqrt{\frac{n}{n_{l}}}\frac{\sum_{i=1}^{n_{l}}({\bf
u}^{\boldsymbol{\pi}}(X_{i}^{(l)};\boldsymbol{d}_{\boldsymbol{q}})-E_{q_{l}}{\bf
u}^{\boldsymbol{\pi}}(X;\boldsymbol{d}_{\boldsymbol{q}}))}{\sqrt{n_{l}}}.$

Since $\Phi_{l}$ is polynomially ergodic of order $m$ and
$E_{q_{l}}|u^{\pi}(X;\boldsymbol{d}_{\boldsymbol{q}})|^{2+\delta}$ is finite
for each $\pi\in\Pi$ where $m>1+2/\delta$, it follows that
$\sum_{i=1}^{n_{l}}({\bf
u}^{\boldsymbol{\pi}}(X_{i}^{(l)};\boldsymbol{d}_{\boldsymbol{q}})-E_{q_{l}}{\bf
u}^{\boldsymbol{\pi}}(X;\boldsymbol{d}_{\boldsymbol{q}}))/\sqrt{n_{l}}\buildrel
d\over{\rightarrow}N(0,T_{l}(\boldsymbol{d}_{\boldsymbol{q}}))$ where
$T_{l}(\boldsymbol{d}_{\boldsymbol{q}})$ is the matrix with elements defined
in (3.16). As $n_{l}/n\rightarrow s_{l}$ and the Markov chains $\Phi_{l}$’s
are independent, it follows that $\sqrt{n}(\hat{\bf
u}(\boldsymbol{\pi};\boldsymbol{d}_{\boldsymbol{q}})-{\bf
u}(\boldsymbol{\pi},q_{1}))\buildrel
d\over{\rightarrow}N(0,\sum_{l=1}^{k}(a_{l}^{2}/s_{l})T_{l}(\boldsymbol{d}_{\boldsymbol{q}}))$.

Next by Taylor series expansion of
$F(\boldsymbol{d})\equiv\hat{u}(\pi;\boldsymbol{d})$ about
$\boldsymbol{d}_{\boldsymbol{q}}$, we have

$\sqrt{n}(F(\boldsymbol{\hat{d}}_{\boldsymbol{q}})-F(\boldsymbol{d}_{\boldsymbol{q}}))=\sqrt{n}\nabla
F(\boldsymbol{d}_{\boldsymbol{q}})^{\top}(\boldsymbol{\hat{d}}_{\boldsymbol{q}}-\boldsymbol{d}_{\boldsymbol{q}})+\frac{\sqrt{n}}{2}(\boldsymbol{\hat{d}}_{\boldsymbol{q}}-\boldsymbol{d}_{\boldsymbol{q}})^{\top}\nabla^{2}F(\boldsymbol{d}^{*})(\boldsymbol{\hat{d}}_{\boldsymbol{q}}-\boldsymbol{d}_{\boldsymbol{q}}),$

where $\boldsymbol{d}^{*}$ is between $\boldsymbol{d}_{\boldsymbol{q}}$ and
$\boldsymbol{\hat{d}}_{\boldsymbol{q}}$. As in Roy et al., (2018), we can then
show that

$\sqrt{n}(\hat{u}(\pi;\boldsymbol{\hat{d}}_{\boldsymbol{q}})-\hat{u}(\pi;\boldsymbol{d}_{\boldsymbol{q}}))=\sqrt{q}c(\pi;\boldsymbol{d}_{\boldsymbol{q}})\sqrt{N}(\boldsymbol{\hat{d}}_{\boldsymbol{q}}-\boldsymbol{d}_{\boldsymbol{q}})+o_{p}(1).$

Accumulating the terms for all $\pi\in\Pi$, we have

$\sqrt{n}(\hat{\bf
u}(\boldsymbol{\pi};\boldsymbol{\hat{d}}_{\boldsymbol{q}})-\hat{\bf
u}(\boldsymbol{\pi};\boldsymbol{d}_{\boldsymbol{q}}))=\sqrt{q}C(\boldsymbol{\pi};\boldsymbol{d}_{\boldsymbol{q}})\sqrt{N}(\boldsymbol{\hat{d}}_{\boldsymbol{q}}-\boldsymbol{d}_{\boldsymbol{q}})+{\bf
o}_{p}(1).$

Thus for constant vectors $t_{1}$ and $t_{2}$ of dimensions $k-1$ and $|\Pi|$
respectively, we have

$\displaystyle
t_{1}^{\top}\sqrt{n}(\boldsymbol{\hat{d}}_{\boldsymbol{q}}-\boldsymbol{d}_{\boldsymbol{q}})+t_{2}^{\top}\sqrt{n}(\hat{\bf
u}(\boldsymbol{\pi};\boldsymbol{\hat{d}}_{\boldsymbol{q}})-{\bf
u}(\boldsymbol{\pi},q_{1}))$
$\displaystyle=\sqrt{r}(t_{1}^{\top}+t_{2}^{\top}C(\boldsymbol{\pi};\boldsymbol{d}_{\boldsymbol{q}}))\sqrt{N}(\boldsymbol{\hat{d}}_{\boldsymbol{q}}-\boldsymbol{d}_{\boldsymbol{q}})+\sum_{l=1}^{k}a_{l}\sqrt{\frac{n}{n_{l}}}\frac{\sum_{i=1}^{n_{l}}t_{2}^{\top}({\bf
u}^{\boldsymbol{\pi}}(X_{i}^{(l)};\boldsymbol{d}_{\boldsymbol{q}})-E_{q_{l}}{\bf
u}^{\boldsymbol{\pi}}(X;\boldsymbol{d}_{\boldsymbol{q}}))}{\sqrt{n_{l}}}+o_{p}(1)$
$\displaystyle\buildrel
d\over{\rightarrow}N(0,r(t_{1}^{\top}+t_{2}^{\top}C(\boldsymbol{\pi};\boldsymbol{d}_{\boldsymbol{q}}))V_{\boldsymbol{q}}(t_{1}+C(\boldsymbol{\pi};\boldsymbol{d}_{\boldsymbol{q}})^{\top}t_{2})+\sum_{l=1}^{k}(a_{l}^{2}/s_{l})t_{2}^{\top}T_{l}(\boldsymbol{d}_{\boldsymbol{q}})t_{2}),$
(E.2)

where the last step follows from the independence of the Markov chains
involved in the two stages. Note that the variance in (E) is the same as

$(t_{1}^{\top},t_{2}^{\top})\left(\begin{array}[]{cc}rV_{\boldsymbol{q}}&\Sigma_{12}\\\
\Sigma_{21}&\Sigma_{22}\\\
\end{array}\right)(t_{1}^{\top},t_{2}^{\top})^{\top}.$

Hence the Cramér-Wold device implies the joint CLT in (3.18). Thus Theorem 4
(a) is proved.

From the proofs of Theorems 1 and 2 (a), we know that
$r\widehat{C}(\boldsymbol{\pi};\boldsymbol{\hat{d}}_{\boldsymbol{q}})\widehat{V}_{\boldsymbol{q}}\widehat{C}(\boldsymbol{\pi};\boldsymbol{\hat{d}}_{\boldsymbol{q}})^{\top}$
is a consistent estimator of
$rC(\boldsymbol{\pi};\boldsymbol{d}_{\boldsymbol{q}})V_{\boldsymbol{q}}C(\boldsymbol{\pi};\boldsymbol{d}_{\boldsymbol{q}})^{\top}$.
If $\boldsymbol{d}_{\boldsymbol{q}}$ is known from the assumptions of Theorem
3 (b) and the results in Vats et al., (2018), we know that
$T_{l}(\boldsymbol{d}_{\boldsymbol{q}})$ is consistently estimated by its SV
estimator $\widehat{T}_{l}(\boldsymbol{d}_{\boldsymbol{q}})$ defined in
(3.17). Then using similar arguments as in the proof of Theorem 2 (a), we can
show that every element of
$\widehat{T}_{l}(\boldsymbol{d}_{\boldsymbol{q}})-\widehat{T}_{l}(\boldsymbol{\hat{d}}_{\boldsymbol{q}})$
converges to zero (a.e.). Hence Theorem 3 (b) is proved.

##  Appendix F Entropy decomposition for generalized IS estimators of means

In this section, we prove a result similar to Theorem 3 for
$\hat{\eta}^{[f]}(\pi;\boldsymbol{\hat{d}})$. Let
$\hat{{\boldsymbol{\eta}}}^{[f]}(\boldsymbol{\pi};\boldsymbol{\hat{d}}_{\boldsymbol{q}})$
be the vector of length $|\Pi|$ consisting of
$\hat{\eta}^{[f]}(\pi;\boldsymbol{\hat{d}}_{\boldsymbol{q}})$’s, $\pi\in\Pi$
in a fixed order. Similarly define $\hat{\bf
v}^{[f]}(\boldsymbol{\pi};\boldsymbol{\hat{d}}_{\boldsymbol{q}})$ and the
vector of true means ${\bf E}_{\boldsymbol{\pi}}f$. Let $p^{*}\equiv|\Pi|$.
Let $E(\boldsymbol{\pi};\boldsymbol{d}_{\boldsymbol{q}})$ be the
$p^{*}\times(k-1)$ matrix with rows $e(\pi;\boldsymbol{d}_{\boldsymbol{q}})$
(defined in Section 3.2 of the paper), $\pi\in\Pi$. Similarly, define
$\widehat{E}(\boldsymbol{\pi};\boldsymbol{d}_{\boldsymbol{q}})$ with rows
$\hat{e}(\pi;\boldsymbol{d}_{\boldsymbol{q}})$, $\pi\in\Pi$. Let ${\bf
v}^{[f],\boldsymbol{\pi}}(x;\boldsymbol{d}_{\boldsymbol{q}})$ be the $p$
dimensional vector consisting of
$v^{[f],\pi}(x;\boldsymbol{d}_{\boldsymbol{q}})$’s defined in (3.10) of the
paper, $\pi\in\Pi$. Define the $2p^{*}\times 2p^{*}$ matrix

$\Lambda_{l}(\boldsymbol{d}_{\boldsymbol{q}})=\left(\begin{array}[]{cc}\Lambda_{l}^{11}(\boldsymbol{d}_{\boldsymbol{q}})&\Lambda_{l}^{12}(\boldsymbol{d}_{\boldsymbol{q}})\\\
\Lambda_{l}^{21}(\boldsymbol{d}_{\boldsymbol{q}})&T_{l}(\boldsymbol{d}_{\boldsymbol{q}})\\\
\end{array}\right),$ (F.1)

where the elements of $\Lambda_{l}^{11}(\boldsymbol{d}_{\boldsymbol{q}})$ are
given by

$\displaystyle\lambda^{11}_{l}(\pi,\pi^{\prime};\boldsymbol{d}_{\boldsymbol{q}})=\mbox{Cov}_{q_{l}}(v^{[f],\pi}(X_{1}^{(l)};\boldsymbol{d}_{\boldsymbol{q}}),v^{[f],\pi^{\prime}}(X_{1}^{(l)};\boldsymbol{d}_{\boldsymbol{q}}))$
$\displaystyle+\sum_{g=1}^{\infty}\mbox{Cov}_{q_{l}}(v^{[f],\pi}(X_{1}^{(l)};\boldsymbol{d}_{\boldsymbol{q}}),v^{[f],\pi^{\prime}}(X_{1+g}^{(l)};\boldsymbol{d}_{\boldsymbol{q}}))$
$\displaystyle+\sum_{g=1}^{\infty}\mbox{Cov}_{q_{l}}(v^{[f],\pi}(X_{1+g}^{(l)};\boldsymbol{d}_{\boldsymbol{q}}),v^{[f],\pi^{\prime}}(X_{1}^{(l)};\boldsymbol{d}_{\boldsymbol{q}})),$

and the elements of $\Lambda_{l}^{12}(\boldsymbol{d}_{\boldsymbol{q}})$ are
given by

$\displaystyle\lambda^{12}_{l}(\pi,\pi^{\prime};\boldsymbol{d}_{\boldsymbol{q}})=\mbox{Cov}_{q_{l}}(v^{[f],\pi}(X_{1}^{(l)};\boldsymbol{d}_{\boldsymbol{q}}),u^{\pi^{\prime}}(X_{1}^{(l)};\boldsymbol{d}_{\boldsymbol{q}}))$
$\displaystyle+\sum_{g=1}^{\infty}\mbox{Cov}_{q_{l}}(v^{[f],\pi}(X_{1}^{(l)};\boldsymbol{d}_{\boldsymbol{q}}),u^{\pi^{\prime}}(X_{1+g}^{(l)};\boldsymbol{d}_{\boldsymbol{q}}))$
$\displaystyle+\sum_{g=1}^{\infty}\mbox{Cov}_{q_{l}}(v^{[f],\pi}(X_{1+g}^{(l)};\boldsymbol{d}_{\boldsymbol{q}}),u^{\pi^{\prime}}(X_{1}^{(l)};\boldsymbol{d}_{\boldsymbol{q}})).$

Also
$\Lambda_{l}^{21}(\boldsymbol{d}_{\boldsymbol{q}})=\Lambda_{l}^{12}(\boldsymbol{d}_{\boldsymbol{q}})^{\top}$
and let
$\Lambda(\boldsymbol{d}_{\boldsymbol{q}})=\sum_{l=1}^{k}(a_{l}^{2}/s_{l})\Lambda_{l}(\boldsymbol{d}_{\boldsymbol{q}})$.
Define a function ${\bf h}:\mathbb{R}^{2p^{*}}\rightarrow\mathbb{R}^{p^{*}}$
where

${\bf
h}(x_{1},\dots,x_{2p^{*}})=\Big{(}\frac{x_{1}}{x_{p^{*}+1}},\frac{x_{2}}{x_{p^{*}+2}},\dots,\frac{x_{p^{*}}}{x_{2p^{*}}}\Big{)}$

with its gradient given by

$\nabla{\bf
h}(x)=\begin{pmatrix}1/x_{p^{*}+1}&0&\ldots&0&-x_{1}/x^{2}_{p^{*}+1}&0&\ldots&0\\\
0&1/x_{p^{*}+2}&\ldots&0&0&-x_{2}/x^{2}_{p^{*}+2}&\ldots&0\\\
\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\ddots&\vdots\\\
0&0&\ldots&1/x_{2p^{*}}&0&0&\ldots&-x_{p^{*}}/x^{2}_{2p^{*}}\end{pmatrix}.$

Define the $p^{*}\times p^{*}$ matrix

$\boldsymbol{\rho}(\boldsymbol{d}_{\boldsymbol{q}})=\nabla{\bf h}({\bf
E}_{\boldsymbol{\pi}}f\odot{\bf u}(\boldsymbol{\pi},q_{1}),{\bf
u}(\boldsymbol{\pi},q_{1}))^{\top}\Lambda(\boldsymbol{d}_{\boldsymbol{q}})\nabla{\bf
h}({\bf E}_{\boldsymbol{\pi}}f\odot{\bf u}(\boldsymbol{\pi},q_{1}),{\bf
u}(\boldsymbol{\pi},q_{1})^{\top},$

where $\odot$ denotes element-wise multiplication. Let

$\widehat{\Lambda}_{l}(\boldsymbol{d}_{\boldsymbol{q}})=\frac{1}{n_{l}}\sum_{j=-(b_{n_{l}}-1)}^{b_{n_{l}}-1}w_{n_{l}}(j)\sum_{i\in
S_{j,n}}\left(\begin{array}[]{c}{\bf
v}^{[f],\boldsymbol{\pi}}(X_{i}^{(l)};\boldsymbol{d}_{\boldsymbol{q}})-\bar{{\bf
v}^{[f]}}(\boldsymbol{d}_{\boldsymbol{q}})\\\ {\bf
u}^{\boldsymbol{\pi}}(X_{i}^{(l)};\boldsymbol{d}_{\boldsymbol{q}})-\bar{{\bf
u}}(\boldsymbol{d}_{\boldsymbol{q}})\\\
\end{array}\right)\left(\begin{array}[]{c}{\bf
v}^{[f],\boldsymbol{\pi}}(X_{i+j}^{(l)};\boldsymbol{d}_{\boldsymbol{q}})-\bar{{\bf
v}^{[f]}}(\boldsymbol{d}_{\boldsymbol{q}})\\\ {\bf
u}^{\boldsymbol{\pi}}(X_{i+j}^{(l)};\boldsymbol{d}_{\boldsymbol{q}})-\bar{{\bf
u}}(\boldsymbol{d}_{\boldsymbol{q}})\\\ \end{array}\right)^{\top},$ (F.2)

where $b_{n_{l}}$’s are the truncation points, $w_{n_{l}}(j)$’s are lag
window, and $\bar{\bf
v}^{[f]}(\boldsymbol{d}_{\boldsymbol{q}})=\sum_{i=1}^{n_{l}}{\bf
v}^{[f],\boldsymbol{\pi}}(X_{i}^{(l)};\boldsymbol{d}_{\boldsymbol{q}})/n_{l}$.
Let
$\widehat{\Lambda}(\boldsymbol{d}_{\boldsymbol{q}})\equiv\sum_{l=1}^{k}(a_{l}^{2}n/n_{l})\widehat{\Lambda}_{l}(\boldsymbol{d}_{\boldsymbol{q}})$.
Finally, let

$\boldsymbol{\hat{\rho}}(\boldsymbol{\hat{d}}_{\boldsymbol{q}})=\nabla{\bf
h}(\hat{\bf
v}^{[f]}(\boldsymbol{\pi};\boldsymbol{\hat{d}}_{\boldsymbol{q}}),\hat{\bf
u}(\boldsymbol{\pi};\boldsymbol{\hat{d}}_{\boldsymbol{q}}))\widehat{\Lambda}(\boldsymbol{\hat{d}}_{\boldsymbol{q}})\nabla{\bf
h}(\hat{\bf
v}^{[f]}(\boldsymbol{\pi};\boldsymbol{\hat{d}}_{\boldsymbol{q}}),\hat{\bf
u}(\boldsymbol{\pi};\boldsymbol{\hat{d}}_{\boldsymbol{q}}))^{\top},$

Theorem 4 Suppose there exists $r\in[0,\infty)$ such that $n/N\rightarrow r$
where $N\equiv\sum_{l=1}^{k}$ and $n=\sum_{l=1}^{k}n_{l}$ are the total sample
sizes for stages 1 and 2 respectively. In addition, let $n_{l}/n\rightarrow
s_{l}$ for $l=1,\cdots,k$.

  1. (a)

Assume that the stage 1 Markov chains are polynomially ergodic of order $m>1$.
Further, assume that the stage 2 Markov chains $\Phi_{1},\ldots,\Phi_{k}$ are
polynomially ergodic of order $m$, and for some $\delta>0$
$E_{q_{l}}|u^{\pi}(X;\boldsymbol{d}_{\boldsymbol{q}})|^{2+\delta}<\infty$ and
$E_{q_{l}}|v^{[f],\pi}(X;\boldsymbol{d}_{\boldsymbol{q}})|^{2+\delta}<\infty$
for each $\pi\in\Pi$ and $l=1,\cdots,k$ where $m>1+2/\delta$. Then as
$n_{1},\ldots,n_{k}\rightarrow\infty$,

$\sqrt{n}\left(\begin{array}[]{c}\boldsymbol{\hat{d}}_{\boldsymbol{q}}-{\boldsymbol{d}}_{\boldsymbol{q}}\\\
\hat{{\boldsymbol{\eta}}}^{[f]}(\boldsymbol{\pi};\boldsymbol{\hat{d}}_{\boldsymbol{q}})-{\bf
E}_{\boldsymbol{\pi}}f\\\ \end{array}\right)\buildrel
d\over{\rightarrow}N\Bigg{(}0,\left(\begin{array}[]{cc}rV_{\boldsymbol{q}}&\Delta_{12}\\\
\Delta_{21}&\Delta_{22}\\\ \end{array}\right)\Bigg{)},$ (F.3)

where
$\Sigma_{21}=rE(\boldsymbol{\pi};\boldsymbol{d}_{\boldsymbol{q}})V_{\boldsymbol{q}}$,
$\Sigma_{12}=\Sigma_{21}^{\top}$, and
$\Sigma_{22}=rE(\boldsymbol{\pi};\boldsymbol{d}_{\boldsymbol{q}})V_{\boldsymbol{q}}E(\boldsymbol{\pi};\boldsymbol{d}_{\boldsymbol{q}})^{\top}+{\boldsymbol{\rho}}(\boldsymbol{d}_{\boldsymbol{q}})$.

  2. (b)

Suppose that the conditions of Theorem 1 hold for the stage 1 Markov chains.
Let $\widehat{V}_{\boldsymbol{q}}$ be the consistent estimator of
$V_{\boldsymbol{q}}$ given in Theorem 1. Assume that the Markov chains
$\Phi_{1},\ldots,\Phi_{k}$ are polynomially ergodic of order
$m\geq(1+\epsilon)(1+2/\delta)$ for some $\epsilon,\delta>0$ such that
$E_{q_{l}}\|{\bf u}^{\boldsymbol{\pi}}(X;\boldsymbol{d})\|^{4+\delta}<\infty$
and $E_{q_{l}}\|{\bf
v}^{[f],\boldsymbol{\pi}}(X;\boldsymbol{d})\|^{4+\delta}<\infty$, ($\|\cdot\|$
denotes the Euclidean norm) for all $l=1,\ldots,k$, and $w_{n_{l}}$ and
$b_{n_{l}}$ satisfy conditions 1-4 in (Vats et al.,, 2018, Theorem 2). Then
$(n/N)\widehat{E}(\boldsymbol{\pi};\boldsymbol{\hat{d}}_{\boldsymbol{q}})\widehat{V}_{\boldsymbol{q}}\widehat{E}(\boldsymbol{\pi};\boldsymbol{\hat{d}}_{\boldsymbol{q}})^{\top}+\boldsymbol{\hat{\rho}}(\boldsymbol{\hat{d}}_{\boldsymbol{q}})$
is a strongly consistent estimator of $\Delta_{22}$ and
$(n/N)\widehat{E}(\boldsymbol{\pi};\boldsymbol{\hat{d}}_{\boldsymbol{q}})\widehat{V}_{\boldsymbol{q}}$
is a consistent estimator of $\Delta_{21}$.

Using similar arguments as in Section 3.3 of the paper, the joint entropy of
$\hat{{\boldsymbol{\eta}}}^{[f]}(\boldsymbol{\pi};\boldsymbol{\hat{d}}_{\boldsymbol{q}})$
and $\boldsymbol{\hat{d}}_{\boldsymbol{q}}$ is sum of the entropy of
$\boldsymbol{\hat{d}}_{\boldsymbol{q}}$, and the conditional entropy of
$\hat{{\boldsymbol{\eta}}}^{[f]}(\boldsymbol{\pi};\boldsymbol{\hat{d}}_{\boldsymbol{q}})$
given $\boldsymbol{\hat{d}}_{\boldsymbol{q}}$. Thus the maximum entropy
selection of skeleton points boils down to choosing $\boldsymbol{q}$ by
maximizing $\log\mbox{det}(\widehat{V_{\boldsymbol{q}}})$.

###### Proof of Theorem 4.

Since the Markov chains used in stage 1 are polynomially ergodic of order
$m>1$, from Roy et al., (2018, Theorem 1), we have
$N^{1/2}(\boldsymbol{\hat{d}}_{\boldsymbol{q}}-\boldsymbol{d}_{\boldsymbol{q}})\buildrel
d\over{\rightarrow}{\cal N}(0,V_{\boldsymbol{q}})$. Since $n/N\rightarrow r$,
it follows that
$\sqrt{n}(\boldsymbol{\hat{d}}_{\boldsymbol{q}}-\boldsymbol{d}_{\boldsymbol{q}})\buildrel
d\over{\rightarrow}{\cal N}(0,rV_{\boldsymbol{q}})$. Following Roy et al.,
(2018, Proof of Theorem 3) we write

$\sqrt{n}(\hat{{\boldsymbol{\eta}}}^{[f]}(\boldsymbol{\pi};\boldsymbol{\hat{d}}_{\boldsymbol{q}})-{\bf
E}_{\boldsymbol{\pi}}f)=\sqrt{n}(\hat{{\boldsymbol{\eta}}}^{[f]}(\boldsymbol{\pi};\boldsymbol{\hat{d}}_{\boldsymbol{q}})-\hat{{\boldsymbol{\eta}}}^{[f]}(\boldsymbol{\pi};\boldsymbol{d}_{\boldsymbol{q}}))+\sqrt{n}(\hat{{\boldsymbol{\eta}}}^{[f]}(\boldsymbol{\pi};\boldsymbol{d}_{\boldsymbol{q}})-{\bf
E}_{\boldsymbol{\pi}}f).$ (F.4)

The 2nd term involves randomness only from the 2nd stage Markov chains. Note
that

$\hat{\bf
v}(\boldsymbol{\pi};\boldsymbol{d}_{\boldsymbol{q}})\buildrel\text{a.s.}\over{\longrightarrow}\sum_{l=1}^{k}a_{l}E_{\pi_{\xi_{l}}}{\bf
v}^{[f],\boldsymbol{\pi}}(X;\boldsymbol{d}_{\boldsymbol{q}})={\bf
E}_{\boldsymbol{\pi}}f\odot{\bf u}(\boldsymbol{\pi},q_{1}).$

Since $\sum_{l=1}^{k}a_{l}E_{q_{l}}{\bf
u}^{\boldsymbol{\pi}}(X;\boldsymbol{d}_{\boldsymbol{q}})={\bf
u}(\boldsymbol{\pi},q_{1})$, we have

$\sqrt{n}\left(\begin{array}[]{c}\hat{\bf
v}(\boldsymbol{\pi};\boldsymbol{d}_{\boldsymbol{q}})-{\bf
E}_{\boldsymbol{\pi}}f\odot{\bf u}(\boldsymbol{\pi},q_{1})\\\ \hat{\bf
u}(\boldsymbol{\pi};\boldsymbol{d}_{\boldsymbol{q}})-{\bf
u}(\boldsymbol{\pi},q_{1})\\\
\end{array}\right)=\sum_{l=1}^{k}a_{l}\sqrt{\frac{n}{n_{l}}}\frac{1}{\sqrt{n_{l}}}\sum_{i=1}^{n_{l}}\left(\begin{array}[]{c}{\bf
v}^{[f],\boldsymbol{\pi}}(X_{i}^{(l)};\boldsymbol{d}_{\boldsymbol{q}})-E_{q_{l}}{\bf
v}^{[f],\boldsymbol{\pi}}(X;\boldsymbol{d}_{\boldsymbol{q}})\\\ {\bf
u}^{\boldsymbol{\pi}}(X_{i}^{(l)};\boldsymbol{d}_{\boldsymbol{q}})-E_{q_{l}}{\bf
u}^{\boldsymbol{\pi}}(X;\boldsymbol{d}_{\boldsymbol{q}})\\\
\end{array}\right).$ (F.5)

Since $\Phi_{l}$ is polynomially ergodic of order $m$ and
$E_{q_{l}}|u^{\pi}(X;\boldsymbol{d}_{\boldsymbol{q}})|^{2+\delta}$ and
$E_{q_{l}}|v^{[f],\pi}(X;\boldsymbol{d}_{\boldsymbol{q}})|^{2+\delta}<\infty$
are finite for each $\pi\in\Pi$ where $m>1+2/\delta$, it follows that

$\frac{1}{\sqrt{n_{l}}}\sum_{i=1}^{n_{l}}\left(\begin{array}[]{c}{\bf
v}^{[f],\boldsymbol{\pi}}(X_{i}^{(l)};\boldsymbol{d}_{\boldsymbol{q}})-E_{q_{l}}{\bf
v}^{[f],\boldsymbol{\pi}}(X;\boldsymbol{d}_{\boldsymbol{q}})\\\ {\bf
u}^{\boldsymbol{\pi}}(X_{i}^{(l)};\boldsymbol{d}_{\boldsymbol{q}})-E_{q_{l}}{\bf
u}^{\boldsymbol{\pi}}(X;\boldsymbol{d}_{\boldsymbol{q}})\\\
\end{array}\right)\buildrel
d\over{\rightarrow}N(0,\Lambda_{l}(\boldsymbol{d}_{\boldsymbol{q}}))$

where $\Lambda_{l}(\boldsymbol{d}_{\boldsymbol{q}})$ is defined in (F.1). As
$n_{l}/n\rightarrow s_{l}$ and the Markov chains $\Phi_{l}$’s are independent,
it follows that

$\sqrt{n}\left(\begin{array}[]{c}\hat{\bf
v}(\boldsymbol{\pi};\boldsymbol{d}_{\boldsymbol{q}})-{\bf
E}_{\boldsymbol{\pi}}f\odot{\bf u}(\boldsymbol{\pi},q_{1})\\\ \hat{\bf
u}(\boldsymbol{\pi},q_{1};\boldsymbol{d}_{\boldsymbol{q}})-{\bf
u}(\boldsymbol{\pi},q_{1})\\\ \end{array}\right)\buildrel
d\over{\rightarrow}N(0,\Lambda(\boldsymbol{d}_{\boldsymbol{q}}))).$

Then applying the delta method to the function ${\bf h}$ we have a CLT for the
estimator
$\hat{{\boldsymbol{\eta}}}^{[f]}(\boldsymbol{\pi};\boldsymbol{d}_{\boldsymbol{q}})$,
that is, we have
$\sqrt{n}(\hat{{\boldsymbol{\eta}}}^{[f]}(\boldsymbol{\pi};\boldsymbol{d}_{\boldsymbol{q}})-{\bf
E}_{\boldsymbol{\pi}}f)\buildrel
d\over{\rightarrow}N(0,\boldsymbol{\rho}(\boldsymbol{d}_{\boldsymbol{q}}))$.

Next by Taylor series expansion of
$L(\boldsymbol{d})=\hat{\eta}^{[f]}(\pi;\boldsymbol{d})$ about
$\boldsymbol{d}_{\boldsymbol{q}}$, we have

$\sqrt{n}(L(\boldsymbol{\hat{d}}_{\boldsymbol{q}})-L(\boldsymbol{d}_{\boldsymbol{q}}))=\sqrt{n}\nabla
L(\boldsymbol{d}_{\boldsymbol{q}})^{\top}(\boldsymbol{\hat{d}}_{\boldsymbol{q}}-\boldsymbol{d}_{\boldsymbol{q}})+\frac{\sqrt{n}}{2}(\boldsymbol{\hat{d}}_{\boldsymbol{q}}-\boldsymbol{d}_{\boldsymbol{q}})^{\top}\nabla^{2}L(\boldsymbol{d}^{*})(\boldsymbol{\hat{d}}_{\boldsymbol{q}}-\boldsymbol{d}_{\boldsymbol{q}}),$

where $\boldsymbol{d}^{*}$ is between $\boldsymbol{d}_{\boldsymbol{q}}$ and
$\boldsymbol{\hat{d}}_{\boldsymbol{q}}$. As in Roy et al., (2018), we can then
show that

$\sqrt{n}(\hat{\eta}^{[f]}(\pi;\boldsymbol{\hat{d}}_{\boldsymbol{q}})-\hat{\eta}^{[f]}(\pi;\boldsymbol{d}_{\boldsymbol{q}}))=\sqrt{r}e(\pi_{\xi};\boldsymbol{d}_{\boldsymbol{q}})\sqrt{N}(\boldsymbol{\hat{d}}_{\boldsymbol{q}}-\boldsymbol{d}_{\boldsymbol{q}})+o_{p}(1).$

Accumulating the terms for all $\pi\in\Pi$, we have

$\sqrt{n}(\hat{{\boldsymbol{\eta}}}^{[f]}(\boldsymbol{\pi};\boldsymbol{\hat{d}}_{\boldsymbol{q}})-\hat{{\boldsymbol{\eta}}}^{[f]}(\boldsymbol{\pi};\boldsymbol{d}_{\boldsymbol{q}}))=\sqrt{r}E(\boldsymbol{\pi};\boldsymbol{d}_{\boldsymbol{q}})\sqrt{N}(\boldsymbol{\hat{d}}_{\boldsymbol{q}}-\boldsymbol{d}_{\boldsymbol{q}})+{\bf
o}_{p}(1).$

Thus for constant vectors $t_{1}$ and $t_{2}$ of dimensions $k-1$ and $p$
respectively, we have

$\displaystyle
t_{1}^{\top}\sqrt{n}(\boldsymbol{\hat{d}}_{\boldsymbol{q}}-\boldsymbol{d}_{\boldsymbol{q}})+t_{2}^{\top}\sqrt{n}(\hat{{\boldsymbol{\eta}}}^{[f]}(\boldsymbol{\pi};\boldsymbol{\hat{d}}_{\boldsymbol{q}})-\hat{{\boldsymbol{\eta}}}^{[f]}(\boldsymbol{\pi};\boldsymbol{d}_{\boldsymbol{q}}))$
$\displaystyle\buildrel
d\over{\rightarrow}N(0,r(t_{1}^{\top}+t_{2}^{\top}E(\boldsymbol{\pi};\boldsymbol{d}_{\boldsymbol{q}}))V_{\boldsymbol{q}}(t_{1}+E(\boldsymbol{\pi};\boldsymbol{d}_{\boldsymbol{q}})^{\top}t_{2})+t_{2}^{\top}{\boldsymbol{\rho}}(\boldsymbol{d}_{\boldsymbol{q}})t_{2}),$
(F.6)

where the last step follows from the independence of the Markov chains
involved in the two stages. Note that the variance in (F) is the same as

$(t_{1}^{\top},t_{2}^{\top})\left(\begin{array}[]{cc}rV_{\boldsymbol{q}}&\Delta_{12}\\\
\Delta_{21}&\Delta_{22}\\\
\end{array}\right)(t_{1}^{\top},t_{2}^{\top})^{\top}.$

Hence the Cramér-Wold device implies the joint CLT in (F.3). Thus Theorem 4
(a) is proved.

From the proofs of Theorem 1 and Theorem 2 (b), we know that
$r\widehat{E}(\boldsymbol{\pi};\boldsymbol{\hat{d}}_{\boldsymbol{q}})\widehat{V}_{\boldsymbol{q}}\widehat{E}(\boldsymbol{\pi};\boldsymbol{\hat{d}}_{\boldsymbol{q}})^{\top}$
is a consistent estimator of
$rE(\boldsymbol{\pi};\boldsymbol{d}_{\boldsymbol{q}})V_{\boldsymbol{q}}E(\boldsymbol{\pi};\boldsymbol{d}_{\boldsymbol{q}})^{\top}$.
Also, $\nabla{\bf h}(\hat{\bf
v}^{[f]}(\boldsymbol{\pi};\boldsymbol{\hat{d}}_{\boldsymbol{q}}),\hat{\bf
u}(\boldsymbol{\pi};\boldsymbol{\hat{d}}_{\boldsymbol{q}}))\buildrel\text{a.s.}\over{\longrightarrow}\nabla{\bf
h}({\bf E}_{\boldsymbol{\pi}}f\odot{\bf u}(\boldsymbol{\pi},q_{1}),{\bf
u}(\boldsymbol{\pi},q_{1}))$. If $\boldsymbol{d}_{\boldsymbol{q}}$ is known
from the assumptions of Theorem 4 (b) and the results in Vats et al., (2018),
we know that $\Lambda_{l}(\boldsymbol{d}_{\boldsymbol{q}})$ is consistently
estimated by its SV estimator
$\widehat{\Lambda}_{l}(\boldsymbol{d}_{\boldsymbol{q}})$ defined in (F.1).
Then using similar arguments as in the proof of Theorem 3, we can show that
every element of
$\widehat{\Lambda}_{l}(\boldsymbol{d}_{\boldsymbol{q}})-\widehat{\Lambda}_{l}(\boldsymbol{\hat{d}}_{\boldsymbol{q}})$
converges to zero (a.e.). Thus
$\boldsymbol{\hat{\rho}}(\boldsymbol{\hat{d}}_{\boldsymbol{q}})$ is a
consistent estimator of $\boldsymbol{\rho}(\boldsymbol{d}_{\boldsymbol{q}})$.
Hence Theorem 5 (b) is proved. ∎

##  Appendix G Comparison of spectral variance and batch means estimators
using a negative binomial GLM example

In this section using a simulation example, we test the consistency of the SV
estimators given in Theorem 2 of the paper, and compare these estimators with
the BM estimators from Roy et al., (2018) and the empirical standard deviation
estimator obtained by replications.

We consider simulated independent count data $y_{i}$, $i=1,\ldots,51$ from a
negative binomial distribution with size parameter $\xi=3$ and mean
$\mu_{i}=\exp(\beta_{0}+\beta_{1}v_{i})$ with $\beta_{0}=0$, $\beta_{1}=1$,
and $v_{i}=-1+(i-1)/25$, $i=1,\ldots,51$. The parameters $\xi$, and
$\beta=(\beta_{0},\beta_{1})$ are then assumed unknown for making inference.
For $\beta_{0}$ and $\beta_{1}$ we use independent $t$ priors with location
parameter 0, scale 10, and degrees of freedom 1. We draw samples from the
posterior density of $\beta$,
$\pi_{\xi}(\beta|\boldsymbol{y})\equiv\ell_{\xi}(\beta|\boldsymbol{y})\pi(\beta)/\theta_{\xi}$
for $\xi\in{\boldsymbol{\xi}}$, with the skeleton set
${\boldsymbol{\xi}}=\\{\xi_{1}=1,\xi_{2}=5\\}$ where
$\ell_{\xi}(\beta|\boldsymbol{y})$ and $\pi(\beta)$ are the likelihood and the
prior density of $\beta$, respectively, and
$\theta_{\xi}\equiv\theta_{\xi}(\boldsymbol{y})\equiv\int_{\mathbb{R}^{2}}\ell_{\xi}(\beta|\boldsymbol{y})\pi(\beta)d\beta$.
Our goal is to estimate $u=\theta_{3}/\theta_{1}$, (the ratio of the
normalizing constants of the posterior densities corresponding to $\xi=3$ and
$\xi=\xi_{1}$) and $E_{\pi_{3}}\beta_{1}$.

At the first stage we estimate the ratio $d_{2}=\theta_{2}/\theta_{1}$ by
maximizing (2.3) of the paper using samples of size $N_{j}$ from
$\pi_{\xi_{j}}$, $j=1,2$. The same samples are used to estimate the asymptotic
standard error of $\hat{d}_{2}$ from Theorem 1 of the paper. At the second
stage we generate new samples of size $n_{j}$ from $\pi_{\xi_{j}}$, $j=1,2$ to
estimate $u$ by (1.3) and $E_{\pi_{3}}\beta_{1}$ by $\hat{\eta}^{[f]}$ defined
in Section 2 as well as their asymptotic variances given in Theorem 2 of the
paper.

We demonstrate the effect of increasing sample size on the asymptotic standard
error estimate. We vary $n_{1}$ according to
$n_{1}\in\\{5^{2},10^{2},15^{2},\ldots,40^{2}\\}$ and set
$N_{1}=N_{2}=9n_{1}=9n_{2}$. We use RStan (Stan Development Team,) for
sampling from $\pi_{\xi}$ after an initial burn-in of 300.

For a given $\boldsymbol{y}$, the estimation procedure is replicated 200 times
using new Markov chain samples each time. For each replication, we calculate
$\hat{u}$, $\hat{\eta}^{[f]}$ and their asymptotic variance estimates. We use
three different methods for estimating the asymptotic variance: the BM, SV
using the Tukey-Hanning window, and SV using the modified Bartlett window. For
the BM method, the stage one samples are split into $\sqrt{N_{j}}$ batches and
similarly the stage two samples are split into $\sqrt{n_{j}}$ batches. Using
the 200 $\hat{u}$ and $\hat{\eta}^{[f]}$ estimates corresponding to the Monte
Carlo sample sizes $n_{1}=40^{2}=n_{2},N_{1}=9n_{1}=N_{2}$, we compute their
empirical standard deviations (SDs), which we treat as the golden standards,
and compare them with 200 asymptotic SD estimates obtained by each of the
above three methods. We then compute the logarithmic score given by
$\log($ker(ESD)) where ESD is the empirical SD and ker$(\cdot)$ is the kernel
density estimate for the asymptotic SD derived from the 200 values. The
density estimate is obtained using the function ‘density’ in R (R Core Team,)
with default values.

The above procedure is repeated for 50 different $\boldsymbol{y}$’s and the
logarithmic score measures are averaged. Table 1 shows these average scores.
For both $u$ and $E_{\pi_{3}}\beta_{1}$ we observe that the score in the
asymptotic SD estimates is increased with increasing sample size, although the
SV method generally results in higher score than the BM methods. We have also
observed that the variability in the asymptotic SD estimates is reduced as the
sample size increases. Thus the accuracy of the asymptotic SD estimates is
improved as the sample size increases.

Table 1:  Average logarithmic score based on kernel density estimate of asymptotic SD’s evaluated at the empirical SD for the simulated negative binomial regression model for: $u$ (top rows); $E_{\pi_{3}}\beta_{1}$ (bottom rows). Higher values are better. The SD estimates are calculated for different Markov chain lengths as shown on the columns and for different methods: BM; SV using the Tukey-Hanning window (Tu); SV using the modified Bartlett window (Ba). $n_{1}$ | 25 | 100 | 225 | 400 | 625 | 900 | 1225 | 1600
---|---|---|---|---|---|---|---|---
BM |  |  $-$0.209
---
1.268
| 0.342
---
1.680
| 0.626
---
1.867
| 0.783
---
2.008
| 0.890
---
2.067
| 0.984
---
2.154
| 1.100
---
2.171
| 1.142
---
2.250
Tu |  |  $-$0.330
---
-1.155
| 0.344
---
1.714
| 0.679
---
1.926
| 0.825
---
2.069
| 0.985
---
2.123
| 1.041
---
2.199
| 1.159
---
2.253
| 1.192
---
2.302
Ba |  |  $-$0.423
---
-1.071
| 0.270
---
1.673
| 0.636
---
1.916
| 0.787
---
2.055
| 0.961
---
2.090
| 1.022
---
2.192
| 1.155
---
2.245
| 1.183
---
2.303

##  Appendix H Algorithms for computing the optimal skeleton set

In many cases, searching over the whole space $Q$ to find the optimal set of
proposal densities is computationally hard. Often, $Q$ comprises of a
parametric family of densities parameterized by $\xi\in\Xi$ so the problem
becomes choosing the skeleton set
${\boldsymbol{\xi}}=\\{\xi_{1},\xi_{2},\ldots,\xi_{k}\\}$ corresponding to the
parameter values of the proposal densities that minimizes an optimality
criterion $\phi({\boldsymbol{\xi}})$.

For convenience we work with a discretized version $\tilde{\Xi}$ of ${\Xi}$
and the skeleton set ${\boldsymbol{\xi}}$ is constrained to be
${\boldsymbol{\xi}}\subset\tilde{\Xi}$. Finding the optimal skeleton set in
this context has been studied in the sampling design and computer experiments
literature (Ko et al.,, 1995; Royle and Nychka,, 1998; Fang et al.,, 2006).
Here we present the two algorithms we used in this paper, a point-swapping
algorithm and a simulated annealing algorithm. Further details for these
algorithms can be found in Royle and Nychka, (1998) and Bélisle, (1992)
respectively.

###  H.1 Point-swapping algorithm

This algorithm performs many iterations, so it is mostly suited in cases where
the optimality criterion is fast to compute. We used the point-swapping
algorithm for computing the space-filling proposal distributions because in
these cases the criterion depends only on pairwise distances of the proposal
distributions which are fast to compute.

The algorithm proceeds as follows:

  1. Initialization: Initialize the skeleton set at ${\boldsymbol{\xi}}^{(0)}$ with $|{\boldsymbol{\xi}}^{(0)}|=k$.

  2. Iterations: Repeat for $i=1,2,\ldots$:

For $j=1,\ldots,k$: Swap the $j$th element of ${\boldsymbol{\xi}}^{(i-1)}$
with the element
$\xi_{j}^{\prime}\in\tilde{\Xi}\setminus{\boldsymbol{\xi}}^{(i-1)}$ such that
the new set
${\boldsymbol{\xi}}^{(i)}={\boldsymbol{\xi}}^{(i-1)}\setminus\\{\xi_{j}\\}\cup\\{\xi_{j}^{\prime}\\}$
produces the biggest drop in the value of the optimality criterion. If no such
$\xi_{j}^{\prime}$ exists, i.e., if
$\phi({\boldsymbol{\xi}}^{(i-1)})<\phi({\boldsymbol{\xi}}^{(i-1)}\setminus\\{\xi_{j}\\}\cup\\{\xi^{\prime}\\})$
for all $\xi^{\prime}\in\tilde{\Xi}\setminus{\boldsymbol{\xi}}^{(i-1)}$, then
${\boldsymbol{\xi}}^{(i)}={\boldsymbol{\xi}}^{(i-1)}$.

  3. Termination: Stop if, after looping over all elements in ${\boldsymbol{\xi}}^{(i)}$, the skeleton set remains unchanged. Return the final skeleton set.

We used the implementation in the R package fields (Nychka et al.,, 2017) to
compute the optimal skeleton set in our paper.

###  H.2 Simulated annealing algorithm

We use simulated annealing to compute the optimal set of proposal
distributions for MNX and ENT. The criterion for these methods is based on the
SV estimate of the asymptotic SE of the multiple IS estimator which is
computed from Monte Carlo samples. In this case, if samples from a particular
proposal distribution in the skeleton set already exist, then they are reused,
otherwise they are generated and stored for a possible future use.

The algorithm proceeds as follows:

  1. Initialization: Initialize the skeleton set at ${\boldsymbol{\xi}}^{(0)}$ with $|{\boldsymbol{\xi}}^{(0)}|=k$ and an initial temperature at $T_{0}$.

  2. Iterations: Repeat for $i=1,2,\ldots$:

    1. (a)

Set $T=T_{0}/\log(\lfloor(i-1)/B\rfloor B+\exp(1))$, where $B$ is a parameter
of the algorithm denoting the number of iterations before the temperature is
lowered.

    2. (b)

Randomly select $\xi\in{\boldsymbol{\xi}}^{(i-1)}$ and
$\xi^{\prime}\in\tilde{\Xi}\setminus{\boldsymbol{\xi}}^{(i-1)}$. Form the
candidate set
${\boldsymbol{\xi}}^{\prime}={\boldsymbol{\xi}}^{(i-1)}\setminus\\{\xi\\}\cup\\{\xi^{\prime}\\}$.

    3. (c)

With probability
$\min\\{1,\exp[(\phi({\boldsymbol{\xi}}^{(i-1)})-\phi({\boldsymbol{\xi}}^{\prime}))/T]\\}$
set ${\boldsymbol{\xi}}^{(i)}={\boldsymbol{\xi}}^{\prime}$, otherwise set
${\boldsymbol{\xi}}^{(i)}={\boldsymbol{\xi}}^{(i-1)}$.

  3. Termination: Stop if $i>i_{\mathrm{max}}$, a predetermined number of iterations. Return the skeleton set found among all iterations which corresponds to the lowest value of the optimality criterion.

##  Appendix I Finney,’s (1947) vasoconstriction data analysis using robit
model

Finney,’s (1947) vasoconstriction data consist of 39 binary responses denoting
the presence or absence, $y$, of vasoconstriction on the subject’s skin after
he or she inhaled air of volume $V$ at rate $R$. We consider a GLM where the
probability of presence for the $i$th subject, $\alpha_{i}$, is modeled using
a robit link function with df $\xi$,
$F_{\xi}^{-1}(\alpha_{i})=\beta_{0}+\beta_{1}\log V_{i}+\beta_{2}\log R_{i}$,
for $i=1,\ldots,39$. Here, $F_{\xi}(\cdot)$ denotes the distribution function
of the standard Student’s $t$ distribution with df $\xi$. As in Roy, (2014),
we consider a Bayesian analysis of the data with robit model. The prior for
$\beta$ is $\beta\sim t_{3}(0,10^{4}(W^{\top}W)^{-1},3)$, where $W$ is the
design matrix.

Roy, (2014) estimates the df parameter $\xi$ by maximizing the marginal
likelihood, that is,
$\hat{\xi}=\operatorname*{argmax}\theta_{\xi}(\boldsymbol{y})$. In particular,
Roy, (2014) uses the generalized IS estimator (1.3) to estimate the (ratios
of) marginal likelihoods, which in turn provides the estimate $\hat{\xi}$. Our
objective is to choose the importance sampling distributions from the family
of posterior densities $\Pi=\\{\pi_{\xi}(\beta|\boldsymbol{y}):\xi>0\\}$ for
the estimation of (1.3). We consider this in section I.1. Whereas in section
I.2, we analyze this problem with proposal densities from the multivariate
normal family.

###  I.1 Selection of proposals for generalized IS

In this section, we consider $Q=\Pi$, thus, choosing proposal distributions is
the same as choosing appropriate $\xi$ values. Because $\xi$ represents the df
parameter, we consider a wide range of points
$\tilde{\Xi}=\\{0.1,0.2,\ldots{},20\\}$ from where we choose the skeleton sets
for the different methods. The reference density corresponds to $\tilde{\xi}$,
where $\tilde{\xi}=10$ is at the middle of this range. For the multiple IS
methods we choose $k=5$ points, one of which must be $\tilde{\xi}$.

The computation is done in two phases. In the first phase we find the optimal
skeleton set, ${\boldsymbol{\xi}}$, for each method. In the second phase we
compare the relative SE of the naive and multiple IS estimators using the
skeleton sets computed in the first phase. The total number of samples used
for each method in the second phase is kept the same. The required Markov
chain samples were generated by Hamiltonian Monte Carlo implemented in the
stan language (Stan Development Team,, 2017). In all calculations involving
the asymptotic variance of the IS estimator, the SV estimate with the Tukey-
Hanning window,

$w_{n}(j)=0.5[1+\cos(\pi|j|/b_{n})]I(|j|<b_{n}),$

was used, where $b_{n}=\sqrt{n}$.

####  I.1.1 Phase 1: Finding the optimal set of proposal densities

  * •

NIS: This is not required because the proposal density is fixed at
${\boldsymbol{\xi}}_{\mathrm{nis}}=\\{10\\}$.

  * •

SFE: This method is based on the Euclidean distance between the parameters.
Therefore, ${\boldsymbol{\xi}}_{\mathrm{sfe}}=\\{2,6,10,14,18\\}$.

  * •

SFS: This method requires the SKLD between two densities corresponding to two
$\xi$ values. Figure 1 shows the logarithm of pairwise SKLD between densities
corresponding to two different values of $\xi$. The SKLD is computed using the
approximation of Section A. It can be seen that the distance is non-Euclidean.
For example, the densities corresponding to $\xi_{1}=1$ and $\xi_{2}=5$ are
further apart than the densities corresponding to $\xi_{1}=16$ and
$\xi_{2}=20$. Using the algorithm of Section H.1 we select
${\boldsymbol{\xi}}_{\mathrm{sfs}}=\\{0.3,1.1,1.9,3.3,10\\}$. It is noted that
the points concentrate more on the low values in $\tilde{\Xi}$.

  * •

SEQ: In this case the optimal set is computed by starting at
${\boldsymbol{\xi}}^{(1)}={\boldsymbol{\xi}}_{\mathrm{nis}}$. Then, given that
at the $i$th iteration, $i=2,\ldots,k$ we are at ${\boldsymbol{\xi}}^{(i-1)}$,
we obtain
${\boldsymbol{\xi}}^{(i)}={\boldsymbol{\xi}}^{(i-1)}\cup\\{\xi^{\prime}\\}$,
where $\xi^{\prime}$ corresponds to the point in
$\tilde{\Xi}\setminus{\boldsymbol{\xi}}^{(i-1)}$ with the highest relative SE.
The relative SE is again computed using 2,000 samples for stage 1 and 2,000
new samples for stage 2, after a burn in of 400. We find
${\boldsymbol{\xi}}_{\mathrm{seq}}=\\{0.1,0.2,0.3,0.7,10\\}$.

  * •

MNX: The optimal set is found by simulated annealing (Section H.2). We start
the simulated annealing algorithm at ${\boldsymbol{\xi}}_{\mathrm{sfs}}$, and
perform $i_{\mathrm{max}}=250$ iterations with $T_{0}=0.1$ and $B=10$. The
optimality criterion is computed as follows. Using 2,000 Markov chain samples
for stage 1 and 2,000 samples for stage 2, we compute
$\hat{\upsilon}_{1}^{2}(\xi):=\hat{c}(\pi;\hat{\boldsymbol{d}})^{\top}\hat{V}\hat{c}(\pi;\hat{\boldsymbol{d}})$
and $\hat{\upsilon}_{2}^{2}(\xi):=\hat{\tau}^{2}(\pi;\hat{\boldsymbol{d}})$
given in Theorem 2(a) for each $\xi\in\tilde{\Xi}$ as well as
$\hat{u}(\xi)=\hat{u}$ given in equation (1.3) of the main paper. Then, given
stage 1 sample size of $N$ and stage 2 sample size of $n$, the relative SE
estimate is given by

$\mathrm{RelSE}(\xi,N,n):=\frac{\hat{\upsilon}_{1}(\xi)/\sqrt{N}+\hat{\upsilon}_{2}(\xi)/\sqrt{n}}{\hat{u}(\xi)}.$
(I.1)

Assuming that the total sample size $M=N+n$ is fixed, the objective becomes
choosing ${\boldsymbol{\xi}}_{\mathrm{mnx}}$ in order to minimize

$\min_{N\in(0,M)}\max_{\xi\in\tilde{\Xi}}\mathrm{RelSE}(\xi,N,M-N).$ (I.2)

We find ${\boldsymbol{\xi}}_{\mathrm{mnx}}=\\{0.1,0.4,1.6,3.3,10\\}$. One can
also impose a constraint that $N$ is at least some number and at most some
other number while finding ${\boldsymbol{\xi}}_{\mathrm{mnx}}$ minimizing
(I.2).

  * •

ENT: The optimal set is found by simulated annealing (Section H.2). We start
the simulated annealing algorithm at ${\boldsymbol{\xi}}_{\mathrm{sfs}}$, and
perform $i_{\mathrm{max}}=250$ iterations with $T_{0}=1$ and $B=10$. The
optimality criterion used for this method is $-\log\mbox{det}(U)$ where $U$ is
the matrix with $(i,j)$th element
$U_{ij}=(\widehat{V}_{\boldsymbol{q}})_{ij}/(\hat{d}_{i}\hat{d}_{j})$. The
estimates $\hat{\boldsymbol{d}}$ and $\widehat{V}_{\boldsymbol{q}}$ are
computed from 2,000 Markov chain samples. In this example, the ENT skeleton
set turns out to be the same as the MNX set, thus
${\boldsymbol{\xi}}_{\mathrm{ent}}={\boldsymbol{\xi}}_{\mathrm{mnx}}=\\{0.1,0.4,1.6,3.3,10\\}$.

Figure 1:  Contour plot of the logarithm of the SKLD between the densities
corresponding to $\xi_{1}$ and $\xi_{2}$ for the vasoconstriction example.

####  I.1.2 Phase 2: Estimation of the ratio of marginal densities

After the skeleton sets are found, we generate a total of $M=50,000$ samples
from the proposal densities, equally divided among all densities in the set.
The total sample size $M$ is generally determined by available computational
resources. Since we have derived consistent SV estimators of the asymptotic
variance of $\hat{u}$, the sample size $M$ can be chosen such that the overall
SE of $\hat{u}$ is smaller than a pre-determined level of accuracy. For NIS we
simply take $M$ Gibbs samples from the density corresponding to
${\boldsymbol{\xi}}_{\mathrm{nis}}$, and for the multiple IS methods we take
$M/k$ samples from each density in the corresponding ${\boldsymbol{\xi}}$ set.
However, the total of $M/k$ samples must be split into stage 1, for estimating
the ratio of marginals within ${\boldsymbol{\xi}}$, $\boldsymbol{d}$, and
stage 2, for estimating the ratio of marginals over the whole set
$\tilde{\Xi}$. To determine the optimal split we use equation (I.1) where
$\hat{\upsilon}_{1}(\xi)$ and $\hat{\upsilon}_{2}(\xi)$ are calculated from
2,000 stage 1 and 2,000 stage 2 Markov chain samples. For MNX, SEQ, and ENT
the existing samples from Phase 1 are reused but for the space-filling
methods, new samples are generated. We take an equal number of samples from
each density, thus we take $\hat{N}/k$ stage 1 samples from each density in
the skeleton set where $\hat{N}$ is the integer that minimizes
$\max_{\xi\in\tilde{\Xi}}\mathrm{RelSE}(\xi,N,M-N)$, and $\hat{n}/k$ stage 2
samples, where $\hat{n}=M-\hat{N}$. This corresponds to stage 1 sample sizes
of 1500, 6500, 8500, 9000, 9000 and stage 2 sample sizes of 8500, 3500, 1500,
1000, 1000 from each density for SFE, SFS, SEQ, MNX, ENT, respectively.

The estimates of the relative SE of (1.3) and the value of the logarithm of
(1.3) across all $\xi$ values in $\tilde{\Xi}$ corresponding to the different
skeleton sets, ${\boldsymbol{\xi}}$, chosen are plotted in Figure 2. It can be
seen that the relative SE is larger when $\xi$ is small. NIS has the lowest
relative SE at $\tilde{\xi}$ but results in much higher standard errors at low
values of $\xi$. Indeed, at $\xi=0.5$, the SV estimate of the relative SE for
NIS method is about four times larger than that for the MNX and ENT methods.
It can also be seen that SFE does not have as good performance as SFS, again
because it avoids sampling from densities corresponding to low values of $\xi$
which are the ones that produce the largest relative SE. SEQ also does not
have good performance because it concentrates all points in a very narrow
region while, evidently, it is more beneficial to spread the points to cover a
wider area as MNX and ENT do.

Figure 2: Vasoconstriction example: The left and right panels show the values
of $\log\hat{u}(\boldsymbol{\hat{d}}_{\boldsymbol{\xi}})$ and estimates of the
relative SE of $\hat{u}(\boldsymbol{\hat{d}}_{\boldsymbol{\xi}})$ across $\xi$
values for different skeleton sets, respectively. Points included in the
skeleton sets are indicated by $+$.

###  I.2 Multiple IS using a mixture of multivariate normal proposals

The proposed methods of choosing reference distributions are applicable to IS
estimators in the situations where $Q$ can be different from $\Pi$. To
demonstrate this, in this section we analyze the vasoconstriction data using
the model and method discussed in Section I.1 with the difference that the
proposal densities are now chosen from the multivariate normal family instead
of the family of the posterior density of $\beta$. Thus we are now able to
draw iid samples from the importance sampling distributions which is not the
case when these are posterior densities of $\beta$ which require MCMC
sampling. Since $\boldsymbol{d}$ is known, the reverse logistic estimation is
not needed here. The space filling, minimax and sequential approaches
developed in the main paper can be used for selecting multivariate normal
proposals as we describe below.

If $n_{l}$ iid samples, $X_{i}^{(l)}$, $i=1,\ldots,n_{l}$, are drawn from the
density $q_{l}(x)$ (a normal density described later), $l=1,\ldots,k$, then
the normalizing constant, $\theta_{\xi}$, of the posterior density
$\pi_{\xi}(x)\equiv\nu_{\xi}(x)/\theta_{\xi}$, is estimated by

$\hat{\theta}_{\xi}=\frac{1}{|n|}\sum_{l=1}^{k}\sum_{i=1}^{n_{l}}\frac{\nu_{\xi}(X_{i}^{(l)})}{\bar{q}(X_{i}^{(l)})},$
(I.3)

where $n=n_{1}+\ldots+n_{k}$ and
$\bar{q}(x)=(n_{1}/n)q_{1}(x)+\ldots+(n_{k}/n)q_{k}(x)$. The variance of this
estimator is estimated by

$\widehat{\mbox{Var}}(\hat{\theta}_{\xi})=\frac{1}{n}\sum_{l=1}^{k}\sum_{i=1}^{n_{l}}\left(\frac{\nu_{\xi}(X_{i}^{(l)})}{\bar{q}(X_{i}^{(l)})}-\hat{\theta}_{\xi}\right)^{2}.$
(I.4)

To choose the proposal densities corresponding to a skeleton set
${\boldsymbol{\xi}}$, for $\xi\in{\boldsymbol{\xi}}$, let
$\tilde{\beta}_{\xi}$ denote the maximizer of
$\ell_{\xi}(\beta|\boldsymbol{y})\pi(\beta)$ and let $\tilde{H}_{\xi}$ denote
the Hessian matrix of $-\log\\{\ell_{\xi}(\beta|\boldsymbol{y})\pi(\beta)\\}$
evaluated at $\tilde{\beta}_{\xi}$. Then the normal approximation to
$\pi_{\xi}(\beta|\boldsymbol{y})$ is taken to be the multivariate normal with
mean $\tilde{\beta}_{\xi}$ and variance $\tilde{H}^{-1}_{\xi}$. Thus, $q_{l}$
is the normal approximation to the posterior density
$\pi_{\xi}(\beta|\boldsymbol{y})$, where $\xi=\xi_{l}$ is a skeleton point.

Unlike in Section I.1, where we used Hamiltonian Monte Carlo to obtain
approximate samples from the proposal (posterior) densities, here we draw iid
samples from the proposal distributions. We use (I.3) to estimate
$\theta_{\xi}$ for all $\xi$ values in the range $\tilde{\Xi}$ identified
previously. We set $k=5$ and generate $n_{l}=10000$ samples from each of the
proposal densities $q_{l},l=1,\dots,5$. We also use naive importance sampling
by drawing $5\times 10000$ samples from the normal approximation to the
posterior for $\beta$ corresponding to
${\boldsymbol{\xi}}_{\mathrm{nis}}=\\{\tilde{\xi}=10\\}$.

For the two space filling methods, the sets
${\boldsymbol{\xi}}_{\mathrm{sfe}}$ and ${\boldsymbol{\xi}}_{\mathrm{sfs}}$
are obtained as described in Section I.1, thus the skeleton sets are the same
as in that section: ${\boldsymbol{\xi}}_{\mathrm{sfe}}=\\{2,6,10,14,18\\}$ and
${\boldsymbol{\xi}}_{\mathrm{sfs}}=\\{0.3,1.1,1.9,3.3,10\\}$.

For SEQ, we aim to select the set ${\boldsymbol{\xi}}_{\mathrm{seq}}$
sequentially, starting with ${\boldsymbol{\xi}}^{(1)}=\\{\tilde{\xi}\\}$. At
the $j$th iteration, we obtain the set
${\boldsymbol{\xi}}^{(j)}={\boldsymbol{\xi}}^{(j-1)}\cup\\{\xi_{j}\\}$. Then,
we draw samples from the normal approximation to the posterior density
corresponding to this $\xi_{j}$, as described in the previous paragraph. Using
these samples, together with existing samples from
${\boldsymbol{\xi}}^{(j-1)}$, we compute, using (I.3) and (I.4), the relative
SE of $\hat{\theta}_{\xi}$ for $\xi\in\tilde{\Xi}$. The value of $\xi$
corresponding to the highest relative SE, denoted by $\xi_{j+1}$, is added to
the set. The final set obtained using this method is
${\boldsymbol{\xi}}_{\mathrm{seq}}=\\{0.6,0.7,0.8,0.9,10\\}$.

For MNX, we run the simulated annealing algorithm with $T_{0}=10$, $B=10$,
$i_{\mathrm{max}}=250$, starting from ${\boldsymbol{\xi}}_{\mathrm{sfs}}$. For
every new density added to the skeleton set we generate 4000 samples from the
corresponding normal approximation, if not already available. These samples,
along with the existing samples from the other densities previously added are
used to calculate the relative SE, using (I.3) and (I.4), for all
$\xi\in\tilde{\Xi}$. The objective of the simulated annealing is to minimize
the maximum relative SE across $\tilde{\Xi}$. The optimal skeleton set is
${\boldsymbol{\xi}}_{\mathrm{mnx}}=\\{0.2,0.3,0.7,2.6,10\\}$.

Plots of the estimator (I.3) using the different methods are shown in the left
panel of Figure 3. It can be seen that NIS and SFE are significantly different
from the other methods. The plots of the relative SE estimates using (I.4) are
shown in the right panel of Figure 3. It can be seen that the NIS, SFE, and
SEQ methods lead to significantly higher relative SE compared to the proposed
methods for low values of $\xi$. For example, at $\xi=0.5$, the relative SE
for NIS, SFE, and SEQ are 2.5, 7, and 3.8 times higher respectively than MNX.

Figure 3: Logarithm of the IS estimates of the marginal density for Finney’s
vasoconstriction data using normal proposal densities (left); and estimates of
the relative SE of the marginal density estimates (right). Points included in
the skeleton sets are indicated by $+$.

##  Appendix J Further details on the autologistic example used in the main
paper

###  J.1 Derivation of the autologistic pmf

We start with the conditional pmf of $x(s_{i})|\boldsymbol{x}_{-s_{i}}$ given
by

$\pi_{i}(x(s_{i})|\boldsymbol{x}_{-s_{i}})=p_{i}^{x(s_{i})}(1-p_{i})^{1-x(s_{i})},$

where, as in the main paper,

$p_{i}=\frac{\exp\\{\mbox{logit}(\kappa)+(\gamma/w)\sum_{s_{j}\in\mbox{nb}_{i}}(x(s_{j})-\kappa)\\}}{1+\exp\\{\mbox{logit}(\kappa)+(\gamma/w)\sum_{s_{j}\in\mbox{nb}_{i}}(x(s_{j})-\kappa)\\}}.$

Letting
$\pi(\boldsymbol{x}|\gamma,\kappa)\propto\exp(A(\boldsymbol{x}|\theta))$,
since only pairwise dependencies are considered, it is known that (Kaiser and
Cressie,, 2000)

$A(\boldsymbol{x}|\theta)=\sum_{1\leq i\leq
m}\log\frac{\pi_{i}(x(s_{i})|\boldsymbol{x}^{*}_{-s_{i}})}{\pi_{i}(x^{*}(s_{i})|\boldsymbol{x}^{*}_{-s_{i}})}+\sum_{1\leq
i<j\leq
m;s_{j}\in\mbox{nb}_{i}}\log\frac{\pi_{i}(x(s_{i})|x(s_{j}),\boldsymbol{x}^{*}_{-\\{s_{i},s_{j}\\}})\pi_{i}(x^{*}(s_{i})|\boldsymbol{x}^{*}_{-s_{i}})}{\pi_{i}(x^{*}(s_{i})|x(s_{j}),\boldsymbol{x}^{*}_{-\\{s_{i},s_{j}\\}})\pi_{i}(x(s_{i})|\boldsymbol{x}^{*}_{-s_{i}})},$

for a suitably chosen $\boldsymbol{x}^{*}$. Choosing $\boldsymbol{x}^{*}={\bf
0}$, simple calculations show that

$A(\boldsymbol{x}|\theta)=(\mbox{logit}(\kappa)-\gamma\kappa)\sum\nolimits_{i=1}^{m}x(s_{i})+\frac{\gamma}{2w}\sum\nolimits_{i=1}^{m}\sum\nolimits_{s_{j}\in\mbox{nb}_{i}}x(s_{i})x(s_{j}).$

###  J.2 Computational details

Here we give more details on how we produced the results at the end of Section
3 in the main paper. The parameter set where we search over is

$\gamma\in\\{-4,-3.92,\ldots,4\\},$

in the case where $\kappa$ is set fixed at $\kappa=0.5$, and

$(\gamma,\kappa)\in\\{-4,-3.2,\ldots,4\\}\times\\{0.1,0.2,\ldots,0.9\\},$

in the case where both parameters are assumed unknown. In both cases we denote
the parameters by $\xi$ and the parameter set by $\tilde{\Xi}$.

The computation is done in two phases. In the first phase we find the optimal
skeleton set, ${\boldsymbol{\xi}}$, for each method. In the second phase we
compare the relative SE of the naive and multiple IS estimators using the
skeleton sets computed in the first phase. The total number of samples used
for each method in the second phase is kept the same. The required Markov
chain samples were generated by Gibbs sampling from the conditional
distribution of each component of $\boldsymbol{x}$ given its neighbors. In all
calculations involving the asymptotic variance of the IS estimator, the SV
estimate with the Tukey-Hanning window,

$w_{n}(j)=0.5[1+\cos(\pi|j|/b_{n})]I(|j|<b_{n}),$

was used, where $b_{n}=\sqrt{n}$.

####  J.2.1 Phase 1: Finding the optimal set of proposal densities