TMB Documentation  v1.9.1
density.hpp
Go to the documentation of this file.
1 // Copyright (C) 2013-2015 Kasper Kristensen
2 // License: GPL-2
3 
14 namespace density {
15 using namespace tmbutils;
16 
17 #define TYPEDEFS(scalartype_) \
18 public: \
19 typedef scalartype_ scalartype; \
20 typedef vector<scalartype> vectortype; \
21 typedef matrix<scalartype> matrixtype; \
22 typedef array<scalartype> arraytype
23 
24 #define VARIANCE_NOT_YET_IMPLEMENTED \
25 private: \
26 vectortype variance(){return vectortype();} \
27 public:
28 
54 #define SIMULATE_NOT_YET_IMPLEMENTED \
55 private: \
56 vectortype sqrt_cov_scale(vectortype x){} \
57 void simulate(vectortype &u){} \
58 vectortype simulate(){} \
59 public:
60 
61 /* Add this macro to all classes with *dynamic* dimension (e.g. AR1) */
62 #define SIMULATE_IMPLEMENTED_UNKNOWN_SIZE \
63 void simulate(vectortype &x) { \
64  rnorm_fill(x); \
65  x = sqrt_cov_scale(x); \
66  x = zero_derivatives(x); \
67 }
68 /* Add this macro to all classes with *fixed* dimension (e.g. MVNORM) */
69 #define SIMULATE_IMPLEMENTED_KNOWN_SIZE(SIZE) \
70 SIMULATE_IMPLEMENTED_UNKNOWN_SIZE \
71 vectortype simulate() { \
72  vectortype x(SIZE); \
73  simulate(x); \
74  return x; \
75 }
76 
77 /* Utility function: The simulators should not track derivatives when
78  running with AD types. A workaround is to zero out the derivatives
79  after simulation (FIXME: there is a minor efficiency loss by
80  tracking the derivatives in the first place...). */
81 template<class arraytype>
82 arraytype zero_derivatives(arraytype x) {
83  for(int i=0; i<x.size(); i++)
84  x(i) = asDouble(x(i));
85  return x;
86 }
87 
88 /* Utility function */
89 template<class arraytype>
90 void rnorm_fill(arraytype &x) {
91  for(int i=0; i<x.size(); i++)
92  x(i) = Rf_rnorm(0., 1.);
93 }
94 
111 template <class scalartype_>
112 class MVNORM_t{
113  TYPEDEFS(scalartype_);
114  matrixtype Q; /* Inverse covariance matrix */
115  scalartype logdetQ; /* log-determinant of Q */
116  matrixtype Sigma; /* Keep for convenience - not used */
117  matrixtype L_Sigma; /* Used by simulate() */
118 public:
119  MVNORM_t(){}
120  MVNORM_t(matrixtype Sigma_, bool use_atomic=true){
121  setSigma(Sigma_, use_atomic);
122  }
123 
134  matrixtype cov(){return Sigma;}
135 
136  /* initializer via covariance matrix */
137  void setSigma(matrixtype Sigma_, bool use_atomic=true){
138  Sigma = Sigma_;
139  scalartype logdetS;
140  if(use_atomic){
141  Q = atomic::matinvpd(Sigma,logdetS);
142  } else {
143  matrixtype I(Sigma.rows(),Sigma.cols());
144  I.setIdentity();
145  Eigen::LDLT<Eigen::Matrix<scalartype,Dynamic,Dynamic> > ldlt(Sigma);
146  Q = ldlt.solve(I);
147  vectortype D = ldlt.vectorD();
148  logdetS = D.log().sum();
149  }
150  logdetQ = -logdetS;
151  }
152  scalartype Quadform(vectortype x){
153  return (x*(vectortype(Q*x))).sum();
154  }
156  scalartype operator()(vectortype x){
157  return -scalartype(.5)*logdetQ + scalartype(.5)*Quadform(x) + x.size()*scalartype(log(sqrt(2.0*M_PI)));
158  }
162  scalartype operator()(vectortype x, vectortype keep){
163  matrix<scalartype> S = Sigma;
164  vector<scalartype> not_keep = scalartype(1.0) - keep;
165  for(int i = 0; i < S.rows(); i++){
166  for(int j = 0; j < S.cols(); j++){
167  S(i,j) = S(i,j) * keep(i) * keep(j);
168  }
169  S(i,i) += not_keep(i) * scalartype(1.0 / (2.0 * M_PI));
170  }
171  return MVNORM_t<scalartype>(S)(x * keep);
172  }
173  arraytype jacobian(arraytype x){
174  arraytype y(x.dim);
175  matrixtype m(x.size()/x.cols(),x.cols());
176  for(int i=0;i<x.size();i++)m(i)=x[i];
177  matrixtype mQ=m*Q;
178  for(int i=0;i<x.size();i++)y[i]=mQ(i);
179  return y;
180  }
181  int ndim(){return 1;}
182  VARIANCE_NOT_YET_IMPLEMENTED
183  vectortype sqrt_cov_scale(vectortype u) {
184  if(L_Sigma.rows() == 0) {
185  Eigen::LLT<Eigen::Matrix<scalartype,Dynamic,Dynamic> > llt(Sigma);
186  L_Sigma = llt.matrixL();
187  }
188  vectortype ans = L_Sigma * u;
189  return ans;
190  }
191  SIMULATE_IMPLEMENTED_KNOWN_SIZE(Sigma.rows());
192 };
193 
207 template <class scalartype>
208 MVNORM_t<scalartype> MVNORM(matrix<scalartype> Sigma, bool use_atomic = true){
209  return MVNORM_t<scalartype>(Sigma, use_atomic);
210 }
211 
259 template <class scalartype_>
260 class UNSTRUCTURED_CORR_t : public MVNORM_t<scalartype_>{
261  TYPEDEFS(scalartype_);
263  UNSTRUCTURED_CORR_t(vectortype x){
264  // (n*n-n)/2=nx ==> n*n-n-2*nx=0 ==> n=(1+sqrt(1+8*nx))/2
265  int nx=x.size();
266  int n=int((1.0+sqrt(1.0+8*nx))/2.0);
267  if((n*n-n)/2!=nx)std::cout << "vector does not specify an UNSTRUCTERED_CORR\n";
268  matrixtype L(n,n);
269  L.setIdentity();
270  int i,j,k=0;
271  for(i=0;i<L.rows();i++){
272  for(j=0;j<L.cols();j++){
273  if(i>j){L(i,j)=x[k];k++;}
274  }
275  }
276  matrixtype llt=L*L.transpose();
277  matrixtype Sigma=llt;
278  for(i=0;i<Sigma.rows();i++){
279  for(j=0;j<Sigma.cols();j++){
280  Sigma(i,j)/=sqrt(llt(i,i)*llt(j,j));
281  }
282  }
283  this->setSigma(Sigma); /* Call MVNORM_t initializer */
284  }
285 };
288 template <class scalartype>
291 }
292 
301 template<class scalartype_>
302 class N01{
303  TYPEDEFS(scalartype_);
304 public:
306  scalartype operator()(scalartype x){
307  return x*x*.5 + log(sqrt(2.0*M_PI));
308  }
309  scalartype operator()(vectortype x){
310  return (x*x*scalartype(.5) + scalartype(log(sqrt(2.0*M_PI))) ).sum() ;
311  }
312  scalartype operator()(arraytype x){
313  return (x*x*scalartype(.5) + scalartype(log(sqrt(2.0*M_PI))) ).sum() ;
314  }
315  arraytype jacobian(arraytype x){return x;}
316  int ndim(){return 1;}
317  VARIANCE_NOT_YET_IMPLEMENTED
318  // FIXME: 1D will not suffice for e.g. SEPARABLE(N01, OTHER);
319  vectortype simulate() {
320  vectortype x(1);
321  x[0] = rnorm(0.0, 1.0);
322  return x;
323  }
324  // Inplace simulate
325  void simulate(vectortype &x) {
326  rnorm_fill(x);
327  }
328  vectortype sqrt_cov_scale(vectortype u) {
329  return u;
330  }
331 
332 };
333 
374 template <class distribution>
375 class AR1_t{
376  TYPEDEFS(typename distribution::scalartype);
377 private:
378  scalartype phi;
379  distribution MARGINAL;
380 public:
381  AR1_t(){/*phi=phi_;MARGINAL=f_;*/}
382  AR1_t(scalartype phi_, distribution f_) : phi(phi_), MARGINAL(f_) {}
384  scalartype operator()(vectortype x){
385  scalartype value;
386  value=scalartype(0);
387  int n=x.rows();
388  int m=x.size()/n;
389  scalartype sigma=sqrt(scalartype(1)-phi*phi); /* Steady-state standard deviation */
390  value+=MARGINAL(x(0)); /* E.g. x0 ~ N(0,1) */
391  for(int i=1;i<n;i++)value+=MARGINAL((x(i)-x(i-1)*phi)/sigma);/* x(i)-phi*x(i-1) ~ N(0,sigma^2) */
392  value+=scalartype((n-1)*m)*log(sigma);
393  return value;
394  }
395 
396  /* Copied vector version - apply over outermost dimension */
397  scalartype operator()(arraytype x){
398  scalartype value;
399  value=scalartype(0);
400  int n=x.cols();
401  int m=x.size()/n;
402  scalartype sigma=sqrt(scalartype(1)-phi*phi); /* Steady-state standard deviation */
403  value+=MARGINAL(x.col(0)); /* E.g. x0 ~ N(0,1) */
404  for(int i=1;i<n;i++)value+=MARGINAL((x.col(i)-x.col(i-1)*phi)/sigma);/* x(i)-phi*x(i-1) ~ N(0,sigma^2) */
405  value+=scalartype((n-1)*m)*log(sigma);
406  return value;
407  }
408 
409 
410  arraytype jacobian(arraytype x){
411  scalartype sigma=sqrt(scalartype(1)-phi*phi); /* Steady-state standard deviation */
412  int n=x.cols();
413  arraytype y(x.dim);
414  y.setZero();
415  y.col(0) = y.col(0) + MARGINAL.jacobian(x.col(0));
416  for(int i=1;i<n;i++){
417  //MARGINAL((x(i)-x(i-1)*phi)/sigma);
418  y.col(i-1) = y.col(i-1) - (phi/sigma) * MARGINAL.jacobian((x.col(i)-x.col(i-1)*phi)/sigma);
419  y.col(i) = y.col(i) + MARGINAL.jacobian((x.col(i)-x.col(i-1)*phi)/sigma)/sigma;
420  }
421  return y;
422  }
423  int ndim(){return 1;}
424  VARIANCE_NOT_YET_IMPLEMENTED
425 
426  /* Simulation */
427  arraytype sqrt_cov_scale(arraytype u) {
428  if(u.dim.size() == 1) {
429  // FIXME: This entire class should be reconsidered with this
430  // special case in mind !
431  u.dim.resize(2); u.dim << 1, u.size();
432  }
433  int n = u.cols();
434  arraytype x(u.dim);
435  scalartype sigma = sqrt(1.0 - phi * phi); /* Steady-state standard deviation */
436  x.col(0) = MARGINAL.sqrt_cov_scale(u.col(0));
437  for(int i=1; i<n; i++)
438  x.col(i) = phi * x.col(i-1) + sigma * MARGINAL.sqrt_cov_scale(u.col(i));
439  return x;
440  }
441  vectortype sqrt_cov_scale(vectortype u) {
442  vector<int> dim(2);
443  dim << 1, u.size();
444  arraytype u_array(u, dim);
445  return sqrt_cov_scale(u_array);
446  }
448  void simulate(vectortype &x) {
449  rnorm_fill(x);
450  x = sqrt_cov_scale(x);
451  x = zero_derivatives(x);
452  }
453  void simulate(arraytype &x) {
454  rnorm_fill(x);
455  x = sqrt_cov_scale(x);
456  x = zero_derivatives(x);
457  }
458 
459 };
460 template <class scalartype, class distribution>
461 AR1_t<distribution> AR1(scalartype phi_, distribution f_){
462  return AR1_t<distribution>(phi_, f_);
463 }
464 template <class scalartype>
465 AR1_t<N01<scalartype> > AR1(scalartype phi_){
466  return AR1_t<N01<scalartype> >(phi_, N01<scalartype>());
467 }
468 
469 
494 template <class scalartype_>
495 class ARk_t{
496  TYPEDEFS(scalartype_);
497  //private:
498  int k;
499  vectortype phi; /* [phi1,...,phik] */
500  vectortype gamma; /* [gamma(1),...,gamma(k)] (note gamma(0) is 1) */
501  /* Initial distribution matrices. */
502  matrixtype V0; /* kxk variance */
503  matrixtype Q0; /* kxk precision */
504  matrixtype L0; /* kxk Cholesky Q0 = L0*L0' */
505  /* gamma is found through (I-M)*gamma=phi ... */
506  matrixtype M; /* kxk */
507  matrixtype I; /* kxk */
508  scalartype sigma;/* increment standard deviation */
509  scalartype logdetQ0;
510 public:
511  ARk_t(){/*phi=phi_;MARGINAL=f_;*/}
512  ARk_t(vectortype phi_){
513  phi=phi_;
514  k=phi.size();
515  V0.resize(k,k);Q0.resize(k,k);
516  M.resize(k,k);I.resize(k,k);
517  /* build M-matrix */
518  M.setZero();
519  int d;
520  for(int i=0;i<k;i++){
521  for(int j=0;j<k;j++){
522  d=abs(i-j);
523  if(d!=0){
524  M(i,d-1)+=phi[j];
525  }
526  }
527  }
528  I.setIdentity();
529  gamma=((I-M).inverse())*matrixtype(phi);
530  /* Increment sd */
531  sigma=sqrt(scalartype(1)-(phi*gamma).sum());
532  /* build V0 matrix */
533  for(int i=0;i<k;i++){
534  for(int j=0;j<k;j++){
535  d=abs(i-j);
536  if(d==0)
537  V0(i,j)=scalartype(1);
538  else
539  V0(i,j)=gamma(d-1);
540  }
541  }
542  /* build Q0 matrix */
543  Q0=V0.inverse();
544  /* log determinant */
545  L0=Q0.llt().matrixL(); /* L0 L0' = Q0 */
546  logdetQ0=scalartype(0);
547  for(int i=0;i<k;i++)logdetQ0+=scalartype(2)*log(L0(i,i));
548  }
553  vectortype cov(int n){
554  vectortype rho(n);
555  for(int i=0;i<n;i++){
556  if(i==0){rho(0)=scalartype(1);}
557  else if(i<=k){rho(i)=gamma(i-1);}
558  else { /* youle walker */
559  scalartype tmp=0;
560  for(int j=0;j<k;j++)tmp+=phi[j]*rho[i-1-j];
561  rho(i)=tmp;
562  }
563  }
564  return rho;
565  }
566 
568  scalartype operator()(vectortype x){
569  scalartype value=0;
570  /* Initial distribution. For i = k,...,1 the recursions for
571  solving L' y = u (1-based index notation) are:
572 
573  y(i) = L(i,i)^-1 * ( u(i) - L(i+1,i) * y(i+1) - ... - L(k,i) * y(k) )
574 
575  where u(i) ~ N(0,1).
576  */
577  scalartype mu, sd;
578  int col;
579  for(int i=0; (i<k) & (i<x.size()); i++){
580  mu = scalartype(0);
581  col = k-1-i; /* reversed index */
582  for(int j=col+1; j<k; j++) mu -= L0(j,col) * x(k-1-j);
583  mu /= L0(col, col);
584  sd = scalartype(1) / L0(col, col);
585  value -= dnorm(x[i], mu, sd, true);
586  }
587  scalartype tmp;
588  for(int i=k;i<x.size();i++){
589  tmp=scalartype(0);
590  for(int j=0;j<k;j++){
591  tmp+=phi[j]*x[i-1-j];
592  }
593  value-=dnorm(x[i],tmp,sigma,1);
594  }
595  return value;
596  }
597  arraytype jacobian(arraytype x){
598  arraytype y(x.dim);
599  y.setZero();
600  for(int i=0;i<k;i++)
601  for(int j=0;j<k;j++)
602  y.col(i)=y.col(i)+Q0(i,j)*x.col(j);
603  vectortype v(k+1);
604  v(0)=scalartype(1);
605  for(int i=1;i<=k;i++)v[i]=-phi[i-1];
606  v=v/sigma;
607  for(int i=k;i<x.cols();i++){
608  for(int j1=0;j1<=k;j1++){
609  for(int j2=0;j2<=k;j2++){
610  y.col(i-j1)=y.col(i-j1)+v[j1]*v[j2]*x.col(i-j2);
611  }
612  }
613  }
614  return y;
615  }
616  int ndim(){return 1;}
617  VARIANCE_NOT_YET_IMPLEMENTED
618  /* Simulate - copy paste from evaluation operator */
619  vectortype sqrt_cov_scale(vectortype u){
620  vectortype x(u.size());
621  scalartype mu, sd;
622  int col;
623  for (int i = 0; (i < k) & (i < x.size()); i++) {
624  mu = scalartype(0);
625  col = k - 1 - i; /* reversed index */
626  for (int j = col + 1; j < k; j++)
627  mu -= L0(j, col) * x(k - 1 - j);
628  mu /= L0(col, col);
629  sd = scalartype(1) / L0(col, col);
630  x(i) = sd * u(i) + mu;
631  }
632  scalartype tmp;
633  for (int i = k; i < x.size(); i++) {
634  tmp = scalartype(0);
635  for (int j = 0; j < k; j++) {
636  tmp += phi[j] * x[i - 1 - j];
637  }
638  x(i) = sigma * u(i) + tmp;
639  }
640  return x;
641  }
642  SIMULATE_IMPLEMENTED_UNKNOWN_SIZE
643 };
644 
645 template <class scalartype>
647  return ARk_t<scalartype>(phi);
648 }
649 
650 
668 template <class scalartype_>
669 class contAR2_t{
670  TYPEDEFS(scalartype_);
671 private:
672  typedef Matrix<scalartype,2,2> matrix2x2;
673  typedef Matrix<scalartype,2,1> matrix2x1;
674  typedef Matrix<scalartype,4,4> matrix4x4;
675  typedef Matrix<scalartype,4,1> matrix4x1;
676  scalartype shape,scale,c0,c1;
677  vectortype grid;
678  matrix2x2 A, V0, I;
679  matrix4x4 B, iB; /* B=A %x% I + I %x% A */
681  matrix4x1 vecSigma,iBvecSigma;
682  vector<MVNORM_t<scalartype> > neglogdmvnorm; /* Cache the 2-dim increments */
683  vector<matrix2x2 > expAdt; /* Cache matrix exponential for grid increments */
684 public:
685  contAR2_t(){};
686  contAR2_t(vectortype grid_, scalartype shape_, scalartype scale_=1){
687  shape=shape_;scale=scale_;grid=grid_;
688  c0=scalartype(-1);c1=scalartype(2)*shape_;
689  c0=c0/(scale*scale); c1=c1/scale;
690  A << scalartype(0), scalartype(1), c0, c1;
691  V0 << 1,0,0,-c0;
692  I.setIdentity();
693  B=kronecker(I,A)+kronecker(A,I);
694  iB=B.inverse();
695  expA=matexp<scalartype,2>(A);
696  vecSigma << 0,0,0,scalartype(-2)*c1*V0(1,1);
697  iBvecSigma=iB*vecSigma;
698  /* cache increment distribution N(0,V(dt)) - one for each grid point */
699  neglogdmvnorm.resize(grid.size());
700  neglogdmvnorm[0]=MVNORM_t<scalartype>(V0);
701  for(int i=1;i<grid.size();i++)neglogdmvnorm[i]=MVNORM_t<scalartype>(V(grid(i)-grid(i-1)));
702  /* cache matrix exponential */
703  expAdt.resize(grid.size());
704  expAdt[0]=expA(scalartype(0));
705  for(int i=1;i<grid.size();i++)expAdt[i]=expA(grid(i)-grid(i-1));
706  }
707  /* Simple formula for matrix exponential exp(B*t) */
708  matrix4x4 expB(scalartype t){
709  return kronecker(expA(t),expA(t));
710  }
711  /* Variance as fct. of time when started deterministic */
712  matrix2x2 V(scalartype t){
713  matrix4x1 tmp;
714  tmp=expB(t)*iBvecSigma-iBvecSigma;
715  matrix2x2 ans;
716  for(int i=0;i<4;i++)ans(i)=tmp(i);
717  return ans;
718  }
721  scalartype operator()(vectortype x,vectortype dx){
722  matrix2x1 y, y0;
723  scalartype ans;
724  y0 << x(0), dx(0);
725  ans = neglogdmvnorm[0](y0);
726  for(int i=1;i<grid.size();i++){
727  y0 << x(i-1), dx(i-1);
728  y << x(i), dx(i);
729  ans += neglogdmvnorm[i](y-expAdt[i]*y0);
730  }
731  return ans;
732  }
733  /* Experiment: Implementing matrix-vector multiply - Q*x.
734  To be used when creating separable extensions...
735  think best to assume that input array has dimension (2,ntime,...)
736  so that we can easily extract x(t) and multiply Q*x(t) with a 2x2 matrix */
737  scalartype operator()(vectortype x){ /* x.dim=[2,n] */
738  vector<int> dim(2);
739  dim << 2 , x.size()/2 ;
740  array<scalartype> y(x,dim);
741  y=y.transpose();
742  return this->operator()(y.col(0),y.col(1));
743  }
744  arraytype matmult(matrix2x2 Q,arraytype x){
745  arraytype y(x.dim);
746  y.col(0) = Q(0,0)*x.col(0)+Q(0,1)*x.col(1); /* TODO: can we subassign like this in array class? Hack: we use "y.row" for that */
747  y.col(1) = Q(1,0)*x.col(0)+Q(1,1)*x.col(1);
748  return y;
749  }
750  arraytype jacobian(arraytype x){
751  arraytype y(x.dim);
752  y.setZero();
753  arraytype tmp(y(0).dim);
754  y.col(0) = neglogdmvnorm[0].jacobian(x.col(0)); /* Time zero contrib */
755  for(int i=1;i<grid.size();i++){
756  /* When taking derivative of .5*(x(i)-G*x(i-1))'*Q*(x(i)-G*x(i-1)) [where G=expAdt]
757  we get contributions like:
758  x(i-1): -G'*Q*(x(i)-G*x(i-1))
759  x(i): Q*(x(i)-G*x(i-1))
760  */
761  tmp=neglogdmvnorm[i].jacobian( x.col(i) - matmult(expAdt[i],x.col(i-1)) );
762  y.col(i)=y.col(i)+tmp;
763  y.col(i-1)=y.col(i-1)-matmult(expAdt[i].transpose(), tmp );
764  }
765  return y;
766  }
767  int ndim(){return 2;} /* Number of dimensions this structure occupies in total array */
768  VARIANCE_NOT_YET_IMPLEMENTED
770 };
771 template <class scalartype, class vectortype>
772 contAR2_t<scalartype> contAR2(vectortype grid_, scalartype shape_, scalartype scale_=1){
773  return contAR2_t<scalartype>(grid_, shape_, scale_);
774 }
775 template <class scalartype>
776 contAR2_t<scalartype> contAR2(scalartype shape_, scalartype scale_=1){
777  return contAR2_t<scalartype>(shape_, scale_);
778 }
779 
827 template <class scalartype_>
828 class GMRF_t {
829  TYPEDEFS(scalartype_);
830 private:
831  Eigen::SparseMatrix<scalartype> Q;
832  scalartype logdetQ;
833  int sqdist(vectortype x, vectortype x_) {
834  int ans = 0;
835  int tmp;
836  for(int i=0; i<x.size(); i++){
837  tmp = CppAD::Integer(x[i]) - CppAD::Integer(x_[i]);
838  ans += tmp * tmp;
839  }
840  return ans;
841  }
842 public:
843  GMRF_t(){}
844  GMRF_t(Eigen::SparseMatrix<scalartype> Q_, int order_=1, bool normalize=true){
845  setQ(Q_, order_, normalize);
846  }
847  GMRF_t(arraytype x, vectortype delta, int order_=1, bool normalize=true){
848  int n = x.cols();
849  typedef Eigen::Triplet<scalartype> T;
850  std::vector<T> tripletList;
851  for(int i=0; i<n; i++){
852  for(int j=0; j<n; j++){
853  if (sqdist(x.col(i),x.col(j)) == 1){
854  tripletList.push_back(T(i, j, scalartype(-1)));
855  tripletList.push_back(T(i, i, scalartype( 1)));
856  }
857  }
858  }
859  for(int i=0; i<n; i++) {
860  tripletList.push_back(T(i, i, delta[i]));
861  }
862  Eigen::SparseMatrix<scalartype> Q_(n, n);
863  Q_.setFromTriplets(tripletList.begin(), tripletList.end());
864  setQ(Q_, order_, normalize);
865  }
866  void setQ(Eigen::SparseMatrix<scalartype> Q_, int order=1, bool normalize=true){
867  Q = Q_;
868  if (normalize) {
869 #ifndef TMBAD_FRAMEWORK
870  Eigen::SimplicialLDLT< Eigen::SparseMatrix<scalartype> > ldl(Q);
871  vectortype D = ldl.vectorD();
872  logdetQ = (log(D)).sum();
873 #else
874  logdetQ = newton::log_determinant(Q);
875 #endif
876  } else {
877  logdetQ = 0;
878  }
879  /* Q^order */
880  for(int i=1; i<order; i++){
881  Q = Q * Q_;
882  }
883  logdetQ = scalartype(order) * logdetQ;
884  }
885  /* Quadratic form: x'*Q^order*x */
886  scalartype Quadform(vectortype x){
887  return (x * (Q * x.matrix()).array()).sum();
888  }
889  scalartype operator()(vectortype x){
890  return -scalartype(.5) * logdetQ + scalartype(.5) * Quadform(x) + x.size() * scalartype(log(sqrt(2.0 * M_PI)));
891  }
892  /* jacobian */
893  arraytype jacobian(arraytype x){
894  arraytype y(x.dim);
895  matrixtype m(x.size()/x.cols(),x.cols());
896  for(int i=0; i<x.size(); i++) m(i) = x[i];
897  matrixtype mQ = m * Q;
898  for(int i=0; i<x.size(); i++) y[i] = mQ(i);
899  return y;
900  }
901  int ndim() { return 1; }
902  vectortype variance() {
903  int n = Q.rows();
904  vectortype ans(n);
905  matrixtype C = invertSparseMatrix(Q);
906  for(int i=0; i<n; i++) ans[i] = C(i,i);
907  return ans;
908  }
909  /* Simulation */
910  Eigen::SparseMatrix<scalartype> L;
911  Eigen::PermutationMatrix<Dynamic,Dynamic> Pinv;
912  vectortype sqrt_cov_scale(vectortype u) {
913  if(L.rows() == 0) {
914  Eigen::SimplicialLLT<Eigen::SparseMatrix<scalartype> > solver(Q);
915  L = solver.matrixL();
916  Pinv = solver.permutationPinv();
917  }
918  // L*L^T = P*Q*Pinv => Q^-1 = A*A^T where A:=P^-1*L^T^-1
919  matrixtype x = L.transpose().template triangularView<Eigen::Upper>().solve(u.matrix());
920  x = Pinv * x;
921  return x.vec();
922  }
923  SIMULATE_IMPLEMENTED_KNOWN_SIZE(Q.rows())
924 };
933 template <class scalartype>
934 GMRF_t<scalartype> GMRF(Eigen::SparseMatrix<scalartype> Q, int order, bool normalize=true) {
935  return GMRF_t<scalartype>(Q, order, normalize);
936 }
937 template <class scalartype, class arraytype >
938 GMRF_t<scalartype> GMRF(arraytype x, vector<scalartype> delta, int order=1, bool normalize=true) {
939  return GMRF_t<scalartype>(x, delta, order, normalize);
940 }
941 template <class scalartype, class arraytype >
942 GMRF_t<scalartype> GMRF(arraytype x, scalartype delta, int order=1, bool normalize=true) {
943  vector<scalartype> d(x.cols());
944  for(int i=0; i<d.size(); i++) d[i] = delta;
945  return GMRF_t<scalartype>(x, d, order, normalize);
946 }
947 template <class scalartype>
948 GMRF_t<scalartype> GMRF(Eigen::SparseMatrix<scalartype> Q, bool normalize = true) {
949  return GMRF_t<scalartype>(Q, 1, normalize);
950 }
951 
959 template <class distribution>
960 class SCALE_t{
961  TYPEDEFS(typename distribution::scalartype);
962 private:
963  distribution f;
964  scalartype scale;
965 public:
966  SCALE_t(){}
967  SCALE_t(distribution f_, scalartype scale_){scale=scale_;f=f_;}
969  scalartype operator()(arraytype x){
970  scalartype ans=f(x/scale);
971  ans+=x.size()*log(scale);
972  return ans;
973  }
974  scalartype operator()(vectortype x){
975  scalartype ans=f(x/scale);
976  ans+=x.size()*log(scale);
977  return ans;
978  }
979  arraytype jacobian(arraytype x){
980  return f.jacobian(x/scale)/scale;
981  }
982  int ndim(){return f.ndim();}
983  vectortype variance(){
984  return (scale*scale)*f.variance();
985  }
986  vectortype sqrt_cov_scale(vectortype u){
987  return scale * f.sqrt_cov_scale(u);
988  }
989  SIMULATE_IMPLEMENTED_UNKNOWN_SIZE
990 };
991 template <class scalartype, class distribution>
992 SCALE_t<distribution> SCALE(distribution f_, scalartype scale_){
993  return SCALE_t<distribution>(f_,scale_);
994 }
995 
1018 template <class distribution>
1020  TYPEDEFS(typename distribution::scalartype);
1021 private:
1022  distribution f;
1023  vectortype scale;
1024 public:
1025  VECSCALE_t(){}
1026  VECSCALE_t(distribution f_, vectortype scale_){scale=scale_;f=f_;}
1028  scalartype operator()(arraytype x){
1029  // assert that x.size()==scale.size()
1030  scalartype ans=f(x/scale);
1031  ans+=(log(scale)).sum();
1032  return ans;
1033  }
1034  scalartype operator()(vectortype x){
1035  // assert that x.size()==scale.size()
1036  scalartype ans=f(x/scale);
1037  ans+=(log(scale)).sum();
1038  return ans;
1039  }
1040  arraytype jacobian(arraytype x){
1041  // assert that x.rows()==scale.size()
1042  arraytype y(x);
1043  for(int i=0;i<y.cols();i++)y.col(i)=y.col(i)/scale[i];
1044  y=f.jacobian(y);
1045  for(int i=0;i<y.cols();i++)y.col(i)=y.col(i)/scale[i];
1046  return y;
1047  }
1048  int ndim(){return f.ndim();}
1049  VARIANCE_NOT_YET_IMPLEMENTED
1050  vectortype sqrt_cov_scale(vectortype u){
1051  return scale * f.sqrt_cov_scale(u);
1052  }
1053  SIMULATE_IMPLEMENTED_UNKNOWN_SIZE
1054 };
1056 template <class vectortype, class distribution>
1057 VECSCALE_t<distribution> VECSCALE(distribution f_, vectortype scale_){
1058  return VECSCALE_t<distribution>(f_,scale_);
1059 }
1060 
1061 
1104 //template <class scalartype, class vectortype, class arraytype, class distribution1, class distribution2>
1105 template <class distribution1, class distribution2>
1107  TYPEDEFS(typename distribution1::scalartype);
1108 private:
1109  distribution1 f;
1110  distribution2 g;
1111 public:
1112  SEPARABLE_t(){}
1113  SEPARABLE_t(distribution1 f_, distribution2 g_){f=f_;g=g_;}
1114  /*
1115  Example: x.dim=[n1,n2,n3].
1116  Apply f on outer dimension (n3) and rotate:
1117  [n3,n1,n2]
1118  Apply g on new outer dimension (n2) and rotate back:
1119  [n1,n2,n3]
1120  */
1121  arraytype jacobian(arraytype x){
1122  int n=f.ndim();
1123  x=f.jacobian(x);
1124  x=x.rotate(n);
1125  x=g.jacobian(x);
1126  x=x.rotate(-n);
1127  return x;
1128  }
1129  /* Create zero vector corresponding to the last n dimensions of dimension-vector d */
1130  arraytype zeroVector(vector<int> d, int n){
1131  int m=1;
1132  vector<int> revd=d.reverse();
1133  vector<int> revnewdim(n);
1134  for(int i=0;i<n;i++){m=m*revd[i];revnewdim[i]=revd[i];}
1135  vectortype x(m);
1136  x.setZero();
1137  return arraytype(x,revnewdim.reverse());
1138  }
1139  scalartype operator()(arraytype x){
1140  if(this->ndim() != x.dim.size())std::cout << "Wrong dimension in SEPARABLE_t\n";
1141  /* Calculate quadform */
1142  arraytype y(x.dim);
1143  y=jacobian(x);
1144  y=x*y; /* pointwise */
1145  scalartype q=scalartype(.5)*(y.sum());
1146  /* Add normalizing constant */
1147  int n=f.ndim();
1148  arraytype zf=zeroVector(x.dim,n);
1149  q+=f(zf)*(scalartype(x.size())/scalartype(zf.size()));
1150  x=x.rotate(n);
1151  int m=g.ndim();
1152  arraytype zg=zeroVector(x.dim,m);
1153  q+=g(zg)*(scalartype(x.size())/scalartype(zg.size()));
1154  q-=log(sqrt(2.0*M_PI))*(zf.size()*zg.size());
1155  /* done */
1156  return q;
1157  }
1158  int ndim(){return f.ndim()+g.ndim();}
1159  VARIANCE_NOT_YET_IMPLEMENTED
1160 
1161  /* For parallel accumulation:
1162  ==========================
1163  Copied operator() above and added extra argument "i" to divide the accumulation
1164  in chunks. The evaluation of
1165  operator()(x)
1166  is equivalent to summing up
1167  operator()(x,i)
1168  with i running through the _outer_dimension_ of x.
1169  */
1170  scalartype operator()(arraytype x, int i){
1171  if(this->ndim() != x.dim.size())std::cout << "Wrong dimension in SEPARABLE_t\n";
1172  /* Calculate quadform */
1173  arraytype y(x.dim);
1174  y=jacobian(x);
1175  y=x*y; /* pointwise */
1176  scalartype q=scalartype(.5)*(y.col(i).sum());
1177  /* Add normalizing constant */
1178  if(i==0){
1179  int n=f.ndim();
1180  arraytype zf=zeroVector(x.dim,n);
1181  q+=f(zf)*(scalartype(x.size())/scalartype(zf.size()));
1182  x=x.rotate(n);
1183  int m=g.ndim();
1184  arraytype zg=zeroVector(x.dim,m);
1185  q+=g(zg)*(scalartype(x.size())/scalartype(zg.size()));
1186  q-=log(sqrt(2.0*M_PI))*(zf.size()*zg.size());
1187  }
1188  /* done */
1189  return q;
1190  }
1191 
1192  arraytype sqrt_cov_scale(arraytype u) {
1193  vector<int> u_dim = u.dim;
1194  vector<int> f_dim = u_dim.tail(f.ndim());
1195  vector<int> g_dim = u_dim.head(g.ndim());
1196  int f_size = f_dim.prod();
1197  int g_size = g_dim.prod();
1198  // Collapse f dimension to a single dimension:
1199  vector<int> new_dim(g.ndim() + 1);
1200  new_dim << g_dim, f_size;
1201  u.setdim(new_dim);
1202  for(int i=0; i<f_size; i++) {
1203  u.col(i) = g.sqrt_cov_scale( u.col(i) );
1204  }
1205  u = u.rotate(1); // u.dim = (f_size, g_dim)
1206  // Collapse g dimension to a single dimension:
1207  new_dim.resize(f.ndim() + 1);
1208  new_dim << f_dim, g_size;
1209  u.setdim(new_dim);
1210  for(int i=0; i<g_size; i++) {
1211  u.col(i) = f.sqrt_cov_scale( u.col(i) );
1212  }
1213  u = u.rotate(1);
1214  u.setdim(u_dim);
1215  return u;
1216  }
1217 
1218  void simulate(arraytype &x) {
1219  rnorm_fill(x);
1220  x = sqrt_cov_scale(x);
1221  x = zero_derivatives(x);
1222  }
1223 
1224 };
1225 
1246 template <class distribution1, class distribution2>
1247 SEPARABLE_t<distribution1,distribution2> SEPARABLE(distribution1 f_, distribution2 g_){
1249 }
1250 
1296 template <class distribution>
1297 class PROJ_t{
1298  TYPEDEFS(typename distribution::scalartype);
1299 private:
1300  distribution f;
1301  bool initialized;
1302 public:
1303  vector<int> proj;
1304  vector<int> cproj; /* complementary proj _sorted_ */
1305  int n,nA,nB;
1306  matrixtype Q; /* Full precision */
1307  MVNORM_t<scalartype> dmvnorm; /* mean zero gaussian with covariance Q_BB */
1308  PROJ_t(){}
1309  PROJ_t(distribution f_, vector<int> proj_){
1310  f=f_;
1311  proj=proj_;
1312  initialized=false;
1313  }
1314  void initialize(int n_){
1315  if(!initialized){
1316  n=n_;
1317  nA=proj.size();
1318  nB=n-nA;
1319  cproj.resize(nB);
1320  vector<int> mark(n);
1321  mark.setZero();
1322  for(int i=0;i<nA;i++)mark[proj[i]]=1;
1323  int k=0;
1324  for(int i=0;i<n;i++)if(!mark[i])cproj[k++]=i;
1325  // Full precision
1326  //matrixtype I(n,n);
1327  //I.setIdentity();
1328  //vectortype v(I);
1329 
1330  k=0;
1331  vectortype v(n*n);
1332  for(int i=0;i<n;i++)
1333  for(int j=0;j<n;j++)
1334  v(k++)=scalartype(i==j);
1335 
1336  vector<int> dim(2);
1337  dim << n,n;
1338  arraytype a(v,dim);
1339  a=f.jacobian(a);
1340  Q.resize(n,n);
1341  for(int i=0;i<n*n;i++)Q(i)=a[i];
1342  // Get Q_BB
1343  matrixtype QBB(nB,nB);
1344  for(int i=0;i<nB;i++)
1345  for(int j=0;j<nB;j++)
1346  QBB(i,j)=Q(cproj[i],cproj[j]);
1347  dmvnorm=MVNORM_t<scalartype>(QBB);
1348  }
1349  initialized=true;
1350  }
1351  vectortype projB(vectortype x){
1352  vectortype y(nB);
1353  for(int i=0;i<nB;i++)y[i]=x[cproj[i]];
1354  return y;
1355  }
1356  vectortype setZeroB(vectortype x){
1357  for(int i=0;i<nB;i++)x[cproj[i]]=scalartype(0);
1358  return x;
1359  }
1360  scalartype operator()(vectortype x){
1361  initialize(x.size());
1362  x=setZeroB(x);
1363  vector<int> dim(1);
1364  dim << x.size();
1365  arraytype xa(x,dim);
1366  vectortype y=projB(f.jacobian(xa));
1367  // f(x_A,0) - dmvnorm(y_B) + 2*dmvnorm(0)
1368  return f(xa) - dmvnorm(y) + 2*dmvnorm(y*scalartype(0));
1369  }
1370  /* array versions */
1371  arraytype projB(arraytype x){
1372  vectortype z((x.size()/n)*nB);
1373  vector<int> dim(x.dim);
1374  dim[dim.size()-1]=nB;
1375  arraytype y(z,dim);
1376  for(int i=0;i<nB;i++)y.col(i)=x.col(cproj[i]);
1377  return y;
1378  }
1379  arraytype setZeroB(arraytype x){
1380  for(int i=0;i<nB;i++)x.col(cproj[i])=x.col(0)*scalartype(0);
1381  return x;
1382  }
1383  arraytype jacobian(arraytype x){
1384  initialize(x.dim[x.dim.size()-1]);
1385  arraytype xa=setZeroB(x);
1386  arraytype y=projB(f.jacobian(xa));
1387  // WRONG: ----> return f.jacobian(xa) - dmvnorm.jacobian(y);
1388  // y=P*Q*Z*x so should be (P*Q*Z)' * dmvnorm.jacobian(y).
1389  // Note: only P is not symmetric.
1390  arraytype tmp=f.jacobian(xa);
1391  arraytype tmp0=tmp*scalartype(0);
1392  arraytype tmp2=dmvnorm.jacobian(y);
1393  // apply P'
1394  for(int i=0;i<nB;i++){
1395  tmp0.col(cproj[i])=tmp2.col(i);
1396  }
1397  // apply Q'(=Q)
1398  tmp0=f.jacobian(tmp0);
1399  // apply Z'(=Z)
1400  tmp0=setZeroB(tmp0);
1401  // Done:
1402  return tmp-tmp0;
1403  }
1404  int ndim(){return f.ndim();}
1405  VARIANCE_NOT_YET_IMPLEMENTED
1406  SIMULATE_NOT_YET_IMPLEMENTED
1407 };
1408 
1409 template <class distribution>
1410 PROJ_t<distribution> PROJ(distribution f_, vector<int> i){
1411  return PROJ_t<distribution>(f_,i);
1412 }
1413 
1414 
1415 #undef TYPEDEFS
1416 
1417 } // End namespace
#define SIMULATE_NOT_YET_IMPLEMENTED
Definition: density.hpp:54
array< Type > col(int i)
Extract sub-array with write access Index i refers to the outer-most (i.e. final) dimension...
Definition: tmbutils.hpp:202
vectortype cov(int n)
Covariance extractor. Run Youle-Walker recursions and return a vector of length n representing the au...
Definition: density.hpp:553
Collection of multivariate Gaussian distributions (members listed in density.hpp) ...
Definition: density.hpp:14
Array class used by TMB.
Definition: tmbutils.hpp:23
matrix< Type > jacobian(Functor f, vector< Type > x)
Calculate jacobian of vector function with vector values.
Definition: autodiff.hpp:203
scalartype operator()(vectortype x)
Evaluate the negative log density.
Definition: density.hpp:384
Multivariate normal distribution with unstructered correlation matrix.
Definition: density.hpp:260
Type dnorm(Type x, Type mean, Type sd, int give_log=0)
Probability density function of the normal distribution.
Definition: dnorm.hpp:17
MVNORM_t< scalartype > MVNORM(matrix< scalartype > Sigma, bool use_atomic=true)
Construct object to evaluate multivariate zero-mean normal density with user supplied covariance matr...
Definition: density.hpp:208
SEPARABLE_t< distribution1, distribution2 > SEPARABLE(distribution1 f_, distribution2 g_)
Construct object to evaluate the separable extension of two multivariate zero-mean normal densities...
Definition: density.hpp:1247
scalartype operator()(vectortype x, vectortype dx)
Evaluate the negative log density of the process x with nuisance parameters dx.
Definition: density.hpp:721
void simulate(vectortype &x)
Draw a simulation from the process.
Definition: density.hpp:448
Standardized normal distribution.
Definition: density.hpp:302
array< Type > transpose()
Array transpose (Special case of array permutation)
Definition: tmbutils.hpp:277
Stationary AR1 process.
Definition: density.hpp:375
matrixtype cov()
Covariance matrix extractor.
Definition: density.hpp:134
scalartype operator()(arraytype x)
Evaluate the negative log density.
Definition: density.hpp:969
VECSCALE_t< distribution > VECSCALE(distribution f_, vectortype scale_)
Construct object to evaluate a scaled density. See VECSCALE_t for details.
Definition: density.hpp:1057
scalartype operator()(scalartype x)
Evaluate the negative log density.
Definition: density.hpp:306
Continuous AR(2) process.
Definition: density.hpp:669
UNSTRUCTURED_CORR_t< scalartype > UNSTRUCTURED_CORR(vector< scalartype > x)
Construct object to evaluate the density with unstructured correlation matrix. See UNSTRUCTURED_CORR_...
Definition: density.hpp:289
scalartype operator()(vectortype x)
Evaluate the negative log density.
Definition: density.hpp:156
Apply scale transformation on a density.
Definition: density.hpp:960
Taped sorting of a vector.
Definition: order.hpp:20
Matrix exponential: 2x2 case which can be handled efficiently.
Definition: matexp.hpp:46
Matrix class used by TMB.
Definition: tmbutils.hpp:102
Type sum(Vector< Type > x)
Definition: convenience.hpp:33
Type rnorm(Type mu, Type sigma)
Simulate from a normal distribution.
matrix< Type > invertSparseMatrix(Eigen::SparseMatrix< Type > A)
Definition: spmat.hpp:140
Vector class used by TMB.
Definition: tmbutils.hpp:18
Apply a vector scale transformation on a density.
Definition: density.hpp:1019
scalartype operator()(arraytype x)
Evaluate the negative log density.
Definition: density.hpp:1028
scalartype operator()(vectortype x, vectortype keep)
Evaluate projected negative log density.
Definition: density.hpp:162
Matrix< scalartype, n1 *n3, n2 *n4 > kronecker(Matrix< scalartype, n1, n2 > x, Matrix< scalartype, n3, n4 > y)
Kronecker product of two matrices.
Definition: kronecker.hpp:12
Stationary AR(k) process.
Definition: density.hpp:495
scalartype operator()(vectortype x)
Evaluate the negative log density.
Definition: density.hpp:568
matrix< Type > matinvpd(matrix< Type > x, Type &logdet)
Matrix inverse and determinant.
Projection of multivariate gaussian variable.
Definition: density.hpp:1297
Separable extension of two densitites.
Definition: density.hpp:1106
Utility functions for TMB (automatically included)
Definition: concat.hpp:5
Multivariate normal distribution with user supplied covariance matrix.
Definition: density.hpp:112
Gaussian Markov Random Field.
Definition: density.hpp:828
GMRF_t< scalartype > GMRF(Eigen::SparseMatrix< scalartype > Q, int order, bool normalize=true)
Construct object to evaluate density of Gaussian Markov Random Field (GMRF) for sparse Q...
Definition: density.hpp:934
License: GPL v2