I am trying to evaluate an adaptive MCMC algorithm on a multimodal target density. Among other performance measures, I would like to evaluate the sampler in terms of Effective Sample Size (ESS). The problem is that, for highly multimodal densities, the standard ESS used in MCMC is not useful.

For instance:

effectiveSize( c(rnorm(1e3), rnorm(1e3, 10)) )
#  var1 
# 1.79313 

Here the ESS is so low because the autocorrelation decays very slowly with the lag. This is because the first 10^3 samples come from the first mixture component, the next $10^3$ from the second. However, one would like the ESS to be close to 2000 in this example.

Hence my question is: are there ESS measures for MCMC that work when the target is multimodal? At the moment my solution is to cluster the output chain, and to sum the ESSs calculated on each clustered subchain.


EDIT: I was recently talking to someone about this problem and they pointed out that we should not expect the ESS to be 2000 here. The OP has introduced a lag ordering in the sample that causes the sample to be no longer independent. Even though the sample was obtained randomly, dependence is introduced when you order the sample in such a way. This is because ACF is calculated for observations sequenced in time. Thus, in calculating ESS, it is assumed that the samples are sequenced in time, and under this assumption, there is clear dependence. If it is truly a random sample then the samples should be jumbled as in sample(c(rnorm(1e3), rnorm(1e3, 10)), and then we should expect the effective sample size to be 2000.

(+1) for a great question. I don't have a satisfying answer for you, but I will give you alternatives. The R package mcmcse provides two other methods to estimate Effective Sample Size. Their methods are based on consistent estimators, so in theory should work better.

> library(mcmcse)
> x <- c(rnorm(1e3), rnorm(1e3, 10))
> ess(x)

47 effective samples is still not great. As long as you have such high lags, both coda and mcmcse will give you small ESS. This is because ESS is based on the calculation of lag covariance.

The idea of splitting up the Markov chain is reasonable. Although you might want to split it up on the basis of large jumps of the Adaptive chain. Instead of doing ESS separately in that case, you might want to do a multivariate ESS, so the dependence across the modes is also captured.

> x1 <- x[1:1000]
> x2 <- x[1001:2000]
> multiESS(cbind(x1,x2))
[1] 885.9462

You can then multiply this number by 2. However, this will only work if you split it up equally.

Another alternative is to change the estimation quantity. ESS is different for different functions that you want to estimate. So for example, having obtained an MCMC output, if you want to estimate the mean versus estimating the second moment, the ESS for both of those quantities will be different by definition. So (even though you might want the mean), you can argue that ESS does not work well for multi modal densities, and so for the purposes of calculating ESS, we change the function of interest.

> x <- c(rnorm(1e3), rnorm(1e3, 10))
> y <- x^(50)
> ess(y)
> effectiveSize(y)

So here you calculated the ESS for estimating the 50th moment. I had to play around with it till it gave a good answer. But I believe the moment that would work best depends on the heights and distances between the modes of the density.


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.