Assessing and Visualizing Simultaneous Simulation Error

TL;DR

This paper develops a multivariate CLT for simultaneous estimation of means and quantiles in Monte Carlo experiments, providing a practical visualization tool for assessing the reliability of these estimates.

Contribution

It introduces a multivariate CLT for combined means and quantiles under mixing conditions and a fast algorithm for hyperrectangular confidence regions with simultaneous coverage.

Findings

The method applies to i.i.d. and MCMC samples.

It enables visual assessment of Monte Carlo estimate reliability.

The approach improves interpretation of simulation results.

Abstract

Monte Carlo experiments produce samples in order to estimate features of a given distribution. However, simultaneous estimation of means and quantiles has received little attention, despite being common practice. In this setting we establish a multivariate central limit theorem for any finite combination of sample means and quantiles under the assumption of a strongly mixing process, which includes the standard Monte Carlo and Markov chain Monte Carlo settings. We build on this to provide a fast algorithm for constructing hyperrectangular confidence regions having the desired simultaneous coverage probability and a convenient marginal interpretation. The methods are incorporated into standard ways of visualizing the results of Monte Carlo experiments enabling the practitioner to more easily assess the reliability of the results. We demonstrate the utility of this approach in various Monte Carlo settings including simulation studies based on independent and identically distributed samples and Bayesian analyses using Markov chain Monte Carlo sampling.

Equations

$${\theta }_{g_i}=E_{\pi }\left(g_i(X)\right)={\int }_{\mathsf{X}}{g}_i(x)\pi (dx),$$

$${\theta }_g=({\theta }_{g_1},\dots ,{\theta }_{g_{p_1}}{)}^{\prime }.$$

$${\xi }_{q_i}=F_{h_i}^{-1}(q_i)=\inf \{v:F_{h_i}(v)\ge q_i\}.$$

$${\varphi }_h={\left({\xi }_{q_1},\dots ,{\xi }_{q_{p_2}}\right)}^{\prime }.$$

$$\alpha (n)=\underset{k\ge 1}{\sup }\underset{A\in {\mathcal{F}}_1^k,B\in {\mathcal{F}}_{k+n}^{\infty }}{\sup }|P(A\cap B)-P(A)P(B)|.$$

$${\bar{g}}_n=\frac{1}{n}\sum _{j=1}^{n}g(X_j)\to {\theta }_g.$$

$${\hat{\varphi }}_n={\left({\hat{\xi }}_{q_1},\dots ,{\hat{\xi }}_{q_{p_2}}\right)}^{\prime }.$$

$${\bar{F}}_{h_i}(v)=\frac{1}{n}\sum _{j=1}^{n}I(h_i(X_j)\le v),$$

$${\bar{F}}_h(v)={\left({\bar{F}}_{h_1}(v_1),\dots ,{\bar{F}}_{h_{p_2}}(v_{p_2})\right)}^{\prime }.$$

$$\Lambda =\left(\begin{array}{cc}{I}_{p_1}& 0_{p_1\times p_2}\\ 0_{p_2\times p_1}& A_h\end{array}\right).$$

$$\Sigma =\text{Cov}(Y_1,Y_1)+\sum _{j=2}^{\infty }\left[\text{Cov}(Y_1,Y_j)+\text{Cov}(Y_1,Y_j{)}^{\prime }\right]$$

$$\sqrt{n}\left(\begin{array}{c}{\bar{g}}_n-{\theta }_g\\ {\hat{\varphi }}_n-\varphi \end{array}\right)\stackrel{d}{\to }N\left(0,{\Lambda }^{-1}\Sigma {\Lambda }^{-1}\right).$$

$${\Sigma }_g=\text{Cov}(g(X_1),g(X_1))+\sum _{j=2}^{\infty }\left[\text{Cov}(g(X_1),g(X_j))+\text{Cov}(g(X_1),g(X_j){)}^{\prime }\right].$$

$$\sqrt{n}({\bar{g}}_n-{\theta }_g)\stackrel{d}{\to }{N}_p(0,{\Sigma }_g).$$

$${\Sigma }_h$$

$$+\sum _{j=2}^{\infty }\left[\text{Cov}(I(h(X_1)>{\varphi }_h),I(h(X_j)>{\varphi }_h))+\text{Cov}(I(h(X_1)>{\varphi }_h),I(h(X_j)>{\varphi }_h){)}^{\prime }\right]$$

$$\sqrt{n}\left(\left(\begin{array}{c}{\bar{g}}_n\\ 1-{\bar{F}}_h({\varphi }_h)\end{array}\right)-\left(\begin{array}{c}{\theta }_g\\ 1-Q\end{array}\right)\right)\stackrel{d}{\to }{N}_{p_1+p_2}\left(0,\Sigma =\left(\begin{array}{cc}{\Sigma }_g& {\Sigma }_{gh}\\ {\Sigma }_{hg}& {\Sigma }_h\end{array}\right)\right).$$

$${\hat{\xi }}_{q_i}={\xi }_{q_i}+\frac{\left(1-{\bar{F}}_{h_i}({\xi }_{q_i})\right)-(1-q_i)}{f_{h_i}({\xi }_{q_i})}+r_{n,q_i},$$

$$\left(1-{\bar{F}}_h({\varphi }_h)\right)-\left(1-Q\right)=A_h\left({\hat{\varphi }}_n-{\varphi }_h\right)+A_h{R}_n\stackrel{p}{\to }{A}_h\left({\hat{\varphi }}_n-{\varphi }_h\right).$$

$$\begin{array}{rl}\sqrt{n}\left(\left(\begin{array}{c}{\bar{g}}_n\\ 1-{\bar{F}}_h({\varphi }_h)\end{array}\right)-\left(\begin{array}{c}{\theta }_g\\ 1-Q\end{array}\right)\right)& =\sqrt{n}\left(\begin{array}{c}{\bar{g}}_n-{\theta }_g\\ A_h\left({\hat{\varphi }}_n-{\varphi }_h\right)\end{array}\right)+o_p(1)\\ & =\sqrt{n}\Lambda \left(\begin{array}{c}{\bar{g}}_n-{\theta }_g\\ {\hat{\varphi }}_n-{\varphi }_h\end{array}\right)+o_p(1).\end{array}$$

$$K^n(x,A)=\mathrm{Pr}(X_{n+j}\in A\mid X_j=x)$$

$$\Vert K^n(x,\cdot )-\pi (\cdot )\Vert \le M(x)\psi (n).$$

$$\alpha (n)\le \psi (n){\text{E}}_{\pi }M$$

$$\hat{\Sigma }=\frac{b}{a-1}\sum _{k=0}^{a-1}\left({\bar{S}}_k(b)-\bar{S}\right){\left({\bar{S}}_k(b)-\bar{S}\right)}^{\prime }.$$

$$C_{SI}(z)=\prod _{i=1}^{p_1}\left[{\bar{g}}_{n,i}-z\frac{{\hat{\sigma }}_{g_i}}{n},{\bar{g}}_{n,i}+z\frac{{\hat{\sigma }}_{g_i}}{n}\right]\prod _{j=1}^{p_2}\left[{\hat{\varphi }}_{n,j}-z\frac{{\hat{\gamma }}_j}{n},{\hat{\varphi }}_{n,j}+z\frac{{\hat{\gamma }}_j}{n}\right].$$

$$\pi (x)=.3f_1(x;1,2.5)+.5f_2(x;5,4)+.2f_3(x;11,3).$$

$$Y=Z{\beta }^{\ast }+ϵ.$$

$$Y_j\mid {\theta }_j\stackrel{ind}{\sim }N({\theta }_j,{\sigma }_j^2)$$

$${\theta }_j\stackrel{ind}{\sim }N(\mu ,{\tau }^2),$$

Full text

Keywords. Asymptotic normality, Monte Carlo, Markov chain Monte Carlo, quantile limit theorems, simulation studies, strong mixing, visualizations.

1 Introduction

Monte Carlo experiments are a standard tool to estimate features of a given distribution (see e.g. Brooks et al.,, 2010; Fishman,, 1996). The justification for using Monte Carlo methods is asymptotic. For example, using a sample mean to estimate a mean of the distribution is justified by an appropriate strong law. However, due to variability arising from simulation in finite time, it is essential to address whether the simulation has been run sufficiently long so that the resulting estimates are reliable.

Reporting of results from simulation studies often focuses on estimates from empirical marginal distributions, e.g., by using sample boxplots of estimates of one feature of interest. Popular software such as Stan (Stan Development Team,, 2018) and tidybayes (Kay,, 2018) follow similar procedures for reporting credible intervals or prediction intervals. However, consideration of simulation reliability is rare; see Flegal et al., (2008), Jones et al., (2006), and Koehler et al., (2009) for more discussion. Despite these works from more than a decade ago, the situation has not substantially improved.

When simulation reliability is addressed, it usually amounts to either considering the Monte Carlo sample size (or effective sample size) or, in the best case, univariate Monte Carlo standard errors. In Markov chain Monte Carlo (MCMC) settings, so-called convergence diagnostics are often reported. Neither the Monte Carlo sample size (or effective sample size) nor convergence diagnostics directly assess the reliability of the resulting estimation. Univariate Monte Carlo standard errors ignore the fact that most simulation experiments are aimed at estimating multiple quantities with the same simulated data, forcing the use of conservative multiplicity adjustments such as Bonferroni or Scheffé, leading to inefficient simulation practices.

Although there has been recent work on assessing reliability through Monte Carlo error in simultaneous estimation of several expectations in MCMC experiments (see e.g. Dai and Jones,, 2017; Vats et al.,, 2018, 2019), current best practice does not apply to simultaneous estimation of several quantiles or combinations of means and quantiles, a common goal (e.g. in Bayesian applications). Moreover, current multivariate approaches focus on either a relative error interpretation (comparing the Monte Carlo variation to the target distribution variation) or on minimum volume ellipsoidal confidence regions. The relative error approach applies well in MCMC experiments, but does not translate to standard Monte Carlo, that is, independent and identically distributed (IID) sampling. The use of ellipsoidal regions makes visualization difficult and complicates the interpretation of marginal estimates.

We aim to address the reliability of simultaneous Monte Carlo estimation of any finite combination of means and quantiles under the assumption that the simulated process is strongly mixing, which includes standard Monte Carlo and correlated sampling from MCMC. Our approach will also yield a convenient marginal interpretation that, as we will demonstrate, can be easily incorporated into standard existing visualizations of Monte Carlo estimates, making the assessment of estimation reliability convenient to the practitioner.

The techniques we consider provide simultaneous confidence intervals for any combination of quantiles and expectations with the desired level, say $1-\alpha$. Thus we avoid the use of conservative methods where the overall coverage probability often greatly exceeds $1-\alpha$. Our approach begins by establishing asymptotic (multivariate) normality of any finite collection of sample means and quantiles. The covariance matrix of the asymptotic distribution has a complicated form, but we provide consistent estimators for it. Instead of proceeding down the well-trodden route of using conservative hyperrectangular approaches or minimum volume ellipsoidal regions, we develop a fast algorithm for constructing hyperrectangular regions having the desired coverage, $1-\alpha$ (cf. Wasserman,, 2010; Montiel Olea and Plagborg-Møller,, 2019). Using a hyperrectangular approach makes the marginal interpretation of the resulting intervals straightforward, making it easier to incorporate in the standard plots given by existing software. The following example illustrates some of these properties.

Example 1*.*

Consider a motivating example with the goal of estimating the mean and $(.10,.90)$-quantiles for a three component mixture of normal densities. We use a Metropolis-Hastings (MH) random walk to produce 1000 MCMC samples. The sample statistics of these draws are used to estimate the 3-dimensional quantity of interest and corresponding $3\times 3$ asymptotic covariance matrix. Figure 1 shows $90\%$ simultaneous confidence intervals superimposed on a plot containing an empirical density estimate; notice that the confidence regions displayed around each quantity of interest have different lengths enabling practitioners to visualize simultaneous simulation variability and assess the reliability of the simulation experiment. It clearly illustrates both the variability of the target distribution and the variability in Monte Carlo estimation without overemphasizing point estimates. This concludes our initial discussion of this example.

The remainder is organized as follows. In Section 2, we establish the joint asymptotic distribution of sample means and quantiles for strongly mixing processes. Section 3 develops a univariate optimization method for obtaining the simultaneous confidence intervals described above. Several examples are given in Section 4, which demonstrate the versatility of the approach given here and its potential application in existing software. More specifically, we return to the three component mixture of normals considered above, then turn our attention to two Monte Carlo simulation studies, and finally consider a Bayesian analysis example. In one of the simulation studies we investigate the empirical coverage properties of the methods we propose. Some final remarks are given in Section 5.

2 Joint asymptotic distribution

Suppose $\pi$ is a probability distribution with support $\mathsf{X}\subseteq\mathbb{R}^{d},d\geq 1$. Let $p_{1}$ and $p_{2}$ be nonnegative integers such that $p_{1}+p_{2}\geq 1$. Let $g:\mathsf{X}\to\mathbb{R}^{p_{1}}$ such that $g=(g_{1},\ldots,g_{p_{1}})^{\prime}$ and, if $X\sim\pi$, then set

[TABLE]

which is assumed to exist. Collecting the expectations, set

[TABLE]

Let $h:\mathsf{X}\to\mathbb{R}^{p_{2}}$ such that $h=\left(h_{1},\ldots,h_{p_{2}}\right)^{\prime}$. If $V=h_{i}(X)$ with distribution function $F_{h_{i}}(v)$, which we assume is absolutely continuous with a continuous density $f_{h_{i}}(v)$, then define the $q_{i}$-quantile associated with $F_{h_{i}}$ as

[TABLE]

Collecting the quantiles, set

[TABLE]

Then the $(p_{1}+p_{2})$-vector $\nu=(\theta_{g}^{\prime},\phi_{h}^{\prime})^{\prime}$ is the set of features of $\pi$ that we aim to estimate.

Let $\{X_{n}\}$ be a strictly stationary stochastic process on a probability space $(\Omega,{\mathcal{F}},P)$ and set ${\mathcal{F}}_{k}^{l}=\sigma(X_{k},\ldots,X_{l})$, the $\sigma$-algebra generated by $X_{k},\ldots,X_{l}$. Define the $\alpha$-mixing coefficients for $n=1,2,3,\ldots$ as

[TABLE]

Then $\{X_{n}\}$ is strongly mixing if $\alpha(n)\to 0$ as $n\to\infty$ and this will be assumed for the remainder. Essentially, strong mixing coefficients quantify the rate at which events in the distant future become independent of the past. Both IID processes and positive Harris recurrent Markov chains (for definitions see Meyn and Tweedie,, 1993) are strongly mixing; see Bradley, (1986, 2005) and Ibragimov and Linnik, (1971).

Estimation of $\theta_{g}$ is straightforward with a vector of sample means since there is a version of the strong law (see, e.g., Blum and Hanson,, 1960; Ibragimov and Linnik,, 1971) that ensures, as $n\to\infty$, with probability 1,

[TABLE]

Similarly, estimation of $\phi_{h}$ is straightforward using sample order statistics. That is, let $\hat{\xi}_{q_{i}}=h_{i}(X)_{\lceil nq_{i}\rceil:n}$ be the $\lceil nq_{i}\rceil^{th}$ order statistic of $h_{i}(X)$ and denote the vector of $p_{2}$ estimated quantiles as

[TABLE]

Define empirical distributions for $F_{h_{i}}$ as

[TABLE]

and

[TABLE]

Letting $Q=(q_{1},\ldots,q_{p_{2}})^{\prime}$ and noting that $\bar{F}_{h_{i}}(v)$ is a sample mean we have that, as $n\to\infty$, $\bar{F}_{h}(\phi_{h})\to Q$ with probability 1. This is sufficient to ensure that, as $n\to\infty$, $\hat{\phi}_{n}\to\phi_{h}$ with probability 1 (for more discussion on this point see Doss et al.,, 2014).

While the above asymptotics justify $\hat{\nu}_{n}:=(\bar{g}_{n}^{\prime},\hat{\phi}_{n}^{\prime})^{\prime}$ as a simulation-based estimator of $\nu$, no matter how large $n$, there will be an unknown Monte Carlo error $(\hat{\nu}_{n}-\nu)$. In the following theorem, we present the conditions for an asymptotic sampling distribution of the Monte Carlo error and the explicit form of the covariance matrix.

Theorem 1**.**

Let $\{X_{n}\}$ be a stationary $\alpha$-mixing process. Suppose $F_{h_{i}}(v)$ is absolutely continuous and twice-differentiable with density $f_{h_{i}}(v)$ such that $0<f_{h_{i}}(\xi_{q_{i}})<\infty$ and the first derivative $f^{\prime}_{h_{i}}(v)$ is bounded in some neighborhood of $\xi_{q_{i}}$. In addition, suppose either

1. (a)

there exist $B$ such that $|h_{i}(x)|<B$ for all $i$ and $\alpha(n)=O(n^{-5/2-\epsilon})$ for some $\epsilon>0$  or 2. (b)

$\text{E}_{\pi}\|g\|^{2+\delta}<\infty$* for some $\delta>0$ and $\alpha(n)=O(n^{-\kappa})$ for $\kappa>\max\{3,(2+\delta)/\delta\}$ .*

Let $A_{h}$ be a $p_{2}\times p_{2}$ diagonal matrix with $i^{th}$ diagonal elements $f_{h_{i}}(\xi_{q_{i}})$ and

[TABLE]

If $Y_{j}=\left(g(X_{j}),I(h(X_{j})>\phi_{h})\right)^{\prime}$ and

[TABLE]

is positive definite, then, as $n\to\infty$,

[TABLE]

Proof.

Under condition (a) with $\alpha(n)=O(n^{-5/2-\epsilon})$, $\sum_{n}\alpha(n)<\infty$ and under condition (b) with $\alpha=O(n^{-\kappa})$, $\sum_{n}\alpha(n)^{\delta/(2+\delta)}<\infty$ for $\delta>0$. Define

[TABLE]

Then $\Sigma_{g}$ is the top left $p_{1}\times p_{1}$ principal sub-matrix of $\Sigma$, and thus is positive definite. By Ibragimov, (1962) and Jones, (2004), as $n\to\infty$,

[TABLE]

Similarly,

[TABLE]

is positive definite. Thus, due to a joint central limit theorem, there exists a $p_{1}\times p_{2}$ cross-covariance matrices $\Sigma_{gh}=\Sigma_{hg}^{\prime}$ so that

[TABLE]

Under condition (a) (Yoshihara,, 1995) or condition (b) (Zhang et al.,, 2014) there is a Bahadur quantile representation, i.e.

[TABLE]

where $r_{n,q_{i}}$ is $o_{p}(n^{-1/2})$. Letting $R_{n}=\left(r_{n,q_{1}},\dots,r_{n,q_{p_{2}}}\right)^{\prime}$, we have

[TABLE]

Combining (3) and (4) yields

[TABLE]

∎

Remark 1*.*

Wang et al., (2011) weaken the mixing conditions further, but their Bahadur quantile representation is not applicable for Metropolis-Hastings Markov chains and hence is of limited applicability in MCMC settings.

Theorem 1 can be stated more conveniently for Markov chains. Suppose $\{X_{n}\}$ is a positive Harris recurrent Markov chain with Markov transition kernel

[TABLE]

that gives the dynamics of the chain. Let $\|\cdot\|$ denote the total variation norm. Further, let $M:\mathsf{X}\mapsto\mathbb{R}^{+}$ with $E_{\pi}M<\infty$, and $\psi:\mbox{\msbm N}\mapsto\mathbb{R}^{+}$ be decreasing such that

[TABLE]

Polynomial ergodicity of order $m$ where $m>0$ means (5) holds with $\psi(n)=n^{-m}$. Geometric ergodicity means (5) holds with $\psi(n)=t^{n}$ for some $0<t<1$, and is thus stronger than polynomial ergodicity. The standard techniques (Jones and Hobert,, 2004) for establishing geometric ergodicity ensure $\text{E}_{\pi}M<\infty$, but this is more complicated in the polynomially ergodic setting; see Jones, (2004) for some discussion. It is standard (Chan and Geyer,, 1994; Doss et al.,, 2014; Jones,, 2004) that

[TABLE]

which yields the following result.

Corollary 1**.**

Suppose $F_{h_{i}}(v)$ is absolutely continuous and twice-differentiable with density $f_{h_{i}}(v)$ such that $0<f_{h_{i}}(\xi_{q_{i}})<\infty$ and $f^{\prime}_{h_{i}}(v)$ is bounded in some neighborhood of $\xi_{q_{i}}$. Let $\{X_{n}\}$ be a Harris ergodic Markov chain with stationary distribution $\pi$. If

1. (a)

there exist $B$ such that $|h_{i}(x)|<B$ for all $i$ and $\{X_{n}\}$ is polynomial ergodic of order $5/2+\epsilon$ for some $\epsilon>0$, or 2. (b)

$\{X_{n}\}$* is polynomially ergodic of order $\kappa>\max\{3,(2+\delta)/\delta\}$ where $\delta>0$ is such that $\text{E}_{\pi}\|g\|^{2+\delta}<\infty$, or* 3. (c)

there exist $B$ such that $|h_{i}(x)|<B$ for all $i$ and $\{X_{n}\}$ is geometrically ergodic, or 4. (d)

$\{X_{n}\}$* is geometrically ergodic and $E_{\pi}\|g\|^{2+\delta}<\infty$ where $\delta>0$,*

then, as $n\to\infty$, (2) holds for any initial distribution of the Markov chain.

The conclusion that the claim holds for any initial distribution, even point masses, holds by an argument similar to that of Proposition 17.1.6 of Meyn and Tweedie, (1993), which says that for positive Harris Markov chains, if a CLT holds for one initial distribution, then it holds for every initial distribution. Next, there has been significant work on establishing that Markov chains encountered in statistics are at least polynomially ergodic; see Doss and Hobert, (2010), Ekvall and Jones, (2019), Hobert and Geyer, (1998), Hobert et al., (2005), Jarner and Hansen, (2000), Jarner and Roberts, (2002), Johnson et al., (2013), Johnson and Jones, (2015), Jones et al., (2014), Jones and Hobert, (2004), Khare and Hobert, (2013), Marchev and Hobert, (2004), Roberts and Polson, (1994), Tan et al., (2013), Tan and Hobert, (2009), and Vats, (2017), among many others.

Consider the case where $\{X_{n}\}$ is an IID process from $\pi$. While it is immediate that $\{X_{n}\}$ is strongly mixing, we can do better than Theorem 1 since weaker density and moment conditions are required to establish the Bahadur quantile representation (Ghosh,, 1971). The proof the following result is omitted since it is similar to the proof of Theorem 1.

Theorem 2**.**

Suppose $\{X_{n}\}$ is an IID process from $\pi$ with $E_{\pi}\|g\|^{2}<\infty$. If $0<f_{h_{i}}(\xi_{q_{i}})<\infty$ and $f^{\prime}_{h_{i}}$ is bounded in a neighborhood of $\xi_{q_{i}}$ for all $i$, then, as $n\to\infty$, (2) holds with $\Sigma=\text{Var}_{\pi}(Y_{1})$.

Some special cases of Theorem 2 have been studied. Laplace found the joint asymptotic distribution of the sample mean and the sample median (Stigler,, 1973). Ferguson, (1998) considers a sample mean and an arbitrary quantile (also see Lin et al.,, 1980) and Babu and Rao, (1988) provides an expression for the covariance between two quantiles.

2.1 Variance estimation

Making use of Theorems 1 and 2 requires estimation of $\Lambda^{-1}\Sigma\Lambda^{-1}$. Because $\Lambda$ is diagonal with non-zero diagonals $\Lambda^{-1}$ is readily available. The only potential difficulty is that we need to estimate $f_{h_{i}}(\xi_{q_{i}})$, but it is easy to use kernel density estimators with a Gaussian kernel to estimate $f_{h_{i}}(\hat{\xi}_{q_{i}})$, and hence estimate $\Lambda$. This approach can perform well even for dependent sequences (Doss et al.,, 2014).

The matrix $\Sigma$ requires more attention since IID and strongly mixing sequences yield different structures of $\Sigma$. For IID sampling, set $S_{j}=(g(X_{j}),I(h(X_{j})>\hat{\phi}_{n}))^{\prime}$, $\Sigma$ in Theorem 2 can be consistently estimated with the sample covariance of $\{S_{1},S_{2},\dots,S_{n}\}$. For a stationary strongly mixing sequence, an expression for $\Sigma$ is given by (1) and estimating it has been studied using batch means (Chen and Seila,, 1987; Jones et al.,, 2006; Vats et al.,, 2019), weighted batch means (Liu and Flegal,, 2018), spectral variance (Andrews,, 1991; Priestley,, 1981; Vats et al.,, 2018), initial sequence (Dai and Jones,, 2017), recursive (Chan and Yau,, 2017) and regenerative sampling estimators (Seila,, 1982; Hobert et al.,, 2002).

In most simulation settings the amount of simulated data potentially available (letting $p=p_{1}+p_{2}$, then all of $d$, $n$, and $p$ can be large) forces consideration of the computational effort required to estimate $\Sigma$. Batch means has proven to be the most computationally efficient, there are conditions comparable to those in Corollary 1 for strong consistency, and it often enjoys good finite-sample properties (for more discussion on these points see Chen and Seila,, 1987; Vats et al.,, 2019). Other methods such as spectral variance estimators can also work well, especially when $n$ is not too large. We will restrict our attention to batch means estimators, which are now defined. Let $n=ab$ where $a$ is the number of batches and $b$ is the batch size. The mean for the $k^{th}$ batch is $\bar{S}_{k}(b)=b^{-1}\sum_{t=1}^{b}S_{kb+t}$ and the overall mean is $\bar{S}=a^{-1}\sum_{k=1}^{a}\bar{S}_{k}(b)$. Then the batch means estimator with batch size $b$ is

[TABLE]

We will typically use a batch size equal to $\lfloor\sqrt{n}\rfloor$. If $p$ is large, then batch means may exhibit downward bias leading to noticeable undercoverage of the simultaneous confidence intervals developed below. Other variance estimators such as weighted batch means or a lugsail window function (Vats and Flegal,, 2018) may be used to induce upward bias, but require a larger computational effort.

3 Simultaneous confidence intervals

We consider $p$-dimensional simultaneous confidence regions for $\nu$. Let $\hat{\Lambda}^{-1}\hat{\Sigma}\hat{\Lambda}^{-1}$ be a strongly consistent estimator of $\Lambda^{-1}\Sigma\Lambda^{-1}$. Then Theorem 1 allows for the construction of Wald (ellipsoidal) confidence regions having the desired asymptotic coverage, $1-\alpha$. However, visualizing a $p$-dimensional Wald ellipsoid is difficult beyond two dimensions and is difficult to use in assessing reliability of Monte Carlo experiments. Instead, we consider hyperrectangular confidence regions using the information from the full covariance structure $\Lambda^{-1}\Sigma\Lambda^{-1}$, and not just the diagonals. The proposed technique is similar to the balanced bootstrap (Beran,, 1988; Bruder and Wolf,, 2018; Romano et al.,, 2010) and the adjusted-Wald intervals (Lütkepohl et al.,, 2015) used primarily in multiple testing, but we avoid the use of computationally-intensive naive bootstrap methods.

The basic approach is to consider hyperrectangular regions $C_{LB}\subseteq C_{UB}$ where $C_{LB}$ is so small that it will have coverage no greater than $1-\alpha$ while $C_{UB}$ is so large that it will have coverage at least $1-\alpha$. There is some hyperrectangular region, say $C_{\alpha}$, between these, i.e. $C_{LB}\subseteq C_{\alpha}\subseteq C_{UB}$, which will have coverage of $1-\alpha$. In general, finding $C_{\alpha}$ will be computationally intensive, but we provide a simple solution by considering a particular subclass of regions using the results of Section 2.

Let $\hat{\sigma}_{g_{i}}$ be the $i$th diagonal of $\Sigma_{g}$, $i=1,\dots,p_{1}$ and $\hat{\gamma}_{j}$ be the $j$th diagonal of $\hat{\Lambda}^{-1}\hat{\Sigma}_{h}\hat{\Lambda}^{-1}$, $j=1,\dots,p_{2}$. For $z>0$, consider the hyperrectangular confidence regions of the form

[TABLE]

Setting $z=z_{1-\alpha/2}$ yields uncorrected intervals that simultaneously have coverage no greater than $1-\alpha$ and we set $C_{LB}:=C_{SI}(z_{1-\alpha/2})$. With $z=z_{1-\alpha/2p}$, we get the Bonferroni-corrected hyperrectangular region which has coverage at least $1-\alpha$ and we set $C_{UB}:=C_{SI}(z_{1-\alpha/2p})$. While there are many other conservative approaches, such as the S̆idák correction (Šidák,, 1967), Scheffe’s correction (Scheffe,, 1999), and Wald projections (Lütkepohl et al.,, 2015), we have found the Bonferroni approach to be an effective way to construct $C_{UB}$.

All that is left is to construct $C_{\alpha}$, but this is easily done if we can find $z^{*}$ such that $z_{1-\alpha}\leq z^{*}\leq z_{1-\alpha/2p}$ and $C_{SI}(z^{*})$ has coverage $1-\alpha$, in which case we set $C_{\alpha}=C_{SI}(z^{*})$. The approach to finding $z^{*}$ is illustrated in Figure 2, but more formally, the optimization algorithm works as follows: suppose $U\sim N_{p}(\hat{\nu}_{n},\hat{\Lambda}^{-1}\hat{\Sigma}\hat{\Lambda}^{-1}/n)$, then $\Pr(U\in C_{LB})\leq(1-\alpha)$ and $\Pr(U\in C_{UB})\geq(1-\alpha)$. Since $\Pr(U\in\mathcal{C}_{SI}(z))$ is strictly increasing as $z$ increases, we use the bisection method between $z_{1-\alpha/2}$ and $z_{1-\alpha/2p}$ to find $z^{*}$ such that $\Pr(U\in\mathcal{C}_{SI}(z^{*}))\approx(1-\alpha)$ and set $C_{\alpha}=\mathcal{C}_{SI}(z^{*})$. Thus finding the $p$-dimensional region $C_{SI}(z^{*})$ only requires optimizing over a univariate parameter, $z$. The required multivariate normal probabilities can be quickly and accurately calculated using quasi-Monte Carlo methods (Genz et al.,, 2018; Genz and Bretz,, 2009).

4 Applications

This section considers implementation of the methodology in three examples. The first one expands on the mixture model from Section 1 demonstrating use of the methods in both IID sampling and MCMC settings and investigating the finite-sample properties of the proposed confidence regions. The second application is to a setting where we compare the finite-sample properties of three frequentist statistical models, and the third demonstrates the use of our methods in a Bayesian hierarchical model. In each application, we identify a combination of means and quantiles of interest, obtain $1-\alpha$ level simultaneous confidence intervals, and integrate the intervals into a plot which is standard in the respective applications. The full R code is available as part of the supplementary material.

4.1 Mixture of normals

For $j=1,2,3$, let $f_{j}(x;\mu_{j},\sigma_{j}^{2})$ be a normal density with mean $\mu_{j}$ and variance $\sigma_{j}^{2}$. Set

[TABLE]

We consider simultaneous estimation of the mean, $.10$ quantile, and $.90$ quantile denoted $\mu$, $\xi_{.10}$, and $\xi_{.90}$, respectively, and hence $\nu=(\mu,\xi_{.10},\xi_{.90})$. Notice that we can use numerical integration methods to find that $(\mu,\xi_{.10},\xi_{.90})=(5,.2544116,11.0143117)$ with absolute error less than 2.5e-6.

We consider both standard Monte Carlo (IID samples) and MCMC metods for simulating from $\pi$. Producing IID samples is easy, simply choose a component according to its mixture probability and then simulate from it. In the MCMC setting we use a random walk Metropolis-Hastings sampler with a $N(0,9)$ proposal distribution. In either case we estimate $\nu$ and $\Lambda^{-1}\Sigma\Lambda^{-1}$ as described in Section 2. We then calculate $C_{SI}(z^{\ast})$ at the 90% confidence level and estimate the density to create Figure 3 (the lower left panel appeared in the introduction). The estimates of $\nu$ are represented by purple lines with the blue region around each estimate representing the simultaneous confidence regions. Notice that as the number of Monte Carlo samples increases, the width of the intervals decreases. Also, as should be expected, the width of the intervals based on MCMC samples is wider than those based on IID samples. Finally, note that because the shape of the density is asymmetric, the two quantiles occur at different density values. This affects the value of $\Lambda^{-1}$ and contributes to the different lengths of the error regions around each estimate.

Figure 3 concerned a single Monte Carlo simulation where every interval presented (12 total) contained the true parameter. We now turn our attention to investigating the finite-sample simultaneous coverage probabilities of the intervals $\mathcal{C}_{SI}(z^{*})$. More specifically, we do both the IID and MCMC sampling for 2000 independent replications with $n=10000$. In each replication we calculate simultaneous confidence intervals $\mathcal{C}_{SI}(z^{*})$, uncorrected marginal intervals $C_{LB}$, and simultaneous Bonferroni intervals $C_{UB}$ at the 80% and 90% confidence levels for which we record whether each region contains the true value. Thus, we have six binary outcomes, one for each region and confidence level combination, which are naturally correlated since we are using the same simulated data. We then use our procedure, based on Theorem 1, with overall 95% confidence level to plot the point estimate of the coverage along with simultaneous confidence regions in Figure 4.

From Figure 4 we can see that $C_{LB}$ yields significant undercoverage while failing to ever capture the nominal coverage probability within any of its interval estimates. For $C_{SI}(z^{\ast})$, the confidence intervals contain the nominal level illustrating that the simultaneous confidence intervals yield coverage close to the nominal level. Bonferroni intervals, $C_{UB}$, overestimate the nominal level although this overcoverage is relatively small since the adjustment is based on a small number of quantities. However, estimation procedures of higher dimensionality will correspond to more conservative estimates for the upper bound.

4.2 Comparing methods with boxplots

Finite-sample performance of statistical methodologies is often done via loss function comparisons with existing methods over repeated simulations with results presented in side-by-side boxplots. Our methodology can be used to illustrate the variability in the estimation of the quantiles in the boxplots and provides a tool to assess whether sufficient replications have been used or if observed differences are subject to substantial Monte Carlo error.

Consider a comparison of lasso (Tibshirani,, 1996), ridge (Hoerl and Kennard,, 1970), and ordinary least squares (OLS) regressions. We will sample $n\in\{100,2000\}$ independent vectors of observations $y\in\mathbb{R}^{100}$ from the following model. Let $Z$ be a $100\times 21$ dimensional matrix of covariates, and $\beta^{*}\in\mathbb{R}^{21}$ be the true regression coefficient vector. For $\epsilon\sim N_{100}(0,I_{100})$, our data generating model is

[TABLE]

We set $\beta^{*}$ to be such that the first 11 elements are zero, and the last 10 are random draws from a normal distribution with mean 0 and variance 2. The matrix $Z$ is constructed such that the first column is all 1s, and the rows of $Z_{-1}$, the matrix $Z$ with the first column removed, are drawn from $N_{20}(0,\Omega)$, where the $ij$th entry of $\Omega$ is $.90^{|i-j|}$.

For each of the simulated $y\in\mathbb{R}^{100}$ we fit lasso, ridge, and OLS regressions to estimate the vector of coefficients. Lasso and ridge estimates are obtained using the glmnet package (Friedman et al.,, 2009) with tuning parameters chosen using cross-validation. In each replication, we note the squared estimation error of the estimated coefficient, $\hat{\beta}$, that is, $\|\hat{\beta}-\beta^{*}\|^{2}$.

Figure 5 presents boxplots of the squared estimation error over replications with and without simultaneous confidence intervals. Recall a box in the boxplot has 25%, 50%, and 75% quantiles. To construct simultaneous confidence intervals, we appeal to the 9-dimensional joint asymptotic distribution for IID sequences for these quantiles. The simultaneous confidence intervals immediately indicate that with only 100 Monte Carlo replications, all quantile estimates have large variation and that the Monte Carlo error is swamping any differences between the estimated quantiles of all three methods. This variability is substantially smaller with 2000 Monte Carlo replications where we can see that the Monte Carlo error is sufficiently small so that we can claim that ridge and OLS perform similarly while the lasso is superior. Incorporating our simultaneous confidence regions allows us to assess the reliability of the simulations results in a way that is impossible with the top row in Figure 5.

4.3 Bayesian analyses

Consider Rubin,’s (1981) school data. At eight schools students were put into SAT prep classes or in a control group where they did not receive coaching. The data reported is the estimated effect of prep classes on verbal SAT scores and the standard deviation at each school. We are interested in estimating the coaching effect at each school.

A standard model (Gelman et al.,, 2004) for analyzing this data is as follows. For $j=1,\dots,J$, let

[TABLE]

[TABLE]

with $\sigma^{2}_{j}$ known and priors $f(\mu)\propto 1$ and $f(\tau)\propto 1/\tau$. We are interested in producing interval estimates of the coaching effect for each school, that is the $\theta_{j}$’s. To this end we simulate draws from the joint posterior $\theta,\mu,\tau|y$ with a deterministic scan Gibbs sampler described in the supplementary material and estimate posterior credible intervals for each $\theta_{j}$.

In Figure 6 we plot estimated marginal densities and 80% credible intervals for each of the $\theta_{j}$’s based on $n=10,000$ MCMC samples. We then use our procedure, based on Theorem 1 and Section 3, to include simultaneous 90% confidence intervals for the 16 marginal quantiles. Figure 7 presents the same analysis based on an MCMC sample of $n=100,000$. Comparing Figures 6 and 7 illustrates that accounting for the Monte Carlo error can impact our statistical conclusions about the coaching effects. Consider the left endpoints of the credible regions in the plots for $\theta_{2}$, $\theta_{4}$, and $\theta_{8}$. Specifically, in Figure 6 the left endpoints are indistinguishable from zero when we account for the Monte Carlo error, but this is no longer an issue with Figure 7.

5 Final Remarks

Simultaneous estimation of means and quantiles has received little attention in the simulation literature, despite being common practice. We provide an approach for assessing the quality of estimation for Monte Carlo sampling with a goal of ensuring reliable Monte Carlo experimental results. The approach is based on a multivariate central limit theorem for any finite combination of means and quantiles assuming the underlying process is strongly mixing at the appropriate rate. We emphasize the IID and Markov chain sampling cases. However, the proof technique used indicates that the result will hold more broadly since all that was required is a strong law for sample means, a CLT for sample means, and a Bahadur quantile representation.

Given the multivariate CLT, we provide a fast algorithm for constructing hyperrectangular confidence regions having the desired simultaneous coverage probability and a convenient marginal interpretation. The methods are then incorporated into standard ways of visualizing the results of Monte Carlo experiments enabling the practitioner to more easily assess the reliability of the results.

Acknowledgements

The authors thank Vladimir Minin for helpful conversations about assessing Monte Carlo error in Bayesian analyses. We also thank the anonymous associate editor and referees whose comments helped improve this paper.

SUPPLEMENTARY MATERIAL

PDF:

File supp_visual.pdf contains a description of the Gibbs sampler used in Section 4.3 and additional plots that were omitted to save space.

R Code:

Contains R code that reproduces the simulations and plots along with a file README.txt with directions (.zip file).