Skip to content

Commit

Permalink
basic_ess needed for mcse
Browse files Browse the repository at this point in the history
  • Loading branch information
mitzimorris committed Oct 11, 2024
1 parent 7946efe commit 35bdb07
Show file tree
Hide file tree
Showing 3 changed files with 28 additions and 5 deletions.
27 changes: 25 additions & 2 deletions src/stan/analyze/mcmc/split_rank_normalized_ess.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -103,15 +103,16 @@ double ess(const Eigen::MatrixXd& chains) {
/**
* Computes the split effective sample size (split ESS) using rank based
* diagnostic for a set of per-chain draws. Based on paper
* https://arxiv.org/abs/1903.08008
* https://arxiv.org/abs/1903.08008 Computes bulk ESS over entire sample,
* and tail ESS over the 0.05 and 0.95 quantiles.
*
* When the number of total draws N is odd, the last draw is ignored.
*
* See more details in Stan reference manual section "Potential
* Scale Reduction". http://mc-stan.org/users/documentation
* @param chains matrix of per-chain draws, num_iters X chain
* @return potential scale reduction
* @return pair ESS_bulk, ESS_tail
*/
inline std::pair<double, double> split_rank_normalized_ess(
const Eigen::MatrixXd& chains) {
Expand Down Expand Up @@ -142,6 +143,28 @@ inline std::pair<double, double> split_rank_normalized_ess(
return std::make_pair(ess_bulk, ess_tail);
}

/**
* Computes the split effective sample size (split ESS)
* diagnostic for a set of per-chain draws.
*
* When the number of total draws N is odd, the last draw is ignored.
*
* See more details in Stan reference manual section "Potential
* Scale Reduction". http://mc-stan.org/users/documentation
* @param chains matrix of per-chain draws, num_iters X chain
* @return potential scale reduction
*/
inline double split_basic_ess(
const Eigen::MatrixXd& chains) {
Eigen::MatrixXd split_draws_matrix = split_chains(chains);
if (!is_finite_and_varies(split_draws_matrix)
|| split_draws_matrix.rows() < 4) {
return std::numeric_limits<double>::quiet_NaN();
}
return ess(split_draws_matrix);
}

} // namespace analyze
} // namespace stan

Expand Down
2 changes: 1 addition & 1 deletion src/stan/io/stan_csv_reader.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -377,7 +377,7 @@ class stan_csv_reader {

if (!read_samples(in, data.samples, data.timing)) {
if (out)
*out << "Unable to parse sample" << std::endl;
*out << "no draws found" << std::endl;
}
return data;
}
Expand Down
4 changes: 2 additions & 2 deletions src/stan/mcmc/chainset.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -420,8 +420,8 @@ class chainset {
* @return pair (bulk_ess, tail_ess)
*/
double mcse_mean(const int index) const {
double ess_bulk = analyze::split_rank_normalized_ess(samples(index)).first;
return sd(index) / std::sqrt(ess_bulk);
double ess_basic = analyze::split_basic_ess(samples(index));
return sd(index) / std::sqrt(ess_basic);
}

/**
Expand Down

0 comments on commit 35bdb07

Please sign in to comment.