Changeset 145 for bdm/stat

Timestamp:

08/20/08 15:41:21 (16 years ago)

Author:

smidl

Message:

Oprava dokumentace

Location:

bdm/stat

Files:

• emix.cpp (modified) (1 diff)
• emix.h (modified) (4 diffs)
• libBM.cpp (modified) (4 diffs)
• libBM.h (modified) (7 diffs)
• loggers.h (modified) (1 diff)

bdm/stat/emix.cpp

@@ -18 +18 @@
 
         int i=0;
-        while ( ( w ( i ) <u0 ) && ( i< ( w.length()-1 ) ) ) {i++;}
+        while ( ( cumDist ( i ) <u0 ) && ( i< ( w.length()-1 ) ) ) {i++;}
 
         return Coms ( i )->sample();

bdm/stat/emix.h

@@ -30 +30 @@
 
 */
-class emix : public epdf
-{
+class emix : public epdf {
         protected:
                 //! weights of the components
@@ -39 +38 @@
         public:
                 //!Default constructor
-                emix ( RV &rv ) : epdf ( rv ) {};
+                emix(RV &rv) : epdf(rv) {};
                 //! Set weights \c w and components \c R
-                void set_parameters ( const vec &w, const Array<epdf*> &Coms );
+                void set_parameters(const vec &w, const Array<epdf*> &Coms);
 
                 vec sample() const;
-                vec mean() const
-                {
-                        int i; vec mu=zeros ( rv.count() );
-                        for ( i=0;i<w.length();i++ ) {mu+=w ( i ) *Coms ( i )->mean(); }
+                vec mean() const {
+                        int i; vec mu = zeros(rv.count());
+                        for (i = 0;i < w.length();i++) {mu += w(i) * Coms(i)->mean(); }
                         return mu;
                 }
-                double evalpdflog ( const vec &val ) const
-                {
+                double evalpdflog(const vec &val) const {
                         int i;
-                        double sum=0.0;
-                        for ( i=0;i<w.length();i++ ) {sum+=w ( i ) *Coms ( i )->evalpdflog ( val );}
-                        return log ( sum );
+                        double sum = 0.0;
+                        for (i = 0;i < w.length();i++) {sum += w(i) * Coms(i)->evalpdflog(val);}
+                        return log(sum);
                 };
 
@@ -65 +62 @@
 /*! \brief Chain rule decomposition of epdf
 
+Probability density in the form of Chain-rule decomposition:
+\[
+f(x_1,x_2,x_3) = f(x_1|x_2,x_3)f(x_2,x_3)f(x_3)
+\]
+Note that
+*/
+class eprod: public epdf {
+        protected:
+                //
+                int n;
+                // pointers to epdfs
+                Array<epdf*> epdfs;
+                Array<mpdf*> mpdfs;
+                //
+                Array<ivec> rvinds;
+                Array<ivec> rvcinds;
+                //! Indicate independence of its factors
+                bool independent;
+                //! Indicate internal creation of mpdfs which must be destroyed
+                bool intermpdfs;
+        public:
+                //!Constructor from list of eFacs or list of mFacs
+                eprod(Array<mpdf*> mFacs): epdf(RV()), n(mFacs.length()), epdfs(n), mpdfs(mFacs), rvinds(n), rvcinds(n) {
+                        int i;
+                        intermpdfs = false;
+                        for (i = 0;i < n;i++) {
+                                rv.add(mpdfs(i)->_rv()); //add rv to common rvs.
+                                epdfs(i) = &(mpdfs(i)->_epdf()); // add pointer to epdf
+                        };
+                        independent = true;
+                        //test rvc of mpdfs and fill rvinds
+                        for (i = 0;i < n;i++) {
+                                // find ith rv in common rv
+                                rvinds(i) = mpdfs(i)->_rv().dataind(rv);
+                                // find ith rvc in common rv
+                                rvcinds(i) = mpdfs(i)->_rvc().dataind(rv);
+                                if (rvcinds(i).length()>0) {independent = false;}
+                        }
+
+                };
+                eprod(Array<epdf*> eFacs): epdf(RV()), n(eFacs.length()), epdfs(eFacs), mpdfs(n), rvinds(n), rvcinds(n) {
+                        int i;
+                        intermpdfs = true;
+                        for (i = 0;i < n;i++) {
+                                if (rv.add(eFacs(i)->_rv())) {it_error("Incompatible eFacs.rv!");} //add rv to common rvs.
+                                mpdfs(i) = new mepdf(*(epdfs(i)));
+                                epdfs(i) = &(mpdfs(i)->_epdf()); // add pointer to epdf
+                        };
+                        independent = true;
+                        //test rvc of mpdfs and fill rvinds
+                        for (i = 0;i < n;i++) {
+                                // find ith rv in common rv
+                                rvinds(i) = mpdfs(i)->_rv().dataind(rv);
+                        }
+                };
+
+                double evalpdflog(const vec &val) const {
+                        int i;
+                        double res = 0.0;
+                        for (i = n - 1;i > 0;i++) {
+                                if (rvcinds(i).length() > 0)
+                                        {mpdfs(i)->condition(val(rvcinds(i)));}
+                                // add logarithms
+                                res += epdfs(i)->evalpdflog(val(rvinds(i)));
+                        }
+                }
+                vec sample () const{
+                        vec smp=zeros(rv.count());
+                        for (int i = (n - 1);i >= 0;i--) {
+                                if (rvcinds(i).length() > 0)
+                                        {mpdfs(i)->condition(smp(rvcinds(i)));}
+                                set_subvector(smp,rvinds(i), epdfs(i)->sample());
+                        }                       
+                        return smp;
+                }
+                vec mean() const{
+                        vec tmp(rv.count());
+                        if (independent) {
+                                for (int i=0;i<n;i++) {
+                                        vec pom = epdfs(i)->mean();
+                                        set_subvector(tmp,rvinds(i), pom);
+                                }
+                        }
+                        else {
+                                int N=50*rv.count();
+                                it_warning("eprod.mean() computed by sampling");
+                                tmp = zeros(rv.count());
+                                for (int i=0;i<N;i++){ tmp += sample();}
+                                tmp /=N;
+                        }
+                        return tmp;
+                }
+                ~eprod(){if (intermpdfs) {for (int i=0;i<n;i++){delete mpdfs(i);}}};
+};
+
+/*! \brief Mixture of mpdfs with constant weights, all mpdfs are of equal type
 
 */
-class eprod: public epdf
-{
-        protected:
-                Array<epdf*> epdfs;
-                Array<mpdf*> mpdfs;
-        public:
-
-
-};
-
-/*! \brief Mixture of mpdfs with constant weights
-
-*/
-class mmix : public mpdf
-{
+class mmix : public mpdf {
         protected:
                 //! Component (epdfs)
@@ -89 +168 @@
         public:
                 //!Default constructor
-                mmix ( RV &rv, RV &rvc ) : mpdf ( rv, rvc ), Epdf ( rv ) {ep=&Epdf;};
+                mmix(RV &rv, RV &rvc) : mpdf(rv, rvc), Epdf(rv) {ep = &Epdf;};
                 //! Set weights \c w and components \c R
-                void set_parameters ( const vec &w, const Array<mpdf*> &Coms )
-                {
-                        Array<epdf*> Eps ( Coms.length() );
+                void set_parameters(const vec &w, const Array<mpdf*> &Coms) {
+                        Array<epdf*> Eps(Coms.length());
 
-                        for ( int i=0;i<Coms.length();i++ )
-                        {
-                                Eps ( i ) =& ( Coms ( i )->_epdf() );
+                        for (int i = 0;i < Coms.length();i++) {
+                                Eps(i) = & (Coms(i)->_epdf());
                         }
-                        Epdf.set_parameters ( w,Eps );
+                        Epdf.set_parameters(w, Eps);
                 };
 
-                void condition ( const vec &cond )
-                {
-                        for ( int i=0;i<Coms.length();i++ ) {Coms ( i )->condition ( cond );}
+                void condition(const vec &cond) {
+                        for (int i = 0;i < Coms.length();i++) {Coms(i)->condition(cond);}
                 };
 };

bdm/stat/libBM.cpp

@@ -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 ) {
+bool 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 );
+        ivec ind = rv2.findself ( *this ); //should be -1 all the time
+        ivec index = itpp::find(ind==-1); 
+        
+        if ( index.length() < rv2.len ) { //conflict
+                ids = concat ( ids, rv2.ids(index) );
+                sizes = concat ( sizes, rv2.sizes(index) );
+                times = concat ( times, rv2.times(index) );
+                names = concat ( names, rv2.names(to_Arr(index)) );
+        }
+        else {
+                ids = concat ( ids, rv2.ids );
+                sizes = concat ( sizes, rv2.sizes );
+                times = concat ( times, rv2.times );
+                names = concat ( names, rv2.names );
+        }
+        tsize = sum(sizes);
+        len = ids.length();
+        return (index.length()<rv2.len);
+
 //      return *this;
 };
@@ -58 +71 @@
 }
 
-RV RV::subselect ( ivec ind ) {
+RV RV::subselect ( ivec ind ) const {
         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 ) {
+RV RV::operator() ( ivec ind ) const {
         return RV ( ids ( ind ), names ( to_Arr ( ind ) ), sizes ( ind ), times ( ind ) );
 }
 
-bool RV::equal ( RV rv2 ) const {
-        return ( ids==rv2.ids ) && ( times == rv2.times ) && ( sizes==rv2.sizes );
+bool RV::equal ( const 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() ); 
+        mat X = zeros ( rv.count(), N );
+        for ( int i = 0;i < N;i++ ) X.set_col ( i, this->sample() );
         return X;
 };
@@ -88 +101 @@
 }
 
-ivec RV::indexlist() {
-        ivec indlist ( tsize );
+str RV::tostr() const {
+        ivec idlist ( tsize );
+        ivec tmlist ( 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++ ) {
+                idlist.set_subvector ( pos, pos + sizes ( i ) - 1, ids ( i ) );
+                tmlist.set_subvector ( pos, pos + sizes ( i ) - 1, times ( i ) );
+                pos += sizes ( i );
         }
-        return indlist;
+        return str ( idlist, tmlist );
+}
+
+ivec RV::dataind ( RV rv2 ) const {
+        str str2 = rv2.tostr();
+        ivec res ( 0 );
+        ivec part;
+        int i;
+        for ( i = 0;i < len;i++ ) {
+                part = itpp::find ( ( str2.ids == ids ( i ) ) & ( str2.times == times ( i ) ) );
+                res = concat ( res, part );
+        }
+        return res;
+}
+
+RV RV::subt ( const RV rv2 ) const {
+        ivec res = this->findself ( rv2 ); // nonzeros
+        ivec valid = itpp::find ( res == -1 ); //-1 => value not found => it remains
+        return ( *this ) ( valid ); //keep those that were not found in rv2
+}
+
+ivec RV::findself ( const RV &rv2 ) const {
+        int i, j;
+        ivec tmp = -ones_i ( len );
+        for ( i = 0;i < len;i++ ) {
+                for ( j = 0;j < rv2.length();j++ ) {
+                        if ( ( ids ( i ) == rv2.ids ( j ) ) & ( times ( i ) == rv2.times ( j ) ) ) {
+                                tmp ( i ) = j;
+                                break;
+                        }
+                }
+        }
+        return tmp;
 }
 

bdm/stat/libBM.h

@@ -18 +18 @@
 
 using namespace itpp;
+
+//! Structure of RV (used internally)
+class str{
+public:
+        ivec ids;
+        ivec times;
+        str(ivec ids0, ivec times0):ids(ids0),times(times0){
+                it_assert_debug(times0.length()==ids0.length(),"Incompatible input");
+        };
+};
 
 /*!
@@ -61 +71 @@
         //TODO why not inline and later??
 
-        //! Find indexes of another rv in self
-        ivec find ( RV rv2 );
+        //! Find indexes of self in another rv, \return ivec of the same size as self.
+        ivec findself (const RV &rv2 ) const;
         //! 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 );
-        //! Add (concat) another variable to the current one
-        friend RV concat ( const RV &rv1, const RV &rv2 );
+        bool equal (const RV &rv2 ) const;
+        //! Add (concat) another variable to the current one, \return 0 if all rv2 were added, 1 if rv2 is in conflict
+        bool add ( const RV &rv2 );
         //! Subtract  another variable from the current one
-        RV subt ( RV rv2 );
+        RV subt ( const RV rv2 ) const;
         //! Select only variables at indeces ind
-        RV subselect ( ivec ind );
+        RV subselect ( ivec ind ) const;
         //! Select only variables at indeces ind
-        RV operator() ( ivec ind );
-        //! Generate new \c RV with \c time shifted by delta.
+        RV operator() ( ivec ind ) const;
+        //! Shift \c time shifted by delta.
         void t ( int delta );
-        //! generate a list of indeces, i.e. which
-        ivec indexlist();
+        //! generate \c str from rv, by expanding sizes
+        str tostr() const;
+        //! generate indeces into \param crv data vector that form data vector of self.
+        ivec dataind(RV crv) const;
 
         //!access function
@@ -92 +102 @@
         std::string name ( int at ) {return names ( at );};
 };
+
+//! Concat two random variables
+RV concat ( const RV &rv1, const RV &rv2 );
 
 
@@ -146 +159 @@
         //! Destructor for future use;
         virtual ~epdf() {};
-        //! access function
-        RV _rv() const {return rv;}
+        //! access function, possibly dangerous!
+        RV& _rv() {return rv;}
 };
 
@@ -170 +183 @@
         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.
+        //! Returns \param 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 );
@@ -190 +203 @@
         //! access function
         RV _rvc() {return rvc;}
+        //! access function
+        RV _rv() {return rv;}
         //!access function
         epdf& _epdf() {return *ep;}
@@ -201 +216 @@
 public:
         //!Default constructor
-        mepdf ( const RV &rv, const RV &rvc, epdf* em ) :mpdf ( rv,rvc ) {ep=em;};
+        mepdf (epdf &em ) :mpdf ( em._rv(),RV() ) {ep=&em;};
 };
 

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 );}
+        //! Save values into an itfile named after \c fname.
         void itsave(const char* fname);
 };