Stan Math Library  2.10.0
reverse mode automatic differentiation
simplex_constrain.hpp
Go to the documentation of this file.
1 #ifndef STAN_MATH_PRIM_MAT_FUN_SIMPLEX_CONSTRAIN_HPP
2 #define STAN_MATH_PRIM_MAT_FUN_SIMPLEX_CONSTRAIN_HPP
3 
10 #include <cmath>
11 
12 namespace stan {
13 
14  namespace math {
15 
28  template <typename T>
29  Eigen::Matrix<T, Eigen::Dynamic, 1>
30  simplex_constrain(const Eigen::Matrix<T, Eigen::Dynamic, 1>& y) {
31  // cut & paste simplex_constrain(Eigen::Matrix, T) w/o Jacobian
32  using Eigen::Matrix;
33  using Eigen::Dynamic;
36  using stan::math::logit;
37  using stan::math::log1m;
38  using std::log;
39  typedef typename index_type<Matrix<T, Dynamic, 1> >::type size_type;
40 
41 
42  int Km1 = y.size();
43  Matrix<T, Dynamic, 1> x(Km1 + 1);
44  T stick_len(1.0);
45  for (size_type k = 0; k < Km1; ++k) {
46  T z_k(inv_logit(y(k) - log(Km1 - k)));
47  x(k) = stick_len * z_k;
48  stick_len -= x(k);
49  }
50  x(Km1) = stick_len;
51  return x;
52  }
53 
67  template <typename T>
68  Eigen::Matrix<T, Eigen::Dynamic, 1>
69  simplex_constrain(const Eigen::Matrix<T, Eigen::Dynamic, 1>& y,
70  T& lp) {
71  using Eigen::Dynamic;
72  using Eigen::Matrix;
75  using stan::math::logit;
76  using stan::math::log1m;
78  using std::log;
79 
80  typedef typename index_type<Matrix<T, Dynamic, 1> >::type size_type;
81 
82  int Km1 = y.size(); // K = Km1 + 1
83  Matrix<T, Dynamic, 1> x(Km1 + 1);
84  T stick_len(1.0);
85  for (size_type k = 0; k < Km1; ++k) {
86  double eq_share = -log(Km1 - k); // = logit(1.0/(Km1 + 1 - k));
87  T adj_y_k(y(k) + eq_share);
88  T z_k(inv_logit(adj_y_k));
89  x(k) = stick_len * z_k;
90  lp += log(stick_len);
91  lp -= log1p_exp(-adj_y_k);
92  lp -= log1p_exp(adj_y_k);
93  stick_len -= x(k); // equivalently *= (1 - z_k);
94  }
95  x(Km1) = stick_len; // no Jacobian contrib for last dim
96  return x;
97  }
98 
99  }
100 
101 }
102 
103 #endif
fvar< T > log(const fvar< T > &x)
Definition: log.hpp:15
fvar< T > inv_logit(const fvar< T > &x)
Definition: inv_logit.hpp:15
Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic >::Index size_type
Type for sizes and indexes in an Eigen matrix with double e.
Definition: typedefs.hpp:13
Primary template class for the metaprogram to compute the index type of a container.
Definition: index_type.hpp:19
fvar< T > logit(const fvar< T > &x)
Definition: logit.hpp:17
fvar< T > log1p_exp(const fvar< T > &x)
Definition: log1p_exp.hpp:13
Eigen::Matrix< T, Eigen::Dynamic, 1 > simplex_constrain(const Eigen::Matrix< T, Eigen::Dynamic, 1 > &y)
Return the simplex corresponding to the specified free vector.
fvar< T > log1m(const fvar< T > &x)
Definition: log1m.hpp:16

     [ Stan Home Page ] © 2011–2016, Stan Development Team.