Changeset 1189 for library/bdm
- Timestamp:
- 09/16/10 13:22:11 (14 years ago)
- Location:
- library/bdm
- Files:
-
- 3 modified
Legend:
- Unmodified
- Added
- Removed
-
library/bdm/estim/arx.h
r1176 r1189 221 221 222 222 void bayes(const vec &dt, const vec &psi){ 223 ldmat &Pbeta = est.P;224 ldmat &Palpha = alternative.P;223 // ldmat &Pbeta = est.P; 224 // ldmat &Palpha = alternative.P; 225 225 226 226 -
library/bdm/stat/exp_family.cpp
r1102 r1189 49 49 ChLam.inv ( iChLam ); 50 50 51 eWishartCh Omega; //inverse Wishart, result is R, 52 Omega.set_parameters ( iChLam, nu - 2*nPsi - dimx ); // 2*nPsi is there to match numercial simulations - check if analytically correct 51 eWishartCh Omega; //iWishart, same degrees of freedom as in Wishard in terms of delta 52 // delta = nu - npsi- p -1 53 Omega.set_parameters ( iChLam, nu - nPsi -dimx -1); 53 54 Omega.validate(); 54 55 55 mat OmChi;56 chmat OmChi; 56 57 mat Z ( M.rows(), M.cols() ); 57 58 58 59 mat Mi; 59 mat RChiT;60 chmat Ri; 60 61 mat tmp ( dimension(), n ); 61 62 M=M.T();// ugly hack == decide what to do with M. 62 63 for ( int i = 0; i < n; i++ ) { 63 64 OmChi = Omega.sample_mat(); 64 RChiT = inv ( OmChi ); 65 Z = randn ( M.rows(), M.cols() ); 66 Mi = M + RChiT * Z * inv ( Vz._L().T() * diag ( sqrt ( Vz._D() ) ) ); 67 68 tmp.set_col ( i, concat ( cvectorize ( Mi ), cvectorize ( RChiT*RChiT.T() ) ) ); 65 OmChi.inv(Ri ); 66 if (M._datasize()>0){ 67 Z = randn ( M.rows(), M.cols() ); 68 Mi = M + Ri._Ch().T() * Z * inv ( Vz._L().T() * diag ( sqrt ( Vz._D() ) ) ); 69 } 70 tmp.set_col ( i, concat ( cvectorize ( Mi ), cvectorize ( Ri.to_mat() ) ) ); 69 71 } 70 72 return tmp; … … 79 81 factorize ( M, Vz, Lam ); 80 82 81 chmat Ch; 82 Ch.setCh ( Lam._L() *diag ( sqrt ( Lam._D() ) ) ); 83 chmat Ch(Lam.to_mat()); 83 84 chmat iCh; 84 85 Ch.inv ( iCh ); 85 86 86 87 eWishartCh Omega; //inverse Wishart, result is R, 87 Omega.set_parameters ( iCh, nu - 2*nPsi - dimx); // 2*nPsi is there to match numercial simulations - check if analytically correct88 Omega.set_parameters ( iCh, nu - nPsi - dimx -1 ); // 2*nPsi is there to match numercial simulations - check if analytically correct 88 89 Omega.validate(); 89 90 90 91 chmat Omi; 91 Omi.setCh ( Omega.sample_mat() ); 92 92 Omi= Omega.sample_mat(); 93 Omi.inv ( Ri ); 94 93 95 if (M._datasize()>0) { 96 M = M.T(); 94 97 mat Z = randn ( M.rows(), M.cols() ); 95 Mi = M + Omi._Ch() * Z * inv ( Vz._L() * diag ( sqrt ( Vz._D() ) ) ); 96 } 97 Omi.inv ( Ri ); 98 Mi = M + Ri._Ch().T() * Z * inv ( Vz._L() * diag ( sqrt ( Vz._D() ) ) ); 99 } 98 100 } 99 101 -
library/bdm/stat/exp_family.h
r1175 r1189 388 388 * \brief Gauss-inverse-Wishart density stored in LD form 389 389 390 * For \f$p\f$-variate densities, given rv.count() should be \f$p imes\f$ V.rows(). 390 * For \f$p\f$-variate densities, 391 * 392 * \f$ M,R \sim GiW( V_t, \nu_t) \propto |R|^{0.5\nu}\exp(-1/2 tr(R^{-1}[I,M] V_t [I;M'])) \f$ 391 393 * 394 * Factorizes as: 395 * \f$ M|R \sim N( \hat{M}, R \otimes Vz^{-1}) \propto |R|^{-0.5dim(psi)}\exp(-1/2 tr((M-\hat{M})R^{-1}(M-\hat{M})Vz) \f$ 396 * \f$ R \sim iW( \Lambda,\delta) \propto |R|^{-0.5(\nu - dim(psi))} |\Lambda|^{}\exp(-1/2 tr(R^{-1}\Lambda) \f$ 397 * where in standard notation \f$ |R|^{-0.5(\delta + p +1)}\f$, i.e. 398 * \f$ \delta = \nu-dim(psi) -p-1 \f$ 392 399 */ 393 400 class egiw : public eEF { … … 426 433 427 434 void factorize ( mat &M, ldmat &Vz, ldmat &Lam ) const; 428 //! LS estimate of \f$ 435 //! LS estimate of \f$\theta\f$ 429 436 vec est_theta() const; 430 437 … … 1612 1619 1613 1620 /*! \brief Inverse Wishart density defined on Choleski decomposition 1621 * 1622 * here \f$\Omega \sim |\Sigma|^{-0.5(\delta}|\Omega|^{0.5(\delta -p -1)} \exp(-1/2 tr(\Omega\Sigma^{-1}))\$f 1614 1623 */ 1615 1624 class eWishartCh : public epdf { 1616 1625 protected: 1617 1626 //! Upper-Triagle of Choleski decomposition of \f$ \Psi \f$ 1618 chmat Y;1627 chmat Sigma; 1619 1628 //! dimension of matrix \f$ \Psi \f$ 1620 1629 int p; … … 1624 1633 //! Set internal structures 1625 1634 void set_parameters ( const mat &Y0, const double delta0 ) { 1626 Y= chmat ( Y0 );1635 Sigma = chmat ( Y0 ); 1627 1636 delta = delta0; 1628 p = Y.rows();1637 p = Sigma.rows(); 1629 1638 } 1630 1639 //! Set internal structures 1631 1640 void set_parameters ( const chmat &Y0, const double delta0 ) { 1632 Y= Y0;1641 Sigma = Y0; 1633 1642 delta = delta0; 1634 p = Y.rows();1643 p = Sigma.rows(); 1635 1644 } 1636 1645 … … 1640 1649 } 1641 1650 1642 //! Sample matrix argument 1643 mat sample_mat() const { 1644 mat X = zeros ( p, p ); 1651 //! Sample matrix argument - lower 1652 //! Using Bartlet decomposition: W=L A A^T L^T where A is lower triag and L is choleski factor of Sigma. 1653 chmat sample_mat() const { 1654 mat A_T = zeros ( p, p ); // A transpose 1645 1655 1646 1656 //sample diagonal … … 1648 1658 GamRNG.setup ( 0.5* ( delta - i ) , 0.5 ); // no +1 !! index if from 0 1649 1659 #pragma omp critical 1650 X( i, i ) = sqrt ( GamRNG() );1660 A_T ( i, i ) = sqrt ( GamRNG() ); 1651 1661 } 1652 1662 //do the rest … … 1654 1664 for ( int j = i + 1; j < p; j++ ) { 1655 1665 #pragma omp critical 1656 X( i, j ) = NorRNG.sample();1666 A_T ( i, j ) = NorRNG.sample(); 1657 1667 } 1658 1668 } 1659 return X*Y._Ch();// return upper triangular part of the decomposition 1669 chmat tmp; 1670 tmp._Ch()=A_T*Sigma._Ch(); 1671 return tmp;// return upper triangular part of the decomposition 1660 1672 } 1661 1673 1662 1674 vec sample () const { 1663 return vec ( sample_mat(). _data(), p*p );1675 return vec ( sample_mat().to_mat()._data(), p*p ); 1664 1676 } 1665 1677 … … 1673 1685 //! fast access function y0 will be copied into Y.Ch. 1674 1686 void setY ( const mat &Ch0 ) { 1675 copy_vector ( dim, Ch0._data(), Y._Ch()._data() );1687 copy_vector ( dim, Ch0._data(), Sigma._Ch()._data() ); 1676 1688 } 1677 1689 1678 1690 //! fast access function y0 will be copied into Y.Ch. 1679 1691 void _setY ( const vec &ch0 ) { 1680 copy_vector ( dim, ch0._data(), Y._Ch()._data() );1692 copy_vector ( dim, ch0._data(), Sigma._Ch()._data() ); 1681 1693 } 1682 1694 1683 1695 //! access function 1684 1696 const chmat& getY() const { 1685 return Y;1697 return Sigma; 1686 1698 } 1687 1699 }; … … 1715 1727 vec sample() const { 1716 1728 mat iCh; 1717 iCh = inv ( W.sample_mat() );1729 iCh = inv ( W.sample_mat()._Ch() ); 1718 1730 return vec ( iCh._data(), dim ); 1719 1731 }