Changeset 404 for library/bdm/stat

Timestamp:

07/02/09 22:16:05 (15 years ago)

Author:

smidl

Message:

Change in epdf: evallog returns -inf for points out of support. Merger is aware of it now.

Location:

library/bdm/stat

Files:

• exp_family.cpp (modified) (3 diffs)
• exp_family.h (modified) (4 diffs)
• merger.cpp (modified) (1 diff)
• merger.h (modified) (1 diff)

library/bdm/stat/exp_family.cpp

@@ -24 +24 @@
         if ( dimx==1 ) { //same as the following, just quicker.
                 double r = val ( vend ); //last entry!
-                if (r<0) return -1e+100;
+                if (r<0) return -inf;
                 vec Psi ( nPsi+dimx );
                 Psi ( 0 ) = -1.0;
@@ -36 +36 @@
                 fsqmat R ( reshape ( val ( nPsi*dimx,vend ),dimx,dimx ) );
                 double ldetR=R.logdet();
-                if (ldetR) return -1e+100;
+                if (ldetR) return -inf;
                 mat Tmp=concat_vertical ( -eye ( dimx ),Th );
                 fsqmat iR ( dimx );
@@ -190 +190 @@
         int i;
 
-        if (any(val<=0.)) return -1e100;
-        if (any(beta<=0.)) return -1e100;
+        if (any(val<=0.)) return -inf;
+        if (any(beta<=0.)) return -inf;
         for ( i=0; i<dim; i++ ) {
                 res += ( alpha ( i ) - 1 ) *std::log ( val ( i ) ) - beta ( i ) *val ( i );

library/bdm/stat/exp_family.h

@@ -51 +51 @@
                                 double tmp;
                                 tmp= evallog_nn ( val )-lognc();
-                                it_assert_debug ( std::isfinite ( tmp ),"Infinite value" ); 
+//                              it_assert_debug ( std::isfinite ( tmp ),"Infinite value" ); 
                                 return tmp;}
                         //!Evaluate normalized log-probability for many samples
@@ -280 +280 @@
                         double evallog_nn ( const vec &val ) const
                         {
-                                double tmp; tmp= ( beta-1 ) *log ( val );               it_assert_debug ( std::isfinite ( tmp ),"Infinite value" );
+                                double tmp; tmp= ( beta-1 ) *log ( val );
+//                              it_assert_debug ( std::isfinite ( tmp ),"Infinite value" );
                                 return tmp;
                         };
@@ -290 +291 @@
                                 for ( int i=0;i<beta.length();i++ ) {lgb+=lgamma ( beta ( i ) );}
                                 tmp= lgb-lgamma ( gam );
-                                it_assert_debug ( std::isfinite ( tmp ),"Infinite value" );
+//                              it_assert_debug ( std::isfinite ( tmp ),"Infinite value" );
                                 return tmp;
                         };
@@ -489 +490 @@
 
                         double eval ( const vec &val ) const  {return nk;}
-                        double evallog ( const vec &val ) const  {return lnk;}
+                        double evallog ( const vec &val ) const  {
+                                if (any(val<low) && any(val>high)) {return inf;}
+                                else return lnk;
+                        }
                         vec sample() const
                         {

library/bdm/stat/merger.cpp

@@ -5 +5 @@
 namespace bdm
 {
-        vec merger_base::merge_points ( mat &lW )
-        {
+        vec merger_base::merge_points ( mat &lW ) {
                 int nu=lW.rows();
-                switch ( METHOD )
-                {
+                vec result;
+                ivec indW;
+                bool infexist;
+                switch ( METHOD ) {
                         case ARITHMETIC:
-                                return log ( sum ( exp ( lW ) ) ); //ugly!
+                                result= log ( sum ( exp ( lW ) ) ); //ugly!
                                 break;
                         case GEOMETRIC:
-                                return sum ( lW ) /nu;
+                                result= sum ( lW ) /nu;
                                 break;
                         case LOGNORMAL:
-                                vec mu = sum ( lW ) /nu; //mean of logs
-                                // vec lam = sum ( pow ( lW,2 ) )-nu*pow ( mu,2 ); ======= numerically unsafe!
-                                vec lam = sum ( pow ( lW-outer_product ( ones ( lW.rows() ),mu ),2 ) );
+                                vec sumlW=sum ( lW ) ;
+                                indW=find((sumlW<inf) & (sumlW>-inf));
+                                infexist=(indW.size()<lW.cols());
+                                vec mu;
+                                vec lam;
+                                if (infexist){
+                                        mu = sumlW(indW) /nu; //mean of logs
+                                        // 
+                                        mat validlW=lW.get_cols(indW);
+                                        lam = sum ( pow ( validlW-outer_product ( ones ( validlW.rows() ),mu ),2 ) );
+                                }
+                                else {
+                                        mu = sum ( lW ) /nu; //mean of logs
+                                        lam = sum ( pow ( lW-outer_product ( ones ( lW.rows() ),mu ),2 ) );
+                                }
+                                //
                                 double coef=0.0;
                                 vec sq2bl=sqrt ( 2*beta*lam ); //this term is everywhere
-                                switch ( nu )
-                                {
+                                switch ( nu ) {
                                         case 2:
                                                 coef= ( 1-0.5*sqrt ( ( 4.0*beta-3.0 ) /beta ) );
-                                                return  coef*sq2bl + mu ;
+                                                result =coef*sq2bl + mu ;
                                                 break;
-                                        case 3://Ratio of Bessel
-                                                coef = sqrt ( ( 3*beta-2 ) /3*beta );
-                                                return log ( besselk ( 0,sq2bl*coef ) ) - log ( besselk ( 0,sq2bl ) ) +  mu;
-                                                break;
-                                        case 4:
-                                                break;
-                                        default: // Approximate conditional density
+                                        // case 4: == can be done similar to case 2 - is it worth it???
+                                        default: // see accompanying document merge_lognorm_derivation.lyx
+                                                coef = sqrt(1-(nu+1)/(2*beta*nu));
+                                                result =log(besselk((nu-3)/2, sq2bl*coef))-log(besselk((nu-3)/2, sq2bl)) + mu;
                                                 break;
                                 }
                                 break;
                 }
-                return vec ( 0 );
+                if (infexist){
+                        vec tmp =-inf*ones(lW.cols());
+                        set_subvector(tmp, indW, result);
+                        return tmp;
+                }
+                else {return result;}
         }
 

library/bdm/stat/merger.h

@@ -242 +242 @@
                         vec mea = mean();
 
-                        cout << sum (w) << "," << w*w << endl;
+//                      cout << sum (w) << "," << w*w << endl;
 
                         mat Tmp = zeros (dim, dim);