Stan  1.0
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/constants.hpp>
14 #include <stan/math/matrix.hpp>
15 #include <stan/math/error_handling.hpp>
16 #include <stan/math/matrix_error_handling.hpp>
17 #include <stan/math/special_functions.hpp>
18 // #include <stan/agrad/special_functions.hpp>
19 
20 namespace stan {
21 
22  namespace prob {
23 
24 
25  const double CONSTRAINT_TOLERANCE = 1E-8;
26 
27 
40  template<typename T>
41  bool
42  factor_cov_matrix(Eigen::Array<T,Eigen::Dynamic,1>& CPCs,
43  Eigen::Array<T,Eigen::Dynamic,1>& sds,
44  const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& Sigma) {
45 
46  size_t K = sds.rows();
47 
48  sds = Sigma.diagonal().array();
49  if( (sds <= 0.0).any() ) return false;
50  sds = sds.sqrt();
51 
52  Eigen::DiagonalMatrix<T,Eigen::Dynamic> D(K);
53  D.diagonal() = sds.inverse();
54  sds = sds.log(); // now unbounded
55 
56  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> R = D * Sigma * D;
57  // to hopefully prevent pivoting due to floating point error
58  R.diagonal().setOnes();
59  Eigen::LDLT<Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> > ldlt;
60  ldlt = R.ldlt();
61  if( !ldlt.isPositive() ) return false;
62  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> U = ldlt.matrixU();
63 
64  size_t position = 0;
65  size_t pull = K - 1;
66 
67  Eigen::Array<T,Eigen::Dynamic,1> temp = U.row(0).tail(pull);
68  CPCs.head(pull) = temp;
69 
70  Eigen::Array<T,Eigen::Dynamic,1> acc(K);
71  acc(0) = -0.0;
72  acc.tail(pull) = 1.0 - temp.square();
73  for(size_t i = 1; i < (K - 1); i++) {
74  position += pull;
75  pull--;
76  temp = U.row(i).tail(pull);
77  temp /= sqrt(acc.tail(pull) / acc(i));
78  CPCs.segment(position, pull) = temp;
79  acc.tail(pull) *= 1.0 - temp.square();
80  }
81  CPCs = 0.5 * ( (1.0 + CPCs) / (1.0 - CPCs) ).log(); // now unbounded
82  return true;
83  }
84 
85  // MATRIX TRANSFORMS +/- JACOBIANS
86 
108  template <typename T>
109  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>
110  read_corr_L(const Eigen::Array<T,Eigen::Dynamic,1>& CPCs, // on (-1,1)
111  const size_t K) {
112  Eigen::Array<T,Eigen::Dynamic,1> temp;
113  Eigen::Array<T,Eigen::Dynamic,1> acc(K-1);
114  acc.setOnes();
115  // Cholesky factor of correlation matrix
116  Eigen::Array<T,Eigen::Dynamic,Eigen::Dynamic> L(K,K);
117  L.setZero();
118 
119  size_t position = 0;
120  size_t pull = K - 1;
121 
122  L(0,0) = 1.0;
123  L.col(0).tail(pull) = temp = CPCs.head(pull);
124  acc.tail(pull) = 1.0 - temp.square();
125  for(size_t i = 1; i < (K - 1); i++) {
126  position += pull;
127  pull--;
128  temp = CPCs.segment(position, pull);
129  L(i,i) = sqrt(acc(i-1));
130  L.col(i).tail(pull) = temp * acc.tail(pull).sqrt();
131  acc.tail(pull) *= 1.0 - temp.square();
132  }
133  L(K-1,K-1) = sqrt(acc(K-2));
134  return L.matrix();
135  }
136 
150  template <typename T>
151  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>
152  read_corr_matrix(const Eigen::Array<T,Eigen::Dynamic,1>& CPCs,
153  const size_t K) {
154  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> L
155  = read_corr_L(CPCs, K);
156  using stan::math::multiply_lower_tri_self_transpose;
157  return multiply_lower_tri_self_transpose(L);
158  }
159 
186  template <typename T>
187  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>
188  read_corr_L(const Eigen::Array<T,Eigen::Dynamic,1>& CPCs,
189  const size_t K,
190  T& log_prob) {
191 
192  size_t k = 0;
193  size_t i = 0;
194  T log_1cpc2;
195  double lead = K - 2.0;
196  // no need to abs() because this Jacobian determinant
197  // is strictly positive (and triangular)
198  // skip last row (odd indexing) because it adds nothing by design
199  typedef typename Eigen::Matrix<T,Eigen::Dynamic,1>::size_type size_type;
200  for (size_type j = 0;
201  j < (CPCs.rows() - 1);
202  ++j) {
203  using stan::math::log1m;
204  using stan::math::square;
205  log_1cpc2 = log1m(square(CPCs[j]));
206  // derivative of correlation wrt CPC
207  log_prob += lead / 2.0 * log_1cpc2;
208  i++;
209  if (i > K) {
210  k++;
211  i = k + 1;
212  lead = K - k - 1.0;
213  }
214  }
215  return read_corr_L(CPCs, K);
216  }
217 
236  template <typename T>
237  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>
238  read_corr_matrix(const Eigen::Array<T,Eigen::Dynamic,1>& CPCs,
239  const size_t K,
240  T& log_prob) {
241 
242  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> L
243  = read_corr_L(CPCs, K, log_prob);
244  using stan::math::multiply_lower_tri_self_transpose;
245  return multiply_lower_tri_self_transpose(L);
246  }
247 
258  template <typename T>
259  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>
260  read_cov_L(const Eigen::Array<T,Eigen::Dynamic,1>& CPCs,
261  const Eigen::Array<T,Eigen::Dynamic,1>& sds,
262  T& log_prob) {
263  size_t K = sds.rows();
264  // adjust due to transformation from correlations to covariances
265  log_prob += (sds.log().sum() + stan::math::LOG_2) * K;
266  return sds.matrix().asDiagonal() * read_corr_L(CPCs, K, log_prob);
267  }
268 
278  template <typename T>
279  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>
280  read_cov_matrix(const Eigen::Array<T,Eigen::Dynamic,1>& CPCs,
281  const Eigen::Array<T,Eigen::Dynamic,1>& sds,
282  T& log_prob) {
283 
284  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> L
285  = read_cov_L(CPCs, sds, log_prob);
286  using stan::math::multiply_lower_tri_self_transpose;
287  return multiply_lower_tri_self_transpose(L);
288  }
289 
297  template<typename T>
298  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>
299  read_cov_matrix(const Eigen::Array<T,Eigen::Dynamic,1>& CPCs,
300  const Eigen::Array<T,Eigen::Dynamic,1>& sds) {
301 
302  size_t K = sds.rows();
303  Eigen::DiagonalMatrix<T,Eigen::Dynamic> D(K);
304  D.diagonal() = sds;
305  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> L
306  = D * read_corr_L(CPCs, K);
307  using stan::math::multiply_lower_tri_self_transpose;
308  return multiply_lower_tri_self_transpose(L);
309  }
310 
311 
321  template<typename T>
322  const Eigen::Array<T,Eigen::Dynamic,1>
323  make_nu(const T eta, const size_t K) {
324 
325  Eigen::Array<T,Eigen::Dynamic,1> nu(K * (K - 1) / 2);
326 
327  T alpha = eta + (K - 2.0) / 2.0; // from Lewandowski et. al.
328 
329  // Best (1978) implies nu = 2 * alpha for the dof in a t
330  // distribution that generates a beta variate on (-1,1)
331  T alpha2 = 2.0 * alpha;
332 
333  typedef typename Eigen::Matrix<T,Eigen::Dynamic,1>::size_type size_type;
334  for (size_type j = 0; j < (K - 1); j++) {
335  nu(j) = alpha2;
336  }
337  size_t counter = K - 1;
338  for (size_type i = 1; i < (K - 1); i++) {
339  alpha -= 0.5;
340  alpha2 = 2.0 * alpha;
341  for (size_type j = i + 1; j < K; j++) {
342  nu(counter) = alpha2;
343  counter++;
344  }
345  }
346  return nu;
347  }
348 
349 
350 
351 
352  // IDENTITY
353 
366  template <typename T>
367  inline
368  T identity_constrain(T x) {
369  return x;
370  }
371 
385  template <typename T>
386  inline
387  T identity_constrain(const T x, T& lp) {
388  return x;
389  }
390 
402  template <typename T>
403  inline
404  T identity_free(const T y) {
405  return y;
406  }
407 
408 
409  // POSITIVE
410 
421  template <typename T>
422  inline
423  T positive_constrain(const T x) {
424  return exp(x);
425  }
426 
443  template <typename T>
444  inline
445  T positive_constrain(const T x, T& lp) {
446  lp += x;
447  return exp(x);
448  }
449 
466  template <typename T>
467  inline
468  T positive_free(const T y) {
469  stan::math::check_positive("stan::prob::positive_free(%1%)", y, "Positive variable");
470  return log(y);
471  }
472 
473  // LOWER BOUND
474 
494  template <typename T, typename TL>
495  inline
496  T lb_constrain(const T x, const TL lb) {
497  if (lb == -std::numeric_limits<double>::infinity())
498  return identity_constrain(x);
499  return exp(x) + lb;
500  }
501 
518  template <typename T, typename TL>
519  inline
520  typename boost::math::tools::promote_args<T,TL>::type
521  lb_constrain(const T x, const TL lb, T& lp) {
522  if (lb == -std::numeric_limits<double>::infinity())
523  return identity_constrain(x,lp);
524  lp += x;
525  return exp(x) + lb;
526  }
527 
543  template <typename T, typename TL>
544  inline
545  typename boost::math::tools::promote_args<T,TL>::type
546  lb_free(const T y, const TL lb) {
547  if (lb == -std::numeric_limits<double>::infinity())
548  return identity_free(y);
549  stan::math::check_greater_or_equal("stan::prob::lb_free(%1%)",
550  y, lb, "Lower bounded variable");
551  return log(y - lb);
552  }
553 
554 
555  // UPPER BOUND
556 
576  template <typename T, typename TU>
577  inline
578  typename boost::math::tools::promote_args<T,TU>::type
579  ub_constrain(const T x, const TU ub) {
580  if (ub == std::numeric_limits<double>::infinity())
581  return identity_constrain(x);
582  return ub - exp(x);
583  }
584 
608  template <typename T, typename TU>
609  inline
610  typename boost::math::tools::promote_args<T,TU>::type
611  ub_constrain(const T x, const TU ub, T& lp) {
612  if (ub == std::numeric_limits<double>::infinity())
613  return identity_constrain(x,lp);
614  lp += x;
615  return ub - exp(x);
616  }
617 
640  template <typename T, typename TU>
641  inline
642  typename boost::math::tools::promote_args<T,TU>::type
643  ub_free(const T y, const TU ub) {
644  if (ub == std::numeric_limits<double>::infinity())
645  return identity_free(y);
646  stan::math::check_less_or_equal("stan::prob::ub_free(%1%)",
647  y, ub, "Upper bounded variable");
648  return log(ub - y);
649  }
650 
651 
652  // LOWER & UPPER BOUNDS
653 
681  template <typename T, typename TL, typename TU>
682  inline
683  typename boost::math::tools::promote_args<T,TL,TU>::type
684  lub_constrain(const T x, TL lb, TU ub) {
685  stan::math::validate_less(lb,ub,"lb","ub","lub_constrain/3");
686 
687  if (lb == -std::numeric_limits<double>::infinity())
688  return ub_constrain(x,ub);
689  if (ub == std::numeric_limits<double>::infinity())
690  return lb_constrain(x,lb);
691 
692  T inv_logit_x;
693  if (x > 0) {
694  T exp_minus_x = exp(-x);
695  inv_logit_x = 1.0 / (1.0 + exp_minus_x);
696  // Prevent x from reaching one unless it really really should.
697  if ((x < std::numeric_limits<double>::infinity())
698  && (inv_logit_x == 1))
699  inv_logit_x = 1 - 1e-15;
700  } else {
701  T exp_x = exp(x);
702  inv_logit_x = 1.0 - 1.0 / (1.0 + exp_x);
703  // Prevent x from reaching zero unless it really really should.
704  if ((x > -std::numeric_limits<double>::infinity())
705  && (inv_logit_x== 0))
706  inv_logit_x = 1e-15;
707  }
708  return lb + (ub - lb) * inv_logit_x;
709  }
710 
752  template <typename T, typename TL, typename TU>
753  typename boost::math::tools::promote_args<T,TL,TU>::type
754  lub_constrain(const T x, const TL lb, const TU ub, T& lp) {
755  if (!(lb < ub)) {
756  std::stringstream s;
757  s << "domain error in lub_constrain; lower bound = " << lb
758  << " must be strictly less than upper bound = " << ub;
759  throw std::domain_error(s.str());
760  }
761  if (lb == -std::numeric_limits<double>::infinity())
762  return ub_constrain(x,ub,lp);
763  if (ub == std::numeric_limits<double>::infinity())
764  return lb_constrain(x,lb,lp);
765  T inv_logit_x;
766  if (x > 0) {
767  T exp_minus_x = exp(-x);
768  inv_logit_x = 1.0 / (1.0 + exp_minus_x);
769  lp += log(ub - lb) - x - 2 * log1p(exp_minus_x);
770  // Prevent x from reaching one unless it really really should.
771  if ((x < std::numeric_limits<double>::infinity())
772  && (inv_logit_x == 1))
773  inv_logit_x = 1 - 1e-15;
774  } else {
775  T exp_x = exp(x);
776  inv_logit_x = 1.0 - 1.0 / (1.0 + exp_x);
777  lp += log(ub - lb) + x - 2 * log1p(exp_x);
778  // Prevent x from reaching zero unless it really really should.
779  if ((x > -std::numeric_limits<double>::infinity())
780  && (inv_logit_x== 0))
781  inv_logit_x = 1e-15;
782  }
783  return lb + (ub - lb) * inv_logit_x;
784  }
785 
816  template <typename T, typename TL, typename TU>
817  inline
818  typename boost::math::tools::promote_args<T,TL,TU>::type
819  lub_free(const T y, TL lb, TU ub) {
820  using stan::math::logit;
821  stan::math::check_bounded("stan::prob::lub_free(%1%)",
822  y, lb, ub, "Bounded variable");
823  if (lb == -std::numeric_limits<double>::infinity())
824  return ub_free(y,ub);
825  if (ub == std::numeric_limits<double>::infinity())
826  return lb_free(y,lb);
827  return logit((y - lb) / (ub - lb));
828  }
829 
830 
831  // PROBABILITY
832 
846  template <typename T>
847  inline
848  T prob_constrain(const T x) {
849  using stan::math::inv_logit;
850  return inv_logit(x);
851  }
852 
874  template <typename T>
875  inline
876  T prob_constrain(const T x, T& lp) {
877  using stan::math::inv_logit;
878  using stan::math::log1m;
879  T inv_logit_x = inv_logit(x);
880  lp += log(inv_logit_x) + log1m(inv_logit_x);
881  return inv_logit_x;
882  }
883 
898  template <typename T>
899  inline
900  T prob_free(const T y) {
901  using stan::math::logit;
902  stan::math::check_bounded("stan::prob::prob_free(%1%)",
903  y, 0, 1, "Probability variable");
904  return logit(y);
905  }
906 
907 
908  // CORRELATION
909 
922  template <typename T>
923  inline
924  T corr_constrain(const T x) {
925  return tanh(x);
926  }
927 
940  template <typename T>
941  inline
942  T corr_constrain(const T x, T& lp) {
943  using stan::math::log1m;
944  T tanh_x = tanh(x);
945  lp += log1m(tanh_x * tanh_x);
946  return tanh_x;
947  }
948 
965  template <typename T>
966  inline
967  T corr_free(const T y) {
968  stan::math::check_bounded("stan::prob::lub_free(%1%)",
969  y, -1, 1, "Correlation variable");
970  return atanh(y);
971  }
972 
973 
974  // SIMPLEX
975 
976 
989  template <typename T>
990  Eigen::Matrix<T,Eigen::Dynamic,1>
991  simplex_constrain(const Eigen::Matrix<T,Eigen::Dynamic,1>& y) {
992  // cut & paste simplex_constrain(Eigen::Matrix,T) w/o Jacobian
993  typedef typename Eigen::Matrix<T,Eigen::Dynamic,1>::size_type size_type;
994  using stan::math::logit;
995  using stan::math::inv_logit;
996  using stan::math::log1m;
997  int Km1 = y.size();
998  Eigen::Matrix<T,Eigen::Dynamic,1> x(Km1 + 1);
999  T stick_len(1.0);
1000  for (size_type k = 0; k < Km1; ++k) {
1001  T z_k(inv_logit(y(k) - log(Km1 - k)));
1002  x(k) = stick_len * z_k;
1003  stick_len -= x(k);
1004  }
1005  x(Km1) = stick_len;
1006  return x;
1007  }
1008 
1022  template <typename T>
1023  Eigen::Matrix<T,Eigen::Dynamic,1>
1024  simplex_constrain(const Eigen::Matrix<T,Eigen::Dynamic,1>& y,
1025  T& lp) {
1026  using stan::math::logit;
1027  using stan::math::inv_logit;
1028  using stan::math::log1p_exp;
1029  using stan::math::log1m;
1030  typedef typename Eigen::Matrix<T,Eigen::Dynamic,1>::size_type size_type;
1031  int Km1 = y.size(); // K = Km1 + 1
1032  Eigen::Matrix<T,Eigen::Dynamic,1> x(Km1 + 1);
1033  T stick_len(1.0);
1034  for (size_type k = 0; k < Km1; ++k) {
1035  double eq_share = -log(Km1 - k); // = logit(1.0/(Km1 + 1 - k));
1036  T adj_y_k(y(k) + eq_share);
1037  T z_k(inv_logit(adj_y_k));
1038  x(k) = stick_len * z_k;
1039  lp += log(stick_len);
1040  lp -= log1p_exp(-adj_y_k);
1041  lp -= log1p_exp(adj_y_k);
1042  stick_len -= x(k); // equivalently *= (1 - z_k);
1043  }
1044  x(Km1) = stick_len; // no Jacobian contrib for last dim
1045  return x;
1046  }
1047 
1062  template <typename T>
1063  Eigen::Matrix<T,Eigen::Dynamic,1>
1064  simplex_free(const Eigen::Matrix<T,Eigen::Dynamic,1>& x) {
1065  using stan::math::logit;
1066  typedef typename Eigen::Matrix<T,Eigen::Dynamic,1>::size_type size_type;
1067  stan::math::check_simplex("stan::prob::simplex_free(%1%)", x, "Simplex variable");
1068  int Km1 = x.size() - 1;
1069  Eigen::Matrix<T,Eigen::Dynamic,1> y(Km1);
1070  T stick_len(x(Km1));
1071  for (size_type k = Km1; --k >= 0; ) {
1072  stick_len += x(k);
1073  T z_k(x(k) / stick_len);
1074  y(k) = logit(z_k) + log(Km1 - k);
1075  // log(Km-k) = logit(1.0 / (Km1 + 1 - k));
1076  }
1077  return y;
1078  }
1079 
1080 
1081  // ORDERED
1082 
1092  template <typename T>
1093  Eigen::Matrix<T,Eigen::Dynamic,1>
1094  ordered_constrain(const Eigen::Matrix<T,Eigen::Dynamic,1>& x) {
1095  typedef typename Eigen::Matrix<T,Eigen::Dynamic,1>::size_type size_type;
1096  size_type k = x.size();
1097  Eigen::Matrix<T,Eigen::Dynamic,1> y(k);
1098  if (k == 0)
1099  return y;
1100  y[0] = x[0];
1101  for (size_type i = 1;
1102  i < k;
1103  ++i)
1104  y[i] = y[i-1] + exp(x[i]);
1105  return y;
1106  }
1107 
1120  template <typename T>
1121  inline
1122  Eigen::Matrix<T,Eigen::Dynamic,1>
1123  ordered_constrain(const Eigen::Matrix<T,Eigen::Dynamic,1>& x, T& lp) {
1124  typedef typename Eigen::Matrix<T,Eigen::Dynamic,1>::size_type size_type;
1125  for (size_type i = 1; i < x.size(); ++i)
1126  lp += x(i);
1127  return ordered_constrain(x);
1128  }
1129 
1130 
1131 
1145  template <typename T>
1146  Eigen::Matrix<T,Eigen::Dynamic,1>
1147  ordered_free(const Eigen::Matrix<T,Eigen::Dynamic,1>& y) {
1148  stan::math::check_ordered("stan::prob::ordered_free(%1%)",
1149  y, "Ordered variable");
1150  typedef typename Eigen::Matrix<T,Eigen::Dynamic,1>::size_type size_type;
1151  size_type k = y.size();
1152  Eigen::Matrix<T,Eigen::Dynamic,1> x(k);
1153  if (k == 0)
1154  return x;
1155  x[0] = y[0];
1156  for (size_type i = 1; i < k; ++i)
1157  x[i] = log(y[i] - y[i-1]);
1158  return x;
1159  }
1160 
1161 
1162  // POSITIVE ORDERED
1163 
1173  template <typename T>
1174  Eigen::Matrix<T,Eigen::Dynamic,1>
1175  positive_ordered_constrain(const Eigen::Matrix<T,Eigen::Dynamic,1>& x) {
1176  typedef typename Eigen::Matrix<T,Eigen::Dynamic,1>::size_type size_type;
1177  size_type k = x.size();
1178  Eigen::Matrix<T,Eigen::Dynamic,1> y(k);
1179  if (k == 0)
1180  return y;
1181  y[0] = exp(x[0]);
1182  for (size_type i = 1;
1183  i < k;
1184  ++i)
1185  y[i] = y[i-1] + exp(x[i]);
1186  return y;
1187  }
1188 
1201  template <typename T>
1202  inline
1203  Eigen::Matrix<T,Eigen::Dynamic,1>
1204  positive_ordered_constrain(const Eigen::Matrix<T,Eigen::Dynamic,1>& x, T& lp) {
1205  typedef typename Eigen::Matrix<T,Eigen::Dynamic,1>::size_type size_type;
1206  for (size_type i = 0; i < x.size(); ++i)
1207  lp += x(i);
1208  return positive_ordered_constrain(x);
1209  }
1210 
1211 
1212 
1226  template <typename T>
1227  Eigen::Matrix<T,Eigen::Dynamic,1>
1228  positive_ordered_free(const Eigen::Matrix<T,Eigen::Dynamic,1>& y) {
1229  stan::math::check_positive_ordered("stan::prob::positive_ordered_free(%1%)",
1230  y, "Positive ordered variable");
1231  typedef typename Eigen::Matrix<T,Eigen::Dynamic,1>::size_type size_type;
1232  size_type k = y.size();
1233  Eigen::Matrix<T,Eigen::Dynamic,1> x(k);
1234  if (k == 0)
1235  return x;
1236  x[0] = log(y[0]);
1237  for (size_type i = 1; i < k; ++i)
1238  x[i] = log(y[i] - y[i-1]);
1239  return x;
1240  }
1241 
1242 
1243  // CORRELATION MATRIX
1268  template <typename T>
1269  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>
1270  corr_matrix_constrain(const Eigen::Matrix<T,Eigen::Dynamic,1>& x,
1271  typename Eigen::Matrix<T,Eigen::Dynamic,1>::size_type k) {
1272  typedef typename Eigen::Matrix<T,Eigen::Dynamic,1>::size_type size_type;
1273  size_type k_choose_2 = (k * (k - 1)) / 2;
1274  if (k_choose_2 != x.size())
1275  throw std::invalid_argument ("x is not a valid correlation matrix");
1276  Eigen::Array<T,Eigen::Dynamic,1> cpcs(k_choose_2);
1277  for (size_type i = 0; i < k_choose_2; ++i)
1278  cpcs[i] = corr_constrain(x[i]);
1279  return read_corr_matrix(cpcs,k);
1280  }
1281 
1301  template <typename T>
1302  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>
1303  corr_matrix_constrain(const Eigen::Matrix<T,Eigen::Dynamic,1>& x,
1304  typename Eigen::Matrix<T,Eigen::Dynamic,1>::size_type k,
1305  T& lp) {
1306  typedef typename Eigen::Matrix<T,Eigen::Dynamic,1>::size_type size_type;
1307  size_type k_choose_2 = (k * (k - 1)) / 2;
1308  if (k_choose_2 != x.size())
1309  throw std::invalid_argument ("x is not a valid correlation matrix");
1310 
1311  Eigen::Array<T,Eigen::Dynamic,1> cpcs(k_choose_2);
1312  for (size_type i = 0; i < k_choose_2; ++i)
1313  cpcs[i] = corr_constrain(x[i],lp);
1314  return read_corr_matrix(cpcs,k,lp);
1315  }
1316 
1337  template <typename T>
1338  Eigen::Matrix<T,Eigen::Dynamic,1>
1339  corr_matrix_free(const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& y) {
1340  typedef typename
1341  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>::size_type size_type;
1342  size_type k = y.rows();
1343  if (y.cols() != k)
1344  throw std::domain_error("y is not a square matrix");
1345  if (k == 0)
1346  throw std::domain_error("y has no elements");
1347  size_type k_choose_2 = (k * (k-1)) / 2;
1348  Eigen::Array<T,Eigen::Dynamic,1> x(k_choose_2);
1349  Eigen::Array<T,Eigen::Dynamic,1> sds(k);
1350  bool successful = factor_cov_matrix(x,sds,y);
1351  if (!successful)
1352  throw std::runtime_error("factor_cov_matrix failed on y");
1353  for (size_type i = 0; i < k; ++i) {
1354  // sds on log scale unconstrained
1355  if (fabs(sds[i] - 0.0) >= CONSTRAINT_TOLERANCE) {
1356  std::stringstream s;
1357  s << "all standard deviations must be zero."
1358  << " found log(sd[" << i << "])=" << sds[i] << std::endl;
1359  BOOST_THROW_EXCEPTION(std::runtime_error(s.str()));
1360  }
1361  }
1362  return x.matrix();
1363  }
1364 
1365 
1366  // COVARIANCE MATRIX
1367 
1380  template <typename T>
1381  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>
1382  cov_matrix_constrain(const Eigen::Matrix<T,Eigen::Dynamic,1>& x,
1383  typename Eigen::Matrix<T,Eigen::Dynamic,1>::size_type K) {
1384  using std::exp;
1385  typedef typename Eigen::Matrix<T,Eigen::Dynamic,1>::size_type size_type;
1386  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> L(K,K);
1387  if (x.size() != (K * (K + 1)) / 2)
1388  throw std::domain_error("x.size() != K + (K choose 2)");
1389  int i = 0;
1390  for (size_type m = 0; m < K; ++m) {
1391  for (int n = 0; n < m; ++n)
1392  L(m,n) = x(i++);
1393  L(m,m) = exp(x(i++));
1394  for (size_type n = m + 1; n < K; ++n)
1395  L(m,n) = 0.0;
1396  }
1397  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> M = L * L.transpose();
1398  return L * L.transpose();
1399  }
1400 
1401 
1414  template <typename T>
1415  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>
1416  cov_matrix_constrain(const Eigen::Matrix<T,Eigen::Dynamic,1>& x,
1417  typename Eigen::Matrix<T,Eigen::Dynamic,1>::size_type K,
1418  T& lp) {
1419  using std::exp;
1420  typedef typename Eigen::Matrix<T,Eigen::Dynamic,1>::size_type size_type;
1421  if (x.size() != (K * (K + 1)) / 2)
1422  throw std::domain_error("x.size() != K + (K choose 2)");
1423  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> L(K,K);
1424  int i = 0;
1425  for (size_type m = 0; m < K; ++m) {
1426  for (size_type n = 0; n < m; ++n)
1427  L(m,n) = x(i++);
1428  L(m,m) = exp(x(i++));
1429  for (size_type n = m + 1; n < K; ++n)
1430  L(m,n) = 0.0;
1431  }
1432  // Jacobian for complete transform, including exp() above
1433  lp += (K * stan::math::LOG_2); // needless constant; want propto
1434  for (int k = 0; k < K; ++k)
1435  lp += (K - k + 1) * log(L(k,k)); // only +1 because index from 0
1436  return L * L.transpose();
1437  // return tri_multiply_transpose(L);
1438  }
1439 
1462  template <typename T>
1463  Eigen::Matrix<T,Eigen::Dynamic,1>
1464  cov_matrix_free(const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& y) {
1465  using std::log;
1466  int K = y.rows();
1467  if (y.cols() != K)
1468  throw std::domain_error("y is not a square matrix");
1469  if (K == 0)
1470  throw std::domain_error("y has no elements");
1471  for (int k = 0; k < K; ++k)
1472  if (!(y(k,k) > 0.0))
1473  throw std::domain_error("y has non-positive diagonal");
1474  Eigen::Matrix<T,Eigen::Dynamic,1> x((K * (K + 1)) / 2);
1475  // FIXME: see Eigen LDLT for rank-revealing version -- use that
1476  // even if less efficient?
1477  Eigen::LLT<Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> >
1478  llt(y.rows());
1479  llt.compute(y);
1480  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> L = llt.matrixL();
1481  int i = 0;
1482  for (int m = 0; m < K; ++m) {
1483  for (int n = 0; n < m; ++n)
1484  x(i++) = L(m,n);
1485  x(i++) = log(L(m,m));
1486  }
1487  return x;
1488  }
1489 
1509  template <typename T>
1510  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>
1511  cov_matrix_constrain_lkj(const Eigen::Matrix<T,Eigen::Dynamic,1>& x,
1512  size_t k) {
1513  size_t k_choose_2 = (k * (k - 1)) / 2;
1514  Eigen::Array<T,Eigen::Dynamic,1> cpcs(k_choose_2);
1515  int pos = 0;
1516  for (size_t i = 0; i < k_choose_2; ++i)
1517  cpcs[i] = corr_constrain(x[pos++]);
1518  Eigen::Array<T,Eigen::Dynamic,1> sds(k);
1519  for (size_t i = 0; i < k; ++i)
1520  sds[i] = positive_constrain(x[pos++]);
1521  return read_cov_matrix(cpcs, sds);
1522  }
1523 
1548  template <typename T>
1549  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>
1550  cov_matrix_constrain_lkj(const Eigen::Matrix<T,Eigen::Dynamic,1>& x,
1551  size_t k,
1552  T& lp) {
1553  size_t k_choose_2 = (k * (k - 1)) / 2;
1554  Eigen::Array<T,Eigen::Dynamic,1> cpcs(k_choose_2);
1555  int pos = 0;
1556  for (size_t i = 0; i < k_choose_2; ++i)
1557  cpcs[i] = corr_constrain(x[pos++], lp);
1558  Eigen::Array<T,Eigen::Dynamic,1> sds(k);
1559  for (size_t i = 0; i < k; ++i)
1560  sds[i] = positive_constrain(x[pos++],lp);
1561  return read_cov_matrix(cpcs, sds, lp);
1562  }
1563 
1582  template <typename T>
1583  Eigen::Matrix<T,Eigen::Dynamic,1>
1584  cov_matrix_free_lkj(
1585  const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& y) {
1586  typedef typename
1587  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>::size_type size_type;
1588  size_type k = y.rows();
1589  if (y.cols() != k)
1590  throw std::domain_error("y is not a square matrix");
1591  if (k == 0)
1592  throw std::domain_error("y has no elements");
1593  size_type k_choose_2 = (k * (k-1)) / 2;
1594  Eigen::Array<T,Eigen::Dynamic,1> cpcs(k_choose_2);
1595  Eigen::Array<T,Eigen::Dynamic,1> sds(k);
1596  bool successful = factor_cov_matrix(cpcs,sds,y);
1597  if (!successful)
1598  throw std::runtime_error ("factor_cov_matrix failed on y");
1599  Eigen::Matrix<T,Eigen::Dynamic,1> x(k_choose_2 + k);
1600  size_type pos = 0;
1601  for (size_type i = 0; i < k_choose_2; ++i)
1602  x[pos++] = cpcs[i];
1603  for (size_type i = 0; i < k; ++i)
1604  x[pos++] = sds[i];
1605  return x;
1606  }
1607 
1608  }
1609 
1610 }
1611 
1612 #endif
     [ Stan Home Page ] © 2011–2012, Stan Development Team.