Loading [MathJax]/extensions/TeX/AMSsymbols.js
Automatic Differentiation
 
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages
Loading...
Searching...
No Matches
stochastic_row_constrain.hpp
Go to the documentation of this file.
1#ifndef STAN_MATH_PRIM_CONSTRAINT_STOCHASTIC_ROW_CONSTRAIN_HPP
2#define STAN_MATH_PRIM_CONSTRAINT_STOCHASTIC_ROW_CONSTRAIN_HPP
3
11#include <cmath>
12
13namespace stan {
14namespace math {
15
25template <typename Mat, require_eigen_matrix_dynamic_t<Mat>* = nullptr,
26 require_not_st_var<Mat>* = nullptr>
28 auto&& y_ref = to_ref(y);
29 const Eigen::Index N = y_ref.rows();
30 int Km1 = y_ref.cols();
31 plain_type_t<Mat> x(N, Km1 + 1);
32 using eigen_arr = Eigen::Array<scalar_type_t<Mat>, -1, 1>;
33 eigen_arr stick_len = eigen_arr::Constant(N, 1.0);
34 for (Eigen::Index k = 0; k < Km1; ++k) {
35 auto z_k = inv_logit(y_ref.array().col(k) - log(Km1 - k));
36 x.array().col(k) = stick_len * z_k;
37 stick_len -= x.array().col(k);
38 }
39 x.array().col(Km1) = stick_len;
40 return x;
41}
42
55template <typename Mat, typename Lp,
57 require_not_st_var<Mat>* = nullptr,
59inline plain_type_t<Mat> stochastic_row_constrain(const Mat& y, Lp& lp) {
60 auto&& y_ref = to_ref(y);
61 const Eigen::Index N = y_ref.rows();
62 Eigen::Index Km1 = y_ref.cols();
63 plain_type_t<Mat> x(N, Km1 + 1);
64 Eigen::Array<scalar_type_t<Mat>, -1, 1> stick_len
65 = Eigen::Array<scalar_type_t<Mat>, -1, 1>::Constant(N, 1.0);
66 for (Eigen::Index k = 0; k < Km1; ++k) {
67 const auto eq_share = -log(Km1 - k); // = logit(1.0/(Km1 + 1 - k));
68 auto adj_y_k = (y_ref.array().col(k) + eq_share).eval();
69 auto z_k = inv_logit(adj_y_k);
70 x.array().col(k) = stick_len * z_k;
71 lp += -sum(log1p_exp(adj_y_k)) - sum(log1p_exp(-adj_y_k))
72 + sum(log(stick_len));
73 stick_len -= x.array().col(k); // equivalently *= (1 - z_k);
74 }
75 x.col(Km1).array() = stick_len;
76 return x;
77}
78
89template <typename T, require_std_vector_t<T>* = nullptr>
90inline auto stochastic_row_constrain(const T& y) {
92 y, [](auto&& v) { return stochastic_row_constrain(v); });
93}
94
108template <typename T, typename Lp, require_std_vector_t<T>* = nullptr,
109 require_convertible_t<return_type_t<T>, Lp>* = nullptr>
110inline auto stochastic_row_constrain(const T& y, Lp& lp) {
112 y, [&lp](auto&& v) { return stochastic_row_constrain(v, lp); });
113}
114
133template <bool Jacobian, typename Mat, typename Lp,
135inline plain_type_t<Mat> stochastic_row_constrain(const Mat& y, Lp& lp) {
136 if constexpr (Jacobian) {
137 return stochastic_row_constrain(y, lp);
138 } else {
139 return stochastic_row_constrain(y);
140 }
141}
142
143} // namespace math
144} // namespace stan
145
146#endif
require_t< std::is_convertible< std::decay_t< T >, std::decay_t< S > > > require_convertible_t
Require types T and S satisfies std::is_convertible.
require_t< is_eigen_matrix_dynamic< std::decay_t< T > > > require_eigen_matrix_dynamic_t
Require type satisfies is_eigen_matrix_dynamic.
require_not_t< is_var< scalar_type_t< std::decay_t< T > > > > require_not_st_var
Require scalar type does not satisfy is_var.
Definition is_var.hpp:117
plain_type_t< Mat > stochastic_row_constrain(const Mat &y)
Return a row stochastic matrix.
T eval(T &&arg)
Inputs which have a plain_type equal to the own time are forwarded unmodified (for Eigen expressions ...
Definition eval.hpp:20
fvar< T > log(const fvar< T > &x)
Definition log.hpp:18
fvar< T > log1p_exp(const fvar< T > &x)
Definition log1p_exp.hpp:14
ref_type_t< T && > to_ref(T &&a)
This evaluates expensive Eigen expressions.
Definition to_ref.hpp:17
auto sum(const std::vector< T > &m)
Return the sum of the entries of the specified standard vector.
Definition sum.hpp:23
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 ...