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)