-
-
Notifications
You must be signed in to change notification settings - Fork 369
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
Merge pull request #3312 from stan-dev/feature/3299-improved-ESS-Rhat
Feature/3299 improved ess rhat
- Loading branch information
Showing
18 changed files
with
3,947 additions
and
579 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,43 @@ | ||
| #ifndef STAN_ANALYZE_MCMC_CHECK_CHAINS_HPP | ||
| #define STAN_ANALYZE_MCMC_CHECK_CHAINS_HPP | ||
|
|
||
| #include <stan/math/prim.hpp> | ||
| #include <cmath> | ||
| #include <cstdlib> | ||
| #include <limits> | ||
| #include <utility> | ||
| #include <vector> | ||
|
|
||
| namespace stan { | ||
| namespace analyze { | ||
|
|
||
| /** | ||
| * Checks that values across all matrix columns finite and non-identical. | ||
| * | ||
| * @param chains matrix of draws, one column per chain | ||
| * @return bool true if OK, false otherwise | ||
| */ | ||
| inline bool is_finite_and_varies(const Eigen::MatrixXd chains) { | ||
| size_t num_chains = chains.cols(); | ||
| size_t num_samples = chains.rows(); | ||
| Eigen::VectorXd first_draws = Eigen::VectorXd::Zero(num_chains); | ||
| for (std::size_t i = 0; i < num_chains; ++i) { | ||
| first_draws(i) = chains.col(i)(0); | ||
| for (int j = 0; j < num_samples; ++j) { | ||
| if (!std::isfinite(chains.col(i)(j))) | ||
| return false; | ||
| } | ||
| if (chains.col(i).isApproxToConstant(first_draws(i))) { | ||
| return false; | ||
| } | ||
| } | ||
| if (num_chains > 1 && first_draws.isApproxToConstant(first_draws(0))) { | ||
| return false; | ||
| } | ||
| return true; | ||
| } | ||
|
|
||
| } // namespace analyze | ||
| } // namespace stan | ||
|
|
||
| #endif |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,61 @@ | ||
| #ifndef STAN_ANALYZE_MCMC_RANK_NORMALIZATION_HPP | ||
| #define STAN_ANALYZE_MCMC_RANK_NORMALIZATION_HPP | ||
|
|
||
| #include <stan/math/prim.hpp> | ||
| #include <boost/math/distributions/normal.hpp> | ||
| #include <algorithm> | ||
| #include <cmath> | ||
| #include <vector> | ||
| #include <limits> | ||
|
|
||
| namespace stan { | ||
| namespace analyze { | ||
|
|
||
| /** | ||
| * Computes normalized average ranks for pooled draws. Normal scores computed | ||
| * using inverse normal transformation and a fractional offset. Based on paper | ||
| * https://arxiv.org/abs/1903.08008 | ||
| * | ||
| * @param chains matrix of draws, one column per chain | ||
| * @return normal scores for average ranks of draws | ||
| */ | ||
| inline Eigen::MatrixXd rank_transform(const Eigen::MatrixXd& chains) { | ||
| const Eigen::Index rows = chains.rows(); | ||
| const Eigen::Index cols = chains.cols(); | ||
| const Eigen::Index size = rows * cols; | ||
|
|
||
| std::vector<std::pair<double, int>> value_with_index(size); | ||
| for (Eigen::Index i = 0; i < size; ++i) { | ||
| value_with_index[i] = {chains(i), i}; | ||
| } | ||
|
|
||
| std::sort(value_with_index.begin(), value_with_index.end()); | ||
| Eigen::MatrixXd rank_matrix = Eigen::MatrixXd::Zero(rows, cols); | ||
| // Assigning average ranks | ||
| for (Eigen::Index i = 0; i < size; ++i) { | ||
| // Handle ties by averaging ranks | ||
| Eigen::Index j = i + 1; | ||
| double sum_ranks = j; | ||
| Eigen::Index count = 1; | ||
|
|
||
| while (j < size && value_with_index[j].first == value_with_index[i].first) { | ||
| sum_ranks += j + 1; // Rank starts from 1 | ||
| ++j; | ||
| ++count; | ||
| } | ||
| double avg_rank = sum_ranks / count; | ||
| boost::math::normal_distribution<double> dist; | ||
| for (std::size_t k = i; k < j; ++k) { | ||
| double p = (avg_rank - 0.375) / (size + 0.25); | ||
| const Eigen::Index index = value_with_index[k].second; | ||
| rank_matrix(index) = boost::math::quantile(dist, p); | ||
| } | ||
| i = j - 1; // Skip over tied elements | ||
| } | ||
| return rank_matrix; | ||
| } | ||
|
|
||
| } // namespace analyze | ||
| } // namespace stan | ||
|
|
||
| #endif |
Oops, something went wrong.