I’m trying to compare n_eff computation described in BDA3 and the implemented code, and I think there is a small difference. The comment in the beginning says the implementation matches BDA3. BDA3 (at least version 20150505) describes truncation of the autocorrelation series using Geyer’s initial positive sequence (Geyer, 1992) where truncation “T is the first odd positive integer for which $\widehat{\rho}*{T+1}+\widehat{\rho}*{T+2}$ is negative.”. If I understand the code (stan/src/stan/mcmc/chains.hpp) it’s using “T is the last positive integer for which \widehat{\rho}_{T} is positive.”, that is, it’s not using the sum of odd and even lags for which Geyer had a theoretical argument.

Ping @syclik (he wrote the code)

```
/**
* Returns the effective sample size for the specified parameter
* across all kept samples.
*
* The implementation matches BDA3's effective size description.
*
* Current implementation takes the minimum number of samples
* across chains as the number of samples per chain.
*
* @param VectorXd
* @param Dynamic
* @param samples
*
* @return
*/
double effective_sample_size(const Eigen::Matrix<Eigen::VectorXd,
Dynamic, 1> &samples) const {
int chains = samples.size();
// need to generalize to each jagged samples per chain
int n_samples = samples(0).size();
for (int chain = 1; chain < chains; chain++) {
n_samples = std::min(n_samples,
static_cast<int>(samples(chain).size()));
}
Eigen::Matrix<Eigen::VectorXd, Dynamic, 1> acov(chains);
for (int chain = 0; chain < chains; chain++) {
acov(chain) = autocovariance(samples(chain));
}
Eigen::VectorXd chain_mean(chains);
Eigen::VectorXd chain_var(chains);
for (int chain = 0; chain < chains; chain++) {
double n_kept_samples = num_kept_samples(chain);
chain_mean(chain) = mean(samples(chain));
chain_var(chain) = acov(chain)(0)*n_kept_samples/(n_kept_samples-1);
}
double mean_var = mean(chain_var);
double var_plus = mean_var*(n_samples-1)/n_samples;
if (chains > 1)
var_plus += variance(chain_mean);
Eigen::VectorXd rho_hat_t(n_samples);
rho_hat_t.setZero();
double rho_hat = 0;
int max_t = 0;
for (int t = 1; (t < n_samples && rho_hat >= 0); t++) {
Eigen::VectorXd acov_t(chains);
for (int chain = 0; chain < chains; chain++) {
acov_t(chain) = acov(chain)(t);
}
rho_hat = 1 - (mean_var - mean(acov_t)) / var_plus;
if (rho_hat >= 0)
rho_hat_t(t) = rho_hat;
max_t = t;
}
double ess = chains * n_samples;
if (max_t > 1) {
ess /= 1 + 2 * rho_hat_t.sum();
}
return ess;
}
```