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/math/constants.hpp>
13 #include <stan/math/matrix.hpp>
14 #include <stan/math/error_handling.hpp>
15 #include <stan/math/matrix_error_handling.hpp>
16 #include <stan/math/special_functions.hpp>
17 // #include <stan/agrad/special_functions.hpp>
18 
19 namespace stan {
20 
21  namespace prob {
22 
23 
24  const double CONSTRAINT_TOLERANCE = 1E-8;
25 
26 
39  template<typename T>
40  bool
41  factor_cov_matrix(Eigen::Array<T,Eigen::Dynamic,1>& CPCs,
42  Eigen::Array<T,Eigen::Dynamic,1>& sds,
43  const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& Sigma) {
44 
45  size_t K = sds.rows();
46 
47  sds = Sigma.diagonal().array();
48  if( (sds <= 0.0).any() ) return false;
49  sds = sds.sqrt();
50 
51  Eigen::DiagonalMatrix<T,Eigen::Dynamic> D(K);
52  D.diagonal() = sds.inverse();
53  sds = sds.log(); // now unbounded
54 
55  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> R = D * Sigma * D;
56  // to hopefully prevent pivoting due to floating point error
57  R.diagonal().setOnes();
58  Eigen::LDLT<Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> > ldlt;
59  ldlt = R.ldlt();
60  if( !ldlt.isPositive() ) return false;
61  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> U = ldlt.matrixU();
62 
63  size_t position = 0;
64  size_t pull = K - 1;
65 
66  Eigen::Array<T,Eigen::Dynamic,1> temp = U.row(0).tail(pull);
67  CPCs.head(pull) = temp;
68 
69  Eigen::Array<T,Eigen::Dynamic,1> acc(K);
70  acc(0) = -0.0;
71  acc.tail(pull) = 1.0 - temp.square();
72  for(size_t i = 1; i < (K - 1); i++) {
73  position += pull;
74  pull--;
75  temp = U.row(i).tail(pull);
76  temp /= sqrt(acc.tail(pull) / acc(i));
77  CPCs.segment(position, pull) = temp;
78  acc.tail(pull) *= 1.0 - temp.square();
79  }
80  CPCs = 0.5 * ( (1.0 + CPCs) / (1.0 - CPCs) ).log(); // now unbounded
81  return true;
82  }
83 
84  // MATRIX TRANSFORMS +/- JACOBIANS
85 
107  template <typename T>
108  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>
109  read_corr_L(const Eigen::Array<T,Eigen::Dynamic,1>& CPCs, // on (-1,1)
110  const size_t K) {
111  Eigen::Array<T,Eigen::Dynamic,1> temp;
112  Eigen::Array<T,Eigen::Dynamic,1> acc(K-1);
113  acc.setOnes();
114  // Cholesky factor of correlation matrix
115  Eigen::Array<T,Eigen::Dynamic,Eigen::Dynamic> L(K,K);
116  L.setZero();
117 
118  size_t position = 0;
119  size_t pull = K - 1;
120 
121  L(0,0) = 1.0;
122  L.col(0).tail(pull) = temp = CPCs.head(pull);
123  acc.tail(pull) = 1.0 - temp.square();
124  for(size_t i = 1; i < (K - 1); i++) {
125  position += pull;
126  pull--;
127  temp = CPCs.segment(position, pull);
128  L(i,i) = sqrt(acc(i-1));
129  L.col(i).tail(pull) = temp * acc.tail(pull).sqrt();
130  acc.tail(pull) *= 1.0 - temp.square();
131  }
132  L(K-1,K-1) = sqrt(acc(K-2));
133  return L.matrix();
134  }
135 
149  template <typename T>
150  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>
151  read_corr_matrix(const Eigen::Array<T,Eigen::Dynamic,1>& CPCs,
152  const size_t K) {
153  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> L
154  = read_corr_L(CPCs, K);
155  return L.template triangularView<Eigen::Lower>()
156  * L.matrix().transpose();
157  }
158 
185  template <typename T>
186  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>
187  read_corr_L(const Eigen::Array<T,Eigen::Dynamic,1>& CPCs,
188  const size_t K,
189  T& log_prob) {
190 
191  size_t k = 0;
192  size_t i = 0;
193  T log_1cpc2;
194  double lead = K - 2.0;
195  // no need to abs() because this Jacobian determinant
196  // is strictly positive (and triangular)
197  // skip last row (odd indexing) because it adds nothing by design
198  typedef typename Eigen::Matrix<T,Eigen::Dynamic,1>::size_type size_type;
199  for (size_type j = 0;
200  j < (CPCs.rows() - 1);
201  ++j) {
202  using stan::math::log1m;
203  using stan::math::square;
204  log_1cpc2 = log1m(square(CPCs[j]));
205  // derivative of correlation wrt CPC
206  log_prob += lead / 2.0 * log_1cpc2;
207  i++;
208  if (i > K) {
209  k++;
210  i = k + 1;
211  lead = K - k - 1.0;
212  }
213  }
214  return read_corr_L(CPCs, K);
215  }
216 
235  template <typename T>
236  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>
237  read_corr_matrix(const Eigen::Array<T,Eigen::Dynamic,1>& CPCs,
238  const size_t K,
239  T& log_prob) {
240 
241  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> L
242  = read_corr_L(CPCs, K, log_prob);
243  return L.template triangularView<Eigen::Lower>()
244  * L.matrix().transpose();
245  }
246 
257  template <typename T>
258  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>
259  read_cov_L(const Eigen::Array<T,Eigen::Dynamic,1>& CPCs,
260  const Eigen::Array<T,Eigen::Dynamic,1>& sds,
261  T& log_prob) {
262  size_t K = sds.rows();
263  // adjust due to transformation from correlations to covariances
264  log_prob += (sds.log().sum() + stan::math::LOG_2) * K;
265  return sds.matrix().asDiagonal() * read_corr_L(CPCs, K, log_prob);
266  }
267 
277  template <typename T>
278  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>
279  read_cov_matrix(const Eigen::Array<T,Eigen::Dynamic,1>& CPCs,
280  const Eigen::Array<T,Eigen::Dynamic,1>& sds,
281  T& log_prob) {
282 
283  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> L
284  = read_cov_L(CPCs, sds, log_prob);
285  return L.template triangularView<Eigen::Lower>()
286  * L.matrix().transpose();
287  }
288 
296  template<typename T>
297  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>
298  read_cov_matrix(const Eigen::Array<T,Eigen::Dynamic,1>& CPCs,
299  const Eigen::Array<T,Eigen::Dynamic,1>& sds) {
300 
301  size_t K = sds.rows();
302  Eigen::DiagonalMatrix<T,Eigen::Dynamic> D(K);
303  D.diagonal() = sds;
304  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> L
305  = D * read_corr_L(CPCs, K);
306  return L.template triangularView<Eigen::Lower>()
307  * L.matrix().transpose();
308  }
309 
310 
320  template<typename T>
321  const Eigen::Array<T,Eigen::Dynamic,1>
322  make_nu(const T eta, const size_t K) {
323 
324  Eigen::Array<T,Eigen::Dynamic,1> nu(K * (K - 1) / 2);
325 
326  T alpha = eta + (K - 2.0) / 2.0; // from Lewandowski et. al.
327 
328  // Best (1978) implies nu = 2 * alpha for the dof in a t
329  // distribution that generates a beta variate on (-1,1)
330  T alpha2 = 2.0 * alpha;
331 
332  typedef typename Eigen::Matrix<T,Eigen::Dynamic,1>::size_type size_type;
333  for (size_type j = 0; j < (K - 1); j++) {
334  nu(j) = alpha2;
335  }
336  size_t counter = K - 1;
337  for (size_type i = 1; i < (K - 1); i++) {
338  alpha -= 0.5;
339  alpha2 = 2.0 * alpha;
340  for (size_type j = i + 1; j < K; j++) {
341  nu(counter) = alpha2;
342  counter++;
343  }
344  }
345  return nu;
346  }
347 
348 
349 
350 
351  // IDENTITY
352 
365  template <typename T>
366  inline
367  T identity_constrain(T x) {
368  return x;
369  }
370 
384  template <typename T>
385  inline
386  T identity_constrain(const T x, T& lp) {
387  return x;
388  }
389 
401  template <typename T>
402  inline
403  T identity_free(const T y) {
404  return y;
405  }
406 
407 
408  // POSITIVE
409 
420  template <typename T>
421  inline
422  T positive_constrain(const T x) {
423  return exp(x);
424  }
425 
442  template <typename T>
443  inline
444  T positive_constrain(const T x, T& lp) {
445  lp += x;
446  return exp(x);
447  }
448 
465  template <typename T>
466  inline
467  T positive_free(const T y) {
468  stan::math::check_positive("stan::prob::positive_free(%1%)", y, "Positive variable");
469  return log(y);
470  }
471 
472  // LOWER BOUND
473 
490  template <typename T, typename TL>
491  inline
492  T lb_constrain(const T x, const TL lb) {
493  return exp(x) + lb;
494  }
495 
509  template <typename T, typename TL>
510  inline
511  typename boost::math::tools::promote_args<T,TL>::type
512  lb_constrain(const T x, const TL lb, T& lp) {
513  lp += x;
514  return exp(x) + lb;
515  }
516 
529  template <typename T, typename TL>
530  inline
531  typename boost::math::tools::promote_args<T,TL>::type
532  lb_free(const T y, const TL lb) {
533  stan::math::check_greater_or_equal("stan::prob::lb_free(%1%)",
534  y, lb, "Lower bounded variable");
535  return log(y - lb);
536  }
537 
538 
539  // UPPER BOUND
540 
557  template <typename T, typename TU>
558  inline
559  typename boost::math::tools::promote_args<T,TU>::type
560  ub_constrain(const T x, const TU ub) {
561  return ub - exp(x);
562  }
563 
584  template <typename T, typename TU>
585  inline
586  typename boost::math::tools::promote_args<T,TU>::type
587  ub_constrain(const T x, const TU ub, T& lp) {
588  lp += x;
589  return ub - exp(x);
590  }
591 
611  template <typename T, typename TU>
612  inline
613  typename boost::math::tools::promote_args<T,TU>::type
614  ub_free(const T y, const TU ub) {
615  stan::math::check_less_or_equal("stan::prob::ub_free(%1%)",
616  y, ub, "Upper bounded variable");
617  return log(ub - y);
618  }
619 
620 
621  // LOWER & UPPER BOUNDS
622 
642  template <typename T, typename TL, typename TU>
643  inline
644  typename boost::math::tools::promote_args<T,TL,TU>::type
645  lub_constrain(const T x, TL lb, TU ub) {
646  if (!(lb < ub)) {
647  std::stringstream s;
648  s << "domain error in lub_constrain; lower bound = " << lb
649  << " must be strictly less than upper bound = " << ub;
650  throw std::domain_error(s.str());
651  }
652  T inv_logit_x;
653  if (x > 0) {
654  T exp_minus_x = exp(-x);
655  inv_logit_x = 1.0 / (1.0 + exp_minus_x);
656  // Prevent x from reaching one unless it really really should.
657  if ((x < std::numeric_limits<double>::infinity())
658  && (inv_logit_x == 1))
659  inv_logit_x = 1 - 1e-15;
660  } else {
661  T exp_x = exp(x);
662  inv_logit_x = 1.0 - 1.0 / (1.0 + exp_x);
663  // Prevent x from reaching zero unless it really really should.
664  if ((x > -std::numeric_limits<double>::infinity())
665  && (inv_logit_x== 0))
666  inv_logit_x = 1e-15;
667  }
668  return lb + (ub - lb) * inv_logit_x;
669  }
670 
704  template <typename T, typename TL, typename TU>
705  typename boost::math::tools::promote_args<T,TL,TU>::type
706  lub_constrain(const T x, const TL lb, const TU ub, T& lp) {
707  if (!(lb < ub)) {
708  std::stringstream s;
709  s << "domain error in lub_constrain; lower bound = " << lb
710  << " must be strictly less than upper bound = " << ub;
711  throw std::domain_error(s.str());
712  }
713  T inv_logit_x;
714  if (x > 0) {
715  T exp_minus_x = exp(-x);
716  inv_logit_x = 1.0 / (1.0 + exp_minus_x);
717  lp += log(ub - lb) - x - 2 * log1p(exp_minus_x);
718  // Prevent x from reaching one unless it really really should.
719  if ((x < std::numeric_limits<double>::infinity())
720  && (inv_logit_x == 1))
721  inv_logit_x = 1 - 1e-15;
722  } else {
723  T exp_x = exp(x);
724  inv_logit_x = 1.0 - 1.0 / (1.0 + exp_x);
725  lp += log(ub - lb) + x - 2 * log1p(exp_x);
726  // Prevent x from reaching zero unless it really really should.
727  if ((x > -std::numeric_limits<double>::infinity())
728  && (inv_logit_x== 0))
729  inv_logit_x = 1e-15;
730  }
731  return lb + (ub - lb) * inv_logit_x;
732  }
733 
756  template <typename T, typename TL, typename TU>
757  inline
758  typename boost::math::tools::promote_args<T,TL,TU>::type
759  lub_free(const T y, TL lb, TU ub) {
760  using stan::math::logit;
761  stan::math::check_bounded("stan::prob::lub_free(%1%)",
762  y, lb, ub, "Bounded variable");
763  return logit((y - lb) / (ub - lb));
764  }
765 
766 
767  // PROBABILITY
768 
782  template <typename T>
783  inline
784  T prob_constrain(const T x) {
785  using stan::math::inv_logit;
786  return inv_logit(x);
787  }
788 
810  template <typename T>
811  inline
812  T prob_constrain(const T x, T& lp) {
813  using stan::math::inv_logit;
814  using stan::math::log1m;
815  T inv_logit_x = inv_logit(x);
816  lp += log(inv_logit_x) + log1m(inv_logit_x);
817  return inv_logit_x;
818  }
819 
834  template <typename T>
835  inline
836  T prob_free(const T y) {
837  using stan::math::logit;
838  stan::math::check_bounded("stan::prob::prob_free(%1%)",
839  y, 0, 1, "Probability variable");
840  return logit(y);
841  }
842 
843 
844  // CORRELATION
845 
858  template <typename T>
859  inline
860  T corr_constrain(const T x) {
861  return tanh(x);
862  }
863 
876  template <typename T>
877  inline
878  T corr_constrain(const T x, T& lp) {
879  using stan::math::log1m;
880  T tanh_x = tanh(x);
881  lp += log1m(tanh_x * tanh_x);
882  return tanh_x;
883  }
884 
901  template <typename T>
902  inline
903  T corr_free(const T y) {
904  stan::math::check_bounded("stan::prob::lub_free(%1%)",
905  y, -1, 1, "Correlation variable");
906  return atanh(y);
907  }
908 
909 
910  // SIMPLEX
911 
912 
925  template <typename T>
926  Eigen::Matrix<T,Eigen::Dynamic,1>
927  simplex_constrain(const Eigen::Matrix<T,Eigen::Dynamic,1>& y) {
928  // cut & paste simplex_constrain(Eigen::Matrix,T) w/o Jacobian
929  typedef typename Eigen::Matrix<T,Eigen::Dynamic,1>::size_type size_type;
930  using stan::math::logit;
931  using stan::math::inv_logit;
932  using stan::math::log1m;
933  int Km1 = y.size();
934  Eigen::Matrix<T,Eigen::Dynamic,1> x(Km1 + 1);
935  T stick_len(1.0);
936  for (size_type k = 0; k < Km1; ++k) {
937  T z_k(inv_logit(y(k) - log(Km1 - k)));
938  x(k) = stick_len * z_k;
939  stick_len -= x(k);
940  }
941  x(Km1) = stick_len;
942  return x;
943  }
944 
958  template <typename T>
959  Eigen::Matrix<T,Eigen::Dynamic,1>
960  simplex_constrain(const Eigen::Matrix<T,Eigen::Dynamic,1>& y,
961  T& lp) {
962  using stan::math::logit;
963  using stan::math::inv_logit;
964  using stan::math::log1p_exp;
965  using stan::math::log1m;
966  typedef typename Eigen::Matrix<T,Eigen::Dynamic,1>::size_type size_type;
967  int Km1 = y.size(); // K = Km1 + 1
968  Eigen::Matrix<T,Eigen::Dynamic,1> x(Km1 + 1);
969  T stick_len(1.0);
970  for (size_type k = 0; k < Km1; ++k) {
971  double eq_share = -log(Km1 - k); // = logit(1.0/(Km1 + 1 - k));
972  T adj_y_k(y(k) + eq_share);
973  T z_k(inv_logit(adj_y_k));
974  x(k) = stick_len * z_k;
975  lp += log(stick_len);
976  lp -= log1p_exp(-adj_y_k);
977  lp -= log1p_exp(adj_y_k);
978  stick_len -= x(k); // equivalently *= (1 - z_k);
979  }
980  x(Km1) = stick_len; // no Jacobian contrib for last dim
981  return x;
982  }
983 
998  template <typename T>
999  Eigen::Matrix<T,Eigen::Dynamic,1>
1000  simplex_free(const Eigen::Matrix<T,Eigen::Dynamic,1>& x) {
1001  using stan::math::logit;
1002  typedef typename Eigen::Matrix<T,Eigen::Dynamic,1>::size_type size_type;
1003  stan::math::check_simplex("stan::prob::simplex_free(%1%)", x, "Simplex variable");
1004  int Km1 = x.size() - 1;
1005  Eigen::Matrix<T,Eigen::Dynamic,1> y(Km1);
1006  T stick_len(x(Km1));
1007  for (size_type k = Km1; --k >= 0; ) {
1008  stick_len += x(k);
1009  T z_k(x(k) / stick_len);
1010  y(k) = logit(z_k) + log(Km1 - k);
1011  // log(Km-k) = logit(1.0 / (Km1 + 1 - k));
1012  }
1013  return y;
1014  }
1015 
1016 
1017  // ORDERED
1018 
1028  template <typename T>
1029  Eigen::Matrix<T,Eigen::Dynamic,1>
1030  ordered_constrain(const Eigen::Matrix<T,Eigen::Dynamic,1>& x) {
1031  typedef typename Eigen::Matrix<T,Eigen::Dynamic,1>::size_type size_type;
1032  size_type k = x.size();
1033  Eigen::Matrix<T,Eigen::Dynamic,1> y(k);
1034  if (k == 0)
1035  return y;
1036  y[0] = x[0];
1037  for (size_type i = 1;
1038  i < k;
1039  ++i)
1040  y[i] = y[i-1] + exp(x[i]);
1041  return y;
1042  }
1043 
1056  template <typename T>
1057  inline
1058  Eigen::Matrix<T,Eigen::Dynamic,1>
1059  ordered_constrain(const Eigen::Matrix<T,Eigen::Dynamic,1>& x, T& lp) {
1060  typedef typename Eigen::Matrix<T,Eigen::Dynamic,1>::size_type size_type;
1061  for (size_type i = 1; i < x.size(); ++i)
1062  lp += x(i);
1063  return ordered_constrain(x);
1064  }
1065 
1066 
1067 
1081  template <typename T>
1082  Eigen::Matrix<T,Eigen::Dynamic,1>
1083  ordered_free(const Eigen::Matrix<T,Eigen::Dynamic,1>& y) {
1084  stan::math::check_ordered("stan::prob::ordered_free(%1%)",
1085  y, "Ordered variable");
1086  typedef typename Eigen::Matrix<T,Eigen::Dynamic,1>::size_type size_type;
1087  size_type k = y.size();
1088  Eigen::Matrix<T,Eigen::Dynamic,1> x(k);
1089  if (k == 0)
1090  return x;
1091  x[0] = y[0];
1092  for (size_type i = 1; i < k; ++i)
1093  x[i] = log(y[i] - y[i-1]);
1094  return x;
1095  }
1096 
1097 
1098  // POSITIVE ORDERED
1099 
1109  template <typename T>
1110  Eigen::Matrix<T,Eigen::Dynamic,1>
1111  positive_ordered_constrain(const Eigen::Matrix<T,Eigen::Dynamic,1>& x) {
1112  typedef typename Eigen::Matrix<T,Eigen::Dynamic,1>::size_type size_type;
1113  size_type k = x.size();
1114  Eigen::Matrix<T,Eigen::Dynamic,1> y(k);
1115  if (k == 0)
1116  return y;
1117  y[0] = exp(x[0]);
1118  for (size_type i = 1;
1119  i < k;
1120  ++i)
1121  y[i] = y[i-1] + exp(x[i]);
1122  return y;
1123  }
1124 
1137  template <typename T>
1138  inline
1139  Eigen::Matrix<T,Eigen::Dynamic,1>
1140  positive_ordered_constrain(const Eigen::Matrix<T,Eigen::Dynamic,1>& x, T& lp) {
1141  typedef typename Eigen::Matrix<T,Eigen::Dynamic,1>::size_type size_type;
1142  for (size_type i = 0; i < x.size(); ++i)
1143  lp += x(i);
1144  return positive_ordered_constrain(x);
1145  }
1146 
1147 
1148 
1162  template <typename T>
1163  Eigen::Matrix<T,Eigen::Dynamic,1>
1164  positive_ordered_free(const Eigen::Matrix<T,Eigen::Dynamic,1>& y) {
1165  stan::math::check_positive_ordered("stan::prob::positive_ordered_free(%1%)",
1166  y, "Positive ordered variable");
1167  typedef typename Eigen::Matrix<T,Eigen::Dynamic,1>::size_type size_type;
1168  size_type k = y.size();
1169  Eigen::Matrix<T,Eigen::Dynamic,1> x(k);
1170  if (k == 0)
1171  return x;
1172  x[0] = log(y[0]);
1173  for (size_type i = 1; i < k; ++i)
1174  x[i] = log(y[i] - y[i-1]);
1175  return x;
1176  }
1177 
1178 
1179  // CORRELATION MATRIX
1204  template <typename T>
1205  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>
1206  corr_matrix_constrain(const Eigen::Matrix<T,Eigen::Dynamic,1>& x,
1207  typename Eigen::Matrix<T,Eigen::Dynamic,1>::size_type k) {
1208  typedef typename Eigen::Matrix<T,Eigen::Dynamic,1>::size_type size_type;
1209  size_type k_choose_2 = (k * (k - 1)) / 2;
1210  if (k_choose_2 != x.size())
1211  throw std::invalid_argument ("x is not a valid correlation matrix");
1212  Eigen::Array<T,Eigen::Dynamic,1> cpcs(k_choose_2);
1213  for (size_type i = 0; i < k_choose_2; ++i)
1214  cpcs[i] = corr_constrain(x[i]);
1215  return read_corr_matrix(cpcs,k);
1216  }
1217 
1237  template <typename T>
1238  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>
1239  corr_matrix_constrain(const Eigen::Matrix<T,Eigen::Dynamic,1>& x,
1240  typename Eigen::Matrix<T,Eigen::Dynamic,1>::size_type k,
1241  T& lp) {
1242  typedef typename Eigen::Matrix<T,Eigen::Dynamic,1>::size_type size_type;
1243  size_type k_choose_2 = (k * (k - 1)) / 2;
1244  if (k_choose_2 != x.size())
1245  throw std::invalid_argument ("x is not a valid correlation matrix");
1246 
1247  Eigen::Array<T,Eigen::Dynamic,1> cpcs(k_choose_2);
1248  for (size_type i = 0; i < k_choose_2; ++i)
1249  cpcs[i] = corr_constrain(x[i],lp);
1250  return read_corr_matrix(cpcs,k,lp);
1251  }
1252 
1273  template <typename T>
1274  Eigen::Matrix<T,Eigen::Dynamic,1>
1275  corr_matrix_free(const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& y) {
1276  typedef typename
1277  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>::size_type size_type;
1278  size_type k = y.rows();
1279  if (y.cols() != k)
1280  throw std::domain_error("y is not a square matrix");
1281  if (k == 0)
1282  throw std::domain_error("y has no elements");
1283  size_type k_choose_2 = (k * (k-1)) / 2;
1284  Eigen::Array<T,Eigen::Dynamic,1> x(k_choose_2);
1285  Eigen::Array<T,Eigen::Dynamic,1> sds(k);
1286  bool successful = factor_cov_matrix(x,sds,y);
1287  if (!successful)
1288  throw std::runtime_error("factor_cov_matrix failed on y");
1289  for (size_type i = 0; i < k; ++i) {
1290  // sds on log scale unconstrained
1291  if (fabs(sds[i] - 0.0) >= CONSTRAINT_TOLERANCE) {
1292  std::stringstream s;
1293  s << "all standard deviations must be zero."
1294  << " found log(sd[" << i << "])=" << sds[i] << std::endl;
1295  BOOST_THROW_EXCEPTION(std::runtime_error(s.str()));
1296  }
1297  }
1298  return x.matrix();
1299  }
1300 
1301 
1302  // COVARIANCE MATRIX
1303 
1316  template <typename T>
1317  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>
1318  cov_matrix_constrain(const Eigen::Matrix<T,Eigen::Dynamic,1>& x,
1319  typename Eigen::Matrix<T,Eigen::Dynamic,1>::size_type K) {
1320  using std::exp;
1321  typedef typename Eigen::Matrix<T,Eigen::Dynamic,1>::size_type size_type;
1322  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> L(K,K);
1323  if (x.size() != (K * (K + 1)) / 2)
1324  throw std::domain_error("x.size() != K + (K choose 2)");
1325  int i = 0;
1326  for (size_type m = 0; m < K; ++m) {
1327  for (int n = 0; n < m; ++n)
1328  L(m,n) = x(i++);
1329  L(m,m) = exp(x(i++));
1330  for (size_type n = m + 1; n < K; ++n)
1331  L(m,n) = 0.0;
1332  }
1333  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> M = L * L.transpose();
1334  return L * L.transpose();
1335  }
1336 
1337 
1350  template <typename T>
1351  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>
1352  cov_matrix_constrain(const Eigen::Matrix<T,Eigen::Dynamic,1>& x,
1353  typename Eigen::Matrix<T,Eigen::Dynamic,1>::size_type K,
1354  T& lp) {
1355  using std::exp;
1356  typedef typename Eigen::Matrix<T,Eigen::Dynamic,1>::size_type size_type;
1357  if (x.size() != (K * (K + 1)) / 2)
1358  throw std::domain_error("x.size() != K + (K choose 2)");
1359  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> L(K,K);
1360  int i = 0;
1361  for (size_type m = 0; m < K; ++m) {
1362  for (size_type n = 0; n < m; ++n)
1363  L(m,n) = x(i++);
1364  L(m,m) = exp(x(i++));
1365  for (size_type n = m + 1; n < K; ++n)
1366  L(m,n) = 0.0;
1367  }
1368  // Jacobian for complete transform, including exp() above
1369  lp += (K * stan::math::LOG_2); // needless constant; want propto
1370  for (int k = 0; k < K; ++k)
1371  lp += (K - k + 1) * log(L(k,k)); // only +1 because index from 0
1372  return L * L.transpose();
1373  // return tri_multiply_transpose(L);
1374  }
1375 
1398  template <typename T>
1399  Eigen::Matrix<T,Eigen::Dynamic,1>
1400  cov_matrix_free(const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& y) {
1401  using std::log;
1402  int K = y.rows();
1403  if (y.cols() != K)
1404  throw std::domain_error("y is not a square matrix");
1405  if (K == 0)
1406  throw std::domain_error("y has no elements");
1407  for (int k = 0; k < K; ++k)
1408  if (!(y(k,k) > 0.0))
1409  throw std::domain_error("y has non-positive diagonal");
1410  Eigen::Matrix<T,Eigen::Dynamic,1> x((K * (K + 1)) / 2);
1411  // FIXME: see Eigen LDLT for rank-revealing version -- use that
1412  // even if less efficient?
1413  Eigen::LLT<Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> >
1414  llt(y.rows());
1415  llt.compute(y);
1416  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> L = llt.matrixL();
1417  int i = 0;
1418  for (int m = 0; m < K; ++m) {
1419  for (int n = 0; n < m; ++n)
1420  x(i++) = L(m,n);
1421  x(i++) = log(L(m,m));
1422  }
1423  return x;
1424  }
1425 
1445  template <typename T>
1446  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>
1447  cov_matrix_constrain_lkj(const Eigen::Matrix<T,Eigen::Dynamic,1>& x,
1448  size_t k) {
1449  size_t k_choose_2 = (k * (k - 1)) / 2;
1450  Eigen::Array<T,Eigen::Dynamic,1> cpcs(k_choose_2);
1451  int pos = 0;
1452  for (size_t i = 0; i < k_choose_2; ++i)
1453  cpcs[i] = corr_constrain(x[pos++]);
1454  Eigen::Array<T,Eigen::Dynamic,1> sds(k);
1455  for (size_t i = 0; i < k; ++i)
1456  sds[i] = positive_constrain(x[pos++]);
1457  return read_cov_matrix(cpcs, sds);
1458  }
1459 
1484  template <typename T>
1485  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>
1486  cov_matrix_constrain_lkj(const Eigen::Matrix<T,Eigen::Dynamic,1>& x,
1487  size_t k,
1488  T& lp) {
1489  size_t k_choose_2 = (k * (k - 1)) / 2;
1490  Eigen::Array<T,Eigen::Dynamic,1> cpcs(k_choose_2);
1491  int pos = 0;
1492  for (size_t i = 0; i < k_choose_2; ++i)
1493  cpcs[i] = corr_constrain(x[pos++], lp);
1494  Eigen::Array<T,Eigen::Dynamic,1> sds(k);
1495  for (size_t i = 0; i < k; ++i)
1496  sds[i] = positive_constrain(x[pos++],lp);
1497  return read_cov_matrix(cpcs, sds, lp);
1498  }
1499 
1518  template <typename T>
1519  Eigen::Matrix<T,Eigen::Dynamic,1>
1520  cov_matrix_free_lkj(
1521  const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& y) {
1522  typedef typename
1523  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>::size_type size_type;
1524  size_type k = y.rows();
1525  if (y.cols() != k)
1526  throw std::domain_error("y is not a square matrix");
1527  if (k == 0)
1528  throw std::domain_error("y has no elements");
1529  size_type k_choose_2 = (k * (k-1)) / 2;
1530  Eigen::Array<T,Eigen::Dynamic,1> cpcs(k_choose_2);
1531  Eigen::Array<T,Eigen::Dynamic,1> sds(k);
1532  bool successful = factor_cov_matrix(cpcs,sds,y);
1533  if (!successful)
1534  throw std::runtime_error ("factor_cov_matrix failed on y");
1535  Eigen::Matrix<T,Eigen::Dynamic,1> x(k_choose_2 + k);
1536  size_type pos = 0;
1537  for (size_type i = 0; i < k_choose_2; ++i)
1538  x[pos++] = cpcs[i];
1539  for (size_type i = 0; i < k; ++i)
1540  x[pos++] = sds[i];
1541  return x;
1542  }
1543 
1544  }
1545 
1546 }
1547 
1548 #endif
     [ Stan Home Page ] © 2011–2012, Stan Development Team.