1#ifndef STAN_MATH_PRIM_CONSTRAINT_SUM_TO_ZERO_FREE_HPP
2#define STAN_MATH_PRIM_CONSTRAINT_SUM_TO_ZERO_FREE_HPP
38template <
typename Vec, require_eigen_vector_t<Vec>* =
nullptr>
40 const auto& z_ref =
to_ref(z);
44 const auto N = z.size() - 1;
51 y.coeffRef(N - 1) = -z_ref.coeff(N) *
sqrt(N * (N + 1)) / N;
55 for (
int i = N - 2; i >= 0; --i) {
56 double n =
static_cast<double>(i + 1);
57 auto w = y.coeff(i + 1) /
sqrt((n + 1) * (n + 2));
59 y.coeffRef(i) = (sum_w - z_ref.coeff(i + 1)) *
sqrt(n * (n + 1)) / n;
73template <
typename Mat, require_eigen_matrix_dynamic_t<Mat>* =
nullptr>
75 const auto& z_ref =
to_ref(z);
79 const auto N = z_ref.rows() - 1;
80 const auto M = z_ref.cols() - 1;
87 Eigen::Matrix<value_type_t<Mat>, -1, 1>
beta = Eigen::VectorXd::Zero(N);
89 for (
int j = M - 1; j >= 0; --j) {
92 double a_j =
inv_sqrt((j + 1.0) * (j + 2.0));
93 double b_j = (j + 1.0) * a_j;
95 for (
int i = N - 1; i >= 0; --i) {
96 double a_i =
inv_sqrt((i + 1.0) * (i + 2.0));
97 double b_i = (i + 1.0) * a_i;
99 auto alpha_plus_beta = z_ref.coeff(i, j) +
beta.coeff(i);
101 x.coeffRef(i, j) = (alpha_plus_beta + b_j * ax_previous) / (b_j * b_i);
102 beta.coeffRef(i) += a_j * (b_i * x.coeff(i, j) - ax_previous);
103 ax_previous += a_i * x.coeff(i, j);
117template <
typename T, require_std_vector_t<T>* =
nullptr>
typename value_type< T >::type value_type_t
Helper function for accessing underlying type.
plain_type_t< Vec > sum_to_zero_free(const Vec &z)
Return an unconstrained vector.
void check_sum_to_zero(const char *function, const char *name, const T &theta)
Throw an exception if the specified vector does not sum to 0.
fvar< T > sqrt(const fvar< T > &x)
ref_type_t< T && > to_ref(T &&a)
This evaluates expensive Eigen expressions.
fvar< T > beta(const fvar< T > &x1, const fvar< T > &x2)
Return fvar with the beta function applied to the specified arguments and its gradient.
fvar< T > inv_sqrt(const fvar< T > &x)
typename plain_type< T >::type plain_type_t
The lgamma implementation in stan-math is based on either the reentrant safe lgamma_r implementation ...