Computational Statistics in Data Science

Let denote the determinant. The volume of this ellipsoidal confidence region, which depends on , , and , is given by

(2)

Sometimes a joint sampling distribution may be difficult to obtain, or the limiting variance–covariance matrix may be too complicated to estimate. In such cases, one can consider hyperrectangular confidence regions. Let be a quantile of a standard normal distribution possibly chosen to correct for simultaneous inference. Recall that , let denote the th component of . Further, let denote the th diagonal of . Then

The volume of this hyperrectangular confidence region is

(3)

As more samples are obtained, and converge to 0 so that the variability in the estimator disappears. Sequential stopping rules in Section 5will utilize this feature to terminate simulation.

To construct confidence regions, the asymptotic variance requires estimation. For IID sampling, is estimated by the sample covariance matrix, as discussed in Section 2.3. For MCMC sampling, a rich literature of estimators of is available including spectral variance [14, 15], regeneration‐based [16, 17], and initial sequence estimators [5]18–20]. Considering the size of modern simulation output, we recommend the computationally efficient batch means estimators.

The multivariate batch means estimator considers nonoverlapping batches and constructs a sample covariance matrix from the sample mean vectors of each batch. More formally, let , where is the number of batches, and is the batch sizes. For , define . The batch means estimator of is

Univariate and multivariate batch means estimators have been studied in MCMC and operations research literature [21–26]. Although the batch means estimator has desirable asymptotic properties, it suffers from underestimation in finite samples, particularly for slowly mixing Markov chains. Specifically, let

Then, Vats and Flegal [27] show (ignoring smaller order terms)

When the autocorrelation in the Markov chain is large, or is small, there is significant underestimation in . To combat this issue, Vats and Flegal [27] propose lugsail batch means estimators formed by a linear combination of two batch means estimators with different batch sizes. For and , the lugsail batch means estimator is