Changeset 404 for library/bdm/stat/merger.cpp

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.

Files:

• library/bdm/stat/merger.cpp (modified) (1 diff)

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;}
         }