Multivariate initial sequence estimators in Markov chain Monte Carlo

TL;DR

This paper introduces a new multivariate initial sequence estimator for Markov chain Monte Carlo that accurately estimates the covariance matrix of sample means, improving stability and validity over existing methods.

Contribution

It develops a novel multivariate extension of Geyer's Monte Carlo error estimation method, providing asymptotically valid and stable covariance matrix estimates in MCMC.

Findings

The proposed estimator is asymptotically valid.

Simulation experiments show improved stability.

Finite sample properties are favorable.

Abstract

Markov chain Monte Carlo (MCMC) is a simulation method commonly used for estimating expectations with respect to a given distribution. We consider estimating the covariance matrix of the asymptotic multivariate normal distribution of a vector of sample means. Geyer (1992) developed a Monte Carlo error estimation method for estimating a univariate mean. We propose a novel multivariate version of Geyer's method that provides an asymptotically valid estimator for the covariance matrix and results in stable Monte Carlo estimates. The finite sample properties of the proposed method are investigated via simulation experiments.

Tables

Table 1. Table 1 : Estimated ESS with standard errors. For uIS, only the minimum estimated ESS is reported.

	mK	mIS	mISadj	uIS
Ex1($\times {10}^4$)	5.40 (.002)	5.22 (.001)	5.18 (.001)	3.95 (.002)
Ex2($\times {10}^4$)	4.74 (.007)	3.76 (.002)	3.52 (.003)	1.30 (.001)
Ex3($\times {10}^5$)	8.78 (.000)	8.39 (.000)	8.30 (.001)	7.58 (.001)

Table 2. Table 2 : Average volumes to the p 𝑝 p th ( p = 5 , 8 , 12 𝑝 5 8 12 p=5,8,12 for Ex1, 2, 3) root and standard errors of nominal 90 % percent 90 90\% confidence regions.

	uIS	mK	mIS	mISadj	uIS-Bonferroni
Ex1($\times {10}^{-3}$)	5.53 (.001)	6.31 (.001)	6.41 (.001)	6.44 (.001)	7.82 (.001)
Ex2($\times {10}^{-2}$)	3.51 (.002)	3.95 (.003)	4.43 (.003)	4.58 (.003)	5.33 (.004)
Ex3($\times {10}^{-3}$)	3.84 (.000)	4.78 (.000)	4.89 (.000)	4.92 (.000)	6.16 (.000)

Table 3. Table 3 : Estimated coverage probabilities and standard errors of nominal 90 % percent 90 90\% confidence regions.

	uIS	mK	mIS	mISadj	uIS-Bonferroni
Ex1	.622 (.0108)	.885 (.0071)	.898 (.0068)	.900 (.0067)	.908 (.0065)
Ex2	.386 (.0109)	.660 (.0106)	.845 (.0081)	.881 (.0073)	.862 (.0077)
Ex3	.323 (.0105)	.882 (.0072)	.911 (.0064)	.916 (.0062)	.917 (.0062)

Table 4. Table 4 : Estimated asymptotic generalized variances ( × 10 79 absent superscript 10 79 \times 10^{79} ) of the Monte Carlo errors using the colon cancer dataset.

mK	mIS	mISadj
.044	6.285	77.144

Table 5. Table 5 : Estimated ESS ( × 10 5 absent superscript 10 5 \times 10^{5} ) with Monte Carlo sample size 4 × 10 6 4 superscript 10 6 4\times 10^{6} using the colon cancer dataset. For uIS, only the minimum estimated ESS is reported.

mK	mIS	mISadj	uIS
4.637	4.296	4.134	1.137

Table 6. Table 6 : Volumes to the p 𝑝 p th root ( × 10 − 3 absent superscript 10 3 \times 10^{-3} ) of 90 % percent 90 90\% confidence regions for all components ( p = 65 𝑝 65 p=65 ) and for μ 𝜇 \mu and λ 1 subscript 𝜆 1 \lambda_{1} only ( p = 2 𝑝 2 p=2 ). Computed with Monte Carlo sample size 4 × 10 6 4 superscript 10 6 4\times 10^{6} using the colon cancer dataset.

uIS	mK	mIS	mISadj	uIS-Bonferroni
6.96	8.70	9.04	9.22	13.39

Table 7. (a) 65-dimensional confidence regions.

uIS	mK	mIS	mISadj	uIS-Bonferroni
6.96	8.70	9.04	9.22	13.39

Table 8. (b) Bivariate confidence regions for μ 𝜇 \mu and λ 1 subscript 𝜆 1 \lambda_{1} .

uIS	mK	mIS	mISadj	uIS-Bonferroni
5.60	6.04	6.44	6.56	6.68

Equations

$$\mu ={\mathrm{E}}_F\left\{g(X)\right\}={\int }_{\mathcal{X}}g(x)F(dx).$$

$${\mu }_n=\frac{1}{n}\sum _{i=1}^{n}g(X_i)\to \mu \text{ as }n\to \infty .$$

$$\sqrt{n}({\mu }_n-\mu )⇝{\mathcal{N}}_p(0,\Sigma ).$$

$$\text{ESS}=n{\left(|\Lambda |/|\Sigma |\right)}^{1/p},$$

$$F(dx)P(x,dy)=F(dy)P(y,dx).$$

$$\Sigma =\sum _{t=-\infty }^{+\infty }{\gamma }_t=-{\gamma }_0+\sum _{t=0}^{\infty }({\gamma }_t+{\gamma }_{-t})=-{\gamma }_0+2\sum _{t=0}^{\infty }{\gamma }_t=-{\gamma }_0+2\sum _{i=0}^{\infty }{\Gamma }_i.$$

$${\Sigma }_m=-{\gamma }_0+\sum _{t=0}^{2m+1}({\gamma }_t+{\gamma }_{-t})=-{\gamma }_0+2\sum _{i=0}^m{\Gamma }_i.$$

$${\gamma }_{n,t}={\gamma }_{n,-t}^{\top }=\frac{1}{n}\sum _{i=1}^{n-t}\left\{g(X_i)-{\mu }_n\right\}{\left\{g(X_{i+t})-{\mu }_n\right\}}^{\top }$$

$${\Sigma }_{n,m}=-{\gamma }_{n,0}+\sum _{t=0}^{2m+1}({\gamma }_{n,t}+{\gamma }_{n,-t})=-{\gamma }_{n,0}+2\sum _{i=0}^m{\Gamma }_{n,i}.$$

$${\tilde{\Sigma }}_{n,t_n}={\Sigma }_{n,s_n}+2\sum _{i=s_n+1}^{t_n}{\Gamma }_{n,i}^{+}$$

$${\sigma }_{\mathrm{pos},n}^2=-{\gamma }_{n,0}+2\sum _{i=0}^{t_n}{\Gamma }_{n,i},$$

$$\underset{n\to \infty }{\mathrm{liminf}}{\sigma }_{\mathrm{pos},n}^2\ge {\sigma }^2\text{with probability}1.$$

$$Y_i|X_i,\beta \stackrel{ind}{\sim }\mathrm{Bernoulli}\left\{\frac{1}{1+\exp (-X_i\beta )}\right\}\text{and}\beta \sim {\mathcal{N}}_5(0,4I_5).$$

$$Y_i\mid {\theta }_i,{\gamma }_i\stackrel{ind.}{\sim }\mathcal{N}({\theta }_i,{\gamma }_i^{-1}),{\theta }_i\mid \mu ,{\lambda }_{\theta },{\lambda }_i\stackrel{ind.}{\sim }\mathcal{N}(\mu ,{\lambda }_{\theta }^{-1}{\lambda }_i^{-1}),$$

$$\mu \sim \mathcal{N}(m_0,v_0^{-1}),{\gamma }_i\stackrel{iid}{\sim }\mathcal{G}(a_3,b_3),{\lambda }_{\theta }\sim \mathcal{G}(a_1,b_1),{\lambda }_i\stackrel{iid}{\sim }\mathcal{G}(a_2,b_2),$$

$${\lambda }_{\min }(A)={\lambda }^{(1)}(A)\le {\lambda }^{(2)}(A)\le \cdots \le {\lambda }^{(p-1)}(A)\le {\lambda }^{(p)}(A)={\lambda }_{\max }(A).$$

$${\lambda }^{(k)}(A)+{\lambda }^{(1)}(B)\le {\lambda }^{(k)}(A+B).$$

$${\lambda }^{(k)}(A)>{\lambda }^{(k)}(B)\ge 0.$$

$${\gamma }_t={\int }_{-1}^1{\lambda }^t{\mathcal{E}}_g(d\lambda ).$$

$${\Gamma }_i={\gamma }_{2i}+{\gamma }_{2i+1}={\int }_{-1}^1{\lambda }^{2i}(1+\lambda ){\mathcal{E}}_g(d\lambda )$$

$${\Gamma }_i-{\Gamma }_{i+1}={\int }_{-1}^1{\lambda }^{2i}(1+\lambda {)}^2(1-\lambda ){\mathcal{E}}_g(d\lambda ).$$

$${\Gamma }_0={\int }_{-1}^1(1+\lambda ){\mathcal{E}}_g(d\lambda )=0\iff {\mathcal{E}}_g(\{-1\})=1,$$

$${\Gamma }_0-{\Gamma }_1={\int }_{-1}^1(1+\lambda {)}^2(1-\lambda ){\mathcal{E}}_g(d\lambda )=0\iff {\mathcal{E}}_g(\{-1,1\})=1.$$

$${\Gamma }_i={\int }_{-1}^1{\lambda }^{2i}(1+\lambda ){\mathcal{E}}_g(d\lambda )=0\iff {\mathcal{E}}_g(\{-1,0\})=1,$$

$${\Gamma }_i-{\Gamma }_{i+1}={\int }_{-1}^1{\lambda }^{2i}(1+\lambda {)}^2(1-\lambda ){\mathcal{E}}_g(d\lambda )=0\iff {\mathcal{E}}_g(\{-1,0,1\})=1.$$

$${\gamma }_t={\int }_{-1}^1{\lambda }^t{\mathcal{E}}_g(d\lambda )=0$$

$${\mathcal{E}}_g(\{0\})<1.$$

$${\mathcal{E}}_g(\{-1,1\})=0.$$

$${\gamma }_t^{\ast }={\gamma }_{-t}^{\ast }={\mathrm{cov}}_F\{h(X_i),h(X_{i+t})\}$$

$${\Gamma }_i^{\ast }$$

Full text

††journal: Journal of Multivariate Analysis

1 Introduction

Many distributions encountered in modern applications are intractable in the sense that it is difficult to calculate expectations without resorting to simulation-based methods. If it is difficult to simulate independent realizations from the target distribution, then it is natural to turn to Markov chain Monte Carlo (MCMC). An MCMC experiment consists of generating a realization of an irreducible Markov chain having the distribution of interest as its stationary distribution [22, 25]. The simulated data may then be used to estimate a vector of means associated with the stationary distribution. The reliability of this estimation can be assessed by forming asymptotically valid confidence regions for the means of the stationary distribution [6, 7, 9, 18, 19, 28]. (There is a simliar approach to quantile estimation [3].) The confidence regions are based on estimating the covariance matrix in a multivariate Markov chain central limit theorem (CLT). We propose and study a novel method for estimating this covariance matrix.

Estimating the covariance matrix has been mostly ignored in the MCMC literature until recently. Vats et al. [28] and Vats et al. [29] studied non-overlapping batch means and spectral methods, respectively, and found that these estimators often underestimate the size of the confidence regions and overestimate the effective sample size unless the Monte Carlo sample sizes are enormous. Kosorok [21] proposed an estimator that is closer in spirit to ours than the spectral and batch means methods, but we will see later that it typically overestimates the effective sample size, resulting in overconfidence in the reliability of the simulation. We propose alternative estimators of the covariance matrix that require weaker mixing conditions on the Markov chain and weaker moment conditions on the function of interest than those required by batch means and spectral methods. Specifically, our method applies as long as a Markov chain CLT holds and detailed balance is satisfied, which is not enough to guarantee the asymptotic validity of batch means or spectral methods. We show that the proposed estimators are asymptotically valid and study their empirical performance. The problem we consider will now be described more formally.

Let $F$ be a distribution having support $\mathcal{X}$ and if $p\geq 1$, let $g:\mathcal{X}\to{}^{p}$ be $F$-integrable and set

[TABLE]

Also let $\Phi=\{X_{0},X_{1},X_{2},\ldots\}$ be a Harris ergodic—namely, irreducible, aperiodic and Harris recurrent—Markov chain having invariant distribution $F$. By averaging the function over a realization of $\Phi$, estimation of $\mu$ is straightforward since, with probability 1,

[TABLE]

The Markov chain strong law justifies the use of MCMC but provides no information about the quality of estimation or how large the simulation size $n$ should be. More specifically, additional information is needed to answer either of the following two questions.

Given a pre-specified run length $n$, how reliable is $\mu_{n}$ as an estimate of $\mu$? Specifically, how do we construct a confidence region for $\mu$? 2. 2.

How large should the simulation size $n$ be to ensure a reliable estimate of $\mu$?

We can address these issues through the approximate sampling distribution of the Monte Carlo error, $\mu_{n}-\mu$. A Markov chain CLT exists when there is a positive definite matrix $\Sigma$ such that, as $n\to\infty$,

[TABLE]

See Jones [17] and Roberts and Rosenthal [26] for conditions which ensure a CLT. Notice that, due to the serial correlation inherent to the Markov chain, $\Sigma\neq\mathrm{var}_{F}\{g(X)\}$ except in trivial cases. In Section 3 we propose two new estimators of $\Sigma$. For now, let $\Sigma_{n}$ be a generic positive definite estimator of $\Sigma$.

A confidence region for $\mu$ constructed using $\Sigma_{n}$ forms an ellipsoid in $p$ dimensions oriented along the directions of the eigenvectors of $\Sigma_{n}$. Let $|\cdot|$ denote determinant. One can verify by straightforward calculation that the volume of the confidence region is proportional to $\sqrt{|\Sigma_{n}|}$ and thus depends on the estimated covariance matrix $\Sigma_{n}$ only through the estimate $|\Sigma_{n}|$ of the generalized variance of the Monte Carlo error, $|\Sigma|$. The volume of the confidence region can describe whether the simulation effort is sufficiently large to achieve the desired level of precision in estimation [19, 7, 28].

Another common and intuitively reasonable method for choosing the simulation effort is to simulate until a desired effective sample size (ESS), i.e., the number with the property that $\mu_{n}$ has the same precision as the sample mean obtained by that number of independent and identically distributed (iid) samples, has been achieved [1, 5, 10]. Let $\Lambda=\mathrm{var}_{F}\{g(X)\}$. Vats et al. [28] introduced the following definition of effective sample size

[TABLE]

which is naturally estimated with $n(|\Lambda_{n}|/|\Sigma_{n}|)^{1/p}$ where $\Lambda_{n}$ is an estimator of $\Lambda$, e.g., the usual sample covariance matrix. Vats et al. [28] showed that terminating the simulation based on the effective sample size is equivalent to termination based on a relative confidence region where the Monte Carlo error is compared to size of the uncertainty in the target distribution. The point is that again a common method for assessing the reliability of the simulation is determined by the estimated generalized variance of the Monte Carlo error.

The estimators of $\Sigma$ studied by Kosorok [21], Vats et al. [28], and Vats et al. [29] typically underestimate the generalized variance. We will propose a different method and show that it is asymptotically valid. Specifically, our method provides a consistent overestimate for the asymptotic generalized variance of the Monte Carlo error and therefore will result in a slightly larger simulation effort, leading to a more stable estimation process.

The rest of the paper is organized as follows. In Section 2 we develop notation and background in preparation for the estimation theory. In Section 3 we propose our method and establish its asymptotic validity. In Section 4 we examine the finite sample properties of the proposed method through a variety of examples. We consider a Bayesian logistic regression example of 5 covariates where a symmetric random walk Metropolis–Hastings algorithm is implemented to calculate the posterior mean of the regression coefficient vector, a Bayesian one-way random effects model where we use a random scan Gibbs sampler to estimate the posterior expectation of all 8 parameters, and a reversible multivariate AR(1) process that takes values in 12. We illustrate the use of multivariate methods in a meta-analysis application where the posterior has dimension 65.

2 Notation and background

Recall that $F$ has support $\mathcal{X}$ and let $\mathcal{B}(\mathcal{X})$ be a $\sigma$-algebra. For $n\in\mathbb{N}^{+}=\{1,2,3,\ldots\}$ let $P^{n}(x,dy)$ be the $n$-step Markov transition kernel so that for $x\in\mathcal{X}$, $B\in\mathcal{B}(\mathcal{X})$, and $k\in\mathbb{N}=\{0,1,2,\ldots\}$ we have $P^{n}(x,B)=\Pr(X_{k+n}\in B\mid X_{k}=x)$, where $\Pr$ denotes probability. We assume that $P$ satisfies detailed balance with respect to $F$. That is,

[TABLE]

Metropolis–Hastings algorithms satisfy (3) by construction as do many component-wise Markov chains, such as random scan or random sequence scan algorithms [16]. By integrating both sides of (3) it is easy to see that $F$ is invariant for $P$. Suppose $X_{0}\sim F$, that is the Markov chain is stationary. The assumption of stationarity is not crucial since, for Harris recurrent chains, if a CLT holds under stationarity, it holds for all initial distributions [23, Proposition 17.1.6].

The lag $t$ autocovariance of the process $g(X_{0}),g(X_{1}),g(X_{2}),\ldots$ is defined as $\gamma_{t}=\gamma_{-t}=\mathrm{cov}_{F}\{g(X_{i}),g(X_{i+t})\}.$ Denote the sum of an adjacent pair of autocovariances by $\Gamma_{i}=\gamma_{2i}+\gamma_{2i+1}$ for $i\in\mathbb{N}$ and its smallest eigenvalue by $\xi_{i}$.

We use the shorthand $\infty$ for $+\infty$ unless otherwise specified. If $\sum_{t=0}^{\infty}\gamma_{t}$ converges, the asymptotic covariance matrix in (1) can be written as [20]

[TABLE]

The following propositions will play a significant role in the development of the new estimation method in Section 3.

Proposition 1**.**

The following properties of the sequences $\{\Gamma_{i}:i\in\mathbb{N}\}$ and $\{\xi_{i}:i\in\mathbb{N}\}$ hold.

1. (i)

$\Gamma_{i}$* is positive-definite, for all $i\in\mathbb{N}$.* 2. (ii)

$\Gamma_{i}-\Gamma_{i+1}$* is positive-definite, for all $i\in\mathbb{N}$.* 3. (iii)

$\lim_{i\to\infty}\Gamma_{i}=0$. 4. (iv)

The sequence $\{\xi_{i}:i\in\mathbb{N}\}$ is positive, decreasing, and converges to [math].

Proof.

See Appendix G. ∎

Recall (4) and let the $m$th partial sum be denoted

[TABLE]

Proposition 2**.**

The following properties of the sequence $\{\Sigma_{m}:m\in\mathbb{N}\}$ hold.

1. (i)

There exists a non-negative integer $m_{0}$ such that $\Sigma_{m}$ is positive definite for $m\geq m_{0}$ and not positive definite for $m<m_{0}$. Specifically, when $m_{0}=0$, $\Sigma_{m}$ is positive definite for all $m$. 2. (ii)

The sequence $\{|\Sigma_{m}|:m=m_{0},m_{0}+1,m_{0}+2,\ldots\}$ is positive, increasing, and converges to $|\Sigma|$.

Proof.

See Appendix G. ∎

Remark 1*.*

The value of $m_{0}$ is difficult to calculate explicitly because $\Sigma_{m}$ is usually not available in closed form. However, in Section 4.3 we consider a multivariate AR(1) Markov chain and verify that $m_{0}=0$. In the other examples, we cannot establish $m_{0}=0$ directly, but in our simulations we never observed anything else in 2000 independent replications.

3 Estimation method

A natural estimator of the lagged autocovariance $\gamma_{t}$ is the empirical autocovariance

[TABLE]

where ⊤ denotes transpose. Set $\widetilde{\gamma}_{n,t}=(\gamma_{n,t}+\gamma_{n,-t})/2$ for $t\in\{0,\ldots,n-1\}$ and write the sum of the $i$th ($0\leq i\leq\lfloor n/2-1\rfloor$) adjacent pair as $\Gamma_{n,i}=\widetilde{\gamma}_{n,2i}+\widetilde{\gamma}_{n,2i+1}.$ By construction, $\Gamma_{n,i}$ is symmetric. Let $\xi_{n,i}$ denote its smallest eigenvalue. The empirical estimator of $\Sigma_{m}$ ($0\leq m\leq\lfloor n/2-1\rfloor$) is

[TABLE]

Notice how (6) parallels (5).

3.1 Multivariate initial sequence estimators

We are now in position to formally define the multivariate initial sequence (mIS) estimator. Let $s_{n}$ be the smallest integer such that $\Sigma_{n,s_{n}}$ is positive definite and let $t_{n}$ be the largest integer $m$ ($s_{n}\leq m\leq\lfloor n/2-1\rfloor$) such that $|\Sigma_{n,i}|>|\Sigma_{n,i-1}|$ for all $i\in\{s_{n}+1,\ldots,m\}$. Then the mIS estimator, denoted $\Sigma_{\mathrm{seq},n}$, is defined as $\Sigma_{\mathrm{seq},n}=\Sigma_{n,t_{n}}.$ It is possible that $\Sigma_{n,m}$ fails to be positive definite for all $m\in\{0,\ldots,\lfloor n/2-1\rfloor\}$, and consequently $s_{n}$ does not exist. Fortunately, when $n$ is sufficiently large, we can always find such $s_{n}$.

Theorem 1**.**

With probability 1, $s_{n}$ exists as $n\to\infty$. In particular, with probability 1, $s_{n}\to m_{0}$ as $n\to\infty$.

Proof.

See Appendix H. ∎

Thus mIS is feasible while the following establishes that it is asymptotically valid.

Theorem 2**.**

With probability 1, $\liminf\limits_{n\to\infty}|\Sigma_{\mathrm{seq},n}|\geq|\Sigma|$.

Proof.

See Appendix H. ∎

In the construction of $\Sigma_{\mathrm{seq},n}$ we update $\Sigma_{n,i}$ to $\Sigma_{n,i+1}=\Sigma_{n,i}+2\Gamma_{n,i+1}$. If $\Gamma_{n,i+1}$ has negative eigenvalues, adding $2\Gamma_{n,i+1}$ will squeeze the corresponding confidence region in undesirable directions. A remedy is to force the negative eigenvalues of $\Gamma_{n,i+1}$ to be 0. Suppose $\Gamma_{n,i+1}$ has eigen-decomposition $\Gamma_{n,i+1}=Q^{\top}\Lambda Q$ where $\Lambda=\mathrm{diag}(\lambda_{1},\ldots,\lambda_{p})$. Define the positive part of $\Gamma_{n,i+1}$ as $\Gamma_{n,i+1}^{+}=Q^{\top}\Lambda^{+}Q$, where $\Lambda^{+}=\mathrm{diag}(\max\{\lambda_{1},0\},\ldots,\max\{\lambda_{p},0\})$.

This leads us to define the adjusted multivariate initial sequence (mISadj) estimator. Let $s_{n}$ and $t_{n}$ be as in the definition of mIS and let

[TABLE]

where $\Gamma^{+}_{n,i}$ is the positive part of $\Gamma_{n,i}$. Then the mISadj estimator, denoted $\Sigma_{\mathrm{adj},n}$, is defined as $\Sigma_{\mathrm{adj},n}=\widetilde{\Sigma}_{n,t_{n}}.$ See Figure 1 for a display of the effect of using mISadj over mIS.

By construction, the mISadj estimator is positive definite. The modification adds a positive semi-definite matrix to the mIS estimator, which by Theorem 2 provides a consistent overestimate for the generalized variance, $|\Sigma|$, and therefore the mISadj estimator also has a larger determinant than the asymptotic covariance matrix, $\Sigma$.

Theorem 3**.**

With probability 1, $\liminf\limits_{n\to\infty}|\Sigma_{\mathrm{adj},n}|\geq|\Sigma|$.

3.1.1 Related estimators

The motivation for our approach can be found in Geyer’s [9] univariate initial positive sequence (uIS) estimator. Suppose $\mu$ is one-dimensional and denote the variance of the asymptotic normal distribution $\sigma^{2}$. In this setting Geyer [9] proposed the uIS estimator

[TABLE]

where $t_{n}$ is the largest integer $m$ such that $\Gamma_{n,i}>0$ for all $i\in\{1,\ldots,m\}$. That is, Geyer’s truncation rule is to stop adding in $2\Gamma_{n,i}$ when it causes $\sigma^{2}_{n,i}=-\gamma_{n,0}+2(\Gamma_{n,0}+\cdots+\Gamma_{n,i})$ to decrease. (Figure 2 depicts the behavior of $\sigma^{2}_{n,i}$ and $\Gamma_{n,i}$ for one of the examples we consider later.) The uIS estimator is therefore the first local maximum of the sequence $\{\sigma^{2}_{n,i}:i=0,\ldots,\lfloor n/2-1\rfloor\}$ and thus gives an asymptotic overestimate of $\sigma^{2}$. This is formally stated in his Theorem 3.2:

[TABLE]

Neither mIS nor mISadj is a straightforward generalization of Geyer’s method in that mIS and mISadj coincide but do not reduce to uIS when $\mu$ is one-dimensional. However, this is not essential because the three methods are asymptotically equivalent in univariate settings.

Kosorok [21] proposed an alternative multivariate estimator (mK) which was also motivated by Geyer’s [9] approach. Recall from Proposition 1 that $\{\xi_{i};i\in\mathbb{N}\}$ is positive, decreasing, and converges to [math], where $\xi_{i}$ is the smallest eigenvalue of $\Gamma_{i}$. In mK the truncation point is chosen to be the largest integer $m$ such that $\xi_{n,i}>0$ for all $i\in\{1,\ldots,m\}$. However, this does not ensure that the generalized variance is adequately estimated and often truncates before the sequence $\{|\Sigma_{n,i}|:i=s_{n},\ldots,\lfloor n/2-1\rfloor\}$ reaches the first local maximum, as demonstrated in Figure 3 and 4.

4 Simulation experiments

Our goal is to investigate the finite-sample properties of mIS, mISadj, mK, and uIS through simulation experiments in a variety of examples. In each of the examples, which are described in more detail below, we compare the approaches in terms of effective sample size as well as volume and coverage probability of a joint confidence region.

We describe the simulation examples and the MCMC algorithms used in Section 4.1–4.3. The results of the simulation experiments are given in Section 4.4. We then consider a meta-analysis application in Section 4.5.

4.1 Bayesian logistic regression

For $i\in\{1,\ldots,100\}$, let $X_{i}=(x_{i1},\ldots,x_{i5})$ be the observed covariates for the $i$th observation and $Y_{i}$ be the binary response. We suppose

[TABLE]

This model results in a posterior on $\mathbb{R}^{5}$, denoted $F$. The data we use is provided in the logit dataset in the mcmc R package.

We are interested in estimating the posterior mean of $\beta$, i.e., $\mu=\mathrm{E}_{F}(\beta)$. However, this expectation is intractable and hence we will use a symmetric random walk Metropolis–Hastings algorithm to estimate it. At each step of the Markov chain, the proposal for the next step is $\mathcal{N}_{5}(0,\,0.3^{2}I_{5})$. The standard deviation of $0.3$ ensures that in our application the acceptance rate is about 0.36.

By construction, the Metropolis–Hastings algorithm satisfies detailed balance (3). Vats et al. [28] established that this Markov chain is geometrically ergodic and that the posterior has a moment generating function and hence a CLT as at (1) holds.

4.2 Bayesian one-way random effects model

Suppose for $i\in\{1,\ldots,K\}$,

[TABLE]

[TABLE]

where we assume the $a_{1},a_{2},a_{3},b_{1},b_{2},b_{3}$ and $v_{0}$ are known positive constants while $m_{0}$ is a known scalar. We consider a data set simulated under the settings $K=2$, $a_{1}=a_{2}=b_{1}=b_{2}=0.1$, $a_{3}=b_{3}=1.5$, $m_{0}=0$ and $v_{0}=0.001$. Let $y$ denote all of the data, $\lambda=(\lambda_{1},\ldots,\lambda_{K})^{\top}$, $\xi=(\theta_{1},\ldots,\theta_{K},\mu)^{\top}$, and $\gamma=(\gamma_{1},\ldots,\gamma_{K})^{\top}$. The hierarchy results in a proper posterior density $f(\xi,\lambda_{\theta},\lambda,\gamma\mid y)$ on $\mathbb{R}^{K+1}\times\mathbb{R}^{2K+1}_{+}$. One can verify that the posterior distribution has a finite second moment.

The posterior is intractable in the sense that posterior expectations are not generally available in closed form. We will use a random scan Gibbs sampler having the posterior as its invariant distribution to estimate the posterior expectation of all parameters. Doss and Hobert [4] derived the full conditional densities $f(\lambda_{\theta}\mid\xi,\lambda,\gamma)$, $f(\lambda\mid\xi,\lambda_{\theta},\gamma)$, $f(\gamma\mid\xi,\lambda_{\theta},\lambda)$, and $f(\xi\mid\lambda_{\theta},\lambda,\gamma)$ required to implement random scan Gibbs.

It is well known that the random scan Gibbs sampler kernel is reversible, namely, satisfies detailed balance (3), with respect to the posterior; see e.g., Roberts and Rosenthal [26]. Johnson and Jones [15] established geometric ergodicity of the random scan Gibbs sampler when $2a_{1}+K-2>0$ and $a_{3}>1$. These conditions combined with the second moment condition establish a Markov chain CLT.

4.3 Multivariate AR(1) process

Consider an AR(1) process $\{X_{n};n\in\mathbb{N}\}$ taking values in $\mathbb{R}^{p}$, i.e., $X_{n+1}=AX_{n}+U_{n+1}$, where $U_{n}$’s are iid $\mathbb{R}^{p}$-valued random variables and $A$ is a $p\times p$ matrix.

Ōsawa [24] proved that when $U_{n}$’s follow a normal distribution $\mathcal{N}_{p}(\theta,V)$, then this $\mathbb{R}^{p}$-valued AR(1) process satisfies detailed balance (3) if and only if the matrix $AV$ is symmetric. Suppose further that $\lim_{n\rightarrow\infty}A^{n}=0$, then it has the stationary distribution $\mathcal{N}_{p}[(I-A)^{-1}\theta,(I-A^{2})^{-1}V]$. It is easy to verify that the second moment is finite.

Under stationarity one can derive the lag $t$ autocovariance, $\gamma_{t}=A^{2t}(I-A^{2})^{-1}V$, and hence the covariance matrix, $\Sigma=\{2(I-A^{2})^{-1}-I\}(I-A^{2})^{-1}V$, as in (4). Noticing that $\Sigma$ is finite, and that the Markov chain is reversible with a finite second moment, we establish a Markov chain CLT (1) with mean $\mu=(I-A)^{-1}\theta$ and covariance matrix $\Sigma$ [11, Corollary 6]. Also notice that $\Sigma_{0}=\gamma_{0}+2\gamma_{1}$ is always positive definite, which satisfies the assumption in Remark 1 and hence guarantees the asymptotic properties of our proposed estimation method.

Let us consider the following choices that satisfy the conditions above: $\theta=\textbf{1}_{p}$, $V=I_{p}$, and $A=p^{-1}H_{p}\mathrm{diag}(2^{-1}$, $\ldots,2^{-p})H_{p}^{\top}$, where $H_{p}$ is a Hadamard matrix of order $p$. We set $p=12$ in our simulation study.

4.4 Results

In this section we refer to the setting of Section 4.1 as Example 1, the setting of Section 4.2 as Example 2, and the setting of Section 4.3 as Example 3. For all examples we ran 2000 independent replications of the Markov chain for $10^{6}$ iterations in Examples 1 and 3 and $5\times 10^{5}$ iterations in Example 2, respectively. We will compare the multivariate methods—namely mIS, mISadj, and mK—in the context of estimating the effective sample size. We then turn our attention to the finite-sample properties of the confidence regions produced by the multivariate methods, yielding ellipsoidal regions, and Geyer’s univariate uIS for individual components, yielding cube-shaped regions. To assess coverage probabilities in Examples 1 and 2 we perform an independent run of length $10^{10}$ of the Markov chain in each example and declared the sample average over those $10^{10}$ iterations to be the truth, while in Example 3, the true mean is obtained through the closed form expression derived.

The results concerning effective sample size of the simulation experiments are given in Table 1. Prior to the work of Vats et al. [28] it was standard to report the minimum of the univariate effective sample size calculated component-wise. This leads to a substantial underestimate of the effective sample size as can be seen in Table 1. In contrast, multivariate error estimation yields more accurate evaluation of the effective sample size. We can approximately order the multivariate methods in terms of estimated effective sample size: mK $>$ mIS $>$ mISadj. That is, mK is more optimistic than mIS and mISadj.

We construct $90\%$ confidence regions using the multivariate estimation methods and uIS. Throughout “uIS” and “uIS-Bonferroni” represent the uncorrected and Bonferroni corrected confidence regions generated by uIS, respectively. Let us first examine the volumes of the confidence regions generated by different methods.

The volumes are presented in ascending order from left to right across Table 2. The uncorrected uIS confidence regions are much smaller than the other methods, while the Bonferroni correction considerably enlarges the confidence regions, resulting in bigger volumes than all the multivariate methods.

Recall that the volume of a confidence region depends on the estimated covariance matrix only through the estimated generalized variance of the Monte Carlo error. Therefore, Table 2 compares the estimation of the generalized variance by different multivariate methods. We observe that mK underestimates the generalized variance relatively to mIS. The mISadj method is comparable to mIS in Examples 1 and 3 but clearly overestimates in Example 2.

Table 3 shows the empirical coverage probabilities of the confidence regions produced by different methods. The proposed method, mIS, exceeds mK in both the volume and the coverage of confidence regions, although the coverage rate does not always reach the expectation. The adjustment moderately increases the coverage probability.

The uncorrected uIS regions have a poor coverage. The Bonferroni regions work well in these examples, but in high-dimensional cases the Bonferroni correction can be overly conservative. Overall, multivariate error estimation methods yield better confidence regions.

4.5 A meta-analysis example

Doss and Hobert [4] carried out meta-analyses in order to study the effect of non-steroidal anti-inflammatory drugs (NSAIDs) on the risk of colon cancer. The dataset consists of 21 studies that relate NSAIDs intake and risk of colon cancer; see Harris et al. [12] and Doss and Hobert [4] for details. We apply the Bayesian one-way random effects model described in Section 4.2 to the colon cancer dataset. The posterior $f(\theta_{1},\ldots,\theta_{K},\mu,\lambda_{\theta},\lambda_{1},\ldots,\lambda_{K},\gamma_{1},\ldots,\gamma_{K}\mid y)$ has dimension $p=65$ when $K=21$.

We run a Markov chain for $4\times 10^{6}$ iterations and compute the multivariate estimators—namely mIS, mISadj, and mK—along with Geyer’s uIS for individual components.

The estimated generalized variances are reported in Table 4. The result agrees with our conclusion from the previous simulation study: mISadj is more conservative than mIS; mK clearly underestimates the generalized variance.

Table 5 shows the estimated effective sample sizes. The uIS method results in 65 estimated effective sample sizes, each of which corresponds to a component of the posterior distribution. Only the minimum estimated univariate effective sample size is reported.

An advantage of using multivariate methods like mIS over univariate estimation like uIS is that only multivariate methods capture the cross-correlation between components. This cross-correlation is often significant as seen in Figure 5.

We construct $90\%$ confidence regions using the multivariate estimation methods and uIS. The left panel of Figure 6 shows the cross-sections of the confidence regions that are cut through the center of the confidence regions parallel to the plane spanned by $\mu$ and $\lambda_{1}$. The reader should not be worried that the cross-sectioned ellipsoids appear much larger than the Bonferroni region.

The full 65-dimensional ellipsoid will have a smaller volume than the 65-dimensional Bonferroni region, but this does not have to be the case for cross-sectioned regions. As a comparison, in the right panel of Figure 6 we present bivariate $90\%$ confidence regions for $\mu$ and $\lambda_{1}$ when we ignore the other 63 components. This clearly shows how multivariate estimation methods generate confidence regions that are not so liberal as uIS, yet not so conservative as uIS-Bonferroni.

Table 6 compares the volumes of the confidence regions generated by different methods. The results agree with our conclusion from the previous simulation study: mISadj is slightly more conservative than mIS; mK clearly underestimates the generalized variance. The volumes generated by multivariate estimators are fairly close to each other but the univariate results are far away. Apparently uIS is too liberal while uIS-Bonferroni is too conservative, but the multivariate methods achieve a balance.

4.6 Discussion

The preceding simulation examples and the theory developed indicate that mIS and mISadj perform as they were designed to in that they provide a consistent overestimate of the asymptotic generalized variance of the Monte Carlo error. Compared to standard univariate methods, our estimators adjust for multivariate issues and thus provide more realistic estimates of Monte Carlo effective sample size and slightly larger confidence regions which result in improved performance in terms of coverage probabilities.

Acknowledgments. The authors are grateful to Charles Geyer and Dootika Vats for helpful conversations.

Appendices

G Proofs of Propositions 1 and 2

We begin with some preliminary results which will be useful later.

Lemma 1** (Harville [13], Lemma 18.2.17).**

Let $A_{0},A_{1},A_{2},\ldots$ represent a sequence of $m\times n$ matrices. If the infinite series $\sum_{k=0}^{\infty}A_{k}$ converges, then $\lim_{k\to\infty}A_{k}=0$.

Since the eigenvalues of a Hermitian $p\times p$ matrix $A$ are real, we may (and do) adopt the convention that they are always arranged in algebraically non-decreasing order:

[TABLE]

Lemma 2** (Horn and Johnson [14], Corollary 4.3.15).**

Let $p\times p$ matrices $A$, $B$ be Hermitian and let the respective eigenvalues of $A$, $B$, and $A+B$ be $\{\lambda^{(k)}(A):k=1,\ldots,p\}$, $\{\lambda^{(k)}(B):k=1,\ldots,p\}$, and $\{\lambda^{(k)}(A+B):k=1,\ldots,p\}$, each algebraically ordered as in (A.1). Then, for all $k\in\{1,\ldots,p\}$,

[TABLE]

Lemma 3**.**

Suppose we have two $p\times p$ Hermitian matrices $A$ and $B$. Let the respective eigenvalues of $A$ and $B$ be $\{\lambda^{(k)}(A):k=1,\ldots,p\}$ and $\{\lambda^{(k)}(B):k=1,\ldots,p\}$, each algebraically ordered as in (A.1). If $A-B$ is positive definite, then $\lambda^{(k)}(A)>\lambda^{(k)}(B)$, for all $k\in\{1,\ldots,p\}$. Further, if $A$ and $B$ are both positive semi-definite, then $|A|>|B|$.

Proof.

Applying (A.2) to $B$ and $A-B$, we get for all $k\in\{1,\ldots,p\}$, $\lambda^{(k)}(B)+\lambda^{(1)}(A-B)\leq\lambda^{(k)}(A).$ Since $A-B$ is positive definite, we have $\lambda^{(1)}(A-B)>0$. Therefore, for all $k\in\{1,\ldots,p\}$, $\lambda^{(k)}(A)>\lambda^{(k)}(B)$. When $A$ and $B$ are both positive semi-definite, we further have for all $k\in\{1,\ldots,p\}$,

[TABLE]

Since the determinant is equal to the product of all eigenvalues, we take product of (A.3) for all $k\in\{1,\ldots,p\}$ and obtain $|A|>|B|$. ∎

Lemma 4** (Vats et al. [29], Theorem 2).**

Let $\Sigma_{n}$ be a strongly consistent estimator of $\Sigma$. Let the respective eigenvalues of $\Sigma_{n}$ and $\Sigma$ be $\{\lambda^{(k)}(\Sigma_{n}):k=1,\ldots,p\}$ and $\{\lambda^{(k)}(\Sigma):k=1,\ldots,p\}$, each algebraically ordered as in (A.1). Then $\lambda^{(k)}(\Sigma_{n})\to\lambda^{(k)}(\Sigma)~{}\text{with probability}~{}1$ as $n\to\infty$ for all $k\in\{1,\ldots,p\}$.

Corollary 1**.**

Let $\Sigma_{n}$ be a strongly consistent estimator of $\Sigma$, then $|\Sigma_{n}|\to|\Sigma|~{}\text{with probability}~{}1$ as $n\to\infty$.

G.1 Proof of Proposition 1

We begin with the univariate case so $g:\mathcal{X}\to\mathbb{R}$. Let $\mathcal{E}$ be the spectral decomposition measure associated with transition kernel $P$ and $\mathcal{E}_{g}$ be the induced spectral measure for $g$. Details on the spectral decomposition measure can be found in Rudin [27], Chan and Geyer [2], and Häggström and Rosenthal [11]. Specifically, for all $t\in\mathbb{N}$,

[TABLE]

It follows that for all $i\in\mathbb{N}$,

[TABLE]

and

[TABLE]

Therefore, $\Gamma_{i}$ and $\Gamma_{i}-\Gamma_{i+1}$ must be non-negative. To prove Proposition 1(i) that $\Gamma_{i}>0$ and Proposition 1(ii) that $\Gamma_{i}-\Gamma_{i+1}>0$, we need to show that neither $\Gamma_{i}$ nor $\Gamma_{i}-\Gamma_{i+1}$ can be zero. For $i=0$,

[TABLE]

and

[TABLE]

For $i\in\mathbb{N}^{+}$,

[TABLE]

and

[TABLE]

By (A.5)–(A.8), for an arbitrary $i\in\mathbb{N}$, a necessary condition for each of $\Gamma_{i}=0$ and $\Gamma_{i}-\Gamma_{i+1}=0$ is $\mathcal{E}_{g}(\{-1,0,1\})=1$. We now show that $\mathcal{E}_{g}(\{-1,0,1\})=1$ cannot hold under our assumptions, so that both $\Gamma_{i}$ and $\Gamma_{i}-\Gamma_{i+1}$ are non-zero, which completes the proof of Proposition 1(i)–(ii).

If $\mathcal{E}_{g}$ is a point mass at 0, then (A.4) yields

[TABLE]

for all $t\in\mathbb{N}^{+}$, which is trivial. Therefore, without loss of generality, we assume

[TABLE]

Häggström and Rosenthal [11] showed that when $P$ is irreducible and aperiodic,

[TABLE]

It follows from (A.9) and (A.10) that $\mathcal{E}_{g}(\{-1,0,1\})<1$. By previous arguments, we have proved Proposition 1(i)-(ii). That is, for all $i\in\mathbb{N}$, $\Gamma_{i}>0$ and $\Gamma_{i}-\Gamma_{i+1}>0$.

Proposition 1(iii), namely $\lim_{i\to\infty}\Gamma_{i}=0$, follows from Lemma 1 and the assumption that $\sum_{i=0}^{\infty}\Gamma_{i}$ converges.

Finally, by Proposition 1(i)–(iii), we obtain Proposition 1(iv), i.e., $\{\Gamma_{i}:i\in\mathbb{N}\}$ is positive, decreasing, and converges to [math].

We now turn to the multivariate case so $g:\mathcal{X}\to\mathbb{R}^{p}$ and $p\geq 2$. Set $h=v^{\top}g$ for an arbitrary $v\in\mathbb{R}^{p}$ and $v\neq 0$. Then $h:\mathcal{X}\to\mathbb{R}$ is measurable and square integrable with respect to $F$. Recall that the Markov chain is assumed stationary. For $t\in\mathbb{N}$ define the lag $t$ autocovariance

[TABLE]

and for $i\in\mathbb{N}$ define $\Gamma_{i}^{*}=\gamma_{2i}^{*}+\gamma_{2i+1}^{*}$. Notice that

[TABLE]

By the univariate case considered above, $\Gamma_{i}^{*}>0$. Since $v$ is arbitrary, $\Gamma_{i}$ is positive definite. A similar argument shows that $\Gamma_{i}-\Gamma_{i+1}$ is positive definite. This establishes Proposition 1(i)–(ii).

Use Lemma 1 and notice that $\sum_{i=0}^{\infty}\Gamma_{i}$ converges by assumption. We obtain $\lim_{i\to\infty}\Gamma_{i}=0$. Thus Proposition 1(iii) is proved.

Since $\Gamma_{i}$ is positive definite, $\xi_{i}>0$ for all $i\in\mathbb{N}$. Since $\Gamma_{i}-\Gamma_{i+1}$ is positive definite we obtain from Lemma 3 that $\xi_{i}>\xi_{i+1}$. Hence $\xi_{i}\to 0$ as $i\to\infty$ which establishes Proposition 1(iv). $\Box$

G.2 Proposition 2

For all $m\in\mathbb{N}$ let $\lambda_{m}$ be the smallest eigenvalue of $\Sigma_{m}$. Notice that $\Sigma_{m}-\Sigma_{m-1}=2\Gamma_{m}$ is positive definite by Proposition 1. Then Lemma 3 implies $\lambda_{m}>\lambda_{m-1}$ and hence$\{\lambda_{m}:m\in\mathbb{N}\}$ is monotonically increasing. Since $\{\Sigma_{m}:m\in\mathbb{N}\}$ converges to the asymptotic covariance matrix $\Sigma$, by Lemma 4 we have

[TABLE]

where $\lambda$ is the smallest eigenvalue of $\Sigma$.

If $\lambda_{0}\leq 0$, there exists a positive integer $m_{0}$ such that $\lambda_{m}>0$ for $m\geq m_{0}$ and $\lambda_{m}\leq 0$ for $m<m_{0}$. If $\lambda_{0}>0$, then $\lambda_{m}>0$ for all $m\in\mathbb{N}$. In this case, let $m_{0}=0$. Immediately we have that $\Sigma_{m}$ is positive definite for $m\geq m_{0}$ and not positive definite for $m<m_{0}$. It then follows that for all $m\geq m_{0}$, $|\Sigma_{m}|>0$.

Now let $m>m_{0}$ and notice that $\Sigma_{m}-\Sigma_{m-1}$ is positive definite. Using Lemma 3 we obtain for all $m>m_{0}$, $|\Sigma_{m}|>|\Sigma_{m-1}|.$ By $\lim_{m\to\infty}\Sigma_{m}=\Sigma$ and Corollary 1,

[TABLE]

Therefore, $\{|\Sigma_{m}|:m=m_{0},m_{0}+1,m_{0}+2,\ldots\}$ is positive, increasing, and converges to $|\Sigma|$. $\Box$

H Proofs of Theorems 1 and 2

Lemma 5**.**

For all $t\in\mathbb{N}$, with probability 1, as $n\to\infty$, $\gamma_{n,t}\to\gamma_{t}.$

Proof.

Notice that

[TABLE]

By repeated application of the Markov chain strong law we see that, with probability 1, as $n\to\infty$,

[TABLE]

∎

Corollary 2**.**

For all $m\in\mathbb{N}$, with probability 1, as $n\to\infty$, $\Sigma_{n,m}\to\Sigma_{m}.$

Proof.

This follows immediately from Lemma 5. ∎

Lemma 6**.**

If a sequence of random variables $X_{1},X_{2},\ldots$ converges to $X$ with probability 1, then, for an arbitrary $x\in\mathbb{R}$ such that $\Pr(X=x)=0$,

[TABLE]

and

[TABLE]

Proof.

We only prove the first part. The second part can be shown by a similar argument.

Recall that two events $A$ and $B$ are equal almost surely if both of the events $A\setminus B$ and $B\setminus A$ are null sets [8, p. 13]. Thus we need only show that both $\liminf_{n\rightarrow\infty}\{X_{n}\leq x\}\setminus\{X\leq x\}$ and $\{X\leq x\}\setminus\liminf_{n\rightarrow\infty}\{X_{n}\leq x\}$ are null sets.

Suppose $\omega\in\liminf_{n\rightarrow\infty}\{X_{n}\leq x\}\setminus\{X\leq x\}$. By definition,

[TABLE]

is equivalent to saying that there exists some $n$ such that for all $m\geq n$, $X_{m}(\omega)\leq x$. This implies that

[TABLE]

where the second inequality is due to $\omega\notin\{X\leq x\}$. It follows that

[TABLE]

Thus we have that

[TABLE]

which is a null set because $X_{n}\overset{a.s.}{\rightarrow}X$.

Suppose $\omega\in\{X\leq x\}\setminus\liminf_{n\rightarrow\infty}\{X_{n}\leq x\}$. By definition,

[TABLE]

is equivalent to saying that for all $n$, there exists some $m\geq n$ such that $X_{m}(\omega)>x$. This implies that

[TABLE]

where the second inequality is due to $\omega\in\{X\leq x\}$. It follows that

[TABLE]

Thus we have that

[TABLE]

which is a null set.

So far we have proved (A.11). A similar argument can be used to prove (A.12). ∎

Equipped with the preceding results, we now prove the following lemma in preparation for Theorems 1 and 2.

Recall that $m_{0}$ is a non-negative integer such that $\Sigma_{m}$ is positive definite for $m\geq m_{0}$ and not positive definite for $m<m_{0}$. Also recall that $s_{n}$ is the smallest integer such that $\Sigma_{n,s_{n}}$ is positive definite and that $t_{n}$ is the largest integer $m$ ($s_{n}\leq m\leq\lfloor n/2-1\rfloor$) such that $|\Sigma_{n,i}|>|\Sigma_{n,i-1}|$ for all $i\in\{s_{n}+1,\ldots,m\}$. The smallest eigenvalues of $\Sigma_{m}$ and $\Sigma_{n,m}$ are denoted $\lambda_{m}$ and $\lambda_{n,m}$, respectively.

Lemma 7**.**

Suppose $\liminf_{n\rightarrow\infty}\{\lambda_{n,m_{0}-1}\leq 0\}$ occurs with probability 1. For all $K\geq m_{0}$,

[TABLE]

Proof.

Define $\Delta_{m}=|\Sigma_{m}|-|\Sigma_{m-1}|$ and $\Delta_{n,m}=|\Sigma_{n,m}|-|\Sigma_{n,m-1}|$. Notice that

[TABLE]

where $\lambda_{n,m}$ denotes the smallest eigenvalue of $\Sigma_{n,m}$. Then we write

[TABLE]

By Lemma 4, Corollary 1 and Corollary 2, for all $m$, with probability 1,

[TABLE]

and for all $i$, with probability 1,

[TABLE]

By Proposition 2(i), $\lambda_{m_{0}}>0$ so that $\Pr(\lambda_{m_{0}}=0)=0$ and $\lambda_{m}\leq 0$ for all $m<m_{0}$. In particular, $\lambda_{m}<0$ so that $\Pr(\lambda_{m}=0)=0$ for all $m<m_{0}-1$. By Proposition 2(ii), $\Delta_{i}>0$ so that $\Pr(\Delta_{i}=0)=0$ for $i>m_{0}$. Then by Lemma 6 we have that for all $m<m_{0}-1$,

[TABLE]

[TABLE]

and for $i>m_{0}$

[TABLE]

Notice that (A.16) holds for $m=m_{0}-1$ if $\lambda_{m_{0}-1}<0$. When $\lambda_{m_{0}-1}=0$, (A.16) is true only if $\liminf_{n\rightarrow\infty}\{\lambda_{n,m_{0}-1}\leq 0\}$ occurs with probability 1.

Under the preceding assumption, we continue to write (A.13) as

[TABLE]

By Proposition 2(i), $\lambda_{m_{0}}>0$ and $\lambda_{m}\leq 0$ for all $m<m_{0}$. Then by (A.14) we have for $m<m_{0}$

[TABLE]

and

[TABLE]

By Proposition 2(ii), $\Delta_{i}>0$ for $i>m_{0}$. Then by (A.15) we have

[TABLE]

It follows from (A.18)–(A.20) that

[TABLE]

Then by (A.17) we obtain the result. ∎

Remark 2*.*

Consider the assumption that $\liminf_{n\rightarrow\infty}\{\lambda_{n,m_{0}-1}\leq 0\}$ occurs with probability 1. If $m_{0}=0$, then this assumption is not required for the Lemma; recall Remark 1. In addition, the assumption holds if $\Sigma_{m_{0}-1}$ is not positive semi-definite. Recall from Proposition 2(i) we have that $\Sigma_{m_{0}-1}$ is not positive definite but, of course, it may still be positive semi-definite.

H.1 Theorem 1: Feasibility of the estimation method

Proof.

When $K\geq m_{0}$ and $n>2m_{0}$, $\{s_{n}\text{ exists}\}\supset\{s_{n}=m_{0}\}\supset\{s_{n}=m_{0},t_{n}\geq K\}.$ Then the result follows from Lemma 7. ∎

H.2 Theorem 2: Overestimation for the Asymptotic Generalized Variance of the Monte Carlo Error

Proof.

We need to prove, for all $\epsilon>0$,

[TABLE]

Recall that $\Delta_{i}$ is defined as $\Delta_{i}=|\Sigma_{i}|-|\Sigma_{i-1}|$.

By Proposition 2(ii) that $\lim_{m\to\infty}|\Sigma_{m}|=|\Sigma|$, we can write

[TABLE]

so $\sum\limits_{i=m_{0}+1}^{\infty}\Delta_{i}$ converges; and hence the tail must converge to 0. Therefore, for all $\epsilon>0$, there exists $K_{\epsilon}\geq m_{0}$ such that

[TABLE]

Notice that

[TABLE]

The second step in (A.23) is due to the definition of mIS:

“$s_{n}=m_{0}$ and $t_{n}\geq K_{\epsilon}$ for some $K_{\epsilon}\geq m_{0}$” implies “$|\Sigma_{n,t_{n}}|\geq|\Sigma_{n,K_{\epsilon}}|$”.

It follows directly from (A.23) that

[TABLE]

Therefore, to prove (A.21) it suffices to show

[TABLE]

and

[TABLE]

By the continuity of measure, (A.25) is equivalent to Lemma 7 and thus holds true. Then it remains to prove (A.24).

By Corollary 1 and 2, with probability 1, $|\Sigma_{n,K_{\epsilon}}|\to|\Sigma_{K_{\epsilon}}|\text{ as }n\to\infty,$ which gives

[TABLE]

where $\mathrm{abs}(\cdot)$ denotes absolute value.

When $|\Sigma|-|\Sigma_{K_{\epsilon}}|<\epsilon/2$ as in (A.22), $\mathrm{abs}(|\Sigma_{n,K_{\epsilon}}|-|\Sigma_{K_{\epsilon}}|)<\epsilon/2~{}\text{ implies }|\Sigma_{n,K_{\epsilon}}|>|\Sigma|-\epsilon,$ so

[TABLE]

Putting (A.26) and (A.27) together, we obtain (A.24). ∎

I Confidence region with the univariate approach

We briefly state here the current methods for constructing confidence regions with univariate estimators. Let $\sigma(i)^{2}$ denote the $(i,i)$th entry of $\Sigma$. We treat the problem as $p$ univariate cases, i.e., to estimate $\sigma(i)^{2}$ using univariate samples. Then we construct cube-shaped confidence regions.

Let $\mu_{n}(i)$ be the $i$th component of $\mu_{n}$, and $\sigma_{n}(i)^{2}$ be the estimator for $\sigma(i)^{2}$. The uncorrected confidence region is given by

[TABLE]

with a volume of

[TABLE]

The Bonferroni confidence region for $\mu$ is

[TABLE]

with a volume of

[TABLE]

References