tl;dr I can only generate samples of a random variable $X$ using MCMC. How can I find the error in the estimate of the expected value of $X$ based on this MCMC data?

The problem

I have a "black box" which gives me a sequence of "random" numbers (MCMC sampler). I have recorded $n$ numbers from this black box, which is all the data I can work with (I can't sample additional points). Now I need to estimate the average of the black box based on this data, and get some idea about how precise the estimate is.

In other words, if I simply calculate the average of these $n$ datapoints, how far is it going to be from the true average of the black box?

The difficulty: Consecutive numbers coming form the black box are strongly correlated, however the $i^\text{th}$ and $(i+k)^\text{th}$ numbers are effectively uncorrelated if $k$ is large.

Are there any methods to deal with this type of situation, and find the error in the estimate? Let's assume that $n$ is large enough that I have some uncorreltaed data points.

What I'm doing now

In the following I'll describe what I am doing right now and why I am not satisfied with this method. This method works but it is not ideal.

Let's look at the case when the numbers $x$ are uncorrelated. Let $\langle x \rangle_n$ denote the average of $n$ consecutive numbers, which is an estimate of the expected value of the black box. Then $\operatorname{Var}(\langle x \rangle_n) = \operatorname{Var}(x)/n$.

$\sqrt{\operatorname{Var}(\langle x \rangle_n)}$ can be used as the error of the estimate $\langle x \rangle_n)$ and it's easy to compute based on $\operatorname{Var}(x)$.

However my actual data is correlated. So I partitioned my dataset of $n$ elements into runs of $k$ elements and calculated the product $k \times \operatorname{Var}(\langle x \rangle_k)$ for various $k$. When $k$ is large enough, and the correlations are not important any more, this product will settle at a constant value, which can be used to extrapolate for $\operatorname{Var}(\langle x \rangle_n)$. This is a plot of $k \operatorname{Var}(\langle x \rangle_k)$ for various values of $k$ for four different datasets (four colours) of size $n = 2^{27}$:

enter image description here

You can see that as $k$ becomes large enough, the product $k \operatorname{Var}(\langle x \rangle_k)$ settles on a constant. But as $k$ becomes too large, I can only parition my length-$n$ dataset into a few parts of length $k$, so the estimate of the variance becomes inaccurate. Hence the noise on the right of the curve.

So I need to find the flat part of the curve and discard the noisy part, then estimate the variance of the average of the full dataset ($\langle x \rangle_n$) based on it.

Why my method isn't ideal

The method I describe above is complicated and involves several steps, and several choices in each step. For how many (and what magnitude of) $k$ values should I compute the variances? How do I determine where is the flat part of the curve? How do I decide which part is too noisy to keep? All of these are separate problems on their own with many valid heuristic solutions, which will all influence the result.

Also, this method isn't very efficient computationally.

I am hoping for a simpler and cleaner robust method where the details won't have such a large effect on the end result.

  • $\begingroup$ It seems to me that you are describing an autocorrelated dataset. Perhaps if you go looking with that keyword it will help you find an answer you feel is suitable. An off the cuff approach might be to use the K previous values as predictors of your K+1 value, the residual would have drawn away most of the variance that is associated with the non-independence of values and yield a reasonable estimate of variance without dependence. $\endgroup$ – russellpierce May 20 '14 at 17:40

This is an active area of research. This first question is whether a central limit theorem (CLT) even exists which depends on the mixing properties of your MCMC, e.g. geometric convergence. Typically, this is a nontrivial question.

Provided a CLT exists, the second question is how to obtain a consistent estimator of the variance in the CLT. Partitioning the chain as you have done is called batching and is the approach used in the R package mcmcse. I'd suggest you take a look at that package and the references therein. As an example of using the package you can do

n = 1000
y = rep(0,n)
for (i in 2:n) y[i] = rnorm(1,0.9*y[i-1])

which will return an estimate of the expectation of y, which is known to be zero in this case, as well as an estimate of the Monte Carlo standard error in the estimator which takes into account the correlation in the chain.

This package allows the user to specify the batch size, but defaults to the square root of the total number of samples. Perhaps the references in the package will give an indication of why this is the default or what other batch sizes will work.

| cite | improve this answer | |
  • $\begingroup$ If this approach is aready being used, then I guess it's not so bad. Are there other common approaches? $\endgroup$ – Szabolcs May 21 '14 at 22:03

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.