Changeset 404
- Timestamp:
- 07/02/09 22:16:05 (15 years ago)
- Location:
- library/bdm
- Files:
-
- 7 modified
-
base/bdmbase.h (modified) (2 diffs)
-
itpp_ext.cpp (modified) (3 diffs)
-
itpp_ext.h (modified) (3 diffs)
-
stat/exp_family.cpp (modified) (3 diffs)
-
stat/exp_family.h (modified) (4 diffs)
-
stat/merger.cpp (modified) (1 diff)
-
stat/merger.h (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
library/bdm/base/bdmbase.h
r395 r404 265 265 virtual mat sample_m(int N) const; 266 266 //! Compute log-probability of argument \c val 267 //! In case the argument is out of suport return -Infinity 267 268 virtual double evallog(const vec &val) const {it_error("not implemneted"); return 0.0;}; 268 269 //! Compute log-probability of multiple values argument \c val … … 381 382 this->condition(cond); 382 383 tmp = ep->evallog(dt); 383 it_assert_debug(std::isfinite(tmp), "Infinite value");384 // it_assert_debug(std::isfinite(tmp), "Infinite value"); 384 385 return tmp; 385 386 }; -
library/bdm/itpp_ext.cpp
r343 r404 63 63 64 64 // Gamma 65 66 Gamma_RNG::Gamma_RNG ( double a, double b ) {67 setup ( a,b );68 }69 70 71 bvec operator& ( const bvec &a, const bvec &b ) {72 it_assert_debug ( b.size() ==a.size(), "operator&(): Vectors of different lengths" );73 74 bvec temp ( a.size() );75 for ( int i = 0;i < a.size();i++ ) {76 temp ( i ) = a ( i ) & b ( i );77 }78 return temp;79 }80 81 bvec operator| ( const bvec &a, const bvec &b ) {82 it_assert_debug ( b.size() !=a.size(), "operator&(): Vectors of different lengths" );83 84 bvec temp ( a.size() );85 for ( int i = 0;i < a.size();i++ ) {86 temp ( i ) = a ( i ) | b ( i );87 }88 return temp;89 }90 65 #define log std::log 91 66 #define exp std::exp … … 93 68 #define R_FINITE std::isfinite 94 69 70 bvec operator>(const vec &t1, const vec &t2) { 71 it_assert_debug(t1.length() == t2.length(), "Vec<>::operator>(): different size of vectors"); 72 bvec temp(t1.length()); 73 for (int i = 0; i < t1.length(); i++) 74 temp(i) = (t1[i] > t2[i]); 75 return temp; 76 } 77 78 bvec operator<(const vec &t1, const vec &t2) { 79 it_assert_debug(t1.length() == t2.length(), "Vec<>::operator>(): different size of vectors"); 80 bvec temp(t1.length()); 81 for (int i = 0; i < t1.length(); i++) 82 temp(i) = (t1[i] < t2[i]); 83 return temp; 84 } 85 86 87 bvec operator& ( const bvec &a, const bvec &b ) { 88 it_assert_debug ( b.size() ==a.size(), "operator&(): Vectors of different lengths" ); 89 90 bvec temp ( a.size() ); 91 for ( int i = 0;i < a.size();i++ ) { 92 temp ( i ) = a ( i ) & b ( i ); 93 } 94 return temp; 95 } 96 97 bvec operator| ( const bvec &a, const bvec &b ) { 98 it_assert_debug ( b.size() !=a.size(), "operator&(): Vectors of different lengths" ); 99 100 bvec temp ( a.size() ); 101 for ( int i = 0;i < a.size();i++ ) { 102 temp ( i ) = a ( i ) | b ( i ); 103 } 104 return temp; 105 } 106 107 //#if 0 108 Gamma_RNG::Gamma_RNG ( double a, double b ) { 109 setup ( a,b ); 110 } 95 111 double Gamma_RNG::sample() { 96 112 //A copy of rgamma code from the R package!! … … 290 306 } 291 307 292 308 //#endif 293 309 std::string num2str(double d){ 294 310 char tmp[20];//that should do -
library/bdm/itpp_ext.h
r328 r404 26 26 bvec operator& ( const bvec &a, const bvec &b ); 27 27 bvec operator| ( const bvec &a, const bvec &b ); 28 29 bvec operator>(const vec &t1, const vec &t2); 28 30 29 // template<class Num_T> 31 bvec operator<(const vec &t1, const vec &t2); 32 33 // template<class Num_T> 30 34 // void set_subvector(vec &ov, ivec &iv, const Vec<Num_T> &v); 31 35 … … 38 42 } 39 43 44 const double inf =std::numeric_limits<double>::infinity(); 45 46 //#if 0 40 47 /*! 41 48 \brief Gamma distribution … … 74 81 }; 75 82 bool qr ( const mat &A, mat &R ); 76 83 //#endif 77 84 //! reimplementation of matlab num2str 78 85 std::string num2str(double d); -
library/bdm/stat/exp_family.cpp
r395 r404 24 24 if ( dimx==1 ) { //same as the following, just quicker. 25 25 double r = val ( vend ); //last entry! 26 if (r<0) return - 1e+100;26 if (r<0) return -inf; 27 27 vec Psi ( nPsi+dimx ); 28 28 Psi ( 0 ) = -1.0; … … 36 36 fsqmat R ( reshape ( val ( nPsi*dimx,vend ),dimx,dimx ) ); 37 37 double ldetR=R.logdet(); 38 if (ldetR) return - 1e+100;38 if (ldetR) return -inf; 39 39 mat Tmp=concat_vertical ( -eye ( dimx ),Th ); 40 40 fsqmat iR ( dimx ); … … 190 190 int i; 191 191 192 if (any(val<=0.)) return - 1e100;193 if (any(beta<=0.)) return - 1e100;192 if (any(val<=0.)) return -inf; 193 if (any(beta<=0.)) return -inf; 194 194 for ( i=0; i<dim; i++ ) { 195 195 res += ( alpha ( i ) - 1 ) *std::log ( val ( i ) ) - beta ( i ) *val ( i ); -
library/bdm/stat/exp_family.h
r395 r404 51 51 double tmp; 52 52 tmp= evallog_nn ( val )-lognc(); 53 it_assert_debug ( std::isfinite ( tmp ),"Infinite value" );53 // it_assert_debug ( std::isfinite ( tmp ),"Infinite value" ); 54 54 return tmp;} 55 55 //!Evaluate normalized log-probability for many samples … … 280 280 double evallog_nn ( const vec &val ) const 281 281 { 282 double tmp; tmp= ( beta-1 ) *log ( val ); it_assert_debug ( std::isfinite ( tmp ),"Infinite value" ); 282 double tmp; tmp= ( beta-1 ) *log ( val ); 283 // it_assert_debug ( std::isfinite ( tmp ),"Infinite value" ); 283 284 return tmp; 284 285 }; … … 290 291 for ( int i=0;i<beta.length();i++ ) {lgb+=lgamma ( beta ( i ) );} 291 292 tmp= lgb-lgamma ( gam ); 292 it_assert_debug ( std::isfinite ( tmp ),"Infinite value" );293 // it_assert_debug ( std::isfinite ( tmp ),"Infinite value" ); 293 294 return tmp; 294 295 }; … … 489 490 490 491 double eval ( const vec &val ) const {return nk;} 491 double evallog ( const vec &val ) const {return lnk;} 492 double evallog ( const vec &val ) const { 493 if (any(val<low) && any(val>high)) {return inf;} 494 else return lnk; 495 } 492 496 vec sample() const 493 497 { -
library/bdm/stat/merger.cpp
r399 r404 5 5 namespace bdm 6 6 { 7 vec merger_base::merge_points ( mat &lW ) 8 { 7 vec merger_base::merge_points ( mat &lW ) { 9 8 int nu=lW.rows(); 10 switch ( METHOD ) 11 { 9 vec result; 10 ivec indW; 11 bool infexist; 12 switch ( METHOD ) { 12 13 case ARITHMETIC: 13 re turnlog ( sum ( exp ( lW ) ) ); //ugly!14 result= log ( sum ( exp ( lW ) ) ); //ugly! 14 15 break; 15 16 case GEOMETRIC: 16 re turnsum ( lW ) /nu;17 result= sum ( lW ) /nu; 17 18 break; 18 19 case LOGNORMAL: 19 vec mu = sum ( lW ) /nu; //mean of logs 20 // vec lam = sum ( pow ( lW,2 ) )-nu*pow ( mu,2 ); ======= numerically unsafe! 21 vec lam = sum ( pow ( lW-outer_product ( ones ( lW.rows() ),mu ),2 ) ); 20 vec sumlW=sum ( lW ) ; 21 indW=find((sumlW<inf) & (sumlW>-inf)); 22 infexist=(indW.size()<lW.cols()); 23 vec mu; 24 vec lam; 25 if (infexist){ 26 mu = sumlW(indW) /nu; //mean of logs 27 // 28 mat validlW=lW.get_cols(indW); 29 lam = sum ( pow ( validlW-outer_product ( ones ( validlW.rows() ),mu ),2 ) ); 30 } 31 else { 32 mu = sum ( lW ) /nu; //mean of logs 33 lam = sum ( pow ( lW-outer_product ( ones ( lW.rows() ),mu ),2 ) ); 34 } 35 // 22 36 double coef=0.0; 23 37 vec sq2bl=sqrt ( 2*beta*lam ); //this term is everywhere 24 switch ( nu ) 25 { 38 switch ( nu ) { 26 39 case 2: 27 40 coef= ( 1-0.5*sqrt ( ( 4.0*beta-3.0 ) /beta ) ); 28 re turncoef*sq2bl + mu ;41 result =coef*sq2bl + mu ; 29 42 break; 30 case 3://Ratio of Bessel 31 coef = sqrt ( ( 3*beta-2 ) /3*beta ); 32 return log ( besselk ( 0,sq2bl*coef ) ) - log ( besselk ( 0,sq2bl ) ) + mu; 33 break; 34 case 4: 35 break; 36 default: // Approximate conditional density 43 // case 4: == can be done similar to case 2 - is it worth it??? 44 default: // see accompanying document merge_lognorm_derivation.lyx 45 coef = sqrt(1-(nu+1)/(2*beta*nu)); 46 result =log(besselk((nu-3)/2, sq2bl*coef))-log(besselk((nu-3)/2, sq2bl)) + mu; 37 47 break; 38 48 } 39 49 break; 40 50 } 41 return vec ( 0 ); 51 if (infexist){ 52 vec tmp =-inf*ones(lW.cols()); 53 set_subvector(tmp, indW, result); 54 return tmp; 55 } 56 else {return result;} 42 57 } 43 58 -
library/bdm/stat/merger.h
r399 r404 242 242 vec mea = mean(); 243 243 244 cout << sum (w) << "," << w*w << endl;244 // cout << sum (w) << "," << w*w << endl; 245 245 246 246 mat Tmp = zeros (dim, dim);
