# Stationary Bootstrap:
A Refined Error Estimation for Equilibrium Time Series

Yoshihiko Nishikawa  Laboratoire Charles Coulomb (L2C), Université de
Montpellier, CNRS, 34095 Montpellier, France  Graduate School of Information
Sciences, Tohoku University, Sendai 980-8579, Japan  Jun Takahashi  Institute
of Physics, Chinese Academy of Sciences, Beijing 100190, China  Department of
Physics, Boston University, 590 Commonwealth Avenue, Boston, Massachusetts
02215, USA  Center for Quantum Information and Control, University of New
Mexico, Albuquerque, New Mexico 87131, USA  Takashi Takahashi  Institute for
Physics of Intelligence
Department of Physics Graduate School of Science
The University of Tokyo 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-0033, Japan

###### Abstract

In Markov-chain Monte Carlo simulations, estimating statistical errors or
confidence intervals of numerically obtained values is an essential task. In
this paper, we review several methods for error estimation, such as simple
empirical estimation with multiple independent runs, the blocking method, and
the stationary bootstrap method. We then study their performance when applied
to an actual Monte-Carlo time series. We find that the stationary bootstrap
method gives a reasonable and stable estimation for any quantity using only
one single time series. In contrast, the simple estimation with few
independent runs can be demonstratively erroneous. We further discuss the
potential use of the stationary bootstrap method in numerical simulations.

##  I Introduction

The Monte Carlo method is a general, versatile tool to solve an extremely wide
range of problems in science, including condensed-matter and high-energy
physics Landau and Binder (2014); Berg (2004); Krauth (2006); Tanabashi  _et
al._ (2018), biology Manly (2018), and machine learning Andrieu  _et al._
(2003). The progress of modern computers allows us to perform Monte Carlo
calculation of statistical problems with millions of degrees of freedom within
a reasonable amount of time and estimate various quantities with great
precision. However, because Monte Carlo method estimates quantities from a
finite number of samples, the results inevitably have statistical
uncertainties and fluctuate from simulation to simulation. Accurately
estimating statistical errors, or the confidence intervals, is thus as
essential as computing the quantities themselves in interest. Note that the
confidence interval is the fluctuation regarding the sampling distribution,
not to be confused with the physically relevant distribution of the observable
(see Sec. II for more discussions). Assessing the fluctuation of the sampling
distribution is what is needed for estimating the statistical errors.

When dealing with samples that are independent and identically distributed
(i.i.d.), we know the sampling distribution will converge to a Gaussian
distribution thanks to the central limit theorem. This allows us to use the
standard error of many samples to construct the $68\%$ confidence interval
quite reliably. More sophisticated resampling methods, including the jackknife
and the bootstrap methods Miller (1974); Efron (1979); Efron and Tibshirani
(1994); Wasserman (2004, 2006), have been proposed and rigorous studies have
proven asymptotic exactness in the limit of infinite number of samples for
some setups. Those methods virtually enhance the number of nearly-independent
data, enabling more accurate estimates of the statistical error. With the
advances of these resampling methods, the statistical analysis for i.i.d.
samples have now been established well.

For most of the statistical problems in computational physics, we use the
Markov-chain Monte Carlo (MCMC) method to produce samples. Such samples often
have correlation in the time direction. This correlation complicates their
statistical analysis and the statistical technique for i.i.d. samples
mentioned above is unavailable in such cases. The most straight-forward and
commonly used way to deal with the statistical correlations is to have several
totally independent runs and then calculate the standard error among them as
the estimated statistical error. Heuristic methods such as the blocking method
Flyvbjerg and Petersen (1989); Newman and Barkema (1999); Krauth (2006);
Jonsson (2018) has also been used by physicists.

On the other hand, in the field of modern computational statistics, various
frameworks for estimating the statistics of time-correlated samples, with no
assumption on its sampling distribution have been proposed Carlstein (1986);
Künsch (1989); Liu and Singh (1992); Politis and Romano (1992, 1994a). Despite
these advances, the extensive knowledge in statistics has not been fully
imported to the field of computational physics. We aim to shed light on this
topic by quantitatively analyzing how the modern method from statistics
compares to the traditional approach used commonly in computational physics.

Figure 1: The sampling distribution of the susceptibility $\chi/N$ of actual
Monte Carlo data of the two-dimensional Ising model (see Sec. IV for the
details of our numerical simulation), with the time-series length $n=10^{3}$
(a) and $10^{5}$ (b), respectively. Each inset shows a semi-log plot of each
panel.

In this paper, we focus on the error estimation of stationary time-series data
with a finite correlation time generated by MCMC methods. Applying several
statistical methods to numerically obtained samples, we illustrate how each
method estimates the confidence interval for various number of samples, and
contrast their performances in a realistic setup. In order to be specific,
here we define the statistical errors as confidence intervals of the computed
physical quantities, which indicate high probability of the true value
existing within the interval. This “high probability” is usually set to
$\mathrm{erf}(1/\sqrt{2})\approx 0.68$ corresponding to the value of the
standard error of a Gaussian distribution 111In experimental fields, usually a
95% bar is set.. The standard error indeed gives a reliable estimation of the
$68\%$ confidence interval when the number of samples is large enough and the
sampling distribution is well-approximated by a Gaussian distribution. When
the sampling distribution is strongly non-Gaussian (see Fig. 1), however, the
straight-forwardly calculated standard error will tend to underestimate (be
overconfident in) the confidence interval.

Perhaps the simplest method of error estimation used frequently in the MCMC
community is calculating the standard error from several independent runs
Krauth (2006); Wang (2015); Nishikawa and Hukushima (2016); Anderson  _et al._
(2017). Although this method is straightforward and easy to implement, it
obviously demands more computational resources with a multiplicative factor of
the number of independent runs. In some cases where computational physicists
attempt to reach the state-of-the-art limit with very large system sizes etc.,
it is not the most efficient way. Importantly, when the number of independent
runs is small, this method would underestimate the statistical error.

Another well-used technique is the “blocking” (or “binning”) method Flyvbjerg
and Petersen (1989); Newman and Barkema (1999); Krauth (2006); Jonsson (2018)
where the time series is divided into $K\gg 1$ “blocks” of certain fixed time
periods. Each block can yield a rather noisy estimate of the desired quantity,
and the statistical error is calculated from the standard error among those
$K$ values. By combining two consecutive blocks as one data, we can have
another estimate of the error from the $K/2$ combined “double blocks”. We can
go on with more combining steps, which could be interpreted as renormalization
in the time direction Flyvbjerg and Petersen (1989). The expectation is that
the estimated error will converge to the “correct” value as long as the
renormalization steps $k$ are not too large (ending up giving only
$K/2^{k}\rightarrow O(1)$ data). However, when the time-series length is
finite, it is hard to systematically determine whether the number of blocks
$K$ and renormalization steps $k$ are appropriate. While this method can give
reasonably good estimation of error in practice, it is difficult in general to
find if the estimation is precise since such “convergent renormalization
steps” $k_{\text{opt}}$ may not even exist for a particular finite data set.

The paper is organized as follows. In Sec. II, we summarize some basic
concepts in statistics and comment on the meaning of statistical errors in
general. We also remind the difficulty of evaluating statistical errors in
numerical simulations of physical systems. In Sec. III, we review the
bootstrap method proposed by B. Efron Efron (1979); Efron and Tibshirani
(1994). We then introduce the stationary bootstrap method in Sec. III.2 and
its practical algorithm. In Sec. IV, we detail our numerical procedures for
producing time series and estimating the errors. We then discuss the
performance of the stationary bootstrap method, together with comparison to
the error estimation using independent runs and the blocking method in Sec. V.
We conclude in the last section Sec. VI, with discussing practical guidelines
of the stationary bootstrap method for practical situations.

##  II Error estimation in sampling methods

In MCMC, a time series $x_{1},x_{2},\cdots,x_{n}$ drawn from a distribution
$q$ is produced, where $x_{t}$ is a physical quantity of interest, such as the
magnetization of the Ising model or the microscopic configuration of the
physical system itself, at time step $t$.

Suppose we are interested in estimating a physical quantity $A$. The physical
quantity $A$ is evaluated from the time series as

$\bar{A}_{n}=\bar{A}_{n}(x_{1},x_{2},\cdots,x_{n}).$ (1)

We expect that $\bar{A}_{n}$ converges to $A$ in the large sample limit
$n\to\infty$. However, this estimator $\bar{A}_{n}$ depends on the realization
of the time series, e.g. the seed of a random number generator etc., which
motivates us to evaluate its statistical error.

Again, we stress that the estimation of the confidence interval is about the
sampling distribution

$P_{n}(a)=\left\langle\delta(a-\bar{A}_{n})\right\rangle,$ (2)

where $\delta$ is Dirac’s delta function. This distribution is different from
the distribution of a physical quantity, such as the magnetization and the
energy density. Here, we denote by $\left\langle\cdots\right\rangle$ an
average over the realizations of time series,

$\left\langle\cdots\right\rangle=\int(\cdots)q(x_{1},x_{2},\cdots,x_{n})dx_{1}dx_{2}\cdots
dx_{n}.$ (3)

The sampling distribution clearly depends on $n$ and it converges to a delta
function at $n\to\infty$. On the other hand, the physical distribution is
independent of $n$ and has a non-zero variance in general.

As an example, we consider the case where $x_{i}$ $(i=1,\cdots n)$ is a
physical quantity that is defined for each microscopic configuration of the
system (i.e. sample), such as the magnetization and the energy density. Then
the distribution of that quantity could be estimated as

$P(x)=\left\langle\frac{1}{n}\sum_{i}\delta(x-x_{i})\right\rangle.$ (4)

This converges to a physical distribution (e.g. the Boltzmann distribution),
and is thus a physically measurable quantity. It has a non-zero variance even
in the limit $n\to\infty$ and its fluctuation has an inherent physical
meaning; e.g. the specific heat when $x_{i}$ is the energy density. On the
other hand, the sampling distribution of the sample mean $\sum_{i}x_{i}/n$,

$P_{n}(a)=\left\langle\delta\left(a-\frac{1}{n}\sum_{i}x_{i}\right)\right\rangle,$
(5)

which simply measures the fluctuation among various estimation values. Even
when $P(x)$ is strongly non-Gaussian, e.g. asymmetric or multimodal, the
sampling distribution $P_{n}(a)$ converges to a Gaussian distribution in the
limit $n\to\infty$ if the variance of $P(x)$ is finite, as the central limit
theorem states. These two distributions must not be confused.

For measures of the statistical error, we would primarily focus on the
standard error $\sigma_{A,n}$, i.e. the variance of the sampling distribution
$P_{n}(a)$, and the $68\%$ confidence interval $I_{A,n}$ defined as follows:

$\displaystyle\sigma_{A,n}^{2}$
$\displaystyle=\left\langle\left(\bar{A}_{n}\right)^{2}\right\rangle-\left\langle\bar{A}_{n}\right\rangle^{2},$
(6) $\displaystyle I_{A,n}$ $\displaystyle=[L_{A,n},U_{A,n}],$
$\displaystyle{\rm s.t.~{}}{\rm Prob}\left[\lim_{n\to\infty}\bar{A}_{n}\in
I_{A,n}\right]=0.68.$ (7)

Estimating lower and upper limits of the interval ($L_{A,n}$ and $U_{A,n}$)
from data yields the confidence interval. When $\bar{A}_{n}$ fluctuates around
$\lim_{n\to\infty}\bar{A}_{n}$ with a Gaussian noise,
$[\bar{A}_{n}-\sigma_{A,n},\bar{A}_{n}+\sigma_{A,n}]$ converges to the $68\%$
confidence interval in the $n\rightarrow\infty$ limit. When $n$ is finite,
however, we often find a non-Gaussian fluctuation and the standard error does
not represent the $68\%$ confidence interval.

In general, finding the analytic form of the sampling distribution for a
generic physical quantity is unfeasible. Instead, one could directly estimate
the sampling distribution by running as many MCMC runs as possible, although
it is computationally very demanding. Thus, it would be nice if we can
estimate the sampling distribution Eq. (2) from a single set of samples with
no knowledge on its analytic form. This kind of estimation problem is known as
the nonparametric inference 222The name comes from the fact that when we have
knowledge of the distribution form, then the problem essentially boils down to
estimating parameters for such distributions. in the literature of statistics
Wasserman (2006). For i.i.d. samples, Efron’s bootstrap method Efron (1979)
and its improved versions Efron and Tibshirani (1994) are known to work quite
well. When the samples have correlations, however, Efron’s bootstrap does not
work and some adaptation taking the correlation into account is necessary. For
the case of time correlations, we will demonstrate one of the nonparametric
bootstrap methods indeed gives the statistical properties accurately in a
typical situation.

We must also mention that in general, the convergence of the time series to
the stationary state becomes a serious issue with MCMC. However, since it is
more of a matter of the MCMC algorithms themselves rather than the error
estimation technique, in this paper we will assume that the MCMC has been
well-converged after the initial equilibration.

##  III The Bootstrap method

To resolve the aforementioned difficulties for evaluating statistical errors
from time-correlated data, statisticians have developed various nonparametric
methods. Among the many methods Carlstein (1986); Künsch (1989); Liu and Singh
(1992); Politis and Romano (1992), the stationary bootstrap Politis and Romano
(1994a), based upon Efron’s original bootstrap method for i.i.d. samples Efron
(1979), is one of the most versatile methods.

To explain the stationary bootstrap method, we first describe Efron’s
bootstrap method developed for analyzing i.i.d. samples Efron (1979),
following the standard textbook Efron and Tibshirani (1994) in the next
subsection.

###  III.1 Efron’s bootstrap

To describe Efron’s bootstrap method, let us assume that the samples have no
correlation in this subsection. In this case, the joint distribution
$q(x_{1},x_{2},\cdots,x_{n})$ of the data points $x_{1},x_{2},\dots,x_{n}$ can
be factorized with a probability distribution $q_{0}$ as follows:

$q(x_{1},x_{2},\dots,x_{n})=\prod_{i=1}^{n}q_{0}(x_{i}).$ (8)

If we could generate new samples from $q_{0}$ at a low computational cost,
then we can directly evaluate the sampling distribution in Eq. (2) by
generating as many samples as needed. However, generating new samples is often
computationally demanding, and thus, we will approximate the distribution
$q_{0}$ by some computationally tractable distribution.

The key observation here is that when the samples are independent, the
empirical distribution $q_{\rm emp}$ constructed by many samples approximates
$q_{0}$ quite precisely:

$q_{\rm emp}(x):=\frac{1}{n}\sum_{i=1}^{n}\delta(x-x_{i})\simeq q_{0}(x),\quad
n\gg 1,$ (9)

and that $q_{\rm emp}$ is unbiased in the following sense:

$\langle q_{\rm emp}(x)\rangle=q_{0}(x),\quad\forall n=1,2,\dots.$ (10)

Furthermore, sampling from the empirical distribution is computationally easy
because we can generate a new set of samples just by sampling from the already
obtained data $x_{1},x_{2},\cdots,x_{n}$ with replacement. These observations
suggest that the empirical distribution is a good approximation of $q_{0}$
that can be constructed from a single set of samples. Note that the above
observations do not apply for the case of any form of correlated samples
because the empirical distribution cannot serve as an approximation of the
true generation process of such samples even at $n\to\infty$.

In practice, the bootstrap method numerically computes the average in Eq. (2)
by resampling pseudo data set from the empirical distribution Eq. (9) as
follows:

  1. 1.

Construct the empirical distribution $q_{\rm emp}$ from the data as in Eq.
(9).

  2. 2.

For $b=1,2,\cdots,B$, repeat the following. Here $B$ is a parameter that
specifies how many pseudo data sets are to be sampled.

    1. (a)

Sample $n$ new data point $x_{1}^{*(b)},x_{2}^{*(b)},\cdots,x_{n}^{*(b)}$ from
the empirical distribution $q_{\rm emp}$. This pseudo data set is called the
bootstrap sample.

    2. (b)

Calculate the estimated value of the physical quantity $A$ as in Eq. (1) for
this pseudo data set
$\bar{A}_{n}^{*(b)}=A(x_{1}^{*(b)},x_{2}^{*(b)},\cdots,x_{n}^{*(b)})$.

  3. 3.

Approximate the sampling distribution in Eq. (2) as

$P_{n}(a)\simeq\tilde{P}_{n,B}(a):=\frac{1}{B}\sum_{b=1}^{B}\delta(a-\bar{A}_{n}^{*(b)}).$
(11)

The confidence interval for $\bar{A}_{n}$ can be calculated from the
approximated distribution obtained in Eq. (11). For example, the most naive
approximation for the $68\%$ confidence interval, which is known as the
bootstrap percentile method Efron and Tibshirani (1994); DiCiccio and Efron
(1996), yields $[L^{*},U^{*}]$, where $L^{*}$ and $U^{*}$ represent the $16$
and the $84\%$ quantiles of $\tilde{P}_{n}(a)$, respectively. Despite its
simplicity, this estimation often has a first-order accuracy Hall (1988), i.e.
$|0.68-{\rm Prob}[\lim_{n\to\infty}\bar{A}_{n}\in[L^{*},U^{*}]|=O(n^{-1/2})$
for sufficiently large $n$. Thus, the bootstrap method gives a reasonable
estimate of the confidence interval for any physical quantity without knowing
its analytic form from a single time-series data when the data satisfies the
i.i.d. assumption.

###  III.2 Stationary bootstrap

Whereas Efron’s original bootstrap method Efron (1979); Efron and Tibshirani
(1994) provides systematic, nonparametric estimation for statistical
properties, its applicability is limited to the i.i.d. case. The bootstrap
method has been argued to give a reasonable estimation even when applied to
correlated samples Newman and Barkema (1999), however it actually turns out
that the simple bootstrap method underestimates the statistical error by a few
orders of magnitude, i.e. it fails to correctly estimate the sampling
distribution even in the limit $n\to\infty$ (see Appendix A for further
discussions). Thus we should not use the i.i.d. bootstrap method for error
estimation of correlated samples.

For dependent samples with time correlations, several variants of the
bootstrap method were proposed Carlstein (1986); Künsch (1989); Liu and Singh
(1992); Politis and Romano (1992, 1994a). Those methods build bootstrap
samples by taking blocks of contiguous samples from a time series and piecing
them together. The only parameter for those algorithms is the (average) size
of the blocks. With an appropriate size of blocks, each block can be
considered as an independent short time series and a pseudo time series made
by piecing these blocks together have similar statistical properties to the
original time series in an asymptotic setting. In some limiting cases, the
literature rigorously proves that the methods actually give the exact
statistical properties. Here, we review the stationary bootstrap method
Politis and Romano (1994a) that respects the stationarity of the time series,
unlike other variants of the bootstrap method.

The stationary bootstrap assumes that a time series $x_{1},x_{2},\cdots,x_{n}$
is strictly stationary and weakly dependent, i.e. the time series is time-
translation invariant and its time correlation asymptotically decays to $0$
with time. These conditions are typically met for a well-relaxed MCMC data.
With these assumptions, we impose a ‘periodic boundary condition’ for the time
series, $x_{n+k}=x_{k}$, and build bootstrap samples. The practical procedures
of the stationary bootstrap method is the following:

  1. 0.

Set a parameter $p\in(0,1)$. This parameter controls the typical size of the
blocks sampled in the following step.

  2. 1.

Build $B$ bootstrap pseudo time series

$\boldsymbol{x}^{*(b)}=(x^{*(b)}_{1},x^{*(b)}_{2},\cdots,x^{*(b)}_{n}),\
b=1,2,\cdots,B,$

as follows:

    1. (a)

Choose $t_{1}\in[1,n]$ at random and let $x^{*}_{1}=x_{t_{1}}$.

    2. (b)

With probability $1-p$, let $t_{i+1}=t_{i}+1$, and with probability $p$,
choose $t_{i+1}\in[1,n]$ at random. Then let $x^{*}_{i+1}=x_{t_{i+1}}$. Regard
$x_{n+k}$ as $x_{k}$.

    3. (c)

Repeat (1b) until $i=n-1$.

  3. 2.

For each bootstrap pseudo time series $\boldsymbol{x}^{*(b)}$, calculate any
quantity of interest $\overline{A}^{*(b)}_{n}=A(\boldsymbol{x}^{*(b)})$, e.g.
$\overline{A}_{n}^{*(b)}(\boldsymbol{x}^{*(b)})=\sum x_{i}^{*(b)}/n$ if $A$ is
the sample mean.

  4. 3.

Approximate the sampling distribution of $A$

$P_{n}(a)\simeq
P^{*}_{n,B}(a):=\frac{1}{B}\sum_{b=1}^{B}\delta(a-\overline{A}^{*(b)}_{n}).$
(12)

This enables calculation of any statistical quantity of the distribution, the
most important values being the mean value as the estimator of the quantity
and the standard error as the confidence interval.

In this algorithm, the number $B$ of bootstrap samples should be large enough
to have the estimation almost independent of $B$. We carefully checked the $B$
dependence of the results, and set $B=10^{3}$ throughout this paper.

The probability $p$, which determines the typical length of the blocks of
successively chosen samples, controls the estimated value of the confidence
interval. Choosing an appropriate value for $p$ is essential in the stationary
bootstrap to handle the effect of correlation among the original time-series
data. It was rigorously shown that having the conditions $p\to 0$ and
$np\to\infty$ results in exact estimation in the limit $B,n\to\infty$ Politis
and Romano (1994a), but the optimal value for a given time series with a
finite length is a priori unknown. In practice, we choose $p$ to minimize the
mean-squared error of the coefficient of the estimated standard error of the
simple mean, $n\sigma^{2}_{n,\text{SB}}$, using the stationary bootstrap
method, Lahiri (1999); Politis and White (2004); Nordman (2009); Patton  _et
al._ (2009),

$\displaystyle\text{MSE}(n\sigma_{n,\text{SB}}^{2})$
$\displaystyle:=\left(\left\langle
n\sigma_{n,\text{SB}}^{2}\right\rangle-n\sigma^{2}_{n}\right)^{2}$
$\displaystyle\quad+\left\langle n\sigma_{n,\text{SB}}^{2}-\left\langle
n\sigma_{n,\text{SB}}^{2}\right\rangle\right\rangle^{2}$ (13)
$\displaystyle=G^{2}p^{2}+\frac{D}{pn}+o(p^{2})+o((pn)^{-1}),$ (14)

yielding

$p_{\text{opt}}=\left(\frac{D}{2nG^{2}}\right)^{1/3},$ (15)

where

$\displaystyle\sigma^{2}_{n}=\left\langle\left(\sum
x_{i}/n\right)^{2}\right\rangle-\left\langle\sum x_{i}/n\right\rangle^{2},$
(16) $\displaystyle G=2\sum_{t=0}^{\infty}|t|C(t),$ (17) $\displaystyle
D=2\left(C(0)+2\sum_{t=1}^{\infty}C(t)\right)^{2},$ (18)

and the unnormalized autocorrelation function

$C(t)=\left\langle\left(x_{0}-\left\langle
x_{i}\right\rangle\right)\left(x_{t}-\left\langle
x_{i}\right\rangle\right)\right\rangle.$ (19)

Eq. (15) automatically satisfies the rigorous condition for the correct
estimation, i.e. $p_{\text{opt}}\sim n^{-1/3}\to 0$ and $np\to\infty$ at $n\gg
1$, if $D/G^{2}$ is independent of $n$. Using these formulae, fortunately, the
optimal value of $p$ can also be estimated from the time series we wish to
analyze Politis and White (2004); Patton  _et al._ (2009) (see Appendix B for
a detailed algorithm to find $p_{\text{opt}}$). Note that, when samples are
i.i.d., $p_{\text{opt}}=1$, meaning that the stationary bootstrap becomes
equivalent to the original bootstrap method in this case.

The precise way we create these bootstrap samples are actually crucial. For
example, one may wonder what happens if we take a single block of contiguous
samples whose length is distributed according to the geometric distribution
with mean $\hat{p}^{-1}_{\rm opt}$ as the bootstrap sample
$\boldsymbol{x}^{*(b)}$. In this case, the length of the bootstrap time series
fluctuates with variance $\hat{p}^{-2}_{\rm opt}-\hat{p}^{-1}_{\rm
opt}\simeq\hat{p}^{-2}_{\rm opt}~{}(\hat{p}_{\rm opt}\ll 1)$, and because of
that fluctuation the procedure fails to give the correct statistical
properties.

Another reasonable attempt that actually works is to randomly take blocks with
a fixed length $\hat{p}^{-1}_{\text{opt}}$ of contiguous samples out of the
original time series, again as bootstrap time series. Blocks in this procedure
may overlap with each other. This algorithm is equivalent to the subsampling
method Politis and Romano (1994b), a variant of the resampling methods
applicable to stationary time series and a particular case of the bootstrap
method Lahiri (2003). Regarding that an optimal block length
$\ell^{\text{(block)}}_{\text{opt}}$ for the blocking method is also
$\ell^{\text{(block)}}_{\text{opt}}\sim n^{1/3}$ Wolff (2004), we find the
subsampling method includes the blocking as a special case: If we take only
non-overlapping blocks with length $\ell^{\text{(block)}}_{\text{opt}}$ from
the time series in the subsampling algorithm, that is the optimal blocking.
From this viewpoint, the blocking, as well as the stationary bootstrap,
belongs to the family of resampling methods for stationary time series.

In summary, the stationary bootstrap method gives error estimation and
statistical properties from only one single time series, using a practical
algorithm for finding an optimal value of the parameter $p_{\text{opt}}$. As
well as the original bootstrap method Efron (1979); Efron and Tibshirani
(1994), the stationary bootstrap makes no assumption on the distribution and
allows us to access the full (approximate) distribution of desired quantities.
Thus, more advanced statistical analyses, such as accelerated estimation of
confidence intervals Efron (1987); DiCiccio and Efron (1996); Efron and
Tibshirani (1994) and statistical hypothesis testing, may be combined with the
stationary bootstrap method. The method further provides a systematic way for
estimating statistical properties of more complicated quantities, e.g.
correlation functions and correlation lengths derived from them, and also the
latent heat at a first-order transition that is usually calculated from the
distribution function. A simple Python implementation can be found at
https://github.com/YoshihikoNishikawa/StationaryBootstrap.

##  IV Numerical test

In this section, we apply the methods we have reviewed in the previous
sections to the actual time series obtained by Monte Carlo simulations. We
calculate the squared standard error times the total number of samples,
$Mn\times\sigma^{2}_{n}$, and the cover rate of the estimated confidence
interval as functions of the length of time series $n$. The former quantity
$Mn\times\sigma^{2}_{n}$ should converge to a constant value in the limit
$n\rightarrow\infty$. By regarding the estimation using $10^{4}$ independent
samples as exact, we evaluate the performance of the other methods and the
estimation with fewer independent runs.

###  IV.1 Model

We consider the two-dimensional ferromagnetic Ising model on a square lattice
at the critical temperature. The Hamiltonian of the system is

$H\left(\boldsymbol{s}\right)=-J\sum_{\langle i,j\rangle}s_{i}s_{j},$ (20)

with $J>0$ and $s_{i}\in\\{\pm 1\\}$. The summation runs over the pairs of
neighboring lattice sites. The linear dimension of the lattice $L=32$ (the
total number of spins $N=L^{2}$), and the boundary conditions in both
directions are periodic. Starting from a random configuration in the high-
temperature limit, we equilibrate the system using the single-cluster Wolff
algorithm Wolff (1989) with $10^{4}$ cluster flips. We then use the Metropolis
algorithm with sequential update to measure the magnetization $m$ every
$10^{2}$ Monte Carlo sweep. For each time series $\\{m_{i}\\}_{i=1,\cdots,n}$,
we calculate three commonly studied quantities: the simple average of the
magnetization $\overline{m}=\sum_{i}m_{i}/n$, the magnetic susceptibility
$\chi/N=\overline{m^{2}}-\overline{m}^{2}$, and the Binder cumulant
$U=\overline{\left(m-\overline{m}\right)^{4}}/\overline{\left(m-\overline{m}\right)}^{2}$.
The correlation time $\tau_{m}$ of the dynamics, defined as the time scale of
the asymptotic exponential decay of the autocorrelation function of $m$, is
$1.2\times 10^{3}$ sweeps per spin for our setup. Since we measure the
magnetization once every $10^{2}$ Monte Carlo sweeps,
$n_{m}=\tau_{m}/10^{2}=1.2\times 10^{1}$ is the time-series length
corresponding to the correlation time. We rescale the time-series length $n$
by $n_{m}$ in the following.

As mentioned above, the initial relaxation to the stationary state needs
additional care. Using the Wolff cluster algorithm, which virtually eliminates
the difficulty at the critical temperature, we consider the ideal case where
the time series are sampled from the stationary state in order to study purely
the performance of the error estimation technique.

###  IV.2 Estimation of confidence interval

To estimate the correct statistical properties of the sampling distribution,
we prepare $10^{4}$ independent time series and calculate the standard errors
for $m$, $\chi/N$, and $U$ using the entire data set naively. As a more
realistic situation, we take a small number of independent samples, $M=5$ and
$10$, out of $10^{4}$ time series and calculate the standard errors as well.
Using $10^{4}/M$ standard errors each estimated with $M$ samples, we show the
median value of the distribution of the standard error in the following 333We
show the median because the empirical distribution of the confidence interval,
built with $10^{4}/M$ blocks, is highly asymmetric with a long tail. In this
case the median serves as the best value for a typically obtained result..

For the stationary bootstrap method, we estimate $p_{\text{opt}}$ for each
time series 444We should note that, while we have estimated $p_{\text{opt}}$
using the magnetization and used it for other quantities for simplicity, it is
possible to estimate $p_{\text{opt}}$ for each quantity of interest Lahiri
(1999, 2003), which would further enhance the performance., and calculate the
standard error and the cover rate of the estimated $68\%$ confidence interval
using the bootstrap percentile method, which we introduced in Sec. III.1. We
take $10^{3}$ independent time series to average the standard error and
confidence interval. We find, with the stationary bootstrap method, the mean
and median values have only a small difference, and thus we show the mean
values.

We also apply the blocking method to the magnetization of each time series. In
practice, we need to find a plateau in the standard error
$\sigma^{\text{(block)}}(R)$ as a function of the number of blocking procedure
$R$ for each time series, which corresponds to a fixed point of the
renormalization group Flyvbjerg and Petersen (1989). However, it is hard to
exactly define where the curve enters a plateau region due to the noise, and
in many cases no true plateaus even exist. Although an algorithm for
automatically selecting the optimal number of blocking was proposed in Ref.
Jonsson (2018), we take the value at which the absolute derivative of the
standard error as a function of the number of blocking

$\left|\Delta\sigma^{\text{(block)}}(k)\right|=\left|\frac{\sigma^{\text{(block)}}(k+1)-\sigma^{\text{(block)}}(k-1)}{2}\right|$
(21)

becomes the minimum for each time series, and average it over $10^{3}$ samples
for simplicity. As we will show below, this simple procedure works
sufficiently well and the estimation becomes precise when $n$ becomes large.

##  V Results

###  V.1 Standard error of $m$, $\chi/N$, and $U$

Figure 2: The coefficient of the standard error, $Mn\times\sigma_{n}^{2}$, as
a function of $n$ for (a) the magnetization $m$, (b) susceptibility $\chi/N$,
and (c) Binder cumulant $U$, estimated using $M=10^{4}$, $10$, $5$ independent
runs, the stationary bootstrap, and the blocking method (only for $m$). Inset
of (b) shows an enlarged view. The error bars represent the standard error
calculated with Efron’s bootstrap method. Figure 3: A typical curve of
$n\times\sigma_{n}^{2}$ of the magnetization $m$ obtained by the blocking
method. The time-series length $n=10^{5}$. The dashed line indicates the value
of $n\times\sigma_{n}^{2}$ estimated using $M=10^{4}$ independent runs.

We first discuss the coefficient of the standard error, that is, the squared
standard error multiplied by the total number of samples
$Mn\times\sigma^{2}_{n}$. We compare the values estimated by different
methods; the stationary bootstrap, blocking, and simple calculation using $M$
independent runs (see Fig. 2). At any $n$, the standard error estimated with
$M=10^{4}$ independent runs should be the most precise. The performances of
other methods are evaluated by how much their estimated coefficients are close
to the value of $M=10^{4}$ independent runs at each $n$. Since the standard
error decays as $(Mn)^{-1/2}$ asymptotically, the coefficient converges to a
positive constant value towards $n\to\infty$.

When the number of independent runs $M=10^{4}$, $Mn\times\sigma^{2}_{n}$
converges at $n/n_{m}\gtrsim 10^{1}$ for the magnetization and Binder
cumulant, and at $n/n_{m}\gtrsim 10^{4}$ for the susceptibility, to the value
at $n\to\infty$. Recall that $n_{m}$ corresponds to the correlation time
measured in the sample index. The simple estimation with practical $M=5,10$
gives smaller $Mn\times\sigma^{2}_{n}$ at large $n$ than that with $M=10^{4}$.
A longer simulation seems not to resolve this problem, as
$Mn\times\sigma^{2}_{n}$ remains smaller and almost independent of $n$ even
when $n/n_{m}$ reaches $10^{4}$. This may be somewhat surprising regarding how
common it is to adopt this method in practical situations. We conclude that
the simple estimation with a small $M$ does not capture the statistical
properties of the distribution correctly, and is unreliable. Since we never
know the true value of $Mn\times\sigma^{2}_{n}$ and its dependence on $M$ in
general, we must be cautious of simply using the standard error as the
confidence interval unless $M\gg 1$. We should note that the original
bootstrap would be useful to estimate the standard error more precisely for
the case with small $M$. Also, bootstrap acceleration methods Efron (1987);
DiCiccio and Efron (1996); Efron and Tibshirani (1994) would give a better
estimation of the $68\%$ confidence interval for the small-$M$ cases.

The stationary bootstrap, which uses only $M=1$ time series for estimating
statistical properties, gives $Mn\times\sigma^{2}_{n}$ that deviates from the
simple estimation with $M=10^{4}$ runs when the time series is short. Towards
larger $n$, it shows a converging behavior for any quantity, and eventually
becomes very close to the estimation of $10^{4}$ samples at large $n$, see
Fig. 2. For the magnetization, the blocking method also successfully estimates
$n\sigma^{2}_{n}$ on average at large $n$ (see Fig. 2 (a)). However, even when
the blocking method works correctly at $n/n_{m}\simeq 10^{4}$, it gives a
rather noisy curve for each time series and some do not even have a visible
plateau; see Fig. 3 for a typical curve obtained by the blocking method. This
suggests that estimating the standard error for a given time series with the
blocking method could be difficult.

It is worthwhile to consider the necessary number of samples in relation to
the correlation time of the dynamics $\tau_{m}\simeq 1.2\times 10^{3}$ Monte
Carlo sweeps per spin for the magnetization $m$. The time scale to reach the
limiting coefficient $Mn\times\sigma_{n}^{2}$ is much longer than
$n_{m}=\tau_{m}/10^{2}$ for all methods. This large discrepancy is natural
because approximating a distribution with precision needs many samples even in
the i.i.d. case. Yet, it is worth noting that the total number of samples
needed to reach the limiting $Mn\times\sigma_{n}^{2}$ value varies with
methods: For the simple calculation using $M=10^{4}$ independent runs, the
number is about $4\times 10^{2}\times M\simeq 4\times 10^{6}$, while that for
the stationary bootstrap and the blocking method is about $10^{5}$, which is
significantly smaller than the simple calculation using $M=10^{4}$ independent
runs.

In the next subsection, we directly test whether each method works properly
and gives a reasonable estimate for the confidence interval by calculating the
cover rate.

###  V.2 Cover rate of estimated confidence intervals

Figure 4: The cover rate of estimated $68\%$ confidence intervals for (a) the
magnetization $m$, (b) susceptibility $\chi/N$, and (c) Binder cumulant $U$,
estimated using $M=10$, $5$ independent runs, the stationary bootstrap, and
the blocking method (only for $m$). The dashed line indicates $68\%$. The
error bars represent the standard error calculated with Efron’s bootstrap
method.

In the previous section, we compared the speed of convergence in terms of the
standard error. A more direct way to compare the different methods is to see
the “cover rate” of the estimated confidence intervals. Although this is not a
quantity usually measured for practical purposes, it gives a reasonable
estimate of how reliable the different methods are.

We use the bootstrap percentile method Efron and Tibshirani (1994); DiCiccio
and Efron (1996) introduced in Sec. III.1. This simple method is applicable to
any distribution and allows us to calculate a confidence interval. Whether or
not the estimated confidence interval works properly can be assessed by
calculating the cover rate, that is, the ratio of independent samples having
the true value within the confidence intervals: If the estimated confidence
interval is accurate, the cover rate is equal to $68\%$.

Here, with the stationary bootstrap method, we estimate the $68\%$ confidence
interval. The cover rate is then calculated using $10^{3}$ independent runs.
Since the true values are unknown apart from the simple mean, which is $0$, we
use the simple averages over $10^{4}$ independent runs as the true values. For
comparison, we also calculate the cover rate of the standard error estimated
using $M=5$ and $10$ independent runs and that of the blocking method. For the
standard error of $M$ independent runs, the cover rate is calculated as the
ratio of $10^{4}/M$ sets of runs having the true value in the range of the
$68\%$ confidence interval. We use
$[\bar{A}_{n}-\sigma_{n},\bar{A}_{n}+\sigma_{n}]$ as the $68\%$ confidence
interval for these cases.

Fig. 4 shows the cover rate for the magnetization $m$, susceptibility
$\chi/N$, and the Binder cumulant $U$. When $M=5$ and $10$, the cover rate has
a clear systematic deviation from $68\%$ and does not reach the correct value
even at $n/n_{m}=10^{4}$. The deviation of the cover rate when $M=5$ and $10$
is consistent with the coefficient of the standard error,
$Mn\times\sigma_{n}^{2}$, that is also lower than the correct value, see Fig.
2. This indicates that the estimated statistical values are over-confident
when $M$ is small. Again, the standard error simply calculated using small $M$
is thus unreliable and should not be used in general, while, as noted in the
previous section, the original bootstrap would resolve the issue of this
method. The cover rate of the stationary bootstrap method, on the other hand,
reaches a value very close to $68\%$ at large $n$, meaning that the estimated
confidence interval indeed works correctly for all the physical quantities.
The cover rate of the blocking method for $m$ also reaches $68\%$ at large
$n$, but more slowly compared to the stationary bootstrap method.

Apparently, the confidence interval for the susceptibility $\chi/N$ is the
most difficult to estimate precisely, as it is the slowest to reach $68\%$,
see Fig. 4 (b). The difficulty comes from the strong asymmetric shape of the
sampling distribution of $\chi/N$ for small $n$, as shown in Fig. 1 (a). Other
quantities involving fluctuations are expected to have the same difficulty.
This suggests that, for those having asymmetric distributions, more
sophisticated acceleration methods Efron (1987); DiCiccio and Efron (1996);
Efron and Tibshirani (1994) would give a much better estimate for the
confidence interval. It is an interesting question whether and how much
acceleration methods improve the accuracy Hall (1988) for time-correlated
samples, as in the case of i.i.d. samples.

##  VI Summary and discussion

In summary, we reviewed some methods for estimating statistical properties of
actual MCMC-generated samples that have temporal correlations, with particular
focus on the stationary bootstrap method. We then discussed the connection
between the family of resampling methods (including the stationary bootstrap)
and blocking method, and argued the blocking method can be regarded as a
special case of the resampling methods. Using a realistic bench-marking data
set, we critically assessed the performance of the statistical methods. The
simple way of estimating the standard error using several independent runs was
revealed to be erroneous and may be over-confident when the number of
independent runs $M$ is too small. The stationary bootstrap method, which uses
only one time series, are much more precise than the simple way with a small
$M$.

When $M$ independent time series are readily available, assembling an
$M$-times-longer time series from them and then applying the stationary
bootstrap method to it should be more precise than the simple estimation using
several independent runs. The combined time series could be regarded as a MCMC
data with occasional insertion of a perfect Monte Carlo update where a
configuration is drawn with the equilibrium probability. The violation of the
temporal stationarity is minimal in this case, and so should be the effect to
the accuracy of the stationary bootstrap.

One of the major advantages of the stationary bootstrap is its ability to
access the full sampling distribution and applicability to a large class of
quantities, such as the specific heat, the Binder cumulant, normalized spatial
correlations, and even the latent heat at a first-order transition. More
complicated quantities involving multiple variables, such as the second-moment
correlation length Cooper  _et al._ (1982); Campostrini  _et al._ (2001);
Hasenbusch (2005); Nishikawa and Hukushima (2016) and the bulk and shear
moduli Xu  _et al._ (2012), are within the scope of the method as well Politis
and Romano (1994a). Another potential use of the method is to reliably
estimate the peak locations of measured functions such as the structure
factor.

Figure 5: The inverse of the estimated probability $\hat{p}^{-1}_{\text{opt}}$
as a function of $n/n_{m}$. At $n/n_{m}\gtrsim 10^{3}$,
$\hat{p}^{-1}_{\text{opt}}$ asymptotically grows as $\sim n^{1/3}$ as expected
from Eq. (15).

Finally, it is crucial to have enough samples in the data compared to the
correlation time. In the stationary bootstrap method, this could be checked by
comparing inverse of the estimated optimal probability
$\hat{p}^{-1}_{\text{opt}}$ with the time-series length $n$: Since $p^{-1}\sim
n^{1/3}$ is a necessary condition for the convergence of the stationary
bootstrap Politis and Romano (1994a), we should confirm if
$\hat{p}^{-1}_{\text{opt}}$ actually grows as $\sim n^{1/3}$. In our numerical
test, $\hat{p}^{-1}_{\text{opt}}$ indeed becomes asymptotically $\sim n^{1/3}$
at $n/n_{m}\gtrsim 10^{3}$, see Fig. 5. Another valid sanity check would be to
see if the estimated error scales as $n^{-1/2}$. Although some advanced
acceleration methods Efron (1987); Efron and Tibshirani (1994); DiCiccio and
Efron (1996) would be useful for better estimation of errors and confidence
intervals when the time series is rather short, the best is to simply perform
a much longer simulation and obtain a long enough time series. Otherwise we
get incorrect statistical errors with any method and would reach a wrong
scientific conclusion.

###### Acknowledgements.

The authors would like to thank W. Krauth and K. Hukushima for their critical
reading of the manuscript and useful comments. They also thank Y. Ozeki and A.
W. Sandvik for discussions. JT was supported by the National Science
Foundation (NSF) Focused Research Hub for Theoretical Physics, Center for
Quantum Information and Control (CQuIC) under Grant No. PHY-2116246 and
supported by the U.S. Department of Energy (DOE), Office of Science, National
Quantum Information Science Research Centers, Quantum Systems Accelerator
(QSA).

##  Appendix A Performance of Efron’s original bootstrap method when applied
to time series

Figure 6: (a) The coefficients of the standard error $n\times\sigma_{n}^{2}$
and (b) the cover rates of estimated $68\%$ confidence intervals for the
magnetization $m$ estimated by the the stationary bootstrap and the original
bootstrap method. The dashed line in (b) indicates $68\%$. The error bars
represent the standard error calculated with Efron’s bootstrap method.

In this appendix, we study and discuss the performance of Efron’s original
bootstrap when applied to time-correlated samples. We calculate the
coefficient of the standard error $n\times\sigma_{n}^{2}$ and the cover rate
of estimated $68\%$ confidence intervals for the magnetization $m$ with the
bootstrap percentile method. Fig. 6 shows them as functions of the time-series
length $n$, together with those of the stationary bootstrap method. The
original bootstrap method clearly fails to estimate the standard error: The
coefficient of the standard error $n\times\sigma_{n}^{2}$ differs from that of
the stationary bootstrap method by an order of magnitude, and the cover rate
is far below $68\%$. When applied to time series, the original bootstrap, as
designed for i.i.d. samples, randomizes the order of time-correlated samples.
Hence, resultant bootstrap samples lose the time correlation encoded in the
original samples. The argument made in Ref. Newman and Barkema (1999), saying
the original bootstrap method gives a reasonable error even if applied to
time-correlated samples, is thus incorrect. If one still wishes to use Efron’s
bootstrap method to a time series, the time interval between two successive
samples needs to be potentially very long, such as $2\tau$ with $\tau$ the
correlation time. This strategy makes the number of samples in the time series
fewer, possibly leading to worse statistics.

##  Appendix B Optimal parameter $p_{\text{opt}}$

In this appendix, a (semi)automatic algorithm to find the optimal value
$p_{\text{opt}}$ for a time series is presented. Following Ref. Politis and
White (2004), we calculate the unnormalized autocorrelation function Eq. (19)
from the time series, and numerically estimate Eq. (17) and Eq. (18).

Using the unnormalized autocorrelation function

$\hat{C}(k)=\sum_{i=1}^{N-|k|}(x_{i}-\overline{x})(x_{i+|k|}-\overline{x})/N$
(22)

with $\overline{x}=\sum_{i}x_{i}/n$, we find $T=2\hat{t}$ with $\hat{t}$ the
smallest number satisfying

$|\hat{C}(\hat{t}+k)/\hat{C}(0)|<c\sqrt{\log_{10}N/N}$ (23)

for $k=1,2,\cdots,K_{n}$ with a constant $c>0$ and $K_{n}=o(\log_{10}n)$. The
constant $c$ and $K_{n}$ should be carefully chosen not to underestimate $T$
although any choice of $c>0$ and $0<K_{n}<n$ should work asymptotically
Politis (2003). $c=2$ and $K_{n}=\max\left(5,\sqrt{\log_{10}n}\right)$ are
suggested as a safe choice in Refs. Politis (2003); Politis and White (2004).
Since $\hat{C}(k)$ is calculated only from one single time series, it is noisy
especially at long time. By finding $T$ in this manner, we estimate a rough
upper bound for the correlation time from the rather noisy autocorrelation.
Note that, when $\hat{t}$ and $T$ are comparable to $n$, the time series is
clearly not long enough and we need much more samples for precise estimation.
With $T$ and a window function Politis and Romano (1995); Politis (2003)

$w(k)=\left\\{\begin{array}[]{lll}1&\text{ if }|k|\in[0,1/2],\\\
2(1-|k|)&\text{ if }|k|\in[1/2,1],\\\ 0&\text{otherwise},\end{array}\right.$
(24)

we calculate

$\displaystyle\hat{G}=2\sum_{k=0}^{T}w(k/T)|k|\hat{C}(k),$ (25)
$\displaystyle\hat{D}=2\left(\hat{C}(0)+2\sum_{k=1}^{T}w(k/T)\hat{C}(k)\right)^{2},$
(26)

and finally reach

$\hat{p}_{\text{opt}}=\left(\frac{2\hat{G}^{2}}{\hat{D}}\right)^{-1/3}n^{-1/3},$
(27)

which is an estimation of the optimal parameter $p_{\text{opt}}$.

If an estimated $\hat{p}_{\text{opt}}$ is larger than $1$, the sample
correlation in the time series is negligible and the block size in the
stationary bootstrap method is $1$, meaning that the method is equivalent to
the original bootstrap method for i.i.d. samples in this case.

## References

  * Landau and Binder (2014) D. P. Landau and K. Binder,  _A Guide to Monte Carlo Simulations in Statistical Physics_ (Cambridge University Press, Cambridge, 2014).
  * Berg (2004) B. A. Berg,  _Markov Chain Monte Carlo Simulations and Their Statistical Analysis_ (WORLD SCIENTIFIC, 2004).
  * Krauth (2006) W. Krauth,  _Statistical Mechanics: Algorithms and Computations_ (Oxford University Press, Oxford, 2006).
  * Tanabashi  _et al._ (2018) M. Tanabashi, K. Hagiwara, K. Hikasa, K. Nakamura, Y. Sumino, F. Takahashi, J. Tanaka, K. Agashe, G. Aielli, C. Amsler, M. Antonelli, D. M. Asner, H. Baer, S. Banerjee, R. M. Barnett, _et al._ (Particle Data Group), Phys. Rev. D 98, 030001 (2018).
  * Manly (2018) B. F. Manly,  _Randomization, Bootstrap and Monte Carlo Methods in Biology_ (Chapman and Hall/CRC, 2018).
  * Andrieu  _et al._ (2003) C. Andrieu, N. De Freitas, A. Doucet, and M. I. Jordan, Machine Learning 50, 5 (2003).
  * Miller (1974) R. G. Miller, Biometrika 61, 1 (1974).
  * Efron (1979) B. Efron, The Annals of Statistics 7, 1 (1979).
  * Efron and Tibshirani (1994) B. Efron and R. Tibshirani,  _An Introduction to the Bootstrap_ (Chapman and Hall/CRC, 1994).
  * Wasserman (2004) L. Wasserman,  _All of Statistics_, Springer Texts in Statistics (Springer New York, New York, NY, 2004).
  * Wasserman (2006) L. Wasserman,  _All of Nonparametric Statistics_, Springer Texts in Statistics (Springer New York, New York, NY, 2006).
  * Flyvbjerg and Petersen (1989) H. Flyvbjerg and H. G. Petersen, The Journal of Chemical Physics 91, 461 (1989).
  * Newman and Barkema (1999) M. E. J. Newman and G. T. Barkema,  _Monte Carlo methods in statistical physics_ (Clarendon Press, Oxford, 1999).
  * Jonsson (2018) M. Jonsson, Physical Review E 98, 043304 (2018).
  * Carlstein (1986) E. Carlstein, The Annals of Statistics 14, 1403 (1986).
  * Künsch (1989) H. R. Künsch, The Annals of Statistics 17, 1403 (1989).
  * Liu and Singh (1992) R. Y. Liu and K. Singh, Moving blocks jackknife and bootstrap capture weak dependence, in  _Exploring the Limits of Bootstrap_ (Wiley, New York, NY, 1992) pp. 225–248.
  * Politis and Romano (1992) D. N. Politis and J. P. Romano, A circular block resampling procedure for stationary data, in  _Exploring the Limits of Bootstrap_ (Wiley, New York, NY, 1992) pp. 263–270.
  * Politis and Romano (1994a) D. N. Politis and J. P. Romano, Journal of the American Statistical Association 89, 1303 (1994a).
  * Note (1) In experimental fields, usually a 95% bar is set.
  * Wang (2015) W. Wang, Physical Review E 91, 053303 (2015).
  * Nishikawa and Hukushima (2016) Y. Nishikawa and K. Hukushima, Physical Review B 94, 064428 (2016).
  * Anderson  _et al._ (2017) J. A. Anderson, J. Antonaglia, J. A. Millan, M. Engel, and S. C. Glotzer, Physical Review X 7, 021001 (2017).
  * Note (2) The name comes from the fact that when we have knowledge of the distribution form, then the problem essentially boils down to estimating parameters for such distributions.
  * DiCiccio and Efron (1996) T. J. DiCiccio and B. Efron, Statistical Science 11, 189 (1996).
  * Hall (1988) P. Hall, The Annals of Statistics 16, 927 (1988).
  * Lahiri (1999) S. N. Lahiri, The Annals of Statistics 27, 386 (1999).
  * Politis and White (2004) D. N. Politis and H. White, Econometric Reviews 23, 53 (2004).
  * Nordman (2009) D. J. Nordman, The Annals of Statistics 37, 359 (2009).
  * Patton  _et al._ (2009) A. Patton, D. N. Politis, and H. White, Econometric Reviews 28, 372 (2009).
  * Politis and Romano (1994b) D. N. Politis and J. P. Romano, The Annals of Statistics 22, 1403 (1994b).
  * Lahiri (2003) S. N. Lahiri,  _Resampling Methods for Dependent Data_, Springer Series in Statistics (Springer New York, New York, NY, 2003).
  * Wolff (2004) U. Wolff, Computer Physics Communications 156, 143 (2004).
  * Efron (1987) B. Efron, Journal of the American Statistical Association 82, 171 (1987).
  * Wolff (1989) U. Wolff, Physical Review Letters 62, 361 (1989).
  * Note (3) We show the median because the empirical distribution of the confidence interval, built with $10^{4}/M$ blocks, is highly asymmetric with a long tail. In this case the median serves as the best value for a typically obtained result.
  * Note (4) We should note that, while we have estimated $p_{\text{opt}}$ using the magnetization and used it for other quantities for simplicity, it is possible to estimate $p_{\text{opt}}$ for each quantity of interest Lahiri (1999, 2003), which would further enhance the performance.
  * Cooper  _et al._ (1982) F. Cooper, B. Freedman, and D. Preston, Nuclear Physics B 210, 210 (1982).
  * Campostrini  _et al._ (2001) M. Campostrini, M. Hasenbusch, A. Pelissetto, P. Rossi, and E. Vicari, Physical Review B 63, 214503 (2001).
  * Hasenbusch (2005) M. Hasenbusch, Journal of Physics A: Mathematical and General 38, 5869 (2005).
  * Xu  _et al._ (2012) H. Xu, J. P. Wittmer, P. Polińska, and J. Baschnagel, Physical Review E 86, 046705 (2012).
  * Politis (2003) D. N. Politis, Journal of Nonparametric Statistics 15, 517 (2003).
  * Politis and Romano (1995) D. N. Politis and J. P. Romano, Journal of Time Series Analysis 16, 67 (1995).


$\displaystyle\quad{\color[rgb]{0.25,0.41,0.88}\definecolor[named]{pgfstrokecolor}{rgb}{0.25,0.41,0.88}{\left\\{\exists\langle\varphi\rangle\ldotp(\forall\langle\alpha\rangle\ldotp
0\leq\varphi(x)\leq\alpha(x)\land
0\leq\varphi(y)\leq\alpha(y)\land\varphi(k)\leq\alpha(k)\land\varphi(i)=\alpha(i))\land\varphi(k){-}\varphi(i)\prec
v\right\\}}}$ $\displaystyle\\}$ (IfSync)
$\displaystyle{\color[rgb]{0.25,0.41,0.88}\definecolor[named]{pgfstrokecolor}{rgb}{0.25,0.41,0.88}{\left\\{\exists\langle\varphi\rangle\ldotp(\forall\langle\alpha\rangle\ldotp
0\leq\varphi(x)\leq\alpha(x)\land
0\leq\varphi(y)\leq\alpha(y)\land\varphi(k)\leq\alpha(k)\land\varphi(i)=\alpha(i))\land\varphi(k){-}\varphi(i)\prec
v\right\\}}}$ Figure 17. Second part of the proof: Proving the first premise
of the rule While-$\exists$,
$\forall
v\ldotp\exists\langle\varphi\rangle\ldotp\vdash{\color[rgb]{0.25,0.41,0.88}\definecolor[named]{pgfstrokecolor}{rgb}{0.25,0.41,0.88}{\\{P_{\varphi}\land\varphi(i)<\varphi(k)\land
v=\varphi(k)-\varphi(i)\\}}}\;\mathbf{if}\;(i<k)\;\\{C_{body}\\}\;{\color[rgb]{0.25,0.41,0.88}\definecolor[named]{pgfstrokecolor}{rgb}{0.25,0.41,0.88}{\\{\exists\langle\varphi\rangle\ldotp
P_{\varphi}\land\varphi(k)-\varphi(i)\prec v\\}}}$.
For Cons (1), we simply choose $u=2$. For Cons (2), we notice that
$\varphi(i)<\varphi(k)$ and $\square(i\geq k)$ are inconsistent (this branch
is not taken at this stage), and thus the entailment trivially holds.
$\displaystyle{\color[rgb]{0.25,0.41,0.88}\definecolor[named]{pgfstrokecolor}{rgb}{0.25,0.41,0.88}{\left\\{\forall\langle\alpha\rangle\ldotp
0\leq\varphi(x)\leq\alpha(x)\land 0\leq\varphi(y)\leq\alpha(y)\right\\}}}$
$\displaystyle\mathbf{if}\;(*)\;\\{$
$\displaystyle\quad{\color[rgb]{0.25,0.41,0.88}\definecolor[named]{pgfstrokecolor}{rgb}{0.25,0.41,0.88}{\left\\{\forall\langle\alpha\rangle\ldotp
0\leq\varphi(x)\leq\alpha(x)\land 0\leq\varphi(y)\leq\alpha(y)\right\\}}}$
(Cons)
$\displaystyle\quad{\color[rgb]{0.25,0.41,0.88}\definecolor[named]{pgfstrokecolor}{rgb}{0.25,0.41,0.88}{\left\\{\forall\langle\alpha\rangle\ldotp\alpha(i)<\alpha(k)\Rightarrow\forall
v\ldotp v\geq 2\Rightarrow 0\leq\varphi(x)\leq 2*\alpha(x)+v\land
0\leq\varphi(y)\leq\alpha(y)+\alpha(x)*v\right\\}}}$
$\displaystyle\quad\mathbf{assume}\;\mathit{i<k}\mathbf{;}\;$ (AssumeS)
$\displaystyle\quad{\color[rgb]{0.25,0.41,0.88}\definecolor[named]{pgfstrokecolor}{rgb}{0.25,0.41,0.88}{\left\\{\forall\langle\alpha\rangle\ldotp\forall
v\ldotp v\geq 2\Rightarrow 0\leq\varphi(x)\leq 2*\alpha(x)+v\land
0\leq\varphi(y)\leq\alpha(y)+\alpha(x)*v\right\\}}}$
$\displaystyle\quad\mathit{r}\coloneqq{}\mathit{nonDet()}\mathbf{;}\;$
(HavocS)
$\displaystyle\quad{\color[rgb]{0.25,0.41,0.88}\definecolor[named]{pgfstrokecolor}{rgb}{0.25,0.41,0.88}{\left\\{\forall\langle\alpha\rangle\ldotp\alpha(r)\geq
2\Rightarrow 0\leq\varphi(x)\leq 2*\alpha(x)+\alpha(r)\land
0\leq\varphi(y)\leq\alpha(y)+\alpha(x)*\alpha(r)\right\\}}}$
$\displaystyle\quad\mathbf{assume}\;\mathit{r\geq 2}\mathbf{;}\;$ (AssumeS)
$\displaystyle\quad{\color[rgb]{0.25,0.41,0.88}\definecolor[named]{pgfstrokecolor}{rgb}{0.25,0.41,0.88}{\left\\{\forall\langle\alpha\rangle\ldotp
0\leq\varphi(x)\leq 2*\alpha(x)+\alpha(r)\land
0\leq\varphi(y)\leq\alpha(y)+\alpha(x)*\alpha(r)\right\\}}}$
$\displaystyle\quad\mathit{t}\coloneqq\mathit{x}\mathbf{;}\;$
$\displaystyle\quad\mathit{x}\coloneqq\mathit{2*x+r}\mathbf{;}\;$
$\displaystyle\quad\mathit{y}\coloneqq\mathit{y+t*r}\mathbf{;}\;$
$\displaystyle\quad\mathit{i}\coloneqq\mathit{i+1}$ (AssignS)
$\displaystyle\quad{\color[rgb]{0.25,0.41,0.88}\definecolor[named]{pgfstrokecolor}{rgb}{0.25,0.41,0.88}{\left\\{\forall\langle\alpha\rangle\ldotp
0\leq\varphi(x)\leq\alpha(x)\land 0\leq\varphi(y)\leq\alpha(y)\right\\}}}$
$\displaystyle\\}$ $\displaystyle\mathbf{else}\;\\{$
$\displaystyle\quad{\color[rgb]{0.25,0.41,0.88}\definecolor[named]{pgfstrokecolor}{rgb}{0.25,0.41,0.88}{\left\\{\forall\langle\alpha\rangle\ldotp
0\leq\varphi(x)\leq\alpha(x)\land 0\leq\varphi(y)\leq\alpha(y)\right\\}}}$
(Cons)
$\displaystyle\quad{\color[rgb]{0.25,0.41,0.88}\definecolor[named]{pgfstrokecolor}{rgb}{0.25,0.41,0.88}{\left\\{\forall\langle\alpha\rangle\ldotp\alpha(i)\geq\alpha(k)\Rightarrow
0\leq\varphi(x)\leq\alpha(x)\land 0\leq\varphi(y)\leq\alpha(y)\right\\}}}$
$\displaystyle\quad\mathbf{assume}\;\mathit{i\geq k}\mathbf{;}\;$
$\displaystyle\quad\mathbf{skip}$ (AssumeS, Skip)
$\displaystyle\quad{\color[rgb]{0.25,0.41,0.88}\definecolor[named]{pgfstrokecolor}{rgb}{0.25,0.41,0.88}{\left\\{\forall\langle\alpha\rangle\ldotp
0\leq\varphi(x)\leq\alpha(x)\land 0\leq\varphi(y)\leq\alpha(y)\right\\}}}$
$\displaystyle\\}$ (Choice)
$\displaystyle{\color[rgb]{0.25,0.41,0.88}\definecolor[named]{pgfstrokecolor}{rgb}{0.25,0.41,0.88}{\left\\{(\forall\langle\alpha\rangle\ldotp
0\leq\varphi(x)\leq\alpha(x)\land
0\leq\varphi(y)\leq\alpha(y))\otimes(\forall\langle\alpha\rangle\ldotp
0\leq\varphi(x)\leq\alpha(x)\land 0\leq\varphi(y)\leq\alpha(y))\right\\}}}$
(Cons)
$\displaystyle{\color[rgb]{0.25,0.41,0.88}\definecolor[named]{pgfstrokecolor}{rgb}{0.25,0.41,0.88}{\left\\{\forall\langle\alpha\rangle\ldotp
0\leq\varphi(x)\leq\alpha(x)\land 0\leq\varphi(y)\leq\alpha(y)\right\\}}}$
Figure 18. Third part of the proof. This proof outline shows
$\forall\varphi\ldotp\vdash{\color[rgb]{0.25,0.41,0.88}\definecolor[named]{pgfstrokecolor}{rgb}{0.25,0.41,0.88}{\\{Q_{\varphi}\\}}}\;\mathbf{if}\;(i<k)\;\\{C_{body}\\}\;{\color[rgb]{0.25,0.41,0.88}\definecolor[named]{pgfstrokecolor}{rgb}{0.25,0.41,0.88}{\\{Q_{\varphi}\\}}}$.

##  Appendix H Synchronous Reasoning Over Different Branches

One thesis of this paper is that reasoning about how sets of states are
affected by _one_ program command is powerful enough to reason about any
program hyperproperty, which is supported by our completeness result (Thm.
10).

However, reasoning about (for example) two executions of the same program
sometimes boils down to reasoning about two executions of two _different_ but
similar programs, because of branching. One _a priori_ appeal of relational
program logics over Hyper Hoare Logic is thus the ability to reason about two
different branches _synchronously_.

As an example, imagine that we want to reason about
$C^{\prime}\triangleq(\mathit{x}\coloneqq\mathit{x*2}\mathbf{;}\;C)+C$. Except
for the assignment that happens only in one branch, the two branches are
extremely similar. In a relational program logic, we can exploit this
similarity by first reasoning about the assignment on its own, and then
reasoning about the two remaining branches $C$ and $C$ synchronously, since
they are the same.

On the other hand, with the rule If from Fig. 3, we would have to reason about
the two branches $\mathit{x}\coloneqq\mathit{x*2}\mathbf{;}\;C$ and $C$
separately, even though they are closely related.

This is not a fundamental limitation of Hyper Hoare Logic. We can indeed
enable this kind of synchronous reasoning in Hyper Hoare Logic, by adding
specialized rules, as illustrated by Prop. 1 below.

Let us first define the following notation:

###### Notation 1.

$\displaystyle(A\otimes_{x=1,2}B)(S)\triangleq(A(\\{(l,\sigma)\mid(l,\sigma)\in
S\land l(x)=1\\})\land B(\\{(l,\sigma)\mid(l,\sigma)\in S\land l(x)=2\\}))$

The assertion $A\otimes_{x=1,2}B$ holds in a set $S$ iff the subset of all
states in $S$ such that $l(x)=1$ satisfies $A$, and the subset of all states
in $S$ such that $l(x)=2$ must satisfy $B$. We have proven in Isabelle/HOL the
following synchronous rule:

###### Proposition 0.

Synchronous if rule. If

  1. (1)

$\models{\color[rgb]{0.25,0.41,0.88}\definecolor[named]{pgfstrokecolor}{rgb}{0.25,0.41,0.88}{\\{P\\}}}\;C_{1}\;{\color[rgb]{0.25,0.41,0.88}\definecolor[named]{pgfstrokecolor}{rgb}{0.25,0.41,0.88}{\\{P_{1}\\}}}$

  2. (2)

$\models{\color[rgb]{0.25,0.41,0.88}\definecolor[named]{pgfstrokecolor}{rgb}{0.25,0.41,0.88}{\\{P\\}}}\;C_{2}\;{\color[rgb]{0.25,0.41,0.88}\definecolor[named]{pgfstrokecolor}{rgb}{0.25,0.41,0.88}{\\{P_{2}\\}}}$

  3. (3)

$\models{\color[rgb]{0.25,0.41,0.88}\definecolor[named]{pgfstrokecolor}{rgb}{0.25,0.41,0.88}{\\{P_{1}\otimes_{x=1,2}P_{2}\\}}}\;C\;{\color[rgb]{0.25,0.41,0.88}\definecolor[named]{pgfstrokecolor}{rgb}{0.25,0.41,0.88}{\\{R_{1}\otimes_{x=1,2}R_{2}\\}}}$

  4. (4)

$\models{\color[rgb]{0.25,0.41,0.88}\definecolor[named]{pgfstrokecolor}{rgb}{0.25,0.41,0.88}{\\{R_{1}\\}}}\;C_{1}^{\prime}\;{\color[rgb]{0.25,0.41,0.88}\definecolor[named]{pgfstrokecolor}{rgb}{0.25,0.41,0.88}{\\{Q_{1}\\}}}$

  5. (5)

$\models{\color[rgb]{0.25,0.41,0.88}\definecolor[named]{pgfstrokecolor}{rgb}{0.25,0.41,0.88}{\\{R_{2}\\}}}\;C_{2}^{\prime}\;{\color[rgb]{0.25,0.41,0.88}\definecolor[named]{pgfstrokecolor}{rgb}{0.25,0.41,0.88}{\\{Q_{2}\\}}}$

  6. (6)

$x\notin\mathit{rd}(P_{1})\cup\mathit{rd}(P_{2})\cup\mathit{rd}(R_{1})\cup\mathit{rd}(R_{2})$

Then
$\models{\color[rgb]{0.25,0.41,0.88}\definecolor[named]{pgfstrokecolor}{rgb}{0.25,0.41,0.88}{\\{P\\}}}\;(C_{1}\mathbf{;}\;C\mathbf{;}\;C_{1}^{\prime})+(C_{2}\mathbf{;}\;C\mathbf{;}\;C_{2}^{\prime})\;{\color[rgb]{0.25,0.41,0.88}\definecolor[named]{pgfstrokecolor}{rgb}{0.25,0.41,0.88}{\\{Q_{1}\otimes
Q_{2}\\}}}$.

This proposition shows how to reason synchronously about the program command
$(C_{1}\mathbf{;}\;C\mathbf{;}\;C_{1}^{\prime})+(C_{2}\mathbf{;}\;C\mathbf{;}\;C_{2}^{\prime})$.
Points 1) and 2) show that we can reason independently about the different
parts of the branches $C_{1}$ and $C_{2}$. Point 3) then shows how we can
reason synchronously about the execution of $C$ in both branches. Finally,
points 4) and 5) show how to go back to reasoning independently about each
branch.