Automatic Differentiation
 
Loading...
Searching...
No Matches
simplex_constrain.hpp
Go to the documentation of this file.
1#ifndef STAN_MATH_PRIM_CONSTRAINT_SIMPLEX_CONSTRAIN_HPP
2#define STAN_MATH_PRIM_CONSTRAINT_SIMPLEX_CONSTRAIN_HPP
3
10#include <cmath>
11
12namespace stan {
13namespace math {
14
27template <typename Vec, require_eigen_vector_t<Vec>* = nullptr,
28 require_not_st_var<Vec>* = nullptr>
30 // cut & paste simplex_constrain(Eigen::Matrix, T) w/o Jacobian
31 using std::log;
32 using T = value_type_t<Vec>;
33
34 int Km1 = y.size();
35 plain_type_t<Vec> x(Km1 + 1);
36 T stick_len(1.0);
37 for (Eigen::Index k = 0; k < Km1; ++k) {
38 T z_k = inv_logit(y.coeff(k) - log(Km1 - k));
39 x.coeffRef(k) = stick_len * z_k;
40 stick_len -= x.coeff(k);
41 }
42 x.coeffRef(Km1) = stick_len;
43 return x;
44}
45
59template <typename Vec, require_eigen_vector_t<Vec>* = nullptr,
60 require_not_st_var<Vec>* = nullptr>
63 using Eigen::Dynamic;
64 using Eigen::Matrix;
65 using std::log;
66 using T = value_type_t<Vec>;
67
68 int Km1 = y.size(); // K = Km1 + 1
69 plain_type_t<Vec> x(Km1 + 1);
70 T stick_len(1.0);
71 for (Eigen::Index k = 0; k < Km1; ++k) {
72 double eq_share = -log(Km1 - k); // = logit(1.0/(Km1 + 1 - k));
73 T adj_y_k = y.coeff(k) + eq_share;
74 T z_k = inv_logit(adj_y_k);
75 x.coeffRef(k) = stick_len * z_k;
76 lp += log(stick_len);
77 lp -= log1p_exp(-adj_y_k);
78 lp -= log1p_exp(adj_y_k);
79 stick_len -= x.coeff(k); // equivalently *= (1 - z_k);
80 }
81 x.coeffRef(Km1) = stick_len; // no Jacobian contrib for last dim
82 return x;
83}
84
101template <bool Jacobian, typename Vec, require_not_std_vector_t<Vec>* = nullptr>
103 return_type_t<Vec>& lp) {
104 if (Jacobian) {
105 return simplex_constrain(y, lp);
106 } else {
107 return simplex_constrain(y);
108 }
109}
110
127template <bool Jacobian, typename T, require_std_vector_t<T>* = nullptr>
128inline auto simplex_constrain(const T& y, return_type_t<T>& lp) {
130 y, [&lp](auto&& v) { return simplex_constrain<Jacobian>(v, lp); });
131}
132
133} // namespace math
134} // namespace stan
135
136#endif
typename value_type< T >::type value_type_t
Helper function for accessing underlying type.
typename return_type< Ts... >::type return_type_t
Convenience type for the return type of the specified template parameters.
fvar< T > log(const fvar< T > &x)
Definition log.hpp:15
fvar< T > log1p_exp(const fvar< T > &x)
Definition log1p_exp.hpp:13
plain_type_t< Vec > simplex_constrain(const Vec &y)
Return the simplex corresponding to the specified free vector.
fvar< T > inv_logit(const fvar< T > &x)
Returns the inverse logit function applied to the argument.
Definition inv_logit.hpp:20
typename plain_type< T >::type plain_type_t
The lgamma implementation in stan-math is based on either the reentrant safe lgamma_r implementation ...