Estimating Monte Carlo variance from multiple Markov chains

The naive importance sampling estimator based on the samples from a single importance density can be extremely numerically unstable. We consider multiple distributions importance sampling estimators where samples from more than one probability distributions are combined to consistently estimate means with respect to given target distributions. These generalized importance sampling estimators provide more stable estimators than the naive importance sampling estimators. Importance sampling estimators can also be used in the Markov chain Monte Carlo (MCMC) context, that is, where iid samples are replaced with positive Harris Markov chains with invariant importance distributions. If these Markov chains converge to their respective target distributions at a geometric rate, then under two finite moment conditions a central limit theorem (CLT) holds for the importance sampling estimators. In order to calculate valid asymptotic standard errors, it is required to consistently estimate the asymptotic variance in the CLT. Recently Tan and Doss and Hobert (2015) developed an approach based on regenerative simulation for obtaining consistent estimators of the asymptotic variance. It is well-known that in practice it is often difficult to construct a useful minorization condition that is required in Tan and Doss and Hobert ’s (2015) regenerative simulation method. We provide an alternative estimator for these standard errors based on the easy to implement batch means methods. The multi-chain importance sampling estimators depend on Geyer’s (1994) reverse logistic estimator (of ratios of normalizing constants) which has wide applications, in its own right, in both frequentist and Bayesian inference. We also provide batch means estimator for calculating asymptotically valid standard errors of Geyer’s (1994) reverse logistic estimator. We illustrate the method with an application in Bayesian variable selection in linear regression. In particular, the multi-chain importance sampling estimator is used to perform empirical Bayes variable selection and the batch means estimator is used to obtain standard errors in the large p situation where regenerative method is not applicable.

We present a new method for obtaining confidence intervals in steady-state simulation. In our replicated batch means method, we perform a small number of independent replications to estimate the steady-state mean of the underlying stochastic process. In order to obtain a variance estimator, we further group the observations from these replications into nonoverlapping batches. We show that for large sample sizes, the new variance estimator is less biased than the batch means variance estimator, the variances of the two variance estimators are approximately equal, and the new steady-state mean estimator has a smaller variance than the batch means estimator when there is positive serial correlation between the observations. For small sample sizes, we compare our replicated batch means method with the (standard) batch means and multiple replications methods empirically, and show that the best overall coverage of confidence intervals is obtained by the replicated batch means method with a small number of replications.

Estimating the variance of the sample mean is a classical problem of steady-state simulation output analysis. Traditional batch means estimators require specification of the simulation run length a priori. To our knowledge, the Dynamic Non-overlapping Batch Means (DNBM) estimator is the only existing variance estimator that requires a constant storage space for any sample size. In this paper, we develop the Dynamic Partial-overlapping Batch Means (DPBM) algorithm, that also requires a constant storage space. In terms of the mean squared error, the statistical performance of the DPBM estimators is superior to that of the DNBM estimators.

The Song rule outperforms optimal-batch-size variance estimators in simulation output analysis

A classic problem of stochastic simulation is estimating the variance of point estimators, the prototype estimator being the sample mean from a steady-state autocorrelated process. The traditional batch means (BM) estimator requires knowledge of the sample size a priori. This paper proposes an algorithm to implement certain BM estimators without knowing the sample size in advance. The proposed algorithm is useful when the run length is random or is extremely long in simulation models.

Calculating a Monte Carlo standard error (MCSE) is an important step in the statistical analysis of the simulation output obtained from a Markov chain Monte Carlo experiment. An MCSE is usually based on an estimate of the variance of the asymptotic normal distribution. We consider spectral and batch means methods for estimating this variance. In particular, we establish conditions which guarantee that these estimators are strongly consistent as the simulation effort increases. In addition, for the batch means and overlapping batch means methods we establish conditions ensuring consistency in the mean-square sense which in turn allows us to calculate the optimal batch size up to a constant of proportionality. Finally, we examine the empirical finite-sample properties of spectral variance and batch means estimators and provide recommendations for practitioners.

Studies a modification of Y. Rinott's (1978) two-stage procedure for selecting the normal population with the largest (or smallest) mean. The modification, which is appropriate for use in the simulation environment, uses, in the procedure's first stage, different variance estimators than the usual batch means (BM) variance estimator. In particular, we use variance estimators arising from the method of standardized time series (STS). On the plus side, certain STS estimators have more degrees of freedom than the BM estimator does. On the other hand, STS variance estimators tend to require larger sample sizes than the BM estimator in order to converge to their assumed distributions. These considerations result in trade-offs involving the procedure's achieved probability of correct selection as well as the procedure's expected sample size.

Estimating the variance of the sample mean is a classical problem of stochastic simulation. Traditional batch means estimators require specification of the simulation run length a priori. Dynamic batch means (DBM) is a new approach to implement the traditional batch means in fixed memory by dynamically changing both batch size and number of batches without the knowledge of the simulation run length. This article further improves the DBM by considering small storage requirements and fast computation. The proposed algorithm is useful when the simulation run length is random and extremely long in simulation models.

Run length not required: Optimal-mse dynamic batch means estimators for steady-state simulations

Spaced batch means

Estimating the variance of the sample mean is a classical problem of stochastic simulation. Traditional batch means estimators require specification of the simulation run length a priori. To our knowledge, the dynamic non-overlapping batch means estimator (DNBM) and dynamic partial-overlapping batch means estimator (DPBM) are the only two existing variance estimators requiring a constant storage space for any sample size. The performance of the DPBM is better than that of DNBM in terms of the mse criteria, but the DPBM requires four times more memory than the DNBM. This paper improves the DPBM by developing a computational version of the DPBM that requires the same storage space as the DNBM estimator.

Estimating the variance of the sample mean is a classical problem of stochastic simulation. Traditional batch means estimators require specification of the simulation run length a priori. To our knowledge, the dynamic non-overlapping batch means estimator (DNBM) and dynamic partial-overlapping batch means estimator (DPBM) are the only two existing variance estimators requiring a constant storage space for any sample size. The performance of the DPBM is better than that of DNBM in terms of the mse criteria, but the DPBM requires four times more memory than the DNBM. This paper improves the DPBM by developing a computational version of the DPBM that requires the same storage space as the DNBM estimator.

When using a batch means methodology for estimation of a nonlinear function of a steady-state mean from the output of simulation experiments, it has been shown that a jackknife estimator may reduce the bias and mean squared error (mse) compared to the classical estimator, whereas the average of the classical estimators from the batches (the batch means estimator) has a worse performance from the point of view of bias and mse. In this paper we show that, under reasonable assumptions, the performance of the jackknife, classical and batch means estimators for the estimation of quantiles of the steady-state distribution exhibit similar properties as in the case of the estimation of a nonlinear function of a steady-state mean. We present some experimental results from the simulation of the waiting time in queue for an M/M/1 system under heavy traffic.

The use of Bayesian inference for the analysis of complex statistical models has increased dramatically in recent years, in part due to the increasing availability of computing power. There are a range of techniques available for carrying out Bayesian inference, but the lack of analytic tractability for the vast majority of models of interest means that most of the techniques are numeric, and many are computationally demanding. Indeed, for high-dimensional non-linear models, the only practical methods for analysis are based on Markov chain Monte Carlo (MCMC) techniques, and these are notoriously compute intensive, with some analyses requiring weeks of CPU time on powerful computers. It is clear therefore that the use of parallel computing technology in the context of Bayesian computation is of great interest to many who analyse complex models using Bayesian techniques. Section 18.2 considers the key elements of Bayesian inference, and the notion of graphical representation of the conditional independence structure underlying a statistical model. This turns out to be key to exploiting partitioning of computation in a parallel environment. Section 18.3 looks at the issues surrounding Monte Carlo simulation techniques in a parallel environment, laying the foundations for the examination of parallel MCMC in Section 18.4. Standard pseudo-random number generators are not suitable for use in a parallel setting, so this section examines the underlying reasons and the solution to the problem provided by parallel pseudo-random number generators. Parallel MCMC is the topic of Section 18.4. There are two essentially different strategies which can be used for parallelising an MCMC scheme (though these may be combined in a variety of ways). One is based on running multiple MCMC chains in parallel and the other is based on parallelisation of a single MCMC chain. There are different issues related to the different strategies, and each is appropriate in different situations. Indeed, since MCMC in complex models is somewhat of an art-form anyway, with a range of different possible algorithms and trade-offs even in the context of a non-parallel computing environment, the use of a parallel computer adds an additional layer of complexity to the MCMC algorithm design process. That is, the trade-offs one would adopt for the design of an efficient MCMC algorithm for the analysis of a given statistical algorithm in a non-parallel environment may

Confronting uncertainty in model-based geostatistics using Markov Chain Monte Carlo simulation