Changeset 286 for bdm/stat

Timestamp:

03/05/09 14:03:35 (15 years ago)

Author:

smidl

Message:

make mpdm work again

Location:

bdm/stat

Files:

• emix.h (modified) (12 diffs)
• libBM.h (modified) (5 diffs)
• libEF.h (modified) (2 diffs)

bdm/stat/emix.h

@@ -18 +18 @@
 //#include <std>
 
-namespace bdm{
+namespace bdm {
 
 //this comes first because it is used inside emix!
@@ -47 +47 @@
         //!Default constructor. By default, the given epdf is not copied!
         //! It is assumed that this function will be used only temporarily.
-        mratio ( const epdf* nom0, const RV &rv, bool copy=false ) :mpdf ( ), dl ( rv,rvc,nom0->_rv() ) {
-                ep->set_rv(_rv());
-                set_rvc(nom0->_rv().subt ( rv ) );
-                if ( copy ) {
-//                      nom = nom0->_copy_();
-                        it_error ( "todo" );
-                        destroynom=true;
-                }
-                else {
-                        nom = nom0;
-                        destroynom = false;
-                }
-                it_assert_debug ( rvc.length() >0,"Makes no sense to use this object!" );
+        mratio ( const epdf* nom0, const RV &rv, bool copy=false ) :mpdf ( ), dl ( ) {
+                // adjust rv and rvc
+                rvc = nom0->_rv().subt ( rv );
+                dimc = rvc._dsize();
+                ep = new epdf;
+                ep->set_parameters(rv._dsize());
+                ep->set_rv(rv);
+                
+                //prepare data structures
+                if ( copy ) {it_error ( "todo" ); destroynom=true; }
+                else { nom = nom0; destroynom = false; }
+                it_assert_debug ( rvc.length() >0,"Makes no sense to use this object!" );               
+                
+                // build denominator
                 den = nom->marginal ( rvc );
+                dl.set_connection(rv,rvc,nom0->_rv());
         };
         double evallogcond ( const vec &val, const vec &cond ) {
@@ -67 +69 @@
                 dl.pushup_cond ( nom_val,val,cond );
                 tmp = exp ( nom->evallog ( nom_val ) - den->evallog ( cond ) );
-                it_assert_debug(std::isfinite(tmp),"Infinite value");
+                it_assert_debug ( std::isfinite ( tmp ),"Infinite value" );
                 return tmp;
         }
@@ -96 +98 @@
 public:
         //!Default constructor
-        emix (  ) : epdf ( ) {};
+        emix ( ) : epdf ( ) {};
         //! Set weights \c w and components \c Coms
         //!By default Coms are copied inside. Parameter \c copy can be set to false if Coms live externally. Use method ownComs() if Coms should be destroyed by the destructor.
@@ -103 +105 @@
         vec sample() const;
         vec mean() const {
-                int i; vec mu = zeros ( dim);
+                int i; vec mu = zeros ( dim );
                 for ( i = 0;i < w.length();i++ ) {mu += w ( i ) * Coms ( i )->mean(); }
                 return mu;
@@ -109 +111 @@
         vec variance() const {
                 //non-central moment
-                vec mom2 = zeros(dim);
-                for ( int i = 0;i < w.length();i++ ) {mom2 += w ( i ) * pow(Coms ( i )->mean(),2); }
+                vec mom2 = zeros ( dim );
+                for ( int i = 0;i < w.length();i++ ) {mom2 += w ( i ) * pow ( Coms ( i )->mean(),2 ); }
                 //central moment
-                return mom2-pow(mean(),2);
+                return mom2-pow ( mean(),2 );
         }
         double evallog ( const vec &val ) const {
@@ -118 +120 @@
                 double sum = 0.0;
                 for ( i = 0;i < w.length();i++ ) {sum += w ( i ) * exp ( Coms ( i )->evallog ( val ) );}
-                if (sum==0.0){sum=std::numeric_limits<double>::epsilon();}
+                if ( sum==0.0 ) {sum=std::numeric_limits<double>::epsilon();}
                 double tmp=log ( sum );
-                it_assert_debug(std::isfinite(tmp),"Infinite");
+                it_assert_debug ( std::isfinite ( tmp ),"Infinite" );
                 return tmp;
         };
@@ -151 +153 @@
         //!access function
         epdf* _Coms ( int i ) {return Coms ( i );}
+        void set_rv(const RV &rv){
+                epdf::set_rv(rv);
+                for(int i=0;i<Coms.length();i++){Coms(i)->set_rv(rv);}
+        }
 };
 
@@ -174 +180 @@
         mprod ( Array<mpdf*> mFacs ) : compositepdf ( mFacs ), mpdf (), epdfs ( n ), dls ( n ) {
                 ep=&dummy;
-                RV rv=getrv(true);
-                set_rv(rv);dummy.set_parameters(rv._dsize());
+                RV rv=getrv ( true );
+                set_rv ( rv );dummy.set_parameters ( rv._dsize() );
                 setrvc ( ep->_rv(),rvc );
                 // rv and rvc established = > we can link them with mpdfs
                 for ( int i = 0;i < n;i++ ) {
-                        dls ( i ) = new datalink_m2m ( mpdfs ( i )->_rv(), mpdfs ( i )->_rvc(), _rv(), _rvc() );
+                        dls ( i ) = new datalink_m2m;
+                        dls(i)->set_connection( mpdfs ( i )->_rv(), mpdfs ( i )->_rvc(), _rv(), _rvc() );
                 }
 
@@ -237 +244 @@
         Array<datalink*> dls;
 public:
-        eprod ( const Array<const epdf*> epdfs0 ) : epdf (  ),epdfs ( epdfs0 ),dls ( epdfs.length() ) {
+        eprod () : epdfs ( 0 ),dls ( 0 ) {};
+        void set_parameters ( const Array<const epdf*> &epdfs0, bool named=true ) {
+                epdfs=epdfs0;//.set_length ( epdfs0.length() );
+                dls.set_length ( epdfs.length() );
+
                 bool independent=true;
-                for ( int i=0;i<epdfs.length();i++ ) {
-                        independent=rv.add ( epdfs ( i )->_rv() );
-                        it_assert_debug ( independent==true, "eprod:: given components are not independent ." );
-                }
-                for ( int i=0;i<epdfs.length();i++ ) {
-                        dls ( i ) = new datalink ( epdfs ( i )->_rv() , rv );
+                if ( named ) {
+                        for ( int i=0;i<epdfs.length();i++ ) {
+                                independent=rv.add ( epdfs ( i )->_rv() );
+                                it_assert_debug ( independent==true, "eprod:: given components are not independent." );
+                        }
+                        dim=rv._dsize();
+                }
+                else {
+                        dim =0; for ( int i=0;i<epdfs.length();i++ ) {
+                                dim+=epdfs ( i )->dimension();
+                        }
+                }
+                //
+                int cumdim=0;
+                int dimi=0;
+                int i;
+                for ( i=0;i<epdfs.length();i++ ) {
+                        dls ( i ) = new datalink;
+                        if ( named ) {dls ( i )->set_connection ( epdfs ( i )->_rv() , rv );}
+                        else {
+                                dimi = epdfs ( i )->dimension();
+                                dls ( i )->set_connection ( dimi, dim, linspace ( cumdim,cumdim+dimi-1 ) );
+                                cumdim+=dimi;
+                        }
                 }
         }
@@ -260 +289 @@
                 for ( int i=0;i<epdfs.length();i++ ) {
                         vec pom = epdfs ( i )->mean();
-                        dls ( i )->pushup ( tmp, pow(pom,2) );
-                }
-                return tmp-pow(mean(),2);
+                        dls ( i )->pushup ( tmp, pow ( pom,2 ) );
+                }
+                return tmp-pow ( mean(),2 );
         }
         vec sample() const {
@@ -277 +306 @@
                         tmp+=epdfs ( i )->evallog ( dls ( i )->pushdown ( val ) );
                 }
-                it_assert_debug(std::isfinite(tmp),"Infinite");
+                it_assert_debug ( std::isfinite ( tmp ),"Infinite" );
                 return tmp;
         }

bdm/stat/libBM.h

@@ -261 +261 @@
         //! Lower and upper bounds of \c percentage % quantile, returns mean-2*sigma as default
         virtual void qbounds ( vec &lb, vec &ub, double percentage=0.95 ) const {
-                vec mea=mean(); vec std=sqrt(variance()); 
+                vec mea=mean(); vec std=sqrt ( variance() );
                 lb = mea-2*std; ub=mea+2*std;
         };
@@ -401 +401 @@
                 it_assert_debug ( v2v_up.length() ==downsize,"rv is not fully in rv_up" );
         }
+        //! set connection using indeces
+        void set_connection ( int ds, int us, const ivec &upind ) {
+                downsize = ds;
+                upsize = us;
+                v2v_up= upind;
+
+                it_assert_debug ( v2v_up.length() ==downsize,"rv is not fully in rv_up" );
+        }
         //! Get val for myself from val of "Up"
         vec pushdown ( const vec &val_up ) {
@@ -425 +433 @@
 
 public:
+        datalink_m2e() {};
         //! Constructor
-        datalink_m2e ( const RV &rv,  const RV &rvc, const RV &rv_up ) :
-                        datalink ( rv,rv_up ), condsize ( rvc._dsize() ) {
+        void set_connection ( const RV &rv,  const RV &rvc, const RV &rv_up ) {
+                datalink::set_connection ( rv,rv_up );
+                condsize=  rvc._dsize();
                 //establish v2c connection
                 rvc.dataind ( rv_up, v2c_lo, v2c_up );
@@ -454 +464 @@
 public:
         //! Constructor
-        datalink_m2m ( const RV &rv, const RV &rvc, const RV &rv_up, const RV &rvc_up ) :
-                        datalink_m2e ( rv, rvc, rv_up ) {
+        datalink_m2m() {};
+        void set_connection ( const RV &rv, const RV &rvc, const RV &rv_up, const RV &rvc_up ) {
+                datalink_m2e::set_connection ( rv, rvc, rv_up );
                 //establish c2c connection
                 rvc.dataind ( rvc_up, c2c_lo, c2c_up );
@@ -714 +725 @@
                 if ( opt_L_bounds ) {
                         vec ub,lb;
-                        posterior().qbounds(lb,ub);
+                        posterior().qbounds ( lb,ub );
                         L->logit ( LIDs ( 1 ), lb );
                         L->logit ( LIDs ( 2 ), ub );

bdm/stat/libEF.h

@@ -94 +94 @@
 //      virtual void flatten ( double nu0 ) {it_error ( "Not implemented" );}
 
-        BMEF* _copy_ ( bool changerv=false ) const {it_error ( "function _copy_ not implemented for this BM" ); return NULL;};
+        BMEF* _copy_ () const {it_error ( "function _copy_ not implemented for this BM" ); return NULL;};
 };
 
@@ -707 +707 @@
 class iW : public epdf{
         public:
-                vec sample(){}
+                vec sample(){return vec();}
 };