Changeset 102

Timestamp:

05/12/08 17:34:07 (16 years ago)

Author:

smidl

Message:

corrections of sampling

Location:

bdm/stat

Files:

• libBM.cpp (modified) (4 diffs)
• libBM.h (modified) (11 diffs)
• libEF.cpp (modified) (1 diff)
• libEF.h (modified) (2 diffs)
• loggers.cpp (modified) (1 diff)
• loggers.h (modified) (1 diff)

bdm/stat/libBM.cpp

@@ -7 +7 @@
 using std::cout;
 
-void RV::init( ivec in_ids, Array<std::string> in_names, ivec in_sizes, ivec in_times ) {
+void RV::init ( ivec in_ids, Array<std::string> in_names, ivec in_sizes, ivec in_times ) {
         //
         int i;
@@ -13 +13 @@
         //PRUDENT_MODE
         // All vetors should be of same length
-        if (     ( len != in_ids.length() ) || \
+        if ( ( len != in_ids.length() ) || \
                 ( len != in_names.length() ) || \
                 ( len != in_sizes.length() ) || \
-                ( len != in_times.length() )) {
-                it_error( "RV::RV inconsistent length of input vectors." );
+                ( len != in_times.length() ) ) {
+                it_error ( "RV::RV inconsistent length of input vectors." );
         }
 
@@ -25 +25 @@
         times = in_times;
         tsize = 0;
-        for(i=0;i<len;i++){tsize+=sizes(i);}
+        for ( i=0;i<len;i++ ) {tsize+=sizes ( i );}
 };
 
@@ -32 +32 @@
 }
 
-RV::RV () : tsize(0),len(0){};
+RV::RV () : tsize ( 0 ),len ( 0 ) {};
 
-void RV::add (const RV &rv2) {
+void RV::add ( const RV &rv2 ) {
         // TODO
         tsize+=rv2.tsize;
         len +=rv2.len;
-        ids=concat(ids,rv2.ids);
-        sizes=concat(sizes,rv2.sizes);
-        times=concat(times,rv2.times);
-        names=concat(names,rv2.names);
+        ids=concat ( ids,rv2.ids );
+        sizes=concat ( sizes,rv2.sizes );
+        times=concat ( times,rv2.times );
+        names=concat ( names,rv2.names );
 //      return *this;
 };
 
 RV::RV ( ivec in_ids ) {
-        
+
         len = in_ids.length();
-        Array<std::string> A( len );
+        Array<std::string> A ( len );
         std::string rvstr = "rv";
 
         for ( int i = 0; i < len;i++ ) {
-                A( i ) = rvstr + to_str(i);
+                A ( i ) = rvstr + to_str ( i );
         }
 
-        init ( in_ids, A, ones_i( len ), zeros_i( len ) );
+        init ( in_ids, A, ones_i ( len ), zeros_i ( len ) );
 }
 
-RV RV::subselect(ivec ind){
-        return RV(ids(ind), names(to_Arr(ind)), sizes(ind), times(ind));
+RV RV::subselect ( ivec ind ) {
+        return RV ( ids ( ind ), names ( to_Arr ( ind ) ), sizes ( ind ), times ( ind ) );
 }
 
-void RV::t(int delta){ times +=delta;}
+void RV::t ( int delta ) { times +=delta;}
 
-RV RV::operator()(ivec ind){
-        return RV(ids(ind), names(to_Arr(ind)), sizes(ind), times(ind));
+RV RV::operator() ( ivec ind ) {
+        return RV ( ids ( ind ), names ( to_Arr ( ind ) ), sizes ( ind ), times ( ind ) );
 }
 
-std::ostream &operator<<( std::ostream &os, const RV &rv ) {
+bool RV::equal ( RV rv2 ) const {
+        return ( ids==rv2.ids ) && ( times == rv2.times ) && ( sizes==rv2.sizes );
+}
+
+mat epdf::sampleN ( int N ) const {
+        mat X =zeros ( rv.count(),N );
+        for ( int i=0;i<N;i++ ) X.set_col ( i,this->sample() ); 
+        return X;
+};
+
+
+std::ostream &operator<< ( std::ostream &os, const RV &rv ) {
 
         for ( int i = 0; i < rv.len ;i++ ) {
-                os << rv.ids( i ) << "(" << rv.sizes( i ) << ")" <<  // id(size)=
-                "=" << rv.names( i )  << "_{"  << rv.times( i ) << "}; "; //name_{time}
+                os << rv.ids ( i ) << "(" << rv.sizes ( i ) << ")" <<  // id(size)=
+                "=" << rv.names ( i )  << "_{"  << rv.times ( i ) << "}; "; //name_{time}
         }
         return os;
 }
 
-ivec RV::indexlist(){
-        ivec indlist(tsize);
+ivec RV::indexlist() {
+        ivec indlist ( tsize );
         int i;
         int pos = 0;
-        for(i=0;i<len;i++){
-                indlist.set_subvector(pos,pos+sizes(i)-1, ids(i));
+        for ( i=0;i<len;i++ ) {
+                indlist.set_subvector ( pos,pos+sizes ( i )-1, ids ( i ) );
         }
         return indlist;
 }
 
-RV concat(const RV &rv1, const RV &rv2 ){
+RV concat ( const RV &rv1, const RV &rv2 ) {
         RV pom = rv1;
-        pom.add(rv2); 
+        pom.add ( rv2 );
         return pom;
 }

bdm/stat/libBM.h

@@ -31 +31 @@
         //! len = number of individual rvs
         int len;
-        //! Vector of unique IDs 
+        //! Vector of unique IDs
         ivec ids;
         //! Vector of sizes
@@ -42 +42 @@
 private:
         //! auxiliary function used in constructor
-        void init ( ivec in_ids, Array<std::string> in_names, ivec in_sizes, ivec in_times);
+        void init ( ivec in_ids, Array<std::string> in_names, ivec in_sizes, ivec in_times );
 public:
         //! Full constructor which is called by the others
-        RV ( ivec in_ids, Array<std::string> in_names, ivec in_sizes, ivec in_times);
+        RV ( ivec in_ids, Array<std::string> in_names, ivec in_sizes, ivec in_times );
         //! default constructor
         RV ( ivec ids );
@@ -63 +63 @@
         //! Find indexes of another rv in self
         ivec find ( RV rv2 );
+        //! Compare if \c rv2 is identical to this \c RV
+        bool equal ( RV rv2 ) const;
         //! Add (concat) another variable to the current one
-        void add (const RV &rv2 );
+        void add ( const RV &rv2 );
         //! Add (concat) another variable to the current one
-        friend RV concat (const RV &rv1, const RV &rv2 );
+        friend RV concat ( const RV &rv1, const RV &rv2 );
         //! Subtract  another variable from the current one
         RV subt ( RV rv2 );
@@ -79 +81 @@
 
         //!access function
-        Array<std::string>& _names(){return names;};
-
-        //!access function
-        int id(int at){return ids(at);};
-        //!access function
-        int size(int at){return sizes(at);};
-        //!access function
-        int time(int at){return times(at);};
-        //!access function
-        std::string name(int at){return names(at);};
+        Array<std::string>& _names() {return names;};
+
+        //!access function
+        int id ( int at ) {return ids ( at );};
+        //!access function
+        int size ( int at ) {return sizes ( at );};
+        //!access function
+        int time ( int at ) {return times ( at );};
+        //!access function
+        std::string name ( int at ) {return names ( at );};
 };
 
@@ -100 +102 @@
 public:
         //!default constructor
-        fnc(int dy):dimy(dy){};
+        fnc ( int dy ) :dimy ( dy ) {};
         //! function evaluates numerical value of \f$f(x)\f$ at \f$x=\f$ \c cond
         virtual vec eval ( const vec &cond ) {
                 return vec ( 0 );
-        }; 
+        };
 
         //! access function
@@ -131 +133 @@
         //! Returns a sample, \f$x\f$ from density \f$epdf(rv)\f$
         virtual vec sample () const =0;
+        //! Returns N samples from density \f$epdf(rv)\f$
+        virtual mat sampleN ( int N ) const;
         //! Compute probability of argument \c val
-        virtual double eval ( const vec &val ) const {return exp(this->evalpdflog(val));};
+        virtual double eval ( const vec &val ) const {return exp ( this->evalpdflog ( val ) );};
 
         //! Compute log-probability of argument \c val
         virtual double evalpdflog ( const vec &val ) const =0;
 
-        //! return expected value 
+        //! return expected value
         virtual vec mean() const =0;
-        
+
         //! Destructor for future use;
         virtual ~epdf() {};
@@ -163 +167 @@
 //      virtual fnc moment ( const int order = 1 );
         //! Returns a sample from the density conditioned on \c cond, \f$x \sim epdf(rv|cond)\f$. \param cond is numeric value of \c rv \param ll is a return value of log-likelihood of the sample.
-        virtual vec samplecond ( vec &cond, double &ll ) {this->condition(cond);vec temp= ep->sample();ll=ep->evalpdflog(temp);return temp;};
+        virtual vec samplecond ( vec &cond, double &ll ) {this->condition ( cond );vec temp= ep->sample();ll=ep->evalpdflog ( temp );return temp;};
+        //! Returns N samples from the density conditioned on \c cond, \f$x \sim epdf(rv|cond)\f$. \param cond is numeric value of \c rv \param ll is a return value of log-likelihood of the sample.
+        virtual mat samplecond ( vec &cond, vec &ll, int N ) {
+                this->condition ( cond );
+                mat temp ( rv.count(),N ); vec smp ( rv.count() ); for ( int i=0;i<N;i++ ) {smp=ep->sample() ;temp.set_col ( i, smp );ll ( i ) =ep->evalpdflog ( smp );}
+                return temp;
+        };
         //! Update \c ep so that it represents this mpdf conditioned on \c rvc = cond
         virtual void condition ( const vec &cond ) {};
-        
+
         //! Shortcut for conditioning and evaluation of the internal epdf. In some cases,  this operation can be implemented efficiently.
-        virtual double evalcond (const vec &dt, const vec &cond ) {this->condition(cond);return ep->eval(dt);};
+        virtual double evalcond ( const vec &dt, const vec &cond ) {this->condition ( cond );return ep->eval ( dt );};
 
         //! Destructor for future use;
@@ -176 +186 @@
         mpdf ( const RV &rv0, const RV &rvc0 ) :rv ( rv0 ),rvc ( rvc0 ) {};
         //! access function
-        RV _rvc(){return rvc;}
-        //!access function
-        epdf& _epdf(){return *ep;}
+        RV _rvc() {return rvc;}
+        //!access function
+        epdf& _epdf() {return *ep;}
 };
 
@@ -230 +240 @@
 
         //!Default constructor
-        BM(const RV &rv0) :rv(rv0), ll ( 0 ),evalll ( true ) {//Fixme: test rv 
+        BM ( const RV &rv0 ) :rv ( rv0 ), ll ( 0 ),evalll ( true ) {//Fixme: test rv
         };
 
@@ -240 +250 @@
         void bayes ( mat Dt );
         //! Returns a pointer to the epdf representing posterior density on parameters. Use with care!
-        virtual epdf& _epdf()=0;
+        virtual epdf& _epdf() =0;
 
         //! Destructor for future use;
@@ -267 +277 @@
         virtual void condition ( const vec &val ) =0;
         //! Default constructor
-        BMcond(RV &rv0):rvc(rv0){};
+        BMcond ( RV &rv0 ) :rvc ( rv0 ) {};
         //! Destructor for future use
-        virtual ~BMcond(){};
+        virtual ~BMcond() {};
         //! access function
         const RV& _rvc() const {return rvc;}

bdm/stat/libEF.cpp

@@ -65 +65 @@
 }
 
-mat egamma::sample ( int N ) const {
-        mat Smp ( rv.count(),N );
-        int i,j;
-
-        for ( i=0; i<rv.count(); i++ ) {
-                GamRNG.setup ( alpha ( i ),beta ( i ) );
-
-                for ( j=0; j<N; j++ ) {
-                        Smp ( i,j ) = GamRNG();
-                }
-        }
-
-        return Smp;
-}
+// mat egamma::sample ( int N ) const {
+//      mat Smp ( rv.count(),N );
+//      int i,j;
+// 
+//      for ( i=0; i<rv.count(); i++ ) {
+//              GamRNG.setup ( alpha ( i ),beta ( i ) );
+// 
+//              for ( j=0; j<N; j++ ) {
+//                      Smp ( i,j ) = GamRNG();
+//              }
+//      }
+// 
+//      return Smp;
+// }
 
 double egamma::evalpdflog ( const vec &val ) const {

bdm/stat/libEF.h

@@ -162 +162 @@
         vec sample() const;
         //! TODO: is it used anywhere?
-        mat sample ( int N ) const;
+//      mat sample ( int N ) const;
         double evalpdflog ( const vec &val ) const;
         double lognc () const;
@@ -206 +206 @@
         vec sample() const {
                 vec smp ( rv.count() ); UniRNG.sample_vector ( rv.count(),smp );
-                return low+distance*smp;
+                return low+elem_mult(distance,smp);
         }
         //! set values of \c low and \c high

bdm/stat/loggers.cpp

@@ -7 +7 @@
 #include <io.h>
 #endif
+
+void memlog::itsave(const char* fname){
+        it_file itf(fname);
+        int i;
+        for (i=0; i<entries.length();i++){
+         itf << Name(names(i)) << vectors(i);
+        }
+}
 
 void dirfilelog::init() {

bdm/stat/loggers.h

@@ -84 +84 @@
         void step(bool final=false) {if ( ind<maxlen ) ind++; else it_error ( "memlog::ind is too high;" );}
         void logit ( int id, vec v ) {vectors ( id ).set_row ( ind,v );}
+        void itsave(const char* fname);
 };