Applying multiESS for multiple (dependent) parallel chains I'm using the affine-invariant sampler from emcee to draw samples from a $p$ dimensional posterior, using $M$ parallel chains ($M>10$).
Since my model is p-dimensional with $p>1$, I'm also using the multiESS method to estimate the effective sample size.
My question is: since I have $M$ chains, how should I apply multiESS? Should I apply it over each chain separately? Should I concatenate all the chains and apply it on the single flat chain? Is one method more appropriate than the other?
One thing to keep in mind (I'm not sure if this will affect the decision) is that in the affine-invariant sampler the chains are not independent.
 A: I am using notation from the Goodman and Weare (2010) paper. This is how they describe their estimation problem in Section 3.
There are $L$ walkers in the ensemble, and each at each time point $t$, what is observed is $X(t) = (X_1(t), \dots, X_L(t) )$, for $t = 1, \dots, T_e$. The final estimator for $E[f(x)]$ is,
$$\hat{A}_e = \dfrac{1}{T_e} \sum_{t=1}^{T_e} \left( \dfrac{1}{L} \sum_{k=1}^{L} f(X_k(t)) \right) \,.$$
If the walkers were all independent, then you could calculate multiESS and ess for each walker and then add them all together (as long as there was proper burn-in. Burn-in in this case more important than burn in for a single long run).
However, the walkers are not independent, so they define
$$F(X(t)) = \dfrac{1}{L} \sum_{k=1}^{L} f(X_k(t)) \, \quad \text{ and thus } \quad\hat{A}_e = \dfrac{1}{T_e}\sum_{t=1}^{T_e} F(X_{t})\,  $$
That is, $F(X(t))$ is the average of $f(x)$ obtained from $L$ walkers at time $t$. 
Now here we make the assumption that $F(X(t))$ has a Markovian structure (or atleast a $\phi$-mixing structure). I can't confirm this assumption because I don't quite understand what the sampler is doing in the paper, but they make a similar assumption in Section 3 when they define their $\tau_e$.
Assuming this, $F(X_1), F(X_2), \dots, F(X_{T_e})$ each has variance $\Omega/L$ and exhibit structured correlation so that for the estimator $\hat{A}_e$, so we have the variance
$$\text{Var}(\hat{A}_e) = \dfrac{\Sigma}{T_e} =  \dfrac{L\Sigma}{LT_e}$$
Compare this with the estimator for if all observations were IID (let $\Lambda = \text{Var}_{\pi} f(x)$),
$$\text{Var}(\hat{A}_s) = \dfrac{\Lambda}{LT_e} $$
So the multivariate effective sample size is then,
$$mESS = LT_e \left(\dfrac{\det(\Lambda)}{\det(L\Sigma)} \right)^{1/p}\,. $$
Note that to implement this using multiESS in mcmcse, you must calculate multi.mcse on $F(X(t))$ and not on $f(X(t))$, but to calculate $\Lambda$ you should concatenate all the data, and then calculate the sample covariance. To avoid confusions, I would suggest, don't use the function multiESS, and calculate using the determinant function in R. The code should be similar to the following.
library(mcmcse)
L <- 3
Te <- 1e3
walk1 <- matrix(rnorm(Te*4), nrow = 1000, ncol = 4)
walk2 <- matrix(rnorm(Te*4), nrow = 1000, ncol = 4)
walk3 <- matrix(rnorm(Te*4), nrow = 1000, ncol = 4)

Fx <- (walk1 + walk2 + walk3)/L
Sigma <- mcse.multi(Fx, size = "cuberoot")$cov

concat_full <- rbind(walk1, walk2, walk3)
Lambda <- var(concat_full)

mess <- (L*Te) * (det(Lambda)/det(L*Sigma) )^(1/4)

