Stan  1.3
probability, sampling & optimization
 All Classes Namespaces Files Functions Variables Typedefs Enumerator Friends Macros Pages
transform.hpp
Go to the documentation of this file.
1 #ifndef __STAN__PROB__TRANSFORM_HPP__
2 #define __STAN__PROB__TRANSFORM_HPP__
3 
4 #include <cmath>
5 #include <cstddef>
6 #include <stdexcept>
7 #include <sstream>
8 #include <vector>
9 #include <boost/multi_array.hpp>
10 #include <boost/throw_exception.hpp>
11 #include <boost/math/tools/promotion.hpp>
12 #include <stan/agrad/matrix.hpp>
13 #include <stan/math.hpp>
14 #include <stan/math/matrix.hpp>
15 #include <stan/math/matrix/validate_less.hpp>
16 #include <stan/math/error_handling.hpp>
17 #include <stan/math/matrix_error_handling.hpp>
18 
19 #include <stan/math/matrix/multiply_lower_tri_self_transpose.hpp>
20 
21 namespace stan {
22 
23  namespace prob {
24 
25 
26  const double CONSTRAINT_TOLERANCE = 1E-8;
27 
28 
41  template<typename T>
42  bool
43  factor_cov_matrix(Eigen::Array<T,Eigen::Dynamic,1>& CPCs,
44  Eigen::Array<T,Eigen::Dynamic,1>& sds,
45  const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& Sigma) {
46 
47  size_t K = sds.rows();
48 
49  sds = Sigma.diagonal().array();
50  if( (sds <= 0.0).any() ) return false;
51  sds = sds.sqrt();
52 
53  Eigen::DiagonalMatrix<T,Eigen::Dynamic> D(K);
54  D.diagonal() = sds.inverse();
55  sds = sds.log(); // now unbounded
56 
57  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> R = D * Sigma * D;
58  // to hopefully prevent pivoting due to floating point error
59  R.diagonal().setOnes();
60  Eigen::LDLT<Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> > ldlt;
61  ldlt = R.ldlt();
62  if (!ldlt.isPositive())
63  return false;
64  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> U = ldlt.matrixU();
65 
66  size_t position = 0;
67  size_t pull = K - 1;
68 
69  Eigen::Array<T,1,Eigen::Dynamic> temp = U.row(0).tail(pull);
70 
71  CPCs.head(pull) = temp;
72 
73  Eigen::Array<T,Eigen::Dynamic,1> acc(K);
74  acc(0) = -0.0;
75  acc.tail(pull) = 1.0 - temp.square();
76  for(size_t i = 1; i < (K - 1); i++) {
77  position += pull;
78  pull--;
79  temp = U.row(i).tail(pull);
80  temp /= sqrt(acc.tail(pull) / acc(i));
81  CPCs.segment(position, pull) = temp;
82  acc.tail(pull) *= 1.0 - temp.square();
83  }
84  CPCs = 0.5 * ( (1.0 + CPCs) / (1.0 - CPCs) ).log(); // now unbounded
85  return true;
86  }
87 
88  // MATRIX TRANSFORMS +/- JACOBIANS
89 
111  template <typename T>
112  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>
113  read_corr_L(const Eigen::Array<T,Eigen::Dynamic,1>& CPCs, // on (-1,1)
114  const size_t K) {
115  Eigen::Array<T,Eigen::Dynamic,1> temp;
116  Eigen::Array<T,Eigen::Dynamic,1> acc(K-1);
117  acc.setOnes();
118  // Cholesky factor of correlation matrix
119  Eigen::Array<T,Eigen::Dynamic,Eigen::Dynamic> L(K,K);
120  L.setZero();
121 
122  size_t position = 0;
123  size_t pull = K - 1;
124 
125  L(0,0) = 1.0;
126  L.col(0).tail(pull) = temp = CPCs.head(pull);
127  acc.tail(pull) = 1.0 - temp.square();
128  for(size_t i = 1; i < (K - 1); i++) {
129  position += pull;
130  pull--;
131  temp = CPCs.segment(position, pull);
132  L(i,i) = sqrt(acc(i-1));
133  L.col(i).tail(pull) = temp * acc.tail(pull).sqrt();
134  acc.tail(pull) *= 1.0 - temp.square();
135  }
136  L(K-1,K-1) = sqrt(acc(K-2));
137  return L.matrix();
138  }
139 
153  template <typename T>
154  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>
155  read_corr_matrix(const Eigen::Array<T,Eigen::Dynamic,1>& CPCs,
156  const size_t K) {
157  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> L
158  = read_corr_L(CPCs, K);
159  using stan::math::multiply_lower_tri_self_transpose;
160  return multiply_lower_tri_self_transpose(L);
161  }
162 
189  template <typename T>
190  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>
191  read_corr_L(const Eigen::Array<T,Eigen::Dynamic,1>& CPCs,
192  const size_t K,
193  T& log_prob) {
194 
195  size_t k = 0;
196  size_t i = 0;
197  T log_1cpc2;
198  double lead = K - 2.0;
199  // no need to abs() because this Jacobian determinant
200  // is strictly positive (and triangular)
201  // skip last row (odd indexing) because it adds nothing by design
202  typedef typename Eigen::Matrix<T,Eigen::Dynamic,1>::size_type size_type;
203  for (size_type j = 0;
204  j < (CPCs.rows() - 1);
205  ++j) {
206  using stan::math::log1m;
207  using stan::math::square;
208  log_1cpc2 = log1m(square(CPCs[j]));
209  // derivative of correlation wrt CPC
210  log_prob += lead / 2.0 * log_1cpc2;
211  i++;
212  if (i > K) {
213  k++;
214  i = k + 1;
215  lead = K - k - 1.0;
216  }
217  }
218  return read_corr_L(CPCs, K);
219  }
220 
239  template <typename T>
240  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>
241  read_corr_matrix(const Eigen::Array<T,Eigen::Dynamic,1>& CPCs,
242  const size_t K,
243  T& log_prob) {
244 
245  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> L
246  = read_corr_L(CPCs, K, log_prob);
247  using stan::math::multiply_lower_tri_self_transpose;
248  return multiply_lower_tri_self_transpose(L);
249  }
250 
261  template <typename T>
262  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>
263  read_cov_L(const Eigen::Array<T,Eigen::Dynamic,1>& CPCs,
264  const Eigen::Array<T,Eigen::Dynamic,1>& sds,
265  T& log_prob) {
266  size_t K = sds.rows();
267  // adjust due to transformation from correlations to covariances
268  log_prob += (sds.log().sum() + stan::math::LOG_2) * K;
269  return sds.matrix().asDiagonal() * read_corr_L(CPCs, K, log_prob);
270  }
271 
281  template <typename T>
282  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>
283  read_cov_matrix(const Eigen::Array<T,Eigen::Dynamic,1>& CPCs,
284  const Eigen::Array<T,Eigen::Dynamic,1>& sds,
285  T& log_prob) {
286 
287  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> L
288  = read_cov_L(CPCs, sds, log_prob);
289  using stan::math::multiply_lower_tri_self_transpose;
290  return multiply_lower_tri_self_transpose(L);
291  }
292 
300  template<typename T>
301  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>
302  read_cov_matrix(const Eigen::Array<T,Eigen::Dynamic,1>& CPCs,
303  const Eigen::Array<T,Eigen::Dynamic,1>& sds) {
304 
305  size_t K = sds.rows();
306  Eigen::DiagonalMatrix<T,Eigen::Dynamic> D(K);
307  D.diagonal() = sds;
308  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> L
309  = D * read_corr_L(CPCs, K);
310  using stan::math::multiply_lower_tri_self_transpose;
311  return multiply_lower_tri_self_transpose(L);
312  }
313 
314 
324  template<typename T>
325  const Eigen::Array<T,Eigen::Dynamic,1>
326  make_nu(const T eta, const size_t K) {
327 
328  Eigen::Array<T,Eigen::Dynamic,1> nu(K * (K - 1) / 2);
329 
330  T alpha = eta + (K - 2.0) / 2.0; // from Lewandowski et. al.
331 
332  // Best (1978) implies nu = 2 * alpha for the dof in a t
333  // distribution that generates a beta variate on (-1,1)
334  T alpha2 = 2.0 * alpha;
335 
336  typedef typename Eigen::Matrix<T,Eigen::Dynamic,1>::size_type size_type;
337  for (size_type j = 0; j < (K - 1); j++) {
338  nu(j) = alpha2;
339  }
340  size_t counter = K - 1;
341  for (size_type i = 1; i < (K - 1); i++) {
342  alpha -= 0.5;
343  alpha2 = 2.0 * alpha;
344  for (size_type j = i + 1; j < K; j++) {
345  nu(counter) = alpha2;
346  counter++;
347  }
348  }
349  return nu;
350  }
351 
352 
353 
354 
355  // IDENTITY
356 
369  template <typename T>
370  inline
371  T identity_constrain(T x) {
372  return x;
373  }
374 
388  template <typename T>
389  inline
390  T identity_constrain(const T x, T& /*lp*/) {
391  return x;
392  }
393 
405  template <typename T>
406  inline
407  T identity_free(const T y) {
408  return y;
409  }
410 
411 
412  // POSITIVE
413 
424  template <typename T>
425  inline
426  T positive_constrain(const T x) {
427  return exp(x);
428  }
429 
446  template <typename T>
447  inline
448  T positive_constrain(const T x, T& lp) {
449  lp += x;
450  return exp(x);
451  }
452 
469  template <typename T>
470  inline
471  T positive_free(const T y) {
472  stan::math::check_positive("stan::prob::positive_free(%1%)", y, "Positive variable");
473  return log(y);
474  }
475 
476  // LOWER BOUND
477 
497  template <typename T, typename TL>
498  inline
499  T lb_constrain(const T x, const TL lb) {
500  if (lb == -std::numeric_limits<double>::infinity())
501  return identity_constrain(x);
502  return exp(x) + lb;
503  }
504 
521  template <typename T, typename TL>
522  inline
523  typename boost::math::tools::promote_args<T,TL>::type
524  lb_constrain(const T x, const TL lb, T& lp) {
525  if (lb == -std::numeric_limits<double>::infinity())
526  return identity_constrain(x,lp);
527  lp += x;
528  return exp(x) + lb;
529  }
530 
546  template <typename T, typename TL>
547  inline
548  typename boost::math::tools::promote_args<T,TL>::type
549  lb_free(const T y, const TL lb) {
550  if (lb == -std::numeric_limits<double>::infinity())
551  return identity_free(y);
552  stan::math::check_greater_or_equal("stan::prob::lb_free(%1%)",
553  y, lb, "Lower bounded variable");
554  return log(y - lb);
555  }
556 
557 
558  // UPPER BOUND
559 
579  template <typename T, typename TU>
580  inline
581  typename boost::math::tools::promote_args<T,TU>::type
582  ub_constrain(const T x, const TU ub) {
583  if (ub == std::numeric_limits<double>::infinity())
584  return identity_constrain(x);
585  return ub - exp(x);
586  }
587 
611  template <typename T, typename TU>
612  inline
613  typename boost::math::tools::promote_args<T,TU>::type
614  ub_constrain(const T x, const TU ub, T& lp) {
615  if (ub == std::numeric_limits<double>::infinity())
616  return identity_constrain(x,lp);
617  lp += x;
618  return ub - exp(x);
619  }
620 
643  template <typename T, typename TU>
644  inline
645  typename boost::math::tools::promote_args<T,TU>::type
646  ub_free(const T y, const TU ub) {
647  if (ub == std::numeric_limits<double>::infinity())
648  return identity_free(y);
649  stan::math::check_less_or_equal("stan::prob::ub_free(%1%)",
650  y, ub, "Upper bounded variable");
651  return log(ub - y);
652  }
653 
654 
655  // LOWER & UPPER BOUNDS
656 
684  template <typename T, typename TL, typename TU>
685  inline
686  typename boost::math::tools::promote_args<T,TL,TU>::type
687  lub_constrain(const T x, TL lb, TU ub) {
688  stan::math::validate_less(lb,ub,"lb","ub","lub_constrain/3");
689 
690  if (lb == -std::numeric_limits<double>::infinity())
691  return ub_constrain(x,ub);
692  if (ub == std::numeric_limits<double>::infinity())
693  return lb_constrain(x,lb);
694 
695  T inv_logit_x;
696  if (x > 0) {
697  T exp_minus_x = exp(-x);
698  inv_logit_x = 1.0 / (1.0 + exp_minus_x);
699  // Prevent x from reaching one unless it really really should.
700  if ((x < std::numeric_limits<double>::infinity())
701  && (inv_logit_x == 1))
702  inv_logit_x = 1 - 1e-15;
703  } else {
704  T exp_x = exp(x);
705  inv_logit_x = 1.0 - 1.0 / (1.0 + exp_x);
706  // Prevent x from reaching zero unless it really really should.
707  if ((x > -std::numeric_limits<double>::infinity())
708  && (inv_logit_x== 0))
709  inv_logit_x = 1e-15;
710  }
711  return lb + (ub - lb) * inv_logit_x;
712  }
713 
755  template <typename T, typename TL, typename TU>
756  typename boost::math::tools::promote_args<T,TL,TU>::type
757  lub_constrain(const T x, const TL lb, const TU ub, T& lp) {
758  if (!(lb < ub)) {
759  std::stringstream s;
760  s << "domain error in lub_constrain; lower bound = " << lb
761  << " must be strictly less than upper bound = " << ub;
762  throw std::domain_error(s.str());
763  }
764  if (lb == -std::numeric_limits<double>::infinity())
765  return ub_constrain(x,ub,lp);
766  if (ub == std::numeric_limits<double>::infinity())
767  return lb_constrain(x,lb,lp);
768  T inv_logit_x;
769  if (x > 0) {
770  T exp_minus_x = exp(-x);
771  inv_logit_x = 1.0 / (1.0 + exp_minus_x);
772  lp += log(ub - lb) - x - 2 * log1p(exp_minus_x);
773  // Prevent x from reaching one unless it really really should.
774  if ((x < std::numeric_limits<double>::infinity())
775  && (inv_logit_x == 1))
776  inv_logit_x = 1 - 1e-15;
777  } else {
778  T exp_x = exp(x);
779  inv_logit_x = 1.0 - 1.0 / (1.0 + exp_x);
780  lp += log(ub - lb) + x - 2 * log1p(exp_x);
781  // Prevent x from reaching zero unless it really really should.
782  if ((x > -std::numeric_limits<double>::infinity())
783  && (inv_logit_x== 0))
784  inv_logit_x = 1e-15;
785  }
786  return lb + (ub - lb) * inv_logit_x;
787  }
788 
819  template <typename T, typename TL, typename TU>
820  inline
821  typename boost::math::tools::promote_args<T,TL,TU>::type
822  lub_free(const T y, TL lb, TU ub) {
823  using stan::math::logit;
824  stan::math::check_bounded("stan::prob::lub_free(%1%)",
825  y, lb, ub, "Bounded variable");
826  if (lb == -std::numeric_limits<double>::infinity())
827  return ub_free(y,ub);
828  if (ub == std::numeric_limits<double>::infinity())
829  return lb_free(y,lb);
830  return logit((y - lb) / (ub - lb));
831  }
832 
833 
834  // PROBABILITY
835 
849  template <typename T>
850  inline
851  T prob_constrain(const T x) {
852  using stan::math::inv_logit;
853  return inv_logit(x);
854  }
855 
877  template <typename T>
878  inline
879  T prob_constrain(const T x, T& lp) {
880  using stan::math::inv_logit;
881  using stan::math::log1m;
882  T inv_logit_x = inv_logit(x);
883  lp += log(inv_logit_x) + log1m(inv_logit_x);
884  return inv_logit_x;
885  }
886 
901  template <typename T>
902  inline
903  T prob_free(const T y) {
904  using stan::math::logit;
905  stan::math::check_bounded("stan::prob::prob_free(%1%)",
906  y, 0, 1, "Probability variable");
907  return logit(y);
908  }
909 
910 
911  // CORRELATION
912 
925  template <typename T>
926  inline
927  T corr_constrain(const T x) {
928  return tanh(x);
929  }
930 
943  template <typename T>
944  inline
945  T corr_constrain(const T x, T& lp) {
946  using stan::math::log1m;
947  T tanh_x = tanh(x);
948  lp += log1m(tanh_x * tanh_x);
949  return tanh_x;
950  }
951 
968  template <typename T>
969  inline
970  T corr_free(const T y) {
971  stan::math::check_bounded("stan::prob::lub_free(%1%)",
972  y, -1, 1, "Correlation variable");
973  return atanh(y);
974  }
975 
976 
977  // Unit vector
978 
987  template <typename T>
988  Eigen::Matrix<T,Eigen::Dynamic,1>
989  unit_vector_constrain(const Eigen::Matrix<T,Eigen::Dynamic,1>& y) {
990  typedef typename Eigen::Matrix<T,Eigen::Dynamic,1>::size_type size_type;
991  int Km1 = y.size();
992  Eigen::Matrix<T,Eigen::Dynamic,1> x(Km1 + 1);
993  x(0) = 1.0;
994  const T half_pi = T(M_PI/2.0);
995  for (size_type k = 1; k <= Km1; ++k) {
996  T yk_1 = y(k-1) + half_pi;
997  T sin_yk_1 = sin(yk_1);
998  x(k) = x(k-1)*sin_yk_1;
999  x(k-1) *= cos(yk_1);
1000  }
1001  return x;
1002  }
1003 
1013  template <typename T>
1014  Eigen::Matrix<T,Eigen::Dynamic,1>
1015  unit_vector_constrain(const Eigen::Matrix<T,Eigen::Dynamic,1>& y, T &lp) {
1016  typedef typename Eigen::Matrix<T,Eigen::Dynamic,1>::size_type size_type;
1017  int Km1 = y.size();
1018  Eigen::Matrix<T,Eigen::Dynamic,1> x(Km1 + 1);
1019  x(0) = 1.0;
1020  const T half_pi = T(M_PI/2.0);
1021  for (size_type k = 1; k <= Km1; ++k) {
1022  T yk_1 = y(k-1) + half_pi;
1023  T sin_yk_1 = sin(yk_1);
1024  x(k) = x(k-1)*sin_yk_1;
1025  x(k-1) *= cos(yk_1);
1026  if (k < Km1)
1027  lp += (Km1 - k)*log(fabs(sin_yk_1));
1028  }
1029  return x;
1030  }
1031 
1032  template <typename T>
1033  Eigen::Matrix<T,Eigen::Dynamic,1>
1034  unit_vector_free(const Eigen::Matrix<T,Eigen::Dynamic,1>& x) {
1035  typedef typename Eigen::Matrix<T,Eigen::Dynamic,1>::size_type size_type;
1036  stan::math::check_unit_vector("stan::prob::unit_vector_free(%1%)", x, "Unit vector variable");
1037  int Km1 = x.size() - 1;
1038  Eigen::Matrix<T,Eigen::Dynamic,1> y(Km1);
1039  T sumSq = x(Km1)*x(Km1);
1040  const T half_pi = T(M_PI/2.0);
1041  for (size_type k = Km1; --k >= 0; ) {
1042  y(k) = atan2(sqrt(sumSq),x(k)) - half_pi;
1043  sumSq += x(k)*x(k);
1044  }
1045  return y;
1046  }
1047 
1048  // SIMPLEX
1049 
1050 
1063  template <typename T>
1064  Eigen::Matrix<T,Eigen::Dynamic,1>
1065  simplex_constrain(const Eigen::Matrix<T,Eigen::Dynamic,1>& y) {
1066  // cut & paste simplex_constrain(Eigen::Matrix,T) w/o Jacobian
1067  typedef typename Eigen::Matrix<T,Eigen::Dynamic,1>::size_type size_type;
1068  using stan::math::logit;
1069  using stan::math::inv_logit;
1070  using stan::math::log1m;
1071  int Km1 = y.size();
1072  Eigen::Matrix<T,Eigen::Dynamic,1> x(Km1 + 1);
1073  T stick_len(1.0);
1074  for (size_type k = 0; k < Km1; ++k) {
1075  T z_k(inv_logit(y(k) - log(Km1 - k)));
1076  x(k) = stick_len * z_k;
1077  stick_len -= x(k);
1078  }
1079  x(Km1) = stick_len;
1080  return x;
1081  }
1082 
1096  template <typename T>
1097  Eigen::Matrix<T,Eigen::Dynamic,1>
1098  simplex_constrain(const Eigen::Matrix<T,Eigen::Dynamic,1>& y,
1099  T& lp) {
1100  using stan::math::logit;
1101  using stan::math::inv_logit;
1102  using stan::math::log1p_exp;
1103  using stan::math::log1m;
1104  typedef typename Eigen::Matrix<T,Eigen::Dynamic,1>::size_type size_type;
1105  int Km1 = y.size(); // K = Km1 + 1
1106  Eigen::Matrix<T,Eigen::Dynamic,1> x(Km1 + 1);
1107  T stick_len(1.0);
1108  for (size_type k = 0; k < Km1; ++k) {
1109  double eq_share = -log(Km1 - k); // = logit(1.0/(Km1 + 1 - k));
1110  T adj_y_k(y(k) + eq_share);
1111  T z_k(inv_logit(adj_y_k));
1112  x(k) = stick_len * z_k;
1113  lp += log(stick_len);
1114  lp -= log1p_exp(-adj_y_k);
1115  lp -= log1p_exp(adj_y_k);
1116  stick_len -= x(k); // equivalently *= (1 - z_k);
1117  }
1118  x(Km1) = stick_len; // no Jacobian contrib for last dim
1119  return x;
1120  }
1121 
1136  template <typename T>
1137  Eigen::Matrix<T,Eigen::Dynamic,1>
1138  simplex_free(const Eigen::Matrix<T,Eigen::Dynamic,1>& x) {
1139  using stan::math::logit;
1140  typedef typename Eigen::Matrix<T,Eigen::Dynamic,1>::size_type size_type;
1141  stan::math::check_simplex("stan::prob::simplex_free(%1%)", x, "Simplex variable");
1142  int Km1 = x.size() - 1;
1143  Eigen::Matrix<T,Eigen::Dynamic,1> y(Km1);
1144  T stick_len(x(Km1));
1145  for (size_type k = Km1; --k >= 0; ) {
1146  stick_len += x(k);
1147  T z_k(x(k) / stick_len);
1148  y(k) = logit(z_k) + log(Km1 - k);
1149  // log(Km-k) = logit(1.0 / (Km1 + 1 - k));
1150  }
1151  return y;
1152  }
1153 
1154 
1155  // ORDERED
1156 
1166  template <typename T>
1167  Eigen::Matrix<T,Eigen::Dynamic,1>
1168  ordered_constrain(const Eigen::Matrix<T,Eigen::Dynamic,1>& x) {
1169  typedef typename Eigen::Matrix<T,Eigen::Dynamic,1>::size_type size_type;
1170  size_type k = x.size();
1171  Eigen::Matrix<T,Eigen::Dynamic,1> y(k);
1172  if (k == 0)
1173  return y;
1174  y[0] = x[0];
1175  for (size_type i = 1;
1176  i < k;
1177  ++i)
1178  y[i] = y[i-1] + exp(x[i]);
1179  return y;
1180  }
1181 
1194  template <typename T>
1195  inline
1196  Eigen::Matrix<T,Eigen::Dynamic,1>
1197  ordered_constrain(const Eigen::Matrix<T,Eigen::Dynamic,1>& x, T& lp) {
1198  typedef typename Eigen::Matrix<T,Eigen::Dynamic,1>::size_type size_type;
1199  for (size_type i = 1; i < x.size(); ++i)
1200  lp += x(i);
1201  return ordered_constrain(x);
1202  }
1203 
1204 
1205 
1219  template <typename T>
1220  Eigen::Matrix<T,Eigen::Dynamic,1>
1221  ordered_free(const Eigen::Matrix<T,Eigen::Dynamic,1>& y) {
1222  stan::math::check_ordered("stan::prob::ordered_free(%1%)",
1223  y, "Ordered variable");
1224  typedef typename Eigen::Matrix<T,Eigen::Dynamic,1>::size_type size_type;
1225  size_type k = y.size();
1226  Eigen::Matrix<T,Eigen::Dynamic,1> x(k);
1227  if (k == 0)
1228  return x;
1229  x[0] = y[0];
1230  for (size_type i = 1; i < k; ++i)
1231  x[i] = log(y[i] - y[i-1]);
1232  return x;
1233  }
1234 
1235 
1236  // POSITIVE ORDERED
1237 
1247  template <typename T>
1248  Eigen::Matrix<T,Eigen::Dynamic,1>
1249  positive_ordered_constrain(const Eigen::Matrix<T,Eigen::Dynamic,1>& x) {
1250  typedef typename Eigen::Matrix<T,Eigen::Dynamic,1>::size_type size_type;
1251  size_type k = x.size();
1252  Eigen::Matrix<T,Eigen::Dynamic,1> y(k);
1253  if (k == 0)
1254  return y;
1255  y[0] = exp(x[0]);
1256  for (size_type i = 1;
1257  i < k;
1258  ++i)
1259  y[i] = y[i-1] + exp(x[i]);
1260  return y;
1261  }
1262 
1275  template <typename T>
1276  inline
1277  Eigen::Matrix<T,Eigen::Dynamic,1>
1278  positive_ordered_constrain(const Eigen::Matrix<T,Eigen::Dynamic,1>& x, T& lp) {
1279  typedef typename Eigen::Matrix<T,Eigen::Dynamic,1>::size_type size_type;
1280  for (size_type i = 0; i < x.size(); ++i)
1281  lp += x(i);
1282  return positive_ordered_constrain(x);
1283  }
1284 
1285 
1286 
1300  template <typename T>
1301  Eigen::Matrix<T,Eigen::Dynamic,1>
1302  positive_ordered_free(const Eigen::Matrix<T,Eigen::Dynamic,1>& y) {
1303  stan::math::check_positive_ordered("stan::prob::positive_ordered_free(%1%)",
1304  y, "Positive ordered variable");
1305  typedef typename Eigen::Matrix<T,Eigen::Dynamic,1>::size_type size_type;
1306  size_type k = y.size();
1307  Eigen::Matrix<T,Eigen::Dynamic,1> x(k);
1308  if (k == 0)
1309  return x;
1310  x[0] = log(y[0]);
1311  for (size_type i = 1; i < k; ++i)
1312  x[i] = log(y[i] - y[i-1]);
1313  return x;
1314  }
1315 
1316 
1317  // CORRELATION MATRIX
1342  template <typename T>
1343  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>
1344  corr_matrix_constrain(const Eigen::Matrix<T,Eigen::Dynamic,1>& x,
1345  typename Eigen::Matrix<T,Eigen::Dynamic,1>::size_type k) {
1346  typedef typename Eigen::Matrix<T,Eigen::Dynamic,1>::size_type size_type;
1347  size_type k_choose_2 = (k * (k - 1)) / 2;
1348  if (k_choose_2 != x.size())
1349  throw std::invalid_argument ("x is not a valid correlation matrix");
1350  Eigen::Array<T,Eigen::Dynamic,1> cpcs(k_choose_2);
1351  for (size_type i = 0; i < k_choose_2; ++i)
1352  cpcs[i] = corr_constrain(x[i]);
1353  return read_corr_matrix(cpcs,k);
1354  }
1355 
1375  template <typename T>
1376  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>
1377  corr_matrix_constrain(const Eigen::Matrix<T,Eigen::Dynamic,1>& x,
1378  typename Eigen::Matrix<T,Eigen::Dynamic,1>::size_type k,
1379  T& lp) {
1380  typedef typename Eigen::Matrix<T,Eigen::Dynamic,1>::size_type size_type;
1381  size_type k_choose_2 = (k * (k - 1)) / 2;
1382  if (k_choose_2 != x.size())
1383  throw std::invalid_argument ("x is not a valid correlation matrix");
1384  Eigen::Array<T,Eigen::Dynamic,1> cpcs(k_choose_2);
1385  for (size_type i = 0; i < k_choose_2; ++i)
1386  cpcs[i] = corr_constrain(x[i],lp);
1387  return read_corr_matrix(cpcs,k,lp);
1388  }
1389 
1410  template <typename T>
1411  Eigen::Matrix<T,Eigen::Dynamic,1>
1412  corr_matrix_free(const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& y) {
1413  typedef typename
1414  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>::size_type size_type;
1415  size_type k = y.rows();
1416  if (y.cols() != k)
1417  throw std::domain_error("y is not a square matrix");
1418  if (k == 0)
1419  throw std::domain_error("y has no elements");
1420  size_type k_choose_2 = (k * (k-1)) / 2;
1421  Eigen::Array<T,Eigen::Dynamic,1> x(k_choose_2);
1422  Eigen::Array<T,Eigen::Dynamic,1> sds(k);
1423  bool successful = factor_cov_matrix(x,sds,y);
1424  if (!successful)
1425  throw std::runtime_error("factor_cov_matrix failed on y");
1426  for (size_type i = 0; i < k; ++i) {
1427  // sds on log scale unconstrained
1428  if (fabs(sds[i] - 0.0) >= CONSTRAINT_TOLERANCE) {
1429  std::stringstream s;
1430  s << "all standard deviations must be zero."
1431  << " found log(sd[" << i << "])=" << sds[i] << std::endl;
1432  BOOST_THROW_EXCEPTION(std::runtime_error(s.str()));
1433  }
1434  }
1435  return x.matrix();
1436  }
1437 
1438 
1439  // COVARIANCE MATRIX
1440 
1453  template <typename T>
1454  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>
1455  cov_matrix_constrain(const Eigen::Matrix<T,Eigen::Dynamic,1>& x,
1456  typename Eigen::Matrix<T,Eigen::Dynamic,1>::size_type K) {
1457  using std::exp;
1458  typedef typename Eigen::Matrix<T,Eigen::Dynamic,1>::size_type size_type;
1459  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> L(K,K);
1460  if (x.size() != (K * (K + 1)) / 2)
1461  throw std::domain_error("x.size() != K + (K choose 2)");
1462  int i = 0;
1463  for (size_type m = 0; m < K; ++m) {
1464  for (int n = 0; n < m; ++n)
1465  L(m,n) = x(i++);
1466  L(m,m) = exp(x(i++));
1467  for (size_type n = m + 1; n < K; ++n)
1468  L(m,n) = 0.0;
1469  }
1470  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> M = L * L.transpose();
1471  return L * L.transpose();
1472  }
1473 
1474 
1487  template <typename T>
1488  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>
1489  cov_matrix_constrain(const Eigen::Matrix<T,Eigen::Dynamic,1>& x,
1490  typename Eigen::Matrix<T,Eigen::Dynamic,1>::size_type K,
1491  T& lp) {
1492  using std::exp;
1493  typedef typename Eigen::Matrix<T,Eigen::Dynamic,1>::size_type size_type;
1494  if (x.size() != (K * (K + 1)) / 2)
1495  throw std::domain_error("x.size() != K + (K choose 2)");
1496  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> L(K,K);
1497  int i = 0;
1498  for (size_type m = 0; m < K; ++m) {
1499  for (size_type n = 0; n < m; ++n)
1500  L(m,n) = x(i++);
1501  L(m,m) = exp(x(i++));
1502  for (size_type n = m + 1; n < K; ++n)
1503  L(m,n) = 0.0;
1504  }
1505  // Jacobian for complete transform, including exp() above
1506  lp += (K * stan::math::LOG_2); // needless constant; want propto
1507  for (int k = 0; k < K; ++k)
1508  lp += (K - k + 1) * log(L(k,k)); // only +1 because index from 0
1509  return L * L.transpose();
1510  // return tri_multiply_transpose(L);
1511  }
1512 
1535  template <typename T>
1536  Eigen::Matrix<T,Eigen::Dynamic,1>
1537  cov_matrix_free(const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& y) {
1538  using std::log;
1539  int K = y.rows();
1540  if (y.cols() != K)
1541  throw std::domain_error("y is not a square matrix");
1542  if (K == 0)
1543  throw std::domain_error("y has no elements");
1544  for (int k = 0; k < K; ++k)
1545  if (!(y(k,k) > 0.0))
1546  throw std::domain_error("y has non-positive diagonal");
1547  Eigen::Matrix<T,Eigen::Dynamic,1> x((K * (K + 1)) / 2);
1548  // FIXME: see Eigen LDLT for rank-revealing version -- use that
1549  // even if less efficient?
1550  Eigen::LLT<Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> >
1551  llt(y.rows());
1552  llt.compute(y);
1553  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> L = llt.matrixL();
1554  int i = 0;
1555  for (int m = 0; m < K; ++m) {
1556  for (int n = 0; n < m; ++n)
1557  x(i++) = L(m,n);
1558  x(i++) = log(L(m,m));
1559  }
1560  return x;
1561  }
1562 
1582  template <typename T>
1583  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>
1584  cov_matrix_constrain_lkj(const Eigen::Matrix<T,Eigen::Dynamic,1>& x,
1585  size_t k) {
1586  size_t k_choose_2 = (k * (k - 1)) / 2;
1587  Eigen::Array<T,Eigen::Dynamic,1> cpcs(k_choose_2);
1588  int pos = 0;
1589  for (size_t i = 0; i < k_choose_2; ++i)
1590  cpcs[i] = corr_constrain(x[pos++]);
1591  Eigen::Array<T,Eigen::Dynamic,1> sds(k);
1592  for (size_t i = 0; i < k; ++i)
1593  sds[i] = positive_constrain(x[pos++]);
1594  return read_cov_matrix(cpcs, sds);
1595  }
1596 
1621  template <typename T>
1622  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>
1623  cov_matrix_constrain_lkj(const Eigen::Matrix<T,Eigen::Dynamic,1>& x,
1624  size_t k,
1625  T& lp) {
1626  size_t k_choose_2 = (k * (k - 1)) / 2;
1627  Eigen::Array<T,Eigen::Dynamic,1> cpcs(k_choose_2);
1628  int pos = 0;
1629  for (size_t i = 0; i < k_choose_2; ++i)
1630  cpcs[i] = corr_constrain(x[pos++], lp);
1631  Eigen::Array<T,Eigen::Dynamic,1> sds(k);
1632  for (size_t i = 0; i < k; ++i)
1633  sds[i] = positive_constrain(x[pos++],lp);
1634  return read_cov_matrix(cpcs, sds, lp);
1635  }
1636 
1655  template <typename T>
1656  Eigen::Matrix<T,Eigen::Dynamic,1>
1657  cov_matrix_free_lkj(
1658  const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& y) {
1659  typedef typename
1660  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>::size_type size_type;
1661  size_type k = y.rows();
1662  if (y.cols() != k)
1663  throw std::domain_error("y is not a square matrix");
1664  if (k == 0)
1665  throw std::domain_error("y has no elements");
1666  size_type k_choose_2 = (k * (k-1)) / 2;
1667  Eigen::Array<T,Eigen::Dynamic,1> cpcs(k_choose_2);
1668  Eigen::Array<T,Eigen::Dynamic,1> sds(k);
1669  bool successful = factor_cov_matrix(cpcs,sds,y);
1670  if (!successful)
1671  throw std::runtime_error ("factor_cov_matrix failed on y");
1672  Eigen::Matrix<T,Eigen::Dynamic,1> x(k_choose_2 + k);
1673  size_type pos = 0;
1674  for (size_type i = 0; i < k_choose_2; ++i)
1675  x[pos++] = cpcs[i];
1676  for (size_type i = 0; i < k; ++i)
1677  x[pos++] = sds[i];
1678  return x;
1679  }
1680 
1681  }
1682 
1683 }
1684 
1685 #endif
     [ Stan Home Page ] © 2011–2013, Stan Development Team.