Skip to content

Commit

Permalink
Merge branch 'feature/3299-chainset' of https://github.com/stan-dev/stan
Browse files Browse the repository at this point in the history
 into feature/3299-chainset
  • Loading branch information
mitzimorris committed Oct 16, 2024
2 parents 3534106 + ead22cc commit e38ddfc
Showing 1 changed file with 37 additions and 38 deletions.
75 changes: 37 additions & 38 deletions src/stan/analyze/mcmc/mcse.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -11,44 +11,43 @@
namespace stan {
namespace analyze {


/**
* Computes the mean Monte Carlo error estimate for the central 90% interval.
* See https://arxiv.org/abs/1903.08008, section 4.4.
* Follows implementation in the R posterior package.
*
* @param chains matrix of draws across all chains
* @return mcse
*/
inline double mcse_mean(const Eigen::MatrixXd& chains) {
const Eigen::Index num_draws = chains.rows();
if (chains.rows() < 4
|| !is_finite_and_varies(chains))
return std::numeric_limits<double>::quiet_NaN();

double sd = (chains.array() - chains.mean()).square().sum() / (chains.size() - 1);
return std::sqrt(sd / ess(chains));
}

/**
* Computes the standard deviation of the Monte Carlo error estimate
* https://arxiv.org/abs/1903.08008, section 4.4.
* Follows implementation in the R posterior package:
* https://github.com/stan-dev/posterior/blob/98bf52329d68f3307ac4ecaaea659276ee1de8df/R/convergence.R#L478-L496
*
* @param chains matrix of draws across all chains
* @return mcse
*/
inline double mcse_sd(const Eigen::MatrixXd& chains) {
if (chains.rows() < 4
|| !is_finite_and_varies(chains))
return std::numeric_limits<double>::quiet_NaN();

Eigen::MatrixXd diffs = (chains.array() - chains.mean()).matrix();
double Evar = diffs.array().square().mean();
double varvar = (math::mean(diffs.array().pow(4) - Evar * Evar)) / ess(diffs.array().abs().matrix());
return std::sqrt(varvar / Evar / 4);
}
/**
* Computes the mean Monte Carlo error estimate for the central 90% interval.
* See https://arxiv.org/abs/1903.08008, section 4.4.
* Follows implementation in the R posterior package.
*
* @param chains matrix of draws across all chains
* @return mcse
*/
inline double mcse_mean(const Eigen::MatrixXd& chains) {
const Eigen::Index num_draws = chains.rows();
if (chains.rows() < 4 || !is_finite_and_varies(chains))
return std::numeric_limits<double>::quiet_NaN();

double sd
= (chains.array() - chains.mean()).square().sum() / (chains.size() - 1);
return std::sqrt(sd / ess(chains));
}

/**
* Computes the standard deviation of the Monte Carlo error estimate
* https://arxiv.org/abs/1903.08008, section 4.4.
* Follows implementation in the R posterior package:
* https://github.com/stan-dev/posterior/blob/98bf52329d68f3307ac4ecaaea659276ee1de8df/R/convergence.R#L478-L496
*
* @param chains matrix of draws across all chains
* @return mcse
*/
inline double mcse_sd(const Eigen::MatrixXd& chains) {
if (chains.rows() < 4 || !is_finite_and_varies(chains))
return std::numeric_limits<double>::quiet_NaN();

Eigen::MatrixXd diffs = (chains.array() - chains.mean()).matrix();
double Evar = diffs.array().square().mean();
double varvar = (math::mean(diffs.array().pow(4) - Evar * Evar))
/ ess(diffs.array().abs().matrix());
return std::sqrt(varvar / Evar / 4);
}

} // namespace analyze
} // namespace stan
Expand Down

0 comments on commit e38ddfc

Please sign in to comment.