Batch Means Method Research Articles

Estimating standard errors for importance sampling estimators with 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.

Modified open-access scheduling for new patient evaluations at an academic chronic pain clinic increased patient access to care, but did not materially reduce their mean cancellation rate: A retrospective, observational study

Study objectiveTo determine if open-access scheduling would reduce the cancellation rate for new patient evaluations in a chronic pain clinic by at least 50%. DesignRetrospective, observational study using electronic health records. SettingChronic pain clinic of an academic anesthesia department. PatientsAll patients scheduled for evaluation or follow-up appointments in the chronic pain clinic between April 1, 2014, and December 31, 2015. InterventionsOpen-access scheduling was instituted in April 2015 with appointments offered on a date of the patient's choosing ≥1 business day after calling, with no limit on the daily number of new patients. MeasurementsMean cancellation rates for new patients were compared between the 12-month baseline period prior to and for 7months after the change, following an intervening 2-month washout period. The method of batch means (by month) and the 2-sided Student t-test were used; P<0.01 required for significance. Main resultsThe new patient mean cancellation rate decreased from a baseline of 35.7% by 4.2% (95% confidence interval [CI] 1.4% to 6.9%; P=0.005); however, this failed to reach the 50% reduction target of 17.8%. Appointment lag time decreased by 4.7days (95% CI 2.3 to 7.0days, P<0.001) from 14.1days to 9.4days in the new patient group. More new patients were seen within 1week compared to baseline (50.6% versus 19.1%; P<0.0001). The mean number of new patient visits per month increased from 158.5 to 225.0 (P=0.0004). The cancellation rate and appointment lag times did not decrease for established patient visits, as expected because open-access scheduling was not implemented for this group. ConclusionsAccess to care for new chronic pain patients improved with modified open-access scheduling. However, their mean cancellation rate only decreased from 35.7% to 31.5%, making this a marginally effective strategy to reduce cancellations.

SPSTS: A sequential procedure for estimating the steady-state mean using standardized time series

ABSTRACTThis article presents SPSTS, an automated sequential procedure for computing point and Confidence-Interval (CI) estimators for the steady-state mean of a simulation-generated process subject to user-specified requirements for the CI coverage probability and relative half-length. SPSTS is the first sequential method based on Standardized Time Series (STS) area estimators of the steady-state variance parameter (i.e., the sum of covariances at all lags). Whereas its leading competitors rely on the method of batch means to remove bias due to the initial transient, estimate the variance parameter, and compute the CI, SPSTS relies on the signed areas corresponding to two orthonormal STS area variance estimators for these tasks. In successive stages of SPSTS, standard tests for normality and independence are applied to the signed areas to determine (i) the length of the warm-up period, and (ii) a batch size sufficient to ensure adequate convergence of the associated STS area variance estimators to their limiting chi-squared distributions. SPSTS's performance is compared experimentally with that of recent batch-means methods using selected test problems of varying degrees of difficulty. SPSTS performed comparatively well in terms of its average required sample size as well as the coverage and average half-length of the final CIs.

Economic and Environmental Considerations During Low Fresh Gas Flow Volatile Agent Administration After Change to a Nonreactive Carbon Dioxide Absorbent.

Reducing fresh gas flow (FGF) during general anesthesia reduces costs by decreasing the consumption of volatile anesthetics and attenuates their contribution to greenhouse gas pollution of the environment. The sevoflurane FGF recommendations in the Food and Drug Administration package insert relate to concern over potential toxicity from accumulation in the breathing circuit of compound A, a by-product of the reaction of the volatile agent with legacy carbon dioxide absorbents containing strong alkali such as sodium or potassium hydroxide. Newer, nonreactive absorbents do not produce compound A, making such restrictions moot. We evaluated 4 hypotheses for sevoflurane comparing intervals before and after converting from a legacy absorbent (soda lime) to a nonreactive absorbent (Litholyme): (1) intraoperative FGF would be reduced; (2) sevoflurane consumption per minute of volatile agent administration would be reduced; (3) cost savings due to reduced sevoflurane consumption would (modestly) exceed the incremental cost of the premium absorbent; and (4) residual wastage in discarded sevoflurane bottles would be <1%. Inspired carbon dioxide (PICO2), expired carbon dioxide, oxygen, air, and nitrous oxide FGF, inspired volatile agent concentrations (FiAgent), and liquid volatile agent consumption were extracted from our anesthesia information management system for 8 4 week intervals before and after the absorbent conversion. Anesthesia providers were notified by e-mail and announcements at Grand Rounds about the impending change and were encouraged to reduce their average intraoperative sevoflurane FGF to 1.25 L/min. Personalized e-mail reports were sent every 4 weeks throughout the study period regarding the average intraoperative FGF (i.e., from surgery begin to surgery end) for each agent. Batch means methods were used to compare FGF, volatile agent consumption, net cost savings, and residual sevoflurane left in bottles to be discarded in the trash after filling vaporizers. The time from reaching a PICO2 = 3 mm Hg for 3 minutes until agent exhaustion (PICO2 = 5 mm Hg for 5 minutes) was evaluated. A total of N = 20,235 cases were analyzed (80.2% sevoflurane, 15.1% desflurane, and 4.7% isoflurane). Intraoperative FGF was reduced for cases in which sevoflurane was administered by 435 mL/min (95% confidence interval [CI], 391 to 479 mL/min; P < 10). Hypothesis 1 was accepted. Sevoflurane consumption per minute of administration decreased by 0.039 mL/min (95% CI, 0.029 to 0.049 mL/min; P < 10) after the change to the nonreactive absorbent. Hypothesis 2 was accepted. The difference in mean cost for the sum of the sevoflurane and absorbent purchases for each of the 10 4-week intervals before and after the absorbent switch was -$293 per 4-week interval (95% CI, -$2853 to $2266; P = 0.81). Hypothesis 3 was rejected. The average amount of residual sevoflurane per bottle was 0.67 ± 0.06 mL (95% CI, 0.54 to 0.81 mL per bottle; P < 10 vs 2.5 mL). Hypothesis 4 was accepted. Once the PICO2 reached 3 mm Hg for at least 3 consecutive minutes, the absorbent became exhausted within 95 minutes in most (i.e., >50%) canisters. We showed that an anesthesia department can transition to a premium, nonreactive carbon dioxide absorbent in a manner that is at least cost neutral by reducing FGF below the lower flow limits recommended in the sevoflurane package insert. This was achieved, in part, by electronically monitoring PICO2, automatically notifying the anesthesia technicians when to change the absorbent, and by providing personalized feedback via e-mail to the anesthesia providers.

Robust Statistical Error Estimation of Local Power Tallies in Monte Carlo Calculation of Light Water Reactor

The overlapping batch means method (OBM) has been investigated for robust statistical error estimation of local power tallies in Monte Carlo (MC) reactor core calculation. Originally, a nonoverlapping version was introduced in MC criticality calculation by Gelbard and Prael. However, the issue of batch size optimization was thought of as a lack of robustness. In this work, OBM with asymptotic bias correction was implemented with the batch size of the square root of the number of generations and compared with the orthonormally weighted standardized time series method (OWSTS). Numerical tests were conducted for various positions of the core of a pressurized water reactor. Results obtained indicate that neither OBM nor OWSTS consistently outperforms the other in terms of an overall performance measure incorporating bias and stability. Therefore, OBM with asymptotic bias correction can be an option to statistical error estimation in production MC criticality codes since OWSTS lacks an automated process to determine the number of weighting functions and can output the estimate only at the final generation. It is also shown that OBM with asymptotic bias correction performs equally regardless of the batch size.

Statistical Approaches Applied on a Greek Lip Cancer Dataset

Having at our disposal a dataset of 186 lip cancer cases in Greece, we attempt to interpret them by applying four different statistical methods: Generalized Linear Model (GLM), Markov Chain Monte Carlo (MCMC), Cox proportional hazards model and Bayes factor. The likelihood score equations from GLM exerted estimators with bounded influence, so that the resulting estimators were robust against outliers while maintaining high efficiency in the absence of outliers. Batch means method of estimating the variance of the asymptotic normal distribution, used in MCMC, gave strong consistency when it was applied to our data. A Cox proportional hazards model done with a weighted expectation-maximization gave efficient parameter estimates. Finally, using Bayes factor to the prior distributions for the parameters in compared regression models was proved to be highly sensitive.

Statistical Analysis with Little's Law

The theory supporting Little's Law (L = λW) is now well developed, applying to both limits of averages and expected values of stationary distributions, but applications of Little's Law with actual system data involve measurements over a finite-time interval, which are neither of these. We advocate taking a statistical approach with such measurements. We investigate how estimates of L and λ can be used to estimate W when the waiting times are not observed. We advocate estimating confidence intervals. Given a single sample-path segment, we suggest estimating confidence intervals using the method of batch means, as is often done in stochastic simulation output analysis. We show how to estimate and remove bias due to interval edge effects when the system does not begin and end empty. We illustrate the methods with data from a call center and simulation experiments.

Steady-State Simulation with Replication-Dependent Initial Transients: Analysis and Examples

The replicated batch means (RBM) method for steady-state simulation output analysis generalizes both the independent replications (IR) and batch means (BM) methods. We analyze the performance of RBM in situations where the underlying stochastic process possesses an additive initial transient. Our analysis differs from prior work in that the initial transient is stochastic, and hence the sample paths of the transient process may be replication dependent, and possibly also correlated across replications. We provide asymptotic expressions for the mean and variance of the RBM estimators of the steady-state mean and variance parameter of the stochastic process being simulated. We then use our results to study the performance of RBM as a function of the number of replications, initialization method for the replications, and decay rate of the associated initialization bias. Our results provide guidance on when IR, BM, or a combination thereof is the best choice, and also on effective choices of initial states for the replications.

A variance-estimation-based stopping rule for symbolic dynamic filtering

As an alternative to the batch means (BM) method in the stopping rule for symbolic dynamic filtering, this short paper presents an analytical procedure to estimate the variance parameter and to obtain a lower bound on the length of symbol blocks for constructing probabilistic finite state automata (PFSA). If the modulus of the second largest eigenvalue of the PFSA’s state transition matrix is relatively small or if the symbol block length is not too large, then the performance of the proposed stopping rule is superior to that of the stopping rule based on BM method. The algorithm of the proposed stopping rule is validated on ultrasonic data collected from a fatigue test apparatus for damage detection in the polycrystalline alloy 7075-T6.

N-Skart: A Nonsequential Skewness- and Autoregression-Adjusted Batch-Means Procedure for Simulation Analysis

We discuss N-Skart, a nonsequential procedure designed to deliver a confidence interval (CI) for the steady-state mean of a simulation output process when the user supplies a single simulation-generated time series of arbitrary size and specifies the required coverage probability for a CI based on that data set. N-Skart is a variant of the method of batch means that exploits separate adjustments to the half-length of the CI so as to account for the effects on the distribution of the underlying Student's t-statistic that arise from skewness (nonnormality) and autocorrelation of the batch means. If the sample size is sufficiently large, then N-Skart delivers not only a CI but also a point estimator for the steady-state mean that is approximately free of initialization bias. In an experimental performance evaluation involving a wide range of test processes and sample sizes, N-Skart exhibited close conformance to the user-specified CI coverage probabilities.

Batch means and spectral variance estimators in Markov chain Monte Carlo

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.

On the validity of the batch quantile method for Markov chains

We discuss the asymptotic validity of confidence intervals for quantiles of performance variables when simulating a Markov chain. We show that a batch quantile methodology (similar to the batch means method) can be applied to obtain confidence intervals that are asymptotically valid under mild assumptions.

Gibbs sampling for a Bayesian hierarchical general linear model

We consider a Bayesian hierarchical version of the normal theory general linear model which is practically relevant in the sense that it is general enough to have many applications and it is not straightforward to sample directly from the corresponding posterior distribution. Thus we study a block Gibbs sampler that has the posterior as its invariant distribution. In particular, we establish that the Gibbs sampler converges at a geometric rate. This allows us to establish conditions for a central limit theorem for the ergodic averages used to estimate features of the posterior. Geometric ergodicity is also a key requirement for using batch means methods to consistently estimate the variance of the asymptotic normal distribution. Together, our results give practitioners the tools to be as confident in inferences based on the observations from the Gibbs sampler as they would be with inferences based on random samples from the posterior. Our theoretical results are illustrated with an application to data on the cost of health plans issued by health maintenance organizations.

Analysing steady-state simulation output using vector autoregressive processes with exogenous variables

A simulation study often requires computation of a point estimate and confidence region for the steady-state mean of a stochastic output process. The literature offers a variety of statistical techniques, including replication/deletion, the batch-means method, and spectrum analysis. We present a new multivariate output-analysis technique that is based on the general autoregressive time-series model with exogenous variables to set up a joint confidence region for the steady-state mean. We demonstrate our technique by an extensive computational experiment, and show that it performs at least as well as other output-analysis techniques, without having some of their drawbacks.

Simulation output analysis using integrated paths

This article considers the steady-state simulation output analysis problem for a process that satisfies a functional central limit theorem. We construct an estimator for the time-average variance constant that is based on iterated integrations of the sample path. When the observations are batched, the method generalizes the method of batch means. One advantage of the method is that it can be used without batching the observations; that is, it can allow for the process variance to be estimated at any time as the simulation runs without waiting for a fixed time horizon to complete. When used in conjunction with batching, the method can improve efficiency (the reciprocal of work times mean-squared error) compared with the standard method of batch means. In numerical experiments, efficiency improvement ranged from a factor of 1.5 (for the waiting time sequence in an M/M/1 queueing system with a single integrated path) up to a factor of 14 (for an autoregressive process and 19 integrated paths).

A statistical process control approach to selecting a warm-up period for a discrete-event simulation

The selection of a warm-up period for a discrete-event simulation continues to be problematic. A variety of selection methods have been devised, and are briefly reviewed. It is apparent that no one method can be recommended above any other. A new approach, based upon the principles of statistical process control, is described (SPC method). Because simulation output data are often highly autocorrelated and potentially non-normal, the batch means method is employed in constructing the control chart. The SPC method is tested on seven data sets and encouraging results are obtained concerning its accuracy. The approach is also discussed with respect to its ease of implementation, simplicity, generality of use and requirements for parameter estimation.

Fixed-Width Output Analysis for Markov Chain Monte Carlo

Markov chain Monte Carlo is a method of producing a correlated sample to estimate features of a target distribution through ergodic averages. A fundamental question is when sampling should stop; that is, at what point the ergodic averages are good estimates of the desired quantities. We consider a method that stops the simulation when the width of a confidence interval based on an ergodic average is less than a user-specified value. Hence calculating a Monte Carlo standard error is a critical step in assessing the simulation output. We consider the regenerative simulation and batch means methods of estimating the variance of the asymptotic normal distribution. We give sufficient conditions for the strong consistency of both methods and investigate their finite-sample properties in various examples.

A principal-components approach to assign confidence intervals in steady-state simulation

This paper presents an approach to the assignment of a confidence interval to the mean of a stream of autocorrelated output data from a steady-state simulation run. Based on the principal-components analysis method, the approach is to derive a linear transformation of the data that yields approximate independence of the transformed data. To aid convergence to normality (on which the confidence interval is based) and to keep the dimension of the transformation reasonable, the original output data are batched prior to performing the transformation. The approach is fairly simple to understand and to implement, and experimental results indicate that it may perform better than the non-overlapping batch means method in terms of coverage with any number of batches.

Replicated batch means for steady-state simulations

This paper studies a new steady-state simulation output analysis method called replicated batch means in which a small number of replications are conducted and the observations in these replications are grouped into batches. This paper also introduces and compares methods for selecting the initial state of each replication. More specifically, we show that confidence intervals constructed by the replicated batch means method are valid for large batch sizes and derive expressions for the expected values and variances of the steady-state mean and variance estimators for stationary processes and large sample sizes. We then use these expressions, analytical examples, and numerical experiments to compare the replicated batch means method with the standard batch means and multiple replications methods. The numerical results, which are obtained from an AR(1) process and a small, nearly-decomposable Markov chain, show that the multiple replications method often gives confidence intervals with poorer coverage than the standard and replicated batch means methods and that the replicated batch means method, implemented with good choices of initialization method and number of replications, provides confidence interval coverages that range from being comparable with to being noticeably better than coverages obtained by the standard batch means method. © 2006 Wiley Periodicals, Inc. Naval Research Logistics, 2006

Estimation Methods for Delays in Non-regenerative Discrete-Event Systems

Many long-run delay characteristics of computer, communication, manufacturing, and transportation systems can be specified as time-average limits of sequences of delays in generalized semi-Markov processes (GSMPs). We consider sequences of delays that are determined from the state transitions of the GSMP using the method of start vectors. In this setting, time-average limits typically must be estimated using simulation. Previous work on estimation methods for delays has focused on GSMPs for which there exists a sequence of regeneration points. For such systems, it is often possible to find an explicit sequence of regeneration points or od-regeneration points for the sequence of delays, so that point estimates and confidence intervals for time-average limits can be obtained using the regenerative method for simulation output analysis or one of its variants. This paper is concerned with GSMPS for which this approach is not feasible, either because regeneration points for the GSMP cannot be identified or because regenerations occur too infrequently. We provide conditions on the building blocks of a GSMP and start-vector mechanism under which the sequence of delays is an od-regenerative process and cycle sums have finite moments—these conditions do not require that there exist regeneration points for the GSMP. Although in our setting the od-regeneration points for the sequence of delays usually cannot be determined explicitly, the mere existence of these points implies that time-average limits are well defined and the sequence of delays obeys a multivariate functional central limit theorem. It then follows from results of Glynn, Iglehart, and Muñoz that methods based on standardized time series can be used to obtain strongly consistent point estimates and asymptotic confidence intervals for time-average limits and functions of time-average limits. In particular, the method of batch means is applicable.