Monte Carlo Error Bars for Markov Chain Summaries

One Trajectory, One Replicate

A simulated Markov chain path is usually dependent over time: the state at one step can influence the next state. That dependence does not prevent us from estimating uncertainty for a summary of the path. The key is to regard each independently generated trajectory as one replicate, calculate the summary once for that trajectory, and estimate variability across the replicate summaries. The number of steps in one path is not the number of independent replicates.

Let \((\Omega,\mathcal{F},\mathbb{P})\) be a probability space, and let \((X_n)_{n\ge0}\) be a discrete-time Markov chain on a finite state space \(S\), with transition matrix \(P=(p_{ij})\). Its natural filtration is \(\mathcal{F}_n=\sigma(X_0,\ldots,X_n)\). Fix a finite path length \(m\), and let \(g:S^{m+1}\to\mathbb{R}\) be a summary of the path. For example, \(g\) might be the fraction of steps spent in a chosen state or the total accumulated reward. The summary from one path is the random variable \(Y=g(X_0,\ldots,X_m)\).

Now simulate \(R\) independent trajectories using the same initial-state rule and transition matrix. Write \(Y_r\) for the summary from trajectory \(r\), where \(1\le r\le R\). The trajectories are generated independently of one another and have the same simulation rule, so the summaries \(Y_1,\ldots,Y_R\) are independent and identically distributed. The quantity of interest is their common expectation \(\theta=\mathbb{E}[Y_r]\), estimated by their average.

These formulas estimate uncertainty in \(\bar{Y}\) as an estimate of \(\theta\). They do not measure how well one particular trajectory represents every feature of a chain, and they do not remove bias caused by choosing a path length or initial-state rule that does not match the target quantity. First specify the target, then make sure each replicate is generated under the intended rule.

Why the Replicate Standard Error Works

Proof. Linearity of expectation gives \(\mathbb{E}[\bar{Y}]=R^{-1}\sum_{r=1}^R\mathbb{E}[Y_r]=\theta\). Independence makes the covariance between distinct replicates zero, so \(\operatorname{Var}(\bar{Y})=R^{-2}\sum_{r=1}^R\operatorname{Var}(Y_r)=R^{-2}R\sigma^2=\sigma^2/R\). For the sample variance, use the identity \(\sum_r(Y_r-\bar{Y})^2=\sum_r(Y_r-\theta)^2-R(\bar{Y}-\theta)^2\). Taking expectations, the first term on the right has expectation \(R\sigma^2\); the second has expectation \(R\operatorname{Var}(\bar{Y})=\sigma^2\). Thus \(\mathbb{E}[\sum_r(Y_r-\bar{Y})^2]=(R-1)\sigma^2\), and dividing by \(R-1\) gives \(\mathbb{E}[s^2]=\sigma^2\). Therefore \(s/\sqrt{R}\) estimates \(\sqrt{\operatorname{Var}(\bar{Y})}\). \(\square\)

The theorem explains the role of the number of replicates: when the replicate variance stays fixed, the standard deviation of the average decreases in proportion to \(1/\sqrt{R}\). Increasing the number of steps within a path can change the summary's distribution and variance, but it does not by itself create more independent replicates.

Confidence Intervals for an Expected Summary

A common approximate 95% confidence interval for \(\theta\), when \(R\) is sufficiently large, is \(\bar{Y}\pm1.96\,s/\sqrt{R}\). This uses the standard independent, identically distributed central limit theorem: with finite variance, the standardized replicate mean approaches a standard normal distribution as \(R\) increases. The theorem is a standard result; its proof is beyond this tutorial. Replacing \(\sigma\) by \(s\) gives the usual estimated-standard-error interval. Its coverage is approximate for general replicate distributions, especially when \(R\) is small or the distribution is strongly skewed. If the replicate summaries are normally distributed, the interval using the appropriate Student \(t\) quantile is exact for every \(R\ge2\); normality should not be assumed merely because the chain has many steps.

A confidence interval describes the repeated-simulation performance of the interval procedure: under its assumptions, intervals constructed from repeated sets of independent trajectories cover \(\theta\) at the stated rate (exactly or approximately as specified). It does not mean that 95% of the path summaries lie inside this interval. Nor does one realized interval have a 95% probability of containing a fixed \(\theta\); after the data are observed, it either contains that value or it does not.

Why Steps Within a Path Cannot Be Counted as Replicates

To see the danger of pooling steps, consider the state space \(S=\{0,1\}\) and the transition matrix \(P=\begin{bmatrix}1&0\\0&1\end{bmatrix}\). Each row sums to \(1\). Give each independently generated path an initial state with probability \(\tfrac12\) at each state. The chain remains at that state forever. Although different paths are independent, states within any one path are perfectly dependent.

This example has extreme dependence, but the general lesson is not limited to chains that never move. Adjacent states in a slowly mixing chain can be strongly related. A formula that treats every step as independent can then give a misleadingly small error bar. The replicate method instead needs independence across trajectories; it makes no assumption that the states within a trajectory are independent.

A Small-Replicate Interval and Its Limits

The usual normal approximation requires independent, identically distributed replicate summaries with a finite variance. If a fixed-length summary is computed from a finite-state path and takes bounded values, its variance is finite. Unbounded summaries, such as some stopping-time quantities, require more care. For example, independent summaries with Pareto tail probability \(\mathbb{P}(Y>y)=y^{-3/2}\) for \(y\ge1\) have a finite mean but infinite variance. In that case a finite population standard deviation \(\sigma\) does not exist, so the finite-variance standard error formula and its usual central-limit justification do not apply.

Independence across replicates is also essential to the variance calculation. The identity-chain example shows exactly what can go wrong when dependent states are counted as independent observations. If the simulations share random inputs or are otherwise dependent across runs, the replicate variance formula needs additional justification; it cannot be assumed solely because the outputs are labelled as separate runs.

Practical Checks Before Reporting Error Bars

Before reporting a Monte Carlo standard error or interval, check the target, the unit of replication, and the assumptions behind the calculation.

• Define the target. Say what one path summary measures, its path length, and its initial-state rule.
• Use independent runs. Generate each trajectory independently under the same simulation rule when using the formulas here.
• Calculate one summary per path. Compute \(Y_r\) first; then use the \(R\) values to calculate \(\bar{Y}\), \(s^2\), and \(s/\sqrt{R}\).
• Check the interval's status. Distinguish a large-sample normal approximation from an exact interval under a normal-summary assumption.
• Report the number of replicates. The uncertainty estimate depends on \(R\), not simply on how many states were generated in total.

Common mistake: If a path contains many dependent states, substituting the total number of states for \(R\) in \(s/\sqrt{R}\) treats those states as independent replicates. In the identity-chain example with \(R=2\) and \(L=3\), this uses six observations and predicts variance \(\tfrac1{24}\), although the correct variance is \(\tfrac18\). The correction is to form one path average from each trajectory and use the two independent path averages as the replicates.