simplex_free.hpp

Go to the documentation of this file.

#ifndef STAN_MATH_PRIM_CONSTRAINT_SIMPLEX_FREE_HPP
#define STAN_MATH_PRIM_CONSTRAINT_SIMPLEX_FREE_HPP
 
#include <stan/math/prim/meta.hpp>
#include <stan/math/prim/err.hpp>
#include <stan/math/prim/fun/Eigen.hpp>
#include <stan/math/prim/fun/log.hpp>
#include <stan/math/prim/fun/sqrt.hpp>
#include <stan/math/prim/fun/to_ref.hpp>
#include <cmath>
 
namespace stan {
namespace math {
 
template <typename Vec, require_eigen_vector_t<Vec>* = nullptr>
inline plain_type_t<Vec> simplex_free(const Vec& x) {
  using T = value_type_t<Vec>;
 
  const auto& x_ref = to_ref(x);
  check_simplex("stan::math::simplex_free", "Simplex variable", x_ref);
  Eigen::Index Km1 = x_ref.size() - 1;
  plain_type_t<Vec> y(Km1);
 
  T cumsum = 0.0;
  for (int i = 0; i < Km1; ++i) {
    cumsum += log(x_ref.coeff(i));
    double n = static_cast<double>(i + 1);
    y.coeffRef(i) = (cumsum - n * log(x_ref.coeff(i + 1))) / sqrt(n * (n + 1));
  }
 
  return y;
}
 
template <typename T, require_std_vector_t<T>* = nullptr>
inline auto simplex_free(const T& x) {
  return apply_vector_unary<T>::apply(x,
                                      [](auto&& v) { return simplex_free(v); });
}
 
}  // namespace math
}  // namespace stan
 
#endif