1#ifndef STAN_MATH_PRIM_CONSTRAINT_SUM_TO_ZERO_CONSTRAIN_HPP
2#define STAN_MATH_PRIM_CONSTRAINT_SUM_TO_ZERO_CONSTRAIN_HPP
38template <
typename Vec, require_eigen_col_vector_t<Vec>* =
nullptr,
39 require_not_st_var<Vec>* =
nullptr>
41 const auto N = y.size();
51 for (
int i = N; i > 0; --i) {
52 double n =
static_cast<double>(i);
53 auto w = y_ref.coeff(i - 1) *
inv_sqrt(n * (n + 1));
56 z.coeffRef(i - 1) += sum_w;
57 z.coeffRef(i) -= w * n;
73template <
typename Mat, require_eigen_matrix_dynamic_t<Mat>* =
nullptr,
74 require_not_st_var<Mat>* =
nullptr>
76 const auto N = x.rows();
77 const auto M = x.cols();
85 Eigen::Matrix<value_type_t<Mat>, -1, 1>
beta = Eigen::VectorXd::Zero(N);
87 for (
int j = M - 1; j >= 0; --j) {
90 double a_j =
inv_sqrt((j + 1.0) * (j + 2.0));
91 double b_j = (j + 1.0) * a_j;
93 for (
int i = N - 1; i >= 0; --i) {
94 double a_i =
inv_sqrt((i + 1.0) * (i + 2.0));
95 double b_i = (i + 1.0) * a_i;
97 auto b_i_x = b_i * x_ref.coeff(i, j) - ax_previous;
99 Z.coeffRef(i, j) = (b_j * b_i_x) -
beta.coeff(i);
100 beta.coeffRef(i) += a_j * b_i_x;
102 Z.coeffRef(N, j) -= Z.coeff(i, j);
103 Z.coeffRef(i, M) -= Z.coeff(i, j);
105 ax_previous += a_i * x_ref.coeff(i, j);
107 Z.coeffRef(N, M) -= Z.coeff(N, j);
125template <
typename T,
typename Lp, require_not_st_var<T>* =
nullptr>
140template <
typename T, require_std_vector_t<T>* =
nullptr>
158template <
typename T,
typename Lp, require_std_vector_t<T>* =
nullptr,
159 require_convertible_t<return_type_t<T>, Lp>* =
nullptr>
177template <
bool Jacobian,
typename T,
typename Lp>
typename value_type< T >::type value_type_t
Helper function for accessing underlying type.
ref_type_t< T && > to_ref(T &&a)
This evaluates expensive Eigen expressions.
plain_type_t< Vec > sum_to_zero_constrain(const Vec &y)
Return a vector with sum zero corresponding to the specified free vector.
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 ...