{'fraction_non_alphanumeric': 0.03975155279503106, 'fraction_numerical': 0.01758893280632411, 'mean_word_length': 4.401189387008234, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 10, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A practical sequential stopping rule for high-dimensional MCMC and its application to spatial-temporal Bayesian models

A current challenge for many Bayesian analyses is determining when to terminate high-dimensional Markov chain Monte Carlo simulations. To this end, we propose using an automated sequential stopping procedure that terminates the simulation when the computational uncertainty is small relative to the posterior uncertainty. Such a stopping rule has previously been shown to work well in settings with posteriors of moderate dimension. In this paper, we illustrate its utility in high-dimensional simulations while overcoming some current computational issues. Further, we investigate the relationship between the stopping rule and effective sample size. As examples, we consider two complex Bayesian analyses on spatially and temporally correlated datasets. The first involves a dynamic space-time model on weather station data and the second a spatial variable selection model on fMRI brain imaging data. Our results show the sequential stopping rule is easy to implement, provides uncertainty estimates, and performs well in high-dimensional settings.


Introduction
Markov chain Monte Carlo (MCMC) simulations are commonly employed in a Bayesian context to estimate features of a posterior distribution by constructing a Markov chain with the target as its stationary distribution. A fundamental challenge is determining when to terminate the simulation, especially for the often high-dimensional problems encountered in modern MCMC. For instance, the visual inspection of trace plots and running means (see e.g. Flegal and Jones, 2011) is extremely challenging in high-dimensions. Further, convergence diagnostics (see e.g. Cowles and Carlin, 1996) were designed for problems of at most moderate dimension and can be essentially impossible to implement in high-dimensions.
Given these problems, most practitioners resort to using a fixed-time rule to terminate the simulation. That is, the procedure terminates after n iterations where n is determined heuristically. In this paper, we present a simple sequential stopping rule applicable for high-dimensional MCMC.
As application, we consider the analysis of large spatially and temporally datasets routinely collected by the scientific community. An unprecedented framework to effectively incorporate spatial-temporal associations is to build dependencies in different stages of Bayesian hierarchical models (Banerjee et al., 2004), which often involves the implementation of high-dimensional MCMC. There is considerable literature in this direction, for example, Huerta et al. (2004) develop a time-varying regression model for studying ozone levels; Gelfand et al. (2005) propose spatial process modeling for dynamic data with an application to climate data; Finley et al. (2012) use Gaussian predictive processes to model large space-time data; Woolrich et al. (2004) implement a fully spatio-temporal model for the noise process in fMRI data; Smith and Fahrmeir (2007) and Lee et al. (2011) develop spatial Bayesian variable selection models to brain image study.
With important economic, ecological and public health implications, these analyses require accurate assessment of their inferential uncertainties. However, few of these studies, which often involve thousands of parameters, carefully describe the stopping criterion utilized. Among them, some use convergence diagnostics (see e.g. Gelfand et al., 2005) and some report Monte Carlo standard errors (MCSEs) to assess the quality of estimates (see e.g. Lee et al., 2011). We assume the rest employ a fixed-time stopping rule where n is determined heuristically. Unfortunately, choosing too small an n can lead to inaccurate statistical inference.
As mentioned, many practitioners utilize convergence diagnostics and visual inspections to evaluate if the chain has been run long enough. Unfortunately, these methods are barely tenable in truly high-dimensional settings. For example, as stated in Gössl et al. (2001), "With this high-dimensional data, convergence diagnostics were reduced to a selection of randomly chosen parameter chains".
Instead, we advocate terminating the simulation using a fixed-width stopping rule (FWSR), which are easy to implement and theoretically justified. The main idea is to terminate the simulation when an estimate is sufficiently accurate for the scientific purpose of the analysis. That is, the simulation is terminated the first time a confidence interval width for a desired quantity is sufficiently small. Hence, the the total simulation effort will be random with these procedures.
A FWSR, first studied in MCMC settings by Jones et al. (2006), stops the simulation when the width of a confidence interval is less than a pre-specified value . Further, Jones et al. (2006) and Flegal et al. (2008) show FWSR is superior to using convergence diagnostics as stopping criteria. Flegal and Gong (2013) propose variants known as relative FWSRs that eliminate the need to specify an absolute value . Further, relative FWSRs are more practical in multivariate settings since a single can be used for multiple parameters without apriori knowledge. In particular, they advocate the use of relative standard deviation FWSR for Bayesian computations. In short, this rule terminates the simulation when the computational uncertainty of an estimate is small relative to its posterior uncertainty.
In this article, we further investigate the finite sample properties of the relative standard deviation FWSR and propose necessary modifications for applications in high-dimensional settings. The proposed modifications are driven by computational issues providing significant improvement with minor tradeoffs. As we see later, the benefits of our modifications are a significant reduction in the required computer memory and improved computational efficiency. To our best knowledge, there are no previous attempts to formally address how long to run a MCMC simulation in such high-dimensional settings. Specifically, we extend the previous application (Flegal and Gong, 2013) of the stopping rule for estimating tens of parameters to a spatial Bayesian dynamic model with hundreds of parameters and a more complicated Bayesian fMRI model with thousands of parameters. Further, we compare our results to a convergence diagnostic used as a stopping criterion and show the latter tends to terminate the simulations prematurely. We also establish a connection between the relative standard deviation FWSR and an alternative effective sample size (ESS) calculation.
The two distinct high-dimensional Bayesian hierarchical analyses considered here are (i) the univariate dynamic space-time regression models introduced by Gelfand et al. (2005) applied to weather station data collected over the northeastern United State between 2000 and 2005 (Finley et al., 2012) and (ii) the spatial variable selection models proposed by Lee et al. (2011) applied to the StarPlus fMRI datasets (Carpenter et al., 1999;Keller et al., 2001;Wang and Mitchell, 2002). Both applications clearly demonstrate the potential of the relative standard deviation FWSR in general high-dimensional settings. Moreover, they illustrate the rule is easily implemented in an almost automated fashion while providing uncertainty estimates with confidence.
The rest of the paper is organized as follows. Section 2 formally introduces the relative standard deviation FWSR, proposes modifications for modern applications, and illustrates connections with an ESS calculation. Section 3 details application to two high-dimensional MCMC inferences for Bayesian hierarchical models, where the model, experimental dataset and several comparative studies are summarized. Section 4 concludes with a discussion.
2 Sequential stopping procedure Suppose we want to make inference about a probability distribution π with support X ⊆ R d , d ≥ 1. In general, we denote θ θ θ = (θ 1 , . . . , θ p ) T ∈ R p , p ≥ 1 as a target parameter of interest with respect to π. Note that p can be smaller or larger than d, with large values of either indicating a high-dimensional setting. In this paper, we will restrict our attention to Unfortunately, in most practical settings we cannot calculate θ θ θ analytically and frequently π is such that MCMC is the only viable technique for estimating θ θ θ. The basic MCMC methods entail constructing a time-homogeneous Harris ergodic Markov chain X = X (0) , X (1) , . . . on state space X with σ-algebra B = B(X) and invariant distribution π (Robert and Casella, 2004).
Suppose we simulate X for n iterations, where n is finite. Let Z Z Z(n) := 1 n n−1 i=0 g X (i) = (Z 1 (n), . . . , Z p (n)) T be an estimator of θ θ θ from the observed chain. Under certain regularity conditions (Chan and Geyer, 1994;Jones, 2004;Roberts and Rosenthal, 2004;Tierney, 1994), we can obtain a marginal Markov chain central limit theorem (CLT) for the sampling distribution of an unknown Monte Carlo standard error, Z Z Z(n) − θ θ θ, as n → ∞ where σ 2 i ∈ (0, ∞), i = 1, . . . , p. One could also consider a multivariate Markov chain CLT for Z Z Z(n)−θ θ θ. However, the often high-dimensionality of the associated asymptotic covariance matrix creates additional challenges and extracting useful information from it is a direction of future research.
For i = 1, . . . , p, letσ i (n) denote an estimator of σ i . Then the CLT allows construction of a (1 − δ)100% marginal confidence interval for θ i , i = 1, . . . , p, with width where z δ/2 is a critical value from the standard Normal distribution. The width at (2) allows analysts to report the uncertainty in their estimates and users to assess the practical reliability. Moreover, we can use it to construct sequential FWSRs that terminates the simulation when the w i s fall below a specific value.

A relative fixed-width stopping rule
Suppose is a pre-specified value, then the simplest FWSR terminates the simulation when w i < for all i = 1, . . . , p. Asymptotic validity of such a rule was established by Glynn and Whitt (1992) and first used in MCMC simulations by Jones et al. (2006). Asymptotic validity is important because it ensures the simulation will terminate w.p.1 and the resulting confidence interval will have the right coverage probability (as → 0). Unfortunately, such a rule is difficult to implement in high-dimensional settings without apriori knowledge of the magnitudes of the components in θ θ θ. Further, a single value is unlikely to be suitable across multiple dimensions.
Instead, we advocate the use of what is known as a relative standard deviation FWSR proposed by Flegal and Gong (2013). The main idea is to terminate the simulation when an estimators computational uncertainty is small relative to its posterior uncertainty. As we will illustrate, this is equivalent to terminating the simulation when the ESS is sufficiently large.
First, we require a bit more notation. Let λ 2 i denote the posterior variance associated with θ i . That is, if an i.i.d. sample from π were available, then λ 2 i is the asymptotic variance in the CLT associated with θ i . It is important to note in general that σ 2 i = λ 2 i due to correlation in the Markov chain. For estimation of E π [g(X)], it is easy to verify that λ 2 i is the i-th diagonal element of V ar π [g(X)]. We further supposeλ 2 i (n) is an estimator of λ 2 i , usually Note that exponentiation on a vector is taken element-wise.
A relative standard deviation FWSR terminates the simulation when the length of all the confidence intervals are less than an th fraction of the magnitude of their posterior standard deviations. That is, when w i < λ i for all i = 1, . . . , p. Formally, the time at which the simulation terminates is defined by where p(n) ≥ 0. The role of p(n) is to ensure that the simulation is not terminated prematurely based on a poor estimate of σ 2 i . A reasonable default is p(n) = I(n ≤ n * ) + n −1 (Glynn and Whitt, 1992;Jones et al., 2006), where n * is the desired minimum simulation effort. The user-specified starting value n * is often based on the complexity of the problem at hand and reflects the desired accuracy for the analytical purpose.
Sufficient conditions for asymptotic validity of the relative standard deviation stopping rule are established in Flegal and Gong (2013). In short, we require three primary assumptions for i = 1, . . . , p. First, the limiting process must satisfy a functional central limit theorem (FCLT). Second, the estimator of the associated asymptotic variance must be strongly consistent, that isσ 2 i (n) → σ 2 i w.p.1 as n → ∞. Finally, the estimator of the posterior variance must be strongly consistent, that isλ 2 i (n) → λ 2 i > 0 w.p.1 as n → ∞. While not trivial, one can establish these conditions in many complex practical MCMC settings. An interested reader is directed to Flegal and Gong (2013), Flegal and Jones (2010), and Jones et al. (2006).
The relative standard deviation FWSR is appealing to Bayesian practitioners because it provides a simple, yet informative automated stopping criterion applicable in multivariate settings. First, one only needs to specify a relative and hence no prior knowledge about the magnitude of the parameters is needed. Second, a single will suffice in multivariate settings, whereas the other fixed-width approaches require a vector of values. For small p, where m is a pre-specified gap determined by an estimated simulation effort.

Variance estimation modification
The MCMC community has expended considerable effort in establishing strongly consistent estimators for the asymptotic variance at (1) including batch means (Flegal and Jones, 2010;Jones et al., 2006), spectral variance estimation (Flegal and Jones, 2010) and regenerative simulation (Hobert et al., 2002;Mykland et al., 1995). We consider only non-overlapping batch means (BM) since it is easy to implement and is feasible in high-dimensional settings.
In standard BM the output is broken into a n batches of equal size b n . Suppose the algorithm is run for a total of n = a n b n iterations and define for j = 1, . . . , a n The BM estimate of the asymptotic variance from the CLT at (1) is Jones et al. (2006) establish necessary conditions forσ 2 i (n) → σ 2 i w.p.1 as n → ∞ if the batch size and the number of batches are allowed to increase as the overall length of the simulation increases. Setting b n = n τ and a n = n/b n , the regularity conditions require that X be geometrically ergodic, E π |g| 2+ 1 + 2 < ∞ for some 1 > 0, 2 > 0 and (1 + 1 /2) −1 < τ < 1. A common choice of τ = 1/2 has been shown to work well in applications (Flegal et al., 2008;Jones et al., 2006). We denote the BM estimate with such sampling plan as the usual batch means (uBM) estimate.
Unfortunately, most sampling plans such as uBM require storage of the entire Markov chain to allow recalculations as the batch size b n grows with n. Given a target vector of dimension p, this means a matrix of size p × n will have to be stored in the memory. Clearly, computer memory soon becomes a serious issue, which one can solve by writing parts of the chain out of memory. However, given the frequency that T ( ) is checked and the already computationally intense task of updating the chain, we prefer to simplify the problem.
To this end, we suggest a new sampling plan utilizing less memory while still providing a strongly consistent BM variance estimator. Specifically, one could setb n = sup 2 k : 2 k ≤ n τ , k ∈ Z + andã n = n/b n . Noticeb n is bounded by n τ /2 ≤b n ≤ n τ . Hence, it is easy to establish strong consistency forσ 2 i (n) with such a sampling plan using results in Jones et al. (2006) and Bednorz and Latuszyński (2007). We denote this BM estimate withb n as the alternative batch means (aBM) estimate. Finally, using the new sampling plan with T ( ) requires a standard recursive recalculation ofλ i (n) as n increases. An interested reader is directed to the technique studied by Biesel (1977).
In the proposed sampling plan, we are only storing the current state of the chain and the batch means. Then in some sense, the unit of interest is per batch rather than per iteration. Thus, a natural adjustment to the frequency with which T ( ) should be checked is to examine the criterion every m batches instead of every m iterations. As before, m is pre-specified by the user but is likely much smaller than used previously. The gap between checks is then m batches, or equivalently mb n iterations. Hence, the iterations between checks automatically increases in accordance to the magnitude of the number of iterations without having to estimate the total simulation effort beforehand. In this regard, a similar increase between checks also performs well in general sampling plans.
The two proposed modifications fit naturally with each other and enable implementation of relative standard deviation FWSR in high-dimensional settings. In addition, they yield improvements in memory requirements, computational efficiency measured by clock time, and stopping procedure automation. One drawback of the new sampling plan is not applicable to estimation problems that require storing the entire chain, hence it is best served in expectation estimation.

Connections with effective sample size
Given n iterations in a Markov chain, the ESS measures the size of an i.i.d. sample with the same standard error. One way to define ESS is described in Kass et al. (1998) and Robert and Casella (2004), for i = 1, · · · , p, where ρ k (θ i ) is the autocorrelation of lag k for θ i . This calculation is implemented in many R packages, such as coda (Best et al., 1995) and mcmcse (Flegal and Hughes, 2012).
An alternative approach to define ESS is through the ratio of σ 2 i and λ 2 i . Given the previous strongly consistent estimates, we have as n → ∞ The two ESS calculations produce comparable results in various simulation studies with relatively small correlations. However, our ESS calculation is systematically smaller due to a larger variance estimate when sample is highly correlated. Moreover, the alternative ESS calculation provides additional insights into the relative standard deviation FWSR. From the definition of T ( ), one can easily show that at termination Combining (3) and (4), we have Hence, the relative standard deviation FWSR is equivalent to terminating a simulation when the smallest ESS is above pre-specified level. For instance, setting = 0.05 and δ = 0.05 produces an estimated ESS at least 6147. That is, such a simulation would provide equivalent inference as 6147 i.i.d. samples.

An alternative stopping criterion
Convergence diagnostics (for a review see Cowles and Carlin, 1996) are widely employed by practitioners as a stopping criteria. Particularly, we are interested in the Geweke diagnostic (GD) from Geweke (1992), which we will compare with relative standard deviation FWSR in the next section. Our simulations use the GD implementation from the R package coda (Best et al., 1995). The GD is based on a hypothesis test that the mean estimates of two non-overlapping parts of the Markov chain have converged. As a rule of thumb, Geweke (1992) suggested to take first 0.1 and last 0.5 proportions of the Markov chain. The resulting test statistic is univariate by its nature and the z-score is constructed as follows, Z =x 1 −x 2 ŝ 1 (0)/n 1 +ŝ 2 (0)/n 2 , wherex 1 ,x 2 are the sample average andŝ 1 (0),ŝ 2 (0) are the consistent spectral density estimates at zero frequency for the two parts of the Markov chain, respectively.
The GD requires a single Markov chain, which is close in spirit to the current work. It is also more practical in high-dimensional settings than the popular Gelman-Rubin diagnostic (Brooks and Gelman, 1998;Gelman and Rubin, 1992), which requires parallel chains. Jones et al. (2006) note that the GD is based on a Markov chain CLT and hence does not apply more generally than a FWSR that is based on the calculation of MCSE.

Applications
In this section, we evaluate the performance of the relative standard deviation FWSR and the proposed modifications using modern Bayesian applications. Particularly, we look at the spatial Bayesian dynamic models (Gelfand et al., 2005) with application to the weather station data (Finley et al., 2012) and the spatial Bayesian variable selection models (Lee et al., 2011) on an experiment fMRI dataset (Carpenter et al., 1999). The modified stopping rule is implemented to terminate the MCMC simulations for posterior inference and comparative studies are given to illustrate its several advantages.

Bayesian dynamic space-time model
This application considers the monthly temperature data collected over 356 weather stations in the northeastern Unite States starting in January 2000 to September 2011, which is available in the R package spBayes (Finley and Banerjee, 2013). We fit the univariate Bayesian dynamic space-time regression model proposed by Gelfand et al. (2005) to a subset of the dataset for illustrative purposes. Specially, we only interested in the data recorded from 10 nearby weather stations in the year 2000. Note that the modeling approach is limited to the settings where space is continuous but time is taken to be discrete.
The response y t (s) denotes the temperature at location s and time t. It is modeled through a measurement equation that provides a regression specification with a space-time varying intercept and serially and spatially uncorrelated zero-centered Gaussian disturbances as measurement error t (s). Next a transition equation introduces a p × 1 coefficient vector β β β t , which is a strictly temporal component, and a spatial-temporal component u t (s).
The overall model is given by The GP (0, C t (·, ψ t )) denotes a spatial Gaussian process with covariance function C t (·; ψ t ).
The prior specifications and MCMC schemes follow the spDynLM function in the spBayes package and we use it to sample from the Markov chain. Interested readers are directed to Finley and Banerjee (2013) for details. Particularly, we are interested in estimating the posterior mean of θ θ θ = {τ 2 t , σ 2 t , φ t , β β β t , Σ η , u t (s)}. The subset of data upon which we fit the model produces the parameters of interest θ θ θ where p = 186.
We terminated simulations under several settings and conducted comparative studies of the outputs. Specifically, the relative standard deviation FWSR is implemented using two batch mean estimates (uBM and aBM) and two values (0.1 and 0.05). Since the proposed sampling plan is applied, which leads to the batch size of the form 2 k , k ∈ Z + , we set n * = 2 14 = 16, 384 and added 20 or 21 batches between checks. The variants of the added batches is to ensure an even number of batches when checking the modified stopping rule.
The nominal level of the coverage probabilities was set to 0.95. We compare our results to the GD with a threshold p-value 0.05 checked after 15,000 iterations in a effort to confirm the chain convergence. Among the results from the FWSR, we see that uBM and aBM perform similarly given in terms of the total simulation effort and ESS. However, notice that, aBM reduces the CPU time significantly as the total simulation effort increases. We further compare the 186 estimates obtained from uBM and aBM in terms of the ratios of estimated posterior means and variances under the two settings. Figure 1 illustrates that there are no significant differences in the estimates, although one tradeoff is that the estimatesσ i (n)s from aBM are biased downward slightly. Given the advantages of aBM when p is large where the uBM may not work, we advocate its use for relative standard deviation FWSR in high-dimensional settings. In the rest of the paper, we set aBM as the default sampling plan for T ( ) if not otherwise specified.

Methods Threshold
On the other hand, the total chain length and the ESS makes it clear that GD terminates the simulation too soon. Further, one can compare GD with T ( ) using the ratios w i (n, δ)/λ i (n). Recall that w i (n, δ)s are the width of the resulting confidence intervals. These ratios serve as a measurement of the quality of our estimates. It assesses the computational uncertainty by using MCMC relative to the posterior standard deviation. Figure 2 visualizes the comparison of these ratios. It is clear the ratios from the FWSR T ( ) are significantly smaller and more concentrated than the ones from GD. Notice the cutoffs at 0.05 and 0.1 in Figure 2 are determined by values. These findings agree with Cowles and Carlin (1996) that GD tends to be premature in diagnosing convergence.

Spatial Bayesian variable selection model
This application considers the Bayesian analysis of functional Magnetic Resonance Imaging (fMRI) study. It studies the physiological changes that accompany brain activation via the blood oxygenation level dependent (BOLD) signal contrast. During the course of a typical fMRI experiment, a single patient performs a set of tasks in response to one or several external stimulus while a series of three dimensional brain images are acquired.
Our goal is to detect activated brain regions associated with external stimulus through the image intensities. Imagine that the patient's brain can be divided into tiny voxels on a  (1999). The experiment was designed to investigate brain activities related to high level cognition, i.e. language comprehension and visuospatial processing. It involved eight subjects (four males, four female) with multiple trails. Specifically, we restrict our attention to subject 04847 and trail No. 3. In this trail (approximately twenty-seven seconds), the subject was shown a picture such as ' + * ' for four seconds. Then the picture was replaced by a blank screen for another four seconds. Later, the second stimulus, i.e. a sentence such as "the star sign is above the plus sign", was presented for four seconds, or until the subject pressed a button indicating if the sentence correctly described the picture, whichever came first. A rest period of fifteen seconds was added after the sentence was removed. Snapshots were taken every 0.5 seconds resulting in about 54 images. Data were preprocessed using standard techniques such as slice timing and spatial smoothing (for a review see Lindquist, 2008) (Friston et al., 1995;Smith and Fahrmeir, 2007), Linear combination z T i a v is the baseline trend to remove stimulus-independent effects. β v is the activation amplitude and x v,i is the transformed stimulus (see Figure 3). In many experiments, the external stimulus {s i ; i = 1, . . . , t} alternates activation/inactivation in a 0-1 'boxcar' pattern. However, instead of proceeding in a 0-1 'boxcar' function, the brain produces a fairly fixed, stereotyped blood flow response with delay d v every time a stimulus hits it, where d v is estimated in a preprocessing step. The so-called hemodynamic response function (HRF) is used to characterizes this process. There are several formulations of HRF (see e.g. Friston et al., 1998;Glover, 1999;Gössl et al., 2001). One approach is to use a canonical HRF consisting of a difference of two gamma functions (Lindquist et al., 2009), where α 1 = 6, α 2 = 16, β 1 = β 2 = 1 and c = 1/6. The only unknown parameter, i.e. the amplitude A, is estimated in a preprocessing step. We can transform the orignal 'boxcar' stimulus by a convolution with the HRF, The measurement error is denoted by v,i . Appropriate distributional assumptions about v,i can be made to incorporate temporal correlation and specific priors can be chosen to reflect spatial dependence.
In this article, we apply the spatial Bayesian variable selection models for single subject proposed by Lee et al. (2011). Here we summarize their model formulation and estimation process. An interested reader is directed to their paper for more details.
Denote y y y v = (y v,1 , . . . , y v,t ) T as the BOLD image intensity at time i = 1, . . . , t for voxel v = 1, . . . , N . Let X v be a t × p design matrix of transformed stimulus and β β β v = (β v,1 , . . . , β v,p ) T be a vector of p regression coefficients for each voxel. We formulate a linear regression mode, Notice that the detection of voxel activation is equivalent to the identification of nonzero β β β v s. To this end, we introduce 0/1 binary indicators γ γ γ v = (γ v,1 , . . . , γ v,p ), v = 1, . . . , N , such that β v,j = 0 if γ v,j = 0 and β v,j = 0 if γ v,j = 1. The γ v,j is used to indicate whether the voxel v is activated by input stimulus j. Given γ γ γ v , let β β β v (γ γ γ v ) be the vector of nonzero regression coefficients and X v (γ γ γ v ) be the corresponding design matrix. Then, the model (5) can be rewritten as Further, we assume the independence among σ 2 v and set its prior π σ 2 v ∝ 1/σ 2 v . Zellner's g-prior on β β β v (γ γ γ v ) |γ γ γ v is placed to undertake variable selection or model averaging. The parameter g is adjusted to obtain similar results with those if BIC were used, Define the corresponding sum of squares for posterior inference We incorporate the temporal dependence between observations on a given voxel through the specification of the structure of Λ v . The AR(1) dependence, i.e. Λ v (i, j) = ρ |i−j| v , is an effective compromise between inferential efficacy and computational efficiency. We specify a point mass prior for ρ ρ ρ = (ρ 1 , . . . , ρ N ) at a fixed pointρ ρ ρ using maximum likelihood methods.
We incorporate the spatial dependence, as well as the anatomical information, by using a binary MRF prior, i.e. Ising prior, on γ γ γ v . Let γ γ γ (j) = (γ 1,j , . . . , γ N,j ) T be the vector of indicators for regressor j over all voxels. Then, let w v,k be pre-specified constants that weigh the interaction between voxels v and k and let ν j be parameter to measure the strength of the interaction between voxels for regressor j. We denote v ∼ k, if two voxels v and k are defined as neighbors by the user. In this article, we employ a widely used three-dimensional structure containing the six immediate neighbors: 1 above, 1 below and 4 adjacent. The weight w v,k is set to be the reciprocal of the Euclidean distance between voxel v and k.
Then, the spatial interaction is described as ν j N v=1 v∼k w v,k I (γ v,j = γ k,j ), where I(x) is the usual 0/1 indicator function. A linear "external field" N v=1 α v,j γ v,j is specified to incorporate anatomical prior information, where α v,j is chosen to reflect prior knowledge.
We rewrite the linear model (5) as where α i z z z i s are the baseline signal, β i s are the activation amplitude corresponding to the two tasks "Semantic" and "Symbol", respectively, The binary indicator γ v = {1, 1, γ v,3 , γ v,4 } is used in the variable selection problem described previously. Notice that we assume all α i s nonzero and set ν max = 1.0 as in Lee et al. (2011). Figure 4 visualize the design matrix for this linear model as we described previously.
We implemented the relative standard deviation FWSR in the simulation study. Provided the ultrahigh-dimensional nature of the posterior, memory issues forbid us to use the uBM for T ( ). Again, because of the high-dimensional nature of the problem, we set    Figure 5a and 5b are the activation maps in all eight slices for two tasks "Semantic" and "Symbol", respectively. The red voxels are identified as activated. They gather into small clusters indicating the spatial correlation.
Again, we compare terminating simulations by T ( ) and through the GD. For the GD, we generated 15,000 iterations and used a p-value 0.05 to confirm the chain convergence for each voxel. Note that p = 9398 (= 2 × 4698 + 2) versus the 15, 000 iterations in calculation.
In such a short sample, we note a small proportion of the voxels stays active or inactive throughout the sample, i.e. a sequence of constant 0 or 1. Although we did not take this scenario into the comparison, it may suggest premature termination in some senses. As before, we compare the two approaches in terms of the ratios w i (n, δ)/λ i (n). In this example, we attempted to accurately estimate the true posterior standard deviation, λ i , from a separate simulation with 10 7 iterations. Figure 6 visualizes the comparison. It is clear that the ratios from T ( ) are significantly smaller and more concentrated than the ones from GD. However, T ( ) did require a greater simulation effort (292,864) than did GD (15,000), which might be a concern given the expensive computational cost to update the chain. However, one can adjust to reach a compromise between the accuracy of estimates and the required simulation effort.

Discussion
This paper considers the relative FWSRs in the context of truly high-dimensional MCMC simulations. In our viewpoint, a practical stopping rule should achieve three properties: (1) it is easy to implement in an automated fashion with a few tuning parameters; (2) it attains confidence in resulting estimates; and (3) it is applicable in both low-and high-dimensional settings. With such properties, practitioners can then apply the stopping rule on a routine basis.
We advocate use of the modified relative standard deviation FWSR since it meets all the properties and is especially applicable in high-dimensional Bayesian settings without prior knowledge of the magnitude of the target parameters. It is controlled by one tuning parameter that measures the accuracy of the estimates. Simply put, the estimates are approximately −1 more accurate than their posterior standard deviation. Another way to understand the is through the alternative ESS calculation. Flegal and Gong (2013)  suggest to set the tuning parameter = 0.02 which leads to an ESS of 38,416 (with δ = 0.05).
They show that such a results in confidence intervals with desired coverage probabilities.
However, the value of significantly affect the total simulation effort. For instance, in the Recall the proposed sampling plan aBM summarizes information along the simulation.
On one hand, it eliminates the requirement of storing the entire chain, solves the memory issues and reduces the computational time. On the other hand, it limits the use of such method to more general estimation problems, such as quantile estimation . Certain tradeoffs are necessary in order to overcome obstacles arise from highdimensionality.
One natural extension of the stopping rule is to consider simultaneous multivariate estimation. Flegal and Gong (2013) apply a Bonferonni approach to adjust for multiplicity.
However, the standard Bonferonni approach will not work for a large dimension p, since the individual confidence interval needs to be set to a nominal level of 0.95 1/p . A sophisticated method is required to adjust for multiplicity. One direction of future research is to control the volume of a desired confidence region rather than the width of multiple confidence intervals separately.
Finally, how well a chain mixes and explores the state space is by all means an vital component in MCMC simulations. Especially in high dimensional settings, it is still a challenging problem to date. In the fMRI application, we found a small proportion of the 0-1 type chains stay at a fixed point for a large number of iterations suggesting poor mixing.
Interested readers can refer to Brooks et al. (2010) for more information.