Merge branch 'develop' into fix/init-err-msgs

RELEASE-NOTES.txt

@@ -5,6 +5,17 @@ Note: these are the release notes for the stan-dev/stan repository.
 Further changes may arise at the interface level (stan-dev/{rstan,
 pystan, cmdstan}) and math library level (stan-dev/math).
 
+v2.35.0 (3 June 2024)
+======================================================================
+
+ - The algorithms no longer catch `std::exception` unconditionally. The Math library uses `std::domain_error` for recoverable errors, and these are the ones which are caught. (#3259)
+ - Allow laplace sampling without evaluating `log_prob` for each draw. (#3261)
+ - Allow laplace sampling to save the Hessian as a diagnostic output. (#3261)
+ - Stan's RNG usages now uses a type definition `stan::rng_t` rather than hard coding a specific Boost RNG. (#3263)
+ - Switched the pRNG used by default in the services and tests to be `boost::mixmax`. Note that this means seeds from previous versions will lead to different numerical results in this version. (#3264)
+ - Add a new ranked R-hat diagnostic from [Vehtari](https://arxiv.org/abs/1903.08008). (#3266)
+ - Fixed an issue where Pathfinder would sometimes return more draws than requested. (#3279)
+
 v2.34.1 (23 January 2024)
 ======================================================================
 

makefile

@@ -2,16 +2,12 @@
 # Stan
 # -----------------
 #
-# To customize your build, set make variables in either:
-#    ~/.config/stan/make.local
-#    make/local
-# Variables in make/local is loaded after ~/.config/stan/make.local
+# To customize your build, set make variables in make/local
 
 
 ## 'help' is the default make target.
 help:
 
--include $(HOME)/.config/stan/make.local  # user-defined variables
 -include make/local                       # user-defined variables
 
 MATH ?= lib/stan_math/

src/stan/analyze/mcmc/check_chains.hpp

@@ -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

src/stan/analyze/mcmc/compute_effective_sample_size.hpp

@@ -12,6 +12,8 @@
 namespace stan {
 namespace analyze {
 /**
+ * \deprecated use split_rank_normalized_ess instead
+ *
  * Computes the effective sample size (ESS) for the specified
  * parameter across all kept samples.  The value returned is the
  * minimum of ESS and the number_total_draws *
@@ -138,6 +140,8 @@ inline double compute_effective_sample_size(std::vector<const double*> draws,
 }
 
 /**
+ * \deprecated use split_rank_normalized_ess instead
+ *
  * Computes the effective sample size (ESS) for the specified
  * parameter across all kept samples.  The value returned is the
  * minimum of ESS and the number_total_draws *
@@ -164,6 +168,8 @@ inline double compute_effective_sample_size(std::vector<const double*> draws,
 }
 
 /**
+ * \deprecated use split_rank_normalized_ess instead
+ *
  * Computes the split effective sample size (ESS) for the specified
  * parameter across all kept samples.  The value returned is the
  * minimum of ESS and the number_total_draws *
@@ -199,6 +205,8 @@ inline double compute_split_effective_sample_size(
 }
 
 /**
+ * \deprecated use split_rank_normalized_ess instead
+ *
  * Computes the split effective sample size (ESS) for the specified
  * parameter across all kept samples.  The value returned is the
  * minimum of ESS and the number_total_draws *

src/stan/analyze/mcmc/compute_potential_scale_reduction.hpp

@@ -18,156 +18,8 @@ namespace stan {
 namespace analyze {
 
 /**
- * Computes normalized average ranks for draws. Transforming them to normal
- * scores using inverse normal transformation and a fractional offset. Based on
- * paper https://arxiv.org/abs/1903.08008
- * @param chains stores chains in columns
- * @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 - 3.0 / 8.0) / (size - 2.0 * 3.0 / 8.0 + 1.0);
-      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;
-}
-
-/**
- * Computes square root of marginal posterior variance of the estimand by the
- * weigted average of within-chain variance W and between-chain variance B.
+ * \deprecated use `rhat` instead
  *
- * @param chains stores chains in columns
- * @return square root of ((N-1)/N)W + B/N
- */
-inline double rhat(const Eigen::MatrixXd& chains) {
-  const Eigen::Index num_chains = chains.cols();
-  const Eigen::Index num_draws = chains.rows();
-
-  Eigen::RowVectorXd within_chain_means = chains.colwise().mean();
-  double across_chain_mean = within_chain_means.mean();
-  double between_variance
-      = num_draws
-        * (within_chain_means.array() - across_chain_mean).square().sum()
-        / (num_chains - 1);
-  double within_variance =
-      // Divide each row by chains and get sum of squares for each chain
-      // (getting a vector back)
-      ((chains.rowwise() - within_chain_means)
-           .array()
-           .square()
-           .colwise()
-           // divide each sum of square by num_draws, sum the sum of squares,
-           // and divide by num chains
-           .sum()
-       / (num_draws - 1.0))
-          .sum()
-      / num_chains;
-
-  return sqrt((between_variance / within_variance + num_draws - 1) / num_draws);
-}
-
-/**
- * Computes the potential scale reduction (Rhat) using rank based diagnostic for
- * the specified parameter across all kept samples. Based on paper
- * https://arxiv.org/abs/1903.08008
- *
- * Current implementation assumes draws are stored in contiguous
- * blocks of memory.  Chains are trimmed from the back to match the
- * length of the shortest chain.
- *
- * @param chain_begins stores pointers to arrays of chains
- * @param chain_sizes stores sizes of chains
- * @return potential scale reduction for the specified parameter
- */
-inline std::pair<double, double> compute_potential_scale_reduction_rank(
-    const std::vector<const double*>& chain_begins,
-    const std::vector<size_t>& chain_sizes) {
-  std::vector<const double*> nonzero_chain_begins;
-  std::vector<std::size_t> nonzero_chain_sizes;
-  nonzero_chain_begins.reserve(chain_begins.size());
-  nonzero_chain_sizes.reserve(chain_sizes.size());
-  for (size_t i = 0; i < chain_sizes.size(); ++i) {
-    if (chain_sizes[i]) {
-      nonzero_chain_begins.push_back(chain_begins[i]);
-      nonzero_chain_sizes.push_back(chain_sizes[i]);
-    }
-  }
-  if (!nonzero_chain_sizes.size()) {
-    return {std::numeric_limits<double>::quiet_NaN(),
-            std::numeric_limits<double>::quiet_NaN()};
-  }
-  std::size_t num_nonzero_chains = nonzero_chain_sizes.size();
-  std::size_t min_num_draws = nonzero_chain_sizes[0];
-  for (std::size_t chain = 1; chain < num_nonzero_chains; ++chain) {
-    min_num_draws = std::min(min_num_draws, nonzero_chain_sizes[chain]);
-  }
-
-  // check if chains are constant; all equal to first draw's value
-  bool are_all_const = false;
-  Eigen::VectorXd init_draw = Eigen::VectorXd::Zero(num_nonzero_chains);
-  Eigen::MatrixXd draws_matrix(min_num_draws, num_nonzero_chains);
-
-  for (std::size_t chain = 0; chain < num_nonzero_chains; chain++) {
-    Eigen::Map<const Eigen::Matrix<double, Eigen::Dynamic, 1>> draws(
-        nonzero_chain_begins[chain], nonzero_chain_sizes[chain]);
-
-    for (std::size_t n = 0; n < min_num_draws; n++) {
-      if (!std::isfinite(draws(n))) {
-        return {std::numeric_limits<double>::quiet_NaN(),
-                std::numeric_limits<double>::quiet_NaN()};
-      }
-      draws_matrix(n, chain) = draws(n);
-    }
-
-    init_draw(chain) = draws(0);
-    are_all_const |= !draws.isApproxToConstant(draws(0));
-  }
-  // If all chains are constant then return NaN
-  if (are_all_const && init_draw.isApproxToConstant(init_draw(0))) {
-    return {std::numeric_limits<double>::quiet_NaN(),
-            std::numeric_limits<double>::quiet_NaN()};
-  }
-
-  double rhat_bulk = rhat(rank_transform(draws_matrix));
-  double rhat_tail = rhat(rank_transform(
-      (draws_matrix.array() - math::quantile(draws_matrix.reshaped(), 0.5))
-          .abs()));
-
-  return std::make_pair(rhat_bulk, rhat_tail);
-}
-
-/**
  * Computes the potential scale reduction (Rhat) for the specified
  * parameter across all kept samples.
  *
@@ -248,34 +100,12 @@ inline double compute_potential_scale_reduction(
                        * num_chains / (num_chains - 1);
   double var_within = chain_var.mean();
 
-  // rewrote [(n-1)*W/n + B/n]/W as (n-1+ B/W)/n
   return sqrt((var_between / var_within + num_draws - 1) / num_draws);
 }
 
 /**
- * Computes the potential scale reduction (Rhat) using rank based diagnostic for
- * the specified parameter across all kept samples. Based on paper
- * https://arxiv.org/abs/1903.08008
- *
- * See more details in Stan reference manual section "Potential
- * Scale Reduction". http://mc-stan.org/users/documentation
- *
- * Current implementation assumes draws are stored in contiguous
- * blocks of memory. Chains are trimmed from the back to match the
- * length of the shortest chain. Argument size will be broadcast to
- * same length as draws.
+ * \deprecated use split_rank_normalized_rhat instead
  *
- * @param chain_begins stores pointers to arrays of chains
- * @param size stores sizes of chains
- * @return potential scale reduction for the specified parameter
- */
-inline std::pair<double, double> compute_potential_scale_reduction_rank(
-    const std::vector<const double*>& chain_begins, size_t size) {
-  std::vector<size_t> sizes(chain_begins.size(), size);
-  return compute_potential_scale_reduction_rank(chain_begins, sizes);
-}
-
-/**
  * Computes the potential scale reduction (Rhat) for the specified
  * parameter across all kept samples.
  *
@@ -299,42 +129,8 @@ inline double compute_potential_scale_reduction(
 }
 
 /**
- * Computes the potential scale reduction (Rhat) using rank based diagnostic for
- * the specified parameter across all kept samples. Based on paper
- * https://arxiv.org/abs/1903.08008
- *
- * When the number of total draws N is odd, the (N+1)/2th draw is ignored.
- *
- * See more details in Stan reference manual section "Potential
- * Scale Reduction". http://mc-stan.org/users/documentation
+ * \deprecated use split_rank_normalized_rhat instead
  *
- * Current implementation assumes draws are stored in contiguous
- * blocks of memory. Chains are trimmed from the back to match the
- * length of the shortest chain.
- *
- * @param chain_begins stores pointers to arrays of chains
- * @param chain_sizes stores sizes of chains
- * @return potential scale reduction for the specified parameter
- */
-inline std::pair<double, double> compute_split_potential_scale_reduction_rank(
-    const std::vector<const double*>& chain_begins,
-    const std::vector<size_t>& chain_sizes) {
-  size_t num_chains = chain_sizes.size();
-  size_t num_draws = chain_sizes[0];
-  for (size_t chain = 1; chain < num_chains; ++chain) {
-    num_draws = std::min(num_draws, chain_sizes[chain]);
-  }
-
-  std::vector<const double*> split_draws
-      = split_chains(chain_begins, chain_sizes);
-
-  size_t half = std::floor(num_draws / 2.0);
-  std::vector<size_t> half_sizes(2 * num_chains, half);
-
-  return compute_potential_scale_reduction_rank(split_draws, half_sizes);
-}
-
-/**
  * Computes the split potential scale reduction (Rhat) for the
  * specified parameter across all kept samples.  When the number of
  * total draws N is odd, the (N+1)/2th draw is ignored.
@@ -367,32 +163,8 @@ inline double compute_split_potential_scale_reduction(
 }
 
 /**
- * Computes the potential scale reduction (Rhat) using rank based diagnostic for
- * the specified parameter across all kept samples. Based on paper
- * https://arxiv.org/abs/1903.08008
- *
- * When the number of total draws N is odd, the (N+1)/2th draw is ignored.
- *
- * See more details in Stan reference manual section "Potential
- * Scale Reduction". http://mc-stan.org/users/documentation
- *
- * Current implementation assumes draws are stored in contiguous
- * blocks of memory.  Chains are trimmed from the back to match the
- * length of the shortest chain.  Argument size will be broadcast to
- * same length as draws.
+ * \deprecated use split_rank_normalized_rhat instead
  *
- * @param chain_begins stores pointers to arrays of chains
- * @param size stores sizes of chains
- * @return potential scale reduction for the specified parameter
- */
-inline std::pair<double, double> compute_split_potential_scale_reduction_rank(
-    const std::vector<const double*>& chain_begins, size_t size) {
-  size_t num_chains = chain_begins.size();
-  std::vector<size_t> sizes(num_chains, size);
-  return compute_split_potential_scale_reduction_rank(chain_begins, sizes);
-}
-
-/**
  * Computes the split potential scale reduction (Rhat) for the
  * specified parameter across all kept samples.  When the number of
  * total draws N is odd, the (N+1)/2th draw is ignored.