
%% ================================================================
\section{Problem statement}
\label{sec:prob_stt}

The following problem occurs throughout Markov chain Monte Carlo (MCMC)
experiments.  Let $X_t$ be an \plainidx{identically distributed}, but not
necessarily independent, \plainidx{Markov process}; let $\muX$ and $\sigmaX^2$
be the common mean and variance, respectively.  (We will construct a specific
process $Y_t$ with the same properties.  We reserve the notation $X_t$ for a
general process with these properties.) Given a time-series\index{time series}
realization $X_0, \ldots, X_{N-1}$, the sole desired expressed in this paper is
to to estimate $\muX$, with an \plainidx{error bar} on that estimate.  The
presence of \plainidx{correlations} between the $X_t$'s make this process more
complicated than in the \plainidx{IID} case.

The standard \plainidx{estimator} for $\muX$ is the \plainidx{sample mean},
$\XNbar$.  Given a time-series realization $X_0, \ldots, X_{N-1}$, we compute a
single value of $\XNbar$.  Since the $X_t$'s are random variables, $\XNbar$ is
itself a random variable.  When we conduct $M$ such experiments, we will get
$M$ different values of $\XNbar$.  (We will quantify below the dependence of
the variance, or error bar, of $\XNbar$, upon the \plainidx{autocorrelation} of
the process $X_t$.) Suppose for the sake of discussion that the autocorrelation
is \plainidx{exponential}: $\Corr(X_i,X_j) = \eta^{|i-j|}$ for\index{$\eta$}
some $\eta \in [0,1)$.  Then $\eta=0$ is the IID case, and higher $\eta$'s
correspond to more highly correlated processes.  A few such realizations are
shown in figure \ref{fig:prob_stt_N_M_plots}.

\begin{figure}[!htb]
\begin{center}
\psfragscanon
\includegraphics[scale=0.9]{figures/eta_0.000_N_1000_M_5.eps}
\includegraphics[scale=0.9]{figures/eta_0.900_N_1000_M_5.eps}
\includegraphics[scale=0.9]{figures/eta_0.999_N_1000_M_5.eps}
\caption[Five realizations each of the correlated-uniform Markov process.]
	{Five realizations each of the correlated-uniform Markov process
	$Y_t$ with $\eta=0.0, 0.9, 0.999$.
	\label{fig:prob_stt_N_M_plots}}
\end{center}
\end{figure}

For any process $W_1, \ldots, W_K$, write $m_K(W)$ for the \plainidx{sample
mean} and $s^2_K(W)$ for the \plainidx{sample variance}, the \plainidx{unbiased
estimator} of $\Var(W)$.  Then:

\begin{itemize}

\item \plainidx{$m_N(X_t)$}, which is $\XNbar$, estimates $\muX$.  This is the
\emphidx{sample mean}, taken over $N$ samples.

\item \plainidx{$s^2_N(X_t)$} estimates $\sigmaX^2$.  This is the
\emphidx{sample variance}, taken over $N$ samples.

\item \plainidx{$m_M(\overline{X}_N)$} estimates $\muXNbar$.  This is also
referred to as the \plainidx{sample mean}; it is taken over $MN$ samples.

\item \plainidx{$s^2_M(\overline{X}_N)$} uses $MN$ data points to estimate
$\sigmaXNbar^2$, which is the \emphidx{variance of the sample mean}.

\item In the IID case, the true variance of the sample mean is
$\sigma^2_{\XNbar} = \sigmaX^2/N$; \plainidx{$t^2_N(X_t)$} $= s_N^2(X_t)/N$ is
the \emphidx{naive estimator} of the variance of the sample mean, using $N$
data points.  It is an unbiased estimator only in the IID case.

\item \plainidx{$u^2_N(X_t)$} is the \emphidx{corrected estimator} of
$\sigmaXNbar^2$.  The estimators $t^2_N(X_t)$ and $u^2_N(X_t)$ will be
discussed graphically, numerically, and theoretically below.  The estimated
\emphidx{integrated autocorrelation time} $\hattauint$ will\index{$\hattauint$}
be used to compute $u^2_N(X_t)$ from $t^2_N(X_t)$.

\item \plainidx{$\Var(u^2_N(X_t))$} is the \emphidx{error of the error bar}.
It turns out that $u^2_N(X_t)$ is a \plainidx{rough estimator} for
$\Var(\XNbar)$, and $\Var(u^2_N(X_t))$ increases with $\eta$.  The very name
``error of the error bar'' sounds overwrought; yet, it is a necessary
consideration in MCMC experiments, and must be thought through.

\end{itemize}

The processes \plainidx{$Y_t$} of figure \ref{fig:prob_stt_N_M_plots}, to be
defined explicity in section \ref{sec:corrunif}, have $\muY=1/2$ and
$\sigmaY^2=1/12$, regardless of the autocorrelation exponent $\eta$.  (Note
that $\sqrt{1/12} \approx 0.2887$.) We observe the following behavior from the
aforementioned estimators.  (See figures \ref{fig:results_YNbar} through
\ref{fig:results_errorbar} starting on page \pageref{fig:results_YNbar}, and
table \ref{table:results} on page \pageref{table:results}.)

\begin{itemize}

\item For all $\eta$, $\YNbar$ is unbiased for $\muY$.  Its uncertainty widens
visibly with autocorrelation exponent $\eta$.  This uncertainty is the quantity
of interest.

\item Quantitatively, $s_M(\YNbar)$ gives a good idea of this increasing
uncertainty.  However, $s_M(\YNbar)$ requires $M$ experiments, where $M$ may be
unacceptably large.  If we were always willing to conduct such a large number
of experiments, it would not be necessary to write this paper.  We wish to
estimate the variance of the sample mean \emph{using only one experiment $Y_0,
\ldots, Y_{N-1}$}.  This is the rub.

\item The \plainidx{corrected estimator} $u_N^2(Y_t)$ corresponds roughly with
$s_M(\YNbar)$, and moreover is computed from a single experiment $Y_0, \ldots,
Y_{N-1}$.  The roughness of the approximation of the error bar is acceptable:
it is only an error bar.

\end{itemize}

To summarize, $s^2_M(\XNbar)$ is a \plainidx{multi-experiment estimator} for
the variance of the sample mean; $u^2_N(X_t)$ is a \plainidx{single-experiment
estimator}.  The former is of higher quality, but is more expensive to obtain;
the latter carries its own uncertainty which worsens as the autocorrelation
$\eta$ increases.

Having motivated the problem, we now develop the notation and theory to make
all of these ideas precise.

%% ================================================================
\section{Autocovariance and autocorrelation}

\begin{defn} A \plainidx{Markov process} $X_t$, $t=0,1,2,\ldots$, is
\emphidx{stationary} if the $X_t$'s have a common mean $\muX=\bbE[X_t]$ and
variance $\sigmaX^2=\Var(X_t)$.  \end{defn}

\begin{defn}
Let $\Xt$ be a \plainidx{stationary Markov process} with $\bbE[\Xt] = \muX$
and $\Var(\Xt) = \sigmaX^2$.  The\index{covariance}\index{correlation}
\emphidx{autocovariance} and \emphidx{autocorrelation} of $\Xt$, respectively,
are\index{$C(t),c(t)$}
\begin{align*}
	C(t) &= \Cov(X_0, X_t)
		= \bbE[X_0 X_t] - \bbE[X_0]\bbE[X_t]
		= \bbE[X_0 X_t] - \muX^2
		\\
	c(t) &= \Corr(X_0, X_t)
		= \frac{\bbE[X_0 X_t] - \bbE[X_0]\bbE[X_t]}{\sigma_{X_0}\sigmaX}
		= \frac{\bbE[X_0 X_t] - \muX^2}{\sigmaX^2}.
\end{align*}
\end{defn}

\begin{rem}
In the literature, what we call the autocovariance is often referred to as
autocorrelation.  This incorrect and misleading terminology is, sadly, quite
widespread.
\end{rem}

\begin{rem}
Recall that, as with all correlations, the autocorrelation takes values between
$-1$ and $1$.
\end{rem}

%% ================================================================
\section{The IID uniform process}

Here\index{IID uniform process}\index{uniform process} we recall familiar
\cite{bib:cb,bib:gs} facts about random numbers $U$ which are uniformly
distributed on a closed interval $[a,b]$.  These will be used as building
blocks in section \ref{sec:corrunif}.  Writing the \plainidx{probability
density function} of $U$ as $f_U(x)$, we have
\begin{align*}
	f_U(x) &= \frac{1}{b-a} \cdot 1_{[a,b]}(x) &
	%
	\mu_U &= \frac{1}{b-a} \int_a^b x\,dx
		= \frac{a+b}{2} \\
	%
	\mu_U^2 + \sigma_U^2 =
	\bbE[U^2] &= \frac{1}{b-a} \int_a^b x^2\,dx
		= \frac{a^2+ab+b^2}{3} &
	%
	\sigma_U^2 &=\frac{(b-a)^2}{12}.
\end{align*}
Now consider an IID sequence $\{U_i\}$ of such random variables, indexed by the
integers.  We develop a particularly phrased formula which will simplify the
calculations in section \ref{sec:corrunif}.  Note that if $X_1,X_2$ are IID
with common mean $\mu_X$ and variance $\sigma_X^2$, then $\bbE[X_1^2] = \mu_X^2
+ \sigma_X^2$ whereas $\bbE[X_1 X_2] = \mu_X^2$.  For a sequence of IID
$X_i$'s, including our particular uniform $U_i$'s, this means
\begin{align}
\label{eqn:Edeltaij}
	\bbE[X_i X_j] = \muX^2 + \delta_{ij} \sigmaX^2.
\end{align}

%% ================================================================
\section{The correlated-uniform Markov process}
\label{sec:corrunif}

This paper addresses correlated Markov processes, focusing in particular on
those with exponential autocorrelation.  Here we construct a simple process for
which the mean, variance, and autocorrelation are exactly solvable.  In
particular, the autocorrelation will be controlled by a parameter
$\eta\in[0,1]$, while the mean and variance will be the same as for IID
$U(0,1)$.

\begin{defn}
\label{defn:cump}
Let \plainidx{$U$} be uniformly distributed on $[a,b]$ as in the previous
section, where $a<b$ are left variable for the moment.  Let $0 \le \eta \le 1$
and $a<b$.  The \emphidx{correlated-uniform Markov process} \plainidx{$Y_t$} is
defined by $Y_0 \sim U(a,b)$, and for $t \ge 1$,
\begin{equation}
\label{eqn:Y_t_def}
	\Yt
	\quad = \quad
	\eta Y_{t-1} + (1-\eta) U_t
	\quad = \quad
	\eta^t U_0 + (1-\eta) \sum_{i=1}^t \eta^{t-i} U_i
\end{equation}
where the first equality is a definition and the second equality follows by an
easy induction argument.
\end{defn}

\begin{rem*}
Note that $\eta=0$ is the IID case from the previous section; $\eta=1$ would
give a constant process with zero variance.  The $\eta$ parameter is the
\plainidx{control knob} with which we specify the autocorrelation of the
process, as will be made precise in section \ref{sec:corrunif_stats}.
\end{rem*}

\begin{defn}
\label{defn:stacump}
Closely related to this is the \emphidx{correlated-uniform stationary Markov
process} (or \plainidx{asymptotic process})
\begin{align}
\label{eqn:Yt_stationary}
	%Y_t &= (1-\eta) \sum_{i=0}^\infty \eta^i U_i,
	Y_t &= (1-\eta) \sum_{i=-\infty}^t \eta^{t-i} U_i.
\end{align}
\end{defn}

In practice, we will run the original process for a number of time steps $s$
until $\eta^s \approx 0$, such that the $\eta^s U_0$ term of equation
\eqref{eqn:Y_t_def} dies out, then consider the values of the process only from
that time forward.  In that regime, the process of definition
\ref{defn:stacump} is an approximation to that of definition \ref{defn:cump},
but it is easier to manipulate algebraically.

We seek $a,b$ such that the mean and variance of $\Yt$ do not depend on $\eta$.
The mean is immediate:
\begin{align*}
	\bbE[\Yt]
	&= (1-\eta) \sum_{i=-\infty}^t \eta^{t-i} \bbE[U_i] = \frac{a+b}{2}.
\end{align*}
The variance $\Var(\Yt)$ is a special case of the covariance $\Cov(\Yt, \Ytk)$,
which will be needed below.  Equation \eqref{eqn:Edeltaij} and expressions for
\plainidx{geometric sums} give us
\begin{align*}
	\bbE[\Yt \Ytk]
	&=
		(1-\eta)^2 \eta^{2t+k}
		\sum_{i=-\infty}^{t}   \eta^{-i}
		\sum_{j=-\infty}^{t+k} \eta^{-j} \;\bbE[U_i U_j] \\
	&=
		(1-\eta)^2 \eta^{2t+k}
		\sum_{i=-\infty}^{t}   \eta^{-i}
		\sum_{j=-\infty}^{t+k} \eta^{-j} \;(\mu_U^2 + \delta_{ij} \sigma_U^2) \\
	&=
		\mu_U^2 (1-\eta)^2 \eta^{2t+k}
		\sum_{i=-\infty}^{t}   \eta^{-i}
		\sum_{j=-\infty}^{t+k} \eta^{-j}
		+
		\sigma_U^2 (1-\eta)^2 \eta^{2t+k}
		\sum_{j=-\infty}^{t} \eta^{-2j} \\
	&=
		\mu_U^2
		+
		\sigma_U^2
		\; \eta^k
		\left(\frac{1-\eta}{1+\eta}\right).
\end{align*}
Then
\begin{align*}
	\Var(\Yt)
	&= \sigma_U^2 \; \left(\frac{1-\eta}{1+\eta}\right)
	= \frac{(b-a)^{2}}{12} \; \left(\frac{1-\eta}{1+\eta}\right).
\end{align*}

\begin{figure}[!htb]
\begin{center}
\psfragscanon
\includegraphics[scale=0.9]{figures/corrunif_plot1_eta_0.0.eps}
\includegraphics[scale=0.9]{figures/corrunif_plot1_eta_0.5.eps}
\includegraphics[scale=0.9]{figures/corrunif_plot1_eta_0.9.eps}
\caption[Realizations of $Y_t$ with $\eta=0.0, 0.5, 0.9$.]
	{Realizations of the correlated-uniform Markov process $Y_t$ with
	$\eta=0.0, 0.5, 0.9$.  Burn-in iterates are included.
	\label{fig:corrunif_plot1}}
\end{center}
\end{figure}

Now we may solve for $a$ and $b$ such that $\mu_U$ and $\sigma_U$ are the same
as for IID $U(0,1)$, namely, $1/2$ and $1/12$ respectively.  Solving the pair
of equations
\begin{align*}
	\frac{a+b}{2} &= \frac{1}{2} &
	& \textrm{and} &
	\frac{(b-a)^2}{12}\left(\frac{1-\eta}{1+\eta}\right) &= \frac{1}{12},
\end{align*}
we obtain
\begin{align}
\label{eqn:a_and_b_of_eta}
	a &= \frac{1}{2} \left(1 - \sqrt{\frac{1+\eta}{1-\eta}}\right) &
	& \textrm{and} &
	b &= \frac{1}{2} \left(1 + \sqrt{\frac{1+\eta}{1-\eta}}\right).
\end{align}
Note in particular that for $\eta=0$, the IID case, we recover $a=0,b=1$ as
expected.  Figure \ref{fig:corrunif_plot1} shows some realizations for
$\eta=0.0, 0.5, 0.9$.  For $\eta=0.9$, correlations are clearly visible.  Also
note that there is a \plainidx{burn-in time} required for the process to forget
its initial state $Y_0$.  In this figure, the asymptotic formula of definition
\ref{defn:stacump} appears valid for $t>50$ or so, at which point
$\eta^t=0.9^{50} \approx 0.005 \approx 0$.  This burn-in phenomenon is
discussed in more detail in section \ref{sec:tauint}.

The following is \plainidx{pseudocode} (technically, it is Python\index{Python
language} code, which is largely the same thing) for displaying $N$ steps of
$Y_t$, given the correlation-control parameter $\eta$ and the number
$N_\textrm{therm}$ of burn-in iterates to be discarded:
\begin{verbatim}
    s = sqrt((1+eta)/(1-eta)); a = 0.5 * (1 - s); b = 0.5 * (1 + s)

    Y = random.uniform(a, b) # Burn-in iterates
    for k in range(0, Ntherm):
        U = random.uniform(a, b)
        Y = eta * Y + (1-eta) * U

    for k in range(0, N):    # Iterates to be displayed
        U = random.uniform(a, b)
        Y = eta * Y + (1-eta) * U
        print Y
\end{verbatim}

%% ================================================================
\section{Statistics of the correlated-uniform Markov process}
\label{sec:corrunif_stats}

We now write all statistics of the correlated-uniform Markov process $Y_t$ in
terms of $\eta$.  With $a$ and $b$ in terms of $\eta$ (equation
\eqref{eqn:a_and_b_of_eta}), we have
\begin{align*}
	\mu_U &= \frac{a+b}{2}
		= \frac{1}{2}, &
	\sigma^2_U &= \frac{(b-a)^2}{12}
		= \frac{1}{12} \left(\frac{1+\eta}{1-\eta}\right), \\
	%
	\bbE[\Yt] =
	\muY &= \frac{a+b}{2} = \frac{1}{2}, &
	\Var(\Yt) =
	\sigmaY^2 &= \frac{(b-a)^2}{12} \left(\frac{1-\eta}{1+\eta}\right)
		= \frac{1}{12}; \\
\end{align*}
\begin{align*}
	\bbE[\Yt \Ytk]
		&= \mu_U^2 + \sigma_U^2 \eta^k \left(\frac{1-\eta}{1+\eta}\right)
		= \frac{1}{4} + \frac{\eta^k}{12}, \\
	%
	\Cov(\Yt, \Ytk) &= \bbE[\Yt \Ytk] - \bbE[\Yt] \bbE[\Ytk]
		= \frac{\eta^k}{12}, \\
	%
	\Corr(\Yt, \Ytk)
		&= \frac{\Cov(\Yt, \Ytk)}{\sigmaY \sigma_{\Ytk}}
		= \eta^k.
\end{align*}

\begin{figure}[!htb]
\begin{center}
\psfragscanon
\includegraphics[scale=0.9]{figures/corrunif_histo_eta_0.0.eps}
\includegraphics[scale=0.9]{figures/corrunif_histo_eta_0.5.eps}
\includegraphics[scale=0.9]{figures/corrunif_histo_eta_0.9.eps}
\caption[Histograms of the correlated-uniform Markov process]
	{Histograms of the correlated-uniform Markov process $Y_t$ with $\eta=0.0,
	0.5, 0.9$: $10^6$ iterates, bins of $0.1$ from $-0.7$ to $1.7$.  Burn-in
	iterates have been discarded.
	\label{fig:corrunif_histos}}
\end{center}
\end{figure}
The remaining step needed to completely specify the correlated-uniform Markov
process is to write down the PDF of $\Yt$.  This could be done using
convolutions, since $\Yt$ is a weighted sum (weighted by powers of $\eta$) of
IID uniform random variables.  The algebra is messy, though, and an expression
for the PDF is not needed in this work.  It is sufficient to point out the
following: ($i$) For $\eta=0$, the density is uniform on $[0,1]$.  ($ii$) For
$\eta$ close to 1, which is the case of interest in this work, the density
closely resembles a normal with mean $1/2$ and variance $1/12$.  The support is
compact, so the density cannot be Gaussian, but the support is wide enough to
include substantial tail mass.  See figure \ref{fig:corrunif_histos} for
empirical \plainidx{histograms}.

%% ----------------------------------------------------------------
\section{The variance of the sample mean}
\label{sec:Var_XNbar}

When we use the data from an MCMC simulation to compute the sample mean of a
random variable, the next order of business is to place an \plainidx{error bar}
on that sample mean.

As before, let $X_t$ be a stationary Markov process with common mean $\muX$,
variance $\sigmaX$, and autocorrelation $\Corr(\Xt,\Xtk) = \eta^k$.  Given
$X_0, \ldots, X_{N-1}$, the \plainidx{sample mean} $\XNbar$ is an unbiased
estimator of $\muX$:
\begin{align*}
	\XNbar &= \frac{1}{N} \sumiN X_i.
\end{align*}
By linearity of expectation, $\bbE[\XNbar] = \muX$.  To find the variance of
$\XNbar$, we first need $\bbE[\XNbar^2]$.  This is
\begin{align*}
	\bbE[\XNbar^2] &= \frac{1}{N^2}
		\sumiN \sumjN \bbE[X_i X_j].
\end{align*}
Since
\begin{align*}
	\Corr(X_i,X_j) &= \eta^{|i-j|} = \frac{\bbE[X_i,X_j] - \muX^2}{\sigmaX^2},
\end{align*}
we have
\begin{align}
\label{eqn:E_Xi_Xj}
	\bbE[X_i X_j] = \muX^2 + \sigmaX^2 \eta^{|i-j|}.
\end{align}
Then
\begin{align*}
	\bbE[\XNbar^2] &= \frac{1}{N^2}
		\sumiN \sumjN (\muX^2 + \sigmaX^2 \eta^{|i-j|})
	= \muX^2 + \frac{\sigmaX^2}{N^2}
		\sumiN \sumjN \eta^{|i-j|} \\
	&= \muX^2 + \frac{\sigmaX^2}{N^2} \left[
		\sumiN 1 +
		\sum_{i=0}^{N-2} \eta^{-i} \sum_{j=i+1}^{N-1} \eta^j +
		\sum_{i=1}^{N-1} \eta^i    \sum_{j=0}^{i-1}   \eta^{-j}
		\right].
\end{align*}
Applying geometric-sum\index{geometric sum} formulas and several lines of
algebra, we get
\begin{align*}
	\bbE[\XNbar^2] &= \muX^2 + \frac{\sigmaX^2}{N}
		+ \frac{2 \sigmaX^2 \eta}{N^2(1-\eta)} \left[
			(N-1) - \left(\frac{\eta-\eta^N}{1-\eta}\right)
		\right].
\end{align*}
With $N \approx N-1$ we have
\begin{align*}
	\bbE[\XNbar^2] &\approx \muX^2
		+ \frac{\sigmaX^2}{N}\left(\frac{1+\eta}{1-\eta}\right)
		- \frac{2 \sigmaX^2 \eta^2}{N^2(1-\eta)^2}(1-\eta^{N-1}).
\end{align*}
With $\eta^N \approx 0$ and a bit more algebra we have
\begin{align}
\label{eqn:Var_XNbar_eta_ratio}
	\bbE[\XNbar^2] &\approx \muX^2
		+ \frac{\sigmaX^2}{N}\left(\frac{1+\eta}{1-\eta}\right) &
	&\tand&
	\Var(\XNbar) &\approx
		\frac{\sigmaX^2}{N}\left(\frac{1+\eta}{1-\eta}\right).
\end{align}
Recall that for the IID case ($\eta=0$) we have $\Var(\XNbar) = \sigmaX^2/N$.
This expression recovers that; furthermore, correlations enlarge the error bar
on the sample mean.

%% ================================================================
\section{Estimates of autocorrelation}
\label{sec:c_hat}

Throughout this section, let $\Xt$ be a stationary Markov process with
$\bbE[\Xt] = \muX$, $\Var(\Xt) = \sigmaX^2$, and $\Corr(\Xt,\Xtk) = \eta^{k}$.
(Without loss of generality, take $k \ge 0$.) The simple correlated-uniform
Markov process of section \ref{sec:corrunif} is one example of this; moreover,
an MCMC process on a finite state space may take this form.  (As described in
\cite{bib:berg}, $\eta$ is related to the \plainidx{second dominant eigenvalue}
of the transition matrix of the Markov process.  If the third dominant
eigenvalue is comparable with the second, then the autocorrelation will not
take the simple exponential form described here.)

\begin{rem}
\label{rem:tauexp}
In the literature, one more often sees $\Corr(X_0, X_t) = \exp(-t/\tauexp)$.
Then $\tauexp$ and $\eta$ are put into one-to-one correspondence by
\begin{align*}
	\tauexp &= -1/\log\eta &
	&\tand&
	\eta = \exp(-1/\tauexp).
\end{align*}
\end{rem}

For the correlated-uniform process, the autocorrelation is already known; for a
general MCMC process, one wishes to estimate $\eta$ (or $\tauexp$) from
realization data.  Recall that
\begin{align}
\label{eqn:stacorr}
	\Corr(\Xt,\Xtk) &= \frac{\bbE[\Xt\Xtk]-\bbE[\Xt]\bbE[\Xtk]}
		{\sigma_{\Xt}\sigma_{\Xtk}}
	= \frac{\bbE[X_0 X_k]-\muX^2}{\sigmaX^2}
\end{align}
where the second equality holds by the stationarity of the process, and that we
always have
\begin{align}
\label{ineq:corr_bounds}
	-1 \le \Corr(\Xt,\Xtk) \le 1.
\end{align}
(This holds for the correlation of any pair of random variables.)
Also recall that
\begin{align}
\label{eqn:sigma_Xt}
	\sigmaX^2 = \bbE[X_t^2] - \bbE[X_t]^2.
\end{align}
Recall as well
\cite{bib:cb} that, for $M$ realizations $\Xt^{(0)}, \ldots, \Xt^{(M-1)}$ of
$\Xt$, the \plainidx{unbiased estimator} for the variance of $\Xt$ is
\begin{align}
\label{eqn:var_ube}
	s^2_{X_t} &= \frac{1}{M-1} \left[
		\sumiM (\Xt^{(i)})^2 -
			\frac{1}{M} \left(\sumiM \Xt^{(i)}\right)^2
		\right].
\end{align}

\begin{defn}
Fix $t$ and $k$.  The \emphidx{multi-realization estimator} of the
autocorrelation $\Corr(\Xt,\Xtk)$, requiring $M$ realizations $\Xt^{(0)},
\ldots \Xt^{(M-1)}$ of the process, is a straightforward combination of
equations \eqref{eqn:stacorr}, \eqref{eqn:sigma_Xt}, and \eqref{eqn:var_ube}.
Namely,
\begin{align*}
	\cmhat(t,k) &= \frac
		{\frac{1}{M} \sumiM \left(\Xt^{(i)}\Xtk^{(i)}\right) -
			\frac{1}{M^2}
			\left(\sumiM \Xt^{(i)}\right)
			\left(\sumjM \Xtk^{(j)}\right)}
		{
		\frac{1}{M-1}
		\left[
			\sumiM (\Xt^{(i)})^2 -
				\frac{\left(\sumiM \Xt^{(i)}\right)^2}{M}
			\right]^{1/2}
		\left[
			\sumjM (\Xtk^{(j)})^2 -
				\frac{\left(\sumjM \Xtk^{(j)}\right)^2}{M}
			\right]^{1/2}
		}.
\end{align*}
\end{defn}

\begin{rem*}
Since the process is stationary, one may be tempted to reuse the $\Xt$ variance
estimator for $\Xtk$ --- after all, they estimate the same quantity
$\sigmaX^2$.  In practice, however, doing so tends to produce autocorrelation
estimates which fall (quite far) outside the range $[-1,1]$, violating
inequality \ref{ineq:corr_bounds}.  That  is, the second equality in equation
\eqref{eqn:stacorr} holds theoretically but not at the estimator level.  This
same remark holds for the sliding-window estimator, to be defined next.
\end{rem*}

The difficulty with the multi-realization estimator is that realizations $\Xt$
can be expensive to compute.  Rather than running $M$ processes from $t=0$ up
to some $N$, which takes $O(MN)$ process-generation time, perhaps we can
(carefully) use the stationarity of the process, estimating the autocorrelation
using only a single realization.  This will take only $O(N)$ process-generation
time.

\begin{defn}
Given a single realization $X_0, \ldots, X_{N-1}$, take $k$ from $0,1,2,
\ldots,N-2$.  The \emphidx{sliding-window estimator}\index{window} of the
autocorrelation $\Corr(X_0,X_k)$, is
\begin{align}
\label{eqn:c_hat_k}
	\chat(k) &= \frac
		{\frac{1}{N-k} \sum_{i=0}^{N-k-1} \left(X_i X_{i+k}\right) -
			\frac{1}{(N-k)^2}
			\left(\sum_{i=0}^{N-k-1} X_i \right)
			\left(\sum_{j=0}^{N-k-1} X_{j+k} \right)
		}{
		\left( \frac{1}{N-k-1} \right)
		\left[
			\sum_{i=0}^{N-k-1} X_i^2 -
				\frac{\left(\sum_{i=0}^{N-k-1} X_i \right)^2}{N-k}
			\right]^{1/2}
		\left[
			\sum_{j=0}^{N-k-1} X_{j+k}^2 -
				\frac{\left(\sum_{j=0}^{N-k-1} X_{j+k}\right)^2}{N-k}
			\right]^{1/2}
		}.
\end{align}
\end{defn}

This formula is perhaps intimidating, but is made quite simple with the
aid of the example below, wherein $N=10$ and $k-2$.  Namely:
\begin{itemize}
\item We consider all pairs separated by $k$ time steps:
	$X_0 X_k$, $X_1 X_{k+1}$, \ldots, $X_{N-k-1}X_{N-1}$.
	There are $N-k$ such pairs.
\item The first elements in each pair form a window from $X_0$ to $X_{N-k-1}$.
\item The second elements in each pair form a window from $X_k$ to $X_{N-1}$.
\item We estimate the mean and variance of $X_0$ by the sample mean and
	sample variance over the first window.
\item We estimate the mean and variance of $X_k$ by the sample mean and
	sample variance over the second window.
\item We estimate the cross-moment $\bbE[\Xt\Xtk]$ by the sample mean
	over pair products.
\end{itemize}

\beginex
There are $N=10$ samples, $X_0$ through $X_9$:
$$
	\begin{tabular}{|c|c|c|c|c|c|c|c|c|c|}
	\hline 
	$X_0$&$X_1$&$X_2$&$X_3$&$X_4$&$X_5$&$X_6$&$X_7$&$X_8$&$X_9$\\
	\hline
	\end{tabular}
$$
Picking $k=2$, there are two windows of length $N-k=8$:
$$
	\begin{tabular}{|c|c|c|c|c|c|c|c|c|c|c|c|}
	\hline 
	$X_0$&$X_1$&$X_2$&$X_3$&$X_4$&$X_5$&$X_6$&$X_7$&  {} &  {} \\
	\hline 
  	{} &  {} &$X_2$&$X_3$&$X_4$&$X_5$&$X_6$&$X_7$&$X_8$&$X_9$\\
	\hline
	\end{tabular}
$$
Equation \eqref{eqn:c_hat_k} has five distinct sums:
the sum of $X_0$ through $X_7$,
the sum of squares of $X_0$ through $X_7$,
the sum of $X_2$ through $X_9$,
the sum of squares of $X_2$ through $X_9$,
and the cross sum $X_0 X_2 + \ldots + X_7 X_9$.
\mendex

\begin{figure}[!htb]
\begin{center}
\psfragscanon
\includegraphics[scale=0.6]{figures/autocorr_cmp_eta_0.9_xmax_1000.eps}
\includegraphics[scale=0.6]{figures/autocorr_cmp_eta_0.9_xmax_50.eps}
\caption[Autocorrelation and estimators thereof for $Y_t$ with $\eta=0.9$.]
	{Autocorrelation and estimators thereof for $Y_t$ with $\eta=0.9$.
	Burn-in iterates have been discarded.  The second plot zooms in on the
	first 50 samples of the first plot.
	\label{fig:corrunif_autocorr}}
\end{center}
\end{figure}

\begin{rem}
\label{rem:chat_bias}
One would hope that $\chat(t)$ is an unbiased estimator of $c(t)$.  Finding its
expectation using the definition is intimidating:  we have a ratio of products
of sums of correlated random variables.  Taking an experimental approach
instead, making multiple plots of the form of figure
\ref{fig:corrunif_autocorr}, one finds that $\chat(t)$ does in fact
fractionally underestimate $c(t)$.  This affects the estimated integrated
autocorrelation time, as discussed in remark \ref{rem:hattauint_bias}.
\end{rem}

This estimator has the benefit of making use of all the data in a single
realization.  Its drawback is that, for larger $k$, the sample size $N-k$ is
small.  Thus, the error in the estimator increases for larger $k$.

Figure \ref{fig:corrunif_autocorr} compares estimators against the true value
$c(k) = \Corr(X_0, X_k) = \eta^k$ for $\eta=0.9$.  Here, $N=1000$ time steps
have been used; $M=1000$ realizations for the multi-realization estimator
$\chat_m(k)$.  Note that the decreasing sample size, $N-k$, of the
sliding-window estimator $\chat(k)$ increases the error of this estimator.  For
this reason, $\chat_s(k)$ is also plotted.  This is the same as $\chat_m(k)$,
but with $M=N-k$.  The first plot shows the autocorrelation estimators for
$k=0$ to $998$; the second zooms in on the first $50$ values of $k$.

\begin{rem}
\label{rem:c_hat_wild}
We observe the following:
\begin{itemize}
\item Comparing the full-length and short-length multi-realization estimators
$\chat_m(k)$ vs. $\chat_s(k)$ shows that decreasing sample size does have an
effect for larger $k$.  Nonetheless, the sliding-window estimator $\chat(k)$
shows markedly wilder behavior for larger $k$, which cannot be accounted for by
small-sample-size effects alone.

\item For all three estimators, errors are small when $k$ is small, which is
when the true autocorrelation $c(k) = \eta^k$ is non-negligible.

\item Thus, one should examine estimators of the autocorrelation only for
values of $k$ until the estimators approach zero.  Values past that point are
neither accurate nor needed.
\end{itemize}
\end{rem}

%% ================================================================
\section{Integrated autocorrelation time}
\label{sec:tauint}

Following \cite{bib:berg}, we develop the notion of \plainidx{integrated
autocorrelation time} as follows.  We reconsider the variance of the sample
mean (see section \ref{sec:Var_XNbar}) from a different point of view.  Again,
$X_t$ is a stationary Markov process with common mean $\muX$ and common
variance $\sigmaX^2$.  We have
\begin{align*}
	\Var(\XNbar) &= \bbE[(\XNbar - \muX)^2]
	= \frac{1}{N^2} \sumiN \sumjN
		\bbE[(X_i - \muX)(X_j - \muX)] \\
	&= \frac{1}{N^2} \sumiN \sumjN
		\bbE[X_i X_j - \muX X_i - \muX X_j + \muX^2] \\
	%
	&= \frac{1}{N^2} \sumiN \sumjN \left(\bbE[X_i X_j] - \muX^2 \right)
	= \frac{1}{N^2} \sumiN \sumjN \Cov(X_i, X_j) \\
	%
	&= \frac{1}{N^2} \left[
		\sumiN \Var(X_i) + 2 \sum_{t=1}^{N-1} (N-t) \; \Cov(X_0, X_t)
		\right] \\
	&= \frac{\sigmaX^2}{N}
		+ 2 \sigmaX^2 \sum_{t=1}^{N-1} (N-t) \; \Corr(X_0, X_t) \\
	%
	&= \frac{\sigmaX^2}{N} \left[
		1 + 2 \sum_{t=1}^{N-1} \left(1 - \frac{t}{N}\right) \; \Corr(X_0, X_t)
		\right]
	\approx \frac{\sigmaX^2}{N} \left[
		1 + 2 \sum_{t=1}^{\infty} \Corr(X_0, X_t)
		\right].
\end{align*}
If $X_t$ is IID then we recover the familiar $\Var(\XNbar) = \sigmaX^2/N$;
otherwise we have
\begin{align}
\label{eqn:Var_XNbar_tauint}
	\Var(\XNbar) = \frac{\sigmaX^2}{N} \; \tauint
\end{align}
where $\tauint$ is\index{$\tauint$} the last bracketed expression
above.  Note as well that if $\Corr(X_0,X_k) = \eta^k$, then
\begin{align}
\label{eqn:tauint_fcn_of_eta}
	\tauint &= 1 + 2 \sum_{t=1}^{\infty} \eta^t
		= 1 + \frac{2 \eta}{1-\eta}
		= \frac{1+\eta}{1-\eta}
\end{align}
which is what we would have expected by comparing equations
\eqref{eqn:Var_XNbar_eta_ratio} and \eqref{eqn:Var_XNbar_tauint}.  As a
consequence, when $c(t) = \eta^t$ we have \begin{align}
\label{eqn:tauint_and_eta}
	\tauint &= \frac{1+\eta}{1-\eta} &
	&\tand&
	\eta &= \frac{\tauint - 1}{\tauint + 1}.
\end{align}
Some values are shown for reference in table \ref{table:eta_vs_tauint}.

\begin{table}[!htb]
\begin{center}
\begin{tabular}{|c||c|c|c|c|c|c|c|c|c|c|c|c|}
\hline $\eta$              & 0     & 0.1   & 0.2   & 0.5   & 0.6   & 0.9 & 0.990 & 0.999 \\
\hline $(1+\eta)/(1-\eta)$ & 1     & 1.222 & 1.500 & 3.000 & 4.000 & 19  & 199   & 1999  \\
\hline
\end{tabular}
\caption{$\eta$ vs. $(1+\eta)/(1-\eta)$.
	\label{table:eta_vs_tauint}}
\end{center}
\end{table}

\begin{rem}
\label{rem:tauint_IID}
If the process is IID, i.e. $\eta=0$, then $c(0)=1$, $c(t)=0$ for all
$t \ge 1$, and $\tauint = 1$.
\end{rem}

\begin{defn}
Recall that $s_N^2(X_t)$ (equation \eqref{eqn:var_ube}) estimates $\sigmaX^2$.
Using equation \eqref{eqn:Var_XNbar_tauint}, the \emphidx{naive estimator} and
\emphidx{corrected estimator} of $\Var(\XNbar)$ are
\begin{align}
\label{eqn:t2_and_u2_def}
	t_N^2(X_t) &= \frac{s_N^2(X_t)}{N} &
	&\tand&
	u_N^2(X_t) &= \frac{s_N^2(X_t)}{N} \; \hattauint,
\end{align}
as long as we have an estimator\index{$\hattauint$} $\hattauint$ of $\tauint$.
\end{defn}

%% ================================================================
\section{Estimation of the integrated autocorrelation time}
\label{sec:estimate_tauint}

Recall from remark \ref{rem:c_hat_wild} that $\chat(t)$ is a rather wild
estimator of $c(t)$ at high $t$.  Since
\begin{align*}
	\hattauint &= 1 + 2 \sum_{t=1}^\infty \chat(t)
\end{align*}
is nothing more than a sum of $c(t)$, we can expect it to be ill-behaved as
well.

\begin{figure}[!htb]
\begin{center}
\psfragscanon
\includegraphics[scale=0.8]{figures/tauint_eta_0.9_xmax_1000.eps}
\includegraphics[scale=0.8]{figures/tauint_eta_0.9_xmax_50.eps}
\caption[Estimated and exact integrated autocorrelation times for $Y_t$.]
	{Estimated and exact integrated autocorrelation times for $Y_t$ with
	$\eta=0.9$, using three realizations similar to the one in figure
	\ref{fig:corrunif_autocorr}.  Burn-in iterates have been discarded.  The
	second plot zooms in on the first 50 samples of the first plot.  The
	flat-spot estimator $\hattauint$ of $\tauint$ is found by reading off the
	vertical coordinate of the first turning point of each solid-line plot; the
	true $\tauint$ is the horizontal asymptote of the dotted-line plot.  Two of
	the three turning points yield a $\hattauint$ which is less than the true
	$\tauint$.  This is the general case: we find that $\hattauint$
	underestimates more often than it overestimates.  See also figure
	\ref{fig:results_hattauint} on page \pageref{fig:results_hattauint}.
	\label{fig:corrunif_tauint}}
\end{center}
\end{figure}

\begin{figure}[!htb]
\begin{center}
\psfragscanon
\includegraphics[scale=0.7]{figures/tauints_0.900.eps}
\includegraphics[scale=0.7]{figures/tauints_0.990.eps}
\includegraphics[scale=0.7]{figures/tauints_0.999.eps}
\caption[Estimated integrated autocorrelation times
	for $Y_t$ with varying $\eta$.]
	{Estimated integrated autocorrelation times for $Y_t$ with
	$\eta=0.9,0.99,0.999$, using ten realizations each.
	$N$ is 100,000; burn-in iterates have been discarded.  Recall that true
	$\tauint$ values are $19$, $199$, and $1999$, respectively.  The variation
	in the vertical coordinate of the first flat spot in each plot, which
	increases with $\eta$, gives rise to the error of the error bar on the
	sample mean.
	\label{fig:corrunif_tauints}}
\end{center}
\end{figure}

\begin{defn}
The \emphidx{running-sum estimator} of $\tauint$ is
\begin{align*}
	\hattauint(t) &= 1 + 2 \sum_{k=1}^t \chat(k).
\end{align*}
\end{defn}
The idea is to accumulate the reliable low-$t$ values of $\chat(t)$
until\index{$c(t)$} the sum becomes approximately constant at some $s$, then
stop and declare $\hattauint$ to be $\hattauint(s)$.  This is the
\emphidx{flat-spot estimator} or \emphidx{turning-point estimator} for
$\tauint$.  See figure \ref{fig:corrunif_tauint} for illustration, where $s$ is
approximately 24 for the blue realization and 29 for the red.  From the plots,
we estimate $\tauint \approx 15$; using equation \eqref{eqn:tauint_and_eta}, we
estimate $\eta=(15-1)/(15+1) = 0.875$.  This is reasonable since the data were
obtained with $\eta=0.9$, for which the true $\tauint$ is 19 by equation
\eqref{eqn:tauint_fcn_of_eta}.

It is clear from the figure that estimators $\hattauint$ can vary noticeably
from one realization to the next.  Our estimator for the variance of the sample
mean, i.e. the error bar on the sample mean, is $u_N^2(X_t)$ (equation
\eqref{eqn:t2_and_u2_def}).  Since $\hattauint$ is a factor in $u_N^2(X_t)$,
variations in $\hattauint$ will result in \plainidx{error of the error bar}.
Figure \ref{fig:corrunif_tauints} shows that variations in $\hattauint$
increase with $\eta$.

At present I know of no solution to this problem other than the running of
multiple experiments --- larger $M$, using the notation of section
\ref{sec:prob_stt}.  As long as $\tauint$ is estimated based on a single
experimental result $X_0, \ldots, X_{N-1}$, one must be aware that the error
bars on the sample mean are crude.

\begin{rem}
\label{rem:hattauint_bias}
As was noted in remark \ref{rem:chat_bias}, $\chat(t)$ underestimates $c(t)$.
Since $\hattauint$ is formed from a sum of $\chat(t)$'s, $\hattauint$ is also
a fractional underestimator of $\tauint$, as will be seen in section
\ref{sec:estimate_meanvar}.
\end{rem}

%% ================================================================
\section{Estimation of the variance of the sample mean}
\label{sec:estimate_meanvar}

Given the flat-spot estimator $\hattauint$ of $\tauint$ from section
\ref{sec:estimate_tauint} and the naive estimator of the variance of the sample
mean from equation \eqref{eqn:t2_and_u2_def}, we may now compute the corrected
estimator of the the variance of the sample mean:
\begin{align*}
	u_N^2(X_t) &= \frac{s_N^2(X_t)}{N} \; \hattauint.
\end{align*}

We use the correlated-uniform Markov process to illustrate, since for this
process all quantities have known theoretical values.  As in section
\ref{sec:prob_stt}, we display standard deviations in our plots and tables,
rather than variances:  the units of measurement of the former match those of
the mean, and they correspond visually to variations in the data.

\begin{itemize}

\item The \plainidx{mean} and \plainidx{variance} of $Y_t$ are $\muY = 1/2$ and
$\sigmaY^2 = 1/12$; $\sigmaY \approx 0.289$.  Using $\eta=0.0, 0.9, 0.999$, the
true $\tauint$\index{$\tauint$} is $1, 19, 1999$, respectively.  We conduct
$M=100$ experiments of collecting and analyzing $N=10000$ time-series samples
$Y_0, \ldots, Y_{N-1}$.

\item The true mean is shown in row 1 of table \ref{table:results}.  Estimators
$\YNbar$ are shown in figure \ref{fig:results_YNbar}.  The average of these
over all $M$ experiments is shown in row 2 of table \ref{table:results}.

\item The true naive variance of the sample mean is $\sigmaY^2/N$, with true
naive standard deviation of the sample mean $\sigmaY/\sqrt{N} \approx
0.00289$.  The true corrected variance of the sample mean
is\index{$\sigmaYNbar^2$} $\sigmaYNbar^2 = \tauint \, \sigmaY^2 / N =
1/120000, 19/120000, 1999/120000$.  The true standard deviations of the sample
means are then $\sigmaYNbar$ $\approx$ $0.0028868$, $0.0125831$, $0.1290672$.
These are shown in row 3 of table \ref{table:results}.

\item The \plainidx{multi-experiment estimator} \plainidx{$s_M(\YNbar)$} of
$\sigmaYNbar$ is the sample standard deviation of the $M$ values $\XNbar^{(0)},
\ldots, \XNbar^{(M-1)}$.  These estimators are shown in row 4 of table
\ref{table:results}.  As expected, the multi-experiment estimator is a good
one.

\item Next we turn to \plainidx{single-experiment estimator}s of the variance
of the sample mean.  The estimated naive standard deviation of the sample mean
is $t_N(Y_t) = s_N(Y_t)/\sqrt{N}$.  These are not plotted for each experiment;
their average over all $M$ experiments is shown in row 5 of table
\ref{table:results}.  Note that they match the true variance of the sample mean
only in the IID ($\eta=0$) case.

\item True values of $\tauint$ for each $\eta$ are shown in row 8 of the table.
The flat-spot estimators\index{$\hattauint$} $\hattauint$ for all $M=100$
experiments are shown in figure \ref{fig:results_hattauint}.  Their average and
sample standard deviation over all $M$ experiments are shown in rows 9 and 10.
As discussed in remark \ref{rem:hattauint_bias}, we see that $\hattauint$
fractionally underestimates $\tauint$.

\item Using the $\hattauint$ values, the corrected estimators $u_N(Y_t) =
t_N(Y_t) \sqrt{\hattauint} $ are shown, for all $M=100$ experiments, in figure
\ref{fig:results_uNYt}.  Their averages over all $M$ experiments are shown in
row 6 of table \ref{table:results}.  (Again, the corresponding true values are
in row 2 of the table.) The fractional underestimation of $\hattauint$ carries
over to $u_N(Y_t)$.  One trades the quality of the estimator for the
feasibility of its computation.

\item Standard deviations over $M$ experiments of $u_N(Y_t)$ are shown in row 7
of the table.  Figure \ref{fig:results_errorbar} shows, for $\eta=0.999$, the
$M=100$ values of $\YNbar$ along with their respective $u_N(Y_t)$'s.  These
show the \plainidx{error of the error bar}.

\end{itemize}

\begin{figure}[!htb]
\begin{center}
\psfragscanon
\includegraphics[scale=0.8]{figures/results_YNbar.eps}
\caption[$\YNbar$ over $M=100$ experiments.]
	{$\YNbar$ over $M=100$ experiments, where the true value is $\muX=0.5$.
	Variance of $\YNbar$ increases with autocorrelation factor $\eta$.
	\label{fig:results_YNbar}}
\end{center}
\end{figure}

\begin{figure}[!htb]
\begin{center}
\psfragscanon
\includegraphics[scale=0.8]{figures/results_hattauint.eps}
\caption[$\hattauint(Y_t)$ over $M=100$ experiments, along with true values.]
	{$\hattauint(Y_t)$ over $M=100$ experiments, along with true values.
	Note that $\hattauint(Y_t)$ fractionally underestimates the true
	$\tauint(Y_t)$.
	\label{fig:results_hattauint}}
\end{center}
\end{figure}

\begin{figure}[!htb]
\begin{center}
\psfragscanon
\includegraphics[scale=0.8]{figures/results_uNYt.eps}
\caption
	[$u_N(Y_t)$ over $M=100$ experiments, along with true values.]
	{$u_N(Y_t)$ over $M=100$ experiments, along with true values.
	Note that $u_N(Y_t)$ fractionally underestimates the true standard
	deviation of the sample mean, $\sigma_{\YNbar}=\sigmaY/\sqrt{N}$.
	\label{fig:results_uNYt}}
\end{center}
\end{figure}

\begin{figure}[!htb]
\begin{center}
\psfragscanon
\includegraphics[scale=0.8]{figures/results_errorbar_0.000_sorted.eps}
\includegraphics[scale=0.8]{figures/results_errorbar_0.900_sorted.eps}
\includegraphics[scale=0.8]{figures/results_errorbar_0.999_sorted.eps}

\caption[$\YNbar$ with single-sigma error bars.]
	{$\YNbar$ with single-sigma error bars, $\eta=0, 0.9, 0.999$, $M=100$
	experiments, sorted by increasing $\YNbar$. The magnitude and the variation
	of the error bars both increase with $\eta$.
	\label{fig:results_errorbar}}
\end{center}
\end{figure}

\begin{table}[!htb]
\begin{center}
\begin{tabular}{|p{2.0in}|c||l|l|l|}
\hline
	\multicolumn{1}{|c}{\textbf{Description}} &
	\multicolumn{1}{|c}{$\boldsymbol\eta$} &
	\multicolumn{1}{|c}{\textbf{0}}     &
	\multicolumn{1}{|c}{\textbf{0.9}}   &
	\multicolumn{1}{|c|}{\textbf{0.999}} \\
\hline
%
\hline
	1. True mean &
	$\muYNbar$     & 0.50000 & 0.50000 & 0.50000 \\
\hline
	2. Sample mean &
	$m_M(\YNbar)$  & 0.49948 & 0.49952 & 0.51100 \\
	{} &
	$m_M(\YNbar)$  & 0.49987 & 0.50231 & 0.47341 \\
	{} &
	$m_M(\YNbar)$  & 0.49991 & 0.49895 & 0.47958 \\
\hline
%
\hline
	3. True standard deviation of sample mean &
	$\sigmaYNbar$  & 0.00288 & 0.01258 & 0.12906 \\
\hline
	4. Multi-experiment &
	$s_M(\YNbar)$  & 0.00274 & 0.01166 & 0.10342 \\
	estimator of $\sigmaYNbar$ &
	$s_M(\YNbar)$  & 0.00274 & 0.01167 & 0.12303 \\
	{} &
	$s_M(\YNbar)$  & 0.00298 & 0.00986 & 0.11929 \\
%
\hline
	5. Averaged &
	$m_M(t_N(Y_t))$  & 0.00288 & 0.00287 & 0.00263 \\
	single-experiment naive &
	$m_M(t_N(Y_t))$  & 0.00288 & 0.00288 & 0.00260 \\
	estimators of $\sigmaYNbar$ &
	$m_M(t_N(Y_t))$  & 0.00288 & 0.00287 & 0.00250 \\
\hline
	6. Averaged &
	$m_M(u_N(Y_t))$  & 0.00288 & 0.01279 & 0.09957 \\
	single-experiment corrected &
	$m_M(u_N(Y_t))$  & 0.00289 & 0.01280 & 0.10037 \\
	estimators of $\sigmaYNbar$  &
	$m_M(u_N(Y_t))$  & 0.00289 & 0.01276 & 0.08947 \\
\hline
%
\hline
	7. Sample standard &
	$s_M(u_N(Y_t))$  & 0.00004 & 0.00105 & 0.04402 \\
	deviation of corrected &
	$s_M(u_N(Y_t))$  & 0.00005 & 0.00108 & 0.04665 \\
	estimators of $\sigmaYNbar$ &
	$s_M(u_N(Y_t))$  & 0.00006 & 0.00118 & 0.03851 \\
\hline
%
\hline
	8. True integrated &
	$\tauint$        & 1 & 19 & 1999 \\
	autocorrelation time & & & & \\
\hline
	9. Averages of estimated &
	$m_M(\hattauint)$     & 0.999 & 19.854 & 1442.627 \\
	integrated autocorrelation &
	$m_M(\hattauint)$     & 1.002 & 19.763 & 1500.137 \\
	time &
	$m_M(\hattauint)$     & 1.008 & 19.857 & 1279.173 \\
%
\hline
	10. Standard deviation &
	$s_M(\hattauint)$ & 0.028 & 3.162 &  865.992 \\
	across $M$ experiments &
	$s_M(\hattauint)$ & 0.031 & 3.092 & 1045.842 \\
	of $\hattauint$ &
	$s_M(\hattauint)$ & 0.039 & 3.628 &  758.903 \\
\hline
\end{tabular}
\caption[Statistics for 3 trials of $100$ experiments on $10000$
	samples of $Y_t$.]
	{Statistics for three trials of $M=100$ experiments on $N=10000$
	samples of $Y_t$: $\eta=0.0, 0.9, 0.999$.
	\label{table:results}}
\end{center}
\end{table}

%% ================================================================
\section{Integrated and exponential autocorrelation times}

In remark \ref{rem:tauexp} of section \ref{sec:c_hat}, we noted that if
$\Corr(X_0, X_t) = \eta^t$ for $\eta\in[0,1)$, then we may define an
\plainidx{exponential autocorrelation time} via
\begin{align*}
	\tauexp = -1/\log\eta
\end{align*}
such that $\Corr(X_0, X_t) = \exp(-t/\tauexp)$.  Yet section \ref{sec:tauint}
gave us something similar: the \plainidx{integrated autocorrelation time}
$\tauint$.  In particular, if $\Corr(X_0,X_t) = \eta^t$, then we had
\begin{align*}
	\tauint &= \frac{1+\eta}{1-\eta}.
\end{align*}
Figure \ref{fig:cmp_taus} compares these two.

\begin{figure}[!htb]
\begin{center}
\psfragscanon
\includegraphics[scale=0.8]{figures/cmp_taus.eps}
\includegraphics[scale=0.8]{figures/cmp_taus_log.eps}
\caption[Integrated and exponential autocorrelation times.]
	{Integrated and exponential autocorrelation times as a function of
	$\eta$.
	\label{fig:cmp_taus}}
\end{center}
\end{figure}

%% ================================================================
\section{Batched means}
\label{sec:batched_means}

Introductory statistics tends to deal with the analysis of \plainidx{IID
samples}.  Yet, realization sequences from an MCMC experiment tend to be highly
correlated.  The sample mean estimates the true mean, since expectation is
linear.  But when one wishes to place an accurate error bar on the sample mean,
correlations must be taken into account.

One approach (see for example \cite{bib:berg}, who calls this process
\emphidx{binning}) is to subdivide $X_0, \ldots, X_{N-1}$ into $K=N/B$
batches\index{batched mean} of size $B$.  The $K$ sample means over batches
may be treated as IID samples.  The independence of the $K$ samples means that
the variance of their sample mean will be reduced, but reducing the sample size
from $N$ to $K$ will increase the variance.  We will show that these two
effects cancel:  binning $N$ samples down to $K$ samples does not change the
variance of the sample mean.  (As shown in \cite{bib:berg}, batched means have
their uses: they may be used to construct a method to estimate $\tauint$, as an
alternative to the method of section \ref{sec:estimate_tauint}.)

\begin{defn}
Given $X_0, \ldots, X_{N-1}$ with common mean $\muX$ and variance
$\sigmaX^2$, let $B$ divide $N$ and $K=N/B$.  Then $B$ is the \emphidx{batch
size} and $K$ is the \emphidx{number of batches}.  For $k=0, \ldots, K-1$, the
$k$th batch consists of $X_{kB}, \ldots, X_{(k+1)B-1}$.  The sample mean of the
$k$th batch is
\begin{align*}
	A_k &= \frac{1}{B} \sumiB X_{kB+i}.
\end{align*}
\end{defn}
We now consider the sequence $A_0, \ldots, A_{K-1}$.  We define the
\emphidx{batched mean} to be
\begin{align*}
	\XNBbar &= \frac{1}{K} \sum_{k=0}^{K-1} A_k.
\end{align*}
By linearity of expectation, we immediately have $\bbE[\XNBbar] = \muX$.   We
next inquire about the variance of the batched mean, then compare that to the
variance of the (non-batched) sample mean.

%% ================================================================
\section{Variance and covariance of batches}

To compute $\Var(A_k)$ and $\Corr(A_0, A_k)$, we first need $\bbE[A_k A_\ell]$
for $k=\ell$ and $k\ne\ell$.  In the $k=\ell$ case, the computation is the same
as in section \ref{sec:Var_XNbar}, with $B$ playing the role of $N$.  We have
\begin{align*}
	\bbE[A_k^2] &\approx \muX^2 +
		\frac{\sigmaX^2}{B}\left(\frac{1+\eta}{1-\eta}\right) &
	&\tand&
	\Var(A_k) &\approx
		\frac{\sigmaX^2}{B}\left(\frac{1+\eta}{1-\eta}\right).
\end{align*}
For $k\ne \ell$, without loss of generality assume $k<\ell$.  Using
equation \eqref{eqn:E_Xi_Xj}, we have
\begin{align*}
	\bbE[A_k A_\ell]
	&= \frac{1}{B^2}
		\sumiB \sumjB \bbE[X_{kB+i} X_{\ell B+j}]
	\qquad = \frac{1}{B^2}
		\sumiB \sumjB
			(\muX^2 + \sigmaX^2 \eta^{\ell B+j-kB-i}) \\
	&= \muX^2 + \frac{\sigmaX^2 \eta^{(\ell-k)B}}{B^2}
		\sumiB \eta^{-i} \sumjB \eta^j
	\qquad = \muX^2 + \frac{\sigmaX^2 \eta^{(\ell-k)B}}{B^2}
		\left(\frac{1-\eta^B}{1-\eta}\right)^2.
\end{align*}
If the batch size is chosen so that $\eta^B$ is negligible, then
\begin{align*}
	\bbE[A_k A_\ell]
	&= \muX^2.
\end{align*}

Now we have (for $\eta^B \approx 0$)
\begin{align}
\label{eqn:Var_Ak_Corr_Ak_Aell}
	\Var(A_k)
		%&= \bbE[A_k^2] - \bbE[A_k]^2
		&\approx \frac{\sigmaX^2}{B}\left(\frac{1+\eta}{1-\eta}\right) &
	&\tand&
	\Corr(A_k, A_\ell) &= \frac{\bbE[A_k A_\ell] - \muX^2}{\sigma_{A_0}^2}
	= \delta_{k,\ell}.
\end{align}
This justifies the hope that batches can be constructed to form an IID
sequence.

%% ================================================================
\section{Variance of the batched mean}
\label{sec:cmp_bm_vars}

We now find out what effect batching has on the variance of the sample mean:
batching produces an IID sequence, which will reduce the variance (equation
\eqref{eqn:Var_XNbar_tauint}), yet it reduces the sample size from $N$ down to
$K=N/B$, which by central-limit reasoning should increase the variance.

For the non-batched mean, we have the random variables $X_0, \ldots, X_{N-1}$;
parameters are mean $\muX$, variance $\sigmaX^2$, autocorrelation $\eta^k$, and
(from equation \eqref{eqn:tauint_and_eta}) integrated autocorrelation time
$\tauint = (1+\eta)/(1-\eta)$.  Equation \eqref{eqn:Var_XNbar_eta_ratio} gives
\begin{align}
\label{eqn:Var_XNbar}
	\Var(\XNbar) &= \frac{\sigmaX^2}{N} \; \left(\frac{1+\eta}{1-\eta}\right).
\end{align}

For the batched mean, we batch $X_0, \ldots, X_{N-1}$ into $K$ IID batches of
size $B$.  We have the random variables $A_0, \ldots, A_{K-1}$, with mean
$\muX$, variance $(\sigmaX^2/B) (1+\eta)/(1-\eta)$ (equation
\eqref{eqn:Var_Ak_Corr_Ak_Aell}), autocorrelation $c(k) = \delta_{0,k}$ (since
$A_k$ is IID) and integrated autocorrelation time $\tauint = 1$ (remark
\ref{rem:tauint_IID}).  Then
\begin{align*}
	\Var(\XNBbar)
		&= \frac{\sigmaX^2}{KB}\left(\frac{1+\eta}{1-\eta}\right)
		= \frac{\sigmaX^2}{N}\left(\frac{1+\eta}{1-\eta}\right).
\end{align*}
Thus, to first order in $\eta$ and $N$, as long as $B$ is large enough that
$\eta^B$ is negligibly small, we do not expect batching to change the variance
of the sample mean.

Table \ref{table:batch_N_M} shows some sample results of these calculations for
the correlated-uniform Markov process $Y_t$.  There are $M=100$ experiments of
$N=10000$ samples.  Each experiment was analyzed as-is ($B=1$), as well as with
batch size $B$ = 64, 512, and 4096.  (Recall from table
\ref{table:eta_vs_tauint} that $\eta=0, 0.9, 0.999$ correspond to $\tauint=1,
19, 1999$, respectively.)  Now the process being analyzed is $A_0, \ldots,
A_{K-1}$ where $K=N/B$.  We note the following:

\begin{itemize}

\item For $B=1$ (the original time series), the estimator $u^2$ of the variance
of the sample mean employed was the corrected estimator of equation
\eqref{eqn:t2_and_u2_def}, while for $B=64, 512, 4096$, the $A_k$'s were treated
as if they were IID.  That is, for $B>1$ we set $u^2=t^2$.

\item Batch size does not, of course, affect the sample mean (column 3 of the
table).  Likewise, it does not affect the multi-experiment estimator of the
variance of the sample mean (column 4).

\item The true variance of the sample mean, $\sigma_{\ANbar} = \sqrt{\tauint
\sigma^2_{A_k} / K}$, is shown in column 5.

\item The last two columns show the first two autocorrelation estimates.  These
show that for $\eta=0$ (the IID case), $Y_t$ samples are indeed approximately
IID.  For $\eta=0.9$, the $B=64$ batches are nearly independent, and the
largest batches are quite weakly correlated.  For $\eta=0.999$, where
$\tauint=1999$, batch sizes of 64 and 512 are too small, but batch size 4096 is
large enough to produce weakly correlated batches.

\item For $\eta=0$, the average of the single-experiment estimator $u$ of the
variance of the sample mean is approximately constant with respect to batch
size.  For $\eta=0.9$ and $\eta=0.999$, once the batch size is large enough to
get weakly correlated samples, the estimator $u^2$ on batches agrees with the
estimator $u^2$ on the original time-series data.

\item The multi-realization estimator and the averaged single-realization
estimator of the variance of the sample mean (columns 4 and 6) roughly agree,
for batch sizes large enough that batches are weakly correlated.

\item The error of the error bar (column 7 of the table) is not improved by use
of batched means.

\end{itemize}

\begin{table}[!htb]
\begin{center}
\begin{footnotesize}
\begin{tabular}{|c|c|c|c|c|c|c|c|c|}
\hline {} & {} & {} & {} & {} & {} & {} & {} & {} \\
	   %
       $\eta$ &
	   $B$ &
	   $m_M(\ANbar)$ &
	   $s_M(\ANbar)$ &
	   $\sigma_{\ANbar}$ &
	   $m_M(u_N(A_k))$ &
	   $s_M(u_N(A_k))$ &
	   $\chat_{A_k}(0)$ &
	   $\chat_{A_k}(1)$ \\
	   %
       {} & {} & {} & {} & {} & {} & {} & {} & {} \\
\hline
%
\hline 0.000 &    1 & 0.5001 & 0.00113 & 0.00113 & 0.00113 & 0.000002 & 1.0000 &  0.0001 \\
\hline 0.000 &   64 & 0.5001 & 0.00113 & 0.00113 & 0.00113 & 0.000027 & 0.9990 &  0.0378 \\
\hline 0.000 &  512 & 0.5001 & 0.00113 & 0.00113 & 0.00114 & 0.000073 & 0.9922 &  0.0031 \\
\hline 0.000 & 4096 & 0.5001 & 0.00113 & 0.00113 & 0.00111 & 0.000206 & 0.9375 &  0.2053 \\
\hline
\hline 0.900 &    1 & 0.4996 & 0.00481 & 0.00492 & 0.00490 & 0.000034 & 1.0000 &  0.8982 \\
\hline 0.900 &   64 & 0.4996 & 0.00481 & 0.00492 & 0.00451 & 0.000103 & 0.9990 &  0.1591 \\
\hline 0.900 &  512 & 0.4996 & 0.00481 & 0.00492 & 0.00480 & 0.000280 & 0.9922 & -0.0552 \\
\hline 0.900 & 4096 & 0.4996 & 0.00481 & 0.00492 & 0.00476 & 0.000894 & 0.9375 &  0.1445 \\
\hline
\hline 0.999 &    1 & 0.5112 & 0.04801 & 0.05042 & 0.04937 & 0.004180 & 1.0000 &  0.9987 \\
\hline 0.999 &   64 & 0.5112 & 0.04801 & 0.05042 & 0.00874 & 0.000756 & 0.9990 &  0.9493 \\
\hline 0.999 &  512 & 0.5112 & 0.04801 & 0.05042 & 0.02301 & 0.002276 & 0.9922 &  0.6462 \\
\hline 0.999 & 4096 & 0.5112 & 0.04801 & 0.05042 & 0.04280 & 0.007354 & 0.9375 & -0.0060 \\
\hline
\end{tabular}
\caption[Statistics for $100$ batched experiments on $65536$ samples of
	$Y_t$.]
	{Statistics for $M=100$ batched experiments on $N=65536$ samples of
	$Y_t$: $\eta=0.0, 0.9, 0.999$.
	%xxx insert some conclusions here.
	\label{table:batch_N_M}}
\end{footnotesize}
\end{center}
\end{table}

%% ----------------------------------------------------------------
\section{Conclusions on error bars, autocorrelation, and batched means}
\label{sec:autocorr_conclusions}

Given an MCMC experimental result $X_0, \ldots, X_{N-1}$, we may compute the
sample mean $\XNbar$ and an estimator $s_N^2(X_t)$ of the sample variance.

Batched means improve neither the bias nor the variation of the error bar.  The
variance reduction obtained by (approximate) independence of batches cancels
out the variance increase caused by reduced sample size.

Computing autocorrelations and summing them as described in section
\ref{sec:estimate_tauint}, we may obtain an estimate $\hattauint$ of the
integrated autocorrelation time $\tauint$.  This is used to update the
naive\index{naive estimator} estimated variance of the sample mean $t_N^2(X_t)
= s^2/N$ to the \plainidx{corrected estimator} $u_N^2(X_t) = \hattauint s^2/N$.
With the understanding that $\hattauint$ has itself a noticeable variance and a
fractional underbias, $u_N$ estimates the standard deviation of the sample
mean.
