Changeset 182 for bdm/stat

Timestamp:

10/17/08 10:28:15 (16 years ago)

Author:

smidl

Message:

compilation error!

Location:

bdm/stat

Files:

• emix.cpp (modified) (2 diffs)
• emix.h (modified) (8 diffs)
• libBM.cpp (modified) (2 diffs)
• libBM.h (modified) (4 diffs)
• libEF.cpp (modified) (2 diffs)
• libEF.h (modified) (6 diffs)

bdm/stat/emix.cpp

@@ -31 +31 @@
 }
 
-// mprod::mprod ( Array<mpdf*> mFacs, bool overlap) : mpdf ( RV(), RV() ), n ( mFacs.length() ), epdfs ( n ), mpdfs ( mFacs ), rvinds ( n ), rvcinrv ( n ), rvcsinrvc ( n ) {
+emix* emix::marginal(const RV &rv) const{
+        Array<epdf*> Cn(Coms.length());
+        for(int i=0;i<Coms.length();i++){Cn(i)=Coms(i)->marginal(rv);}
+        emix* tmp = new emix(rv);
+        tmp->set_parameters(w,Coms,false);
+        tmp->ownComs();
+        return tmp;
+}
+
+mratio* emix::condition(const RV &rv) const{
+        return new mratio(this,rv);
+};
+
+// mprod::mprod ( Array<mpdf*> mFacs, bool overlap) : mpdf ( RV(), RV() ), n ( mFacs.length() ), epdfs ( n ), mpdfs ( mFacs ), rvinds ( n ), rvcinrv ( n ), irvcs_rvc ( n ) {
 //              int i;
 //              bool rvaddok;
@@ -54 +67 @@
 //                      rvcinrv ( i ) = mpdfs ( i )->_rvc().dataind ( rv );
 //                      // find ith rvc in common rv
-//                      rvcsinrvc ( i ) = mpdfs ( i )->_rvc().dataind ( rvc );
+//                      irvcs_rvc ( i ) = mpdfs ( i )->_rvc().dataind ( rvc );
 //                      //
 // /*                   if ( rvcinrv ( i ).length() >0 ) {independent = false;}
-//                      if ( rvcsinrvc ( i ).length() >0 ) {independent = false;}*/
+//                      if ( irvcs_rvc ( i ).length() >0 ) {independent = false;}*/
 //              }
 //      };

bdm/stat/emix.h

@@ -19 +19 @@
 
 using namespace itpp;
+
+//this comes first because it is used inside emix!
+
+/*! \brief Class representing ratio of two densities
+which arise e.g. by applying the Bayes rule.
+It represents density in the form:
+\f[
+f(rv|rvc) = \frac{f(rv,rvc)}{f(rvc)}
+\f]
+where \f$ f(rvc) = \int f(rv,rvc) d\ rv \f$.
+
+In particular this type of arise by conditioning of a mixture model.
+
+At present the only supported operation is evalcond().
+ */
+class mratio: public mpdf {
+        protected:
+        //! Nominator in the form of mpdf
+                const epdf* nom;
+        //!Denominator in the form of epdf
+                epdf* den;
+        //!flag for destructor 
+                bool destroynom;
+        public:
+        //!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(rv,RV()){
+                        if (copy){
+//                      nom = nom0->_copy_();
+                                destroynom=true;
+                        } else {
+                                nom = nom0;
+                                destroynom = false;
+                        }
+                        rvc = nom->_rv().subt(rv);
+                        it_assert_debug(rvc.length()>0,"Makes no sense to use this object!");
+                        den = nom->marginal(rvc);
+                };
+                double evalcond(const vec &val, const vec &cond) const {
+                        return exp(nom->evalpdflog(concat(val,cond)) - den->evalpdflog(cond));
+                }
+        //! Object takes ownership of nom and will destroy it
+                void ownnom(){destroynom=true;}
+        //! Default destructor
+                ~mratio(){delete den; if (destroynom) {delete nom;}}
+};
 
 /*!
@@ -41 +87 @@
         //!Default constructor
         emix ( const RV &rv ) : epdf ( rv ) {};
-        //! Set weights \c w and components \c Coms , Coms are not copied!
+        //! Set weights \c w and components \c Coms
+        //!By default Coms are copied inside. \param copy can be set to false if Coms live externally. Use method ownComs() if Coms should be destroyed by the destructor.
         void set_parameters ( const vec &w, const Array<epdf*> &Coms, bool copy=true );
 
@@ -56 +103 @@
                 return log ( sum );
         };
+        vec evalpdflog_m ( const mat &Val ) const {
+                vec x=ones(Val.cols());
+                vec y(Val.cols());
+                for (int i = 0; i < w.length(); i++ ) {
+                        y = w ( i )*exp(Coms ( i )->evalpdflog_m ( Val ) );
+                        elem_mult_inplace(y,x); //result is in x
+                }
+                return log(x);
+        };
+
+        emix* marginal ( const RV &rv ) const;
+        mratio* condition ( const RV &rv ) const; //why not mratio!!
 
 //Access methods
@@ -78 +137 @@
         Array<epdf*> epdfs;
         //! Indeces of rvc in common rvc
-        Array<ivec> rvcsinrvc;
+        Array<ivec> irvcs_rvc;
         //! Indeces of common rvc in rvcs
-        Array<ivec> rvcinrvcs;
+        Array<ivec> irvc_rvcs;
 public:
         /*!\brief Constructor from list of mFacs,
         */
-        mprod ( Array<mpdf*> mFacs ) : compositepdf ( mFacs ), mpdf ( getrv ( true ),RV() ), epdfs ( n ), rvcsinrvc ( n ) {
+        mprod ( Array<mpdf*> mFacs ) : compositepdf ( mFacs ), mpdf ( getrv ( true ),RV() ), epdfs ( n ), irvcs_rvc ( n ), irvc_rvcs ( n ) {
                 setrvc ( rv,rvc );
                 setindices ( rv );
                 for ( int i = 0;i < n;i++ ) {
-                        // find ith rvc in common rv
-                        rvcsinrvc ( i ) = mpdfs ( i )->_rvc().dataind ( rvc );
-                        //and vice-versa
-                        rvcinrvcs ( i ) = rvc.dataind ( mpdfs ( i )->_rvc() );
+                        // establish link between rvc and rvcs
+                        rvc.dataind ( mpdfs ( i )->_rvc(), irvc_rvcs ( i ), irvcs_rvc ( i ) );
                 }
 
@@ -99 +156 @@
         };
 
-        double evalcond ( const vec &val, const vec &cond ) const {
+        double evalcond ( const vec &val, const vec &cond ) {
                 int i;
                 double res = 0.0;
                 vec condi;
-                for ( i = n - 1;i > 0;i++ ) {
-                        if ( ( rvcsinrvc ( i ).length() > 0 ) || ( rvcsinrv ( i ).length() > 0 ) ) {
-                                condi = zeros ( rvcsinrvc ( i ).length() +rvcsinrv ( i ).length() );
+                for ( i = n - 1;i >= 0;i-- ) {
+                        if ( ( irvcs_rvc ( i ).length() > 0 ) || ( irvcs_rv ( i ).length() > 0 ) ) {
+                                condi = zeros ( irvcs_rvc ( i ).length() +irvcs_rv ( i ).length() );
                                 //copy part of the rvc into condi
-                                set_subvector ( condi, rvcinrvcs ( i ), get_vec ( cond,rvcsinrvc ( i ) ) );
+                                set_subvector ( condi, irvcs_rvc ( i ), get_vec ( cond,irvc_rvcs ( i ) ) );
                                 //copy part of the rv into condi
-                                set_subvector ( condi, rvinrvcs ( i ), get_vec ( val,rvcsinrv ( i ) ) );
+                                set_subvector ( condi, irvcs_rv ( i ), get_vec ( val,irv_rvcs ( i ) ) );
                                 mpdfs ( i )->condition ( condi );
                         }
@@ -115 +172 @@
                         res += epdfs ( i )->evalpdflog ( val ( rvsinrv ( i ) ) );
                 }
-                return res;
+                return exp ( res );
         }
         vec samplecond ( const vec &cond, double &ll ) {
@@ -123 +180 @@
                 ll = 0;
                 for ( int i = ( n - 1 );i >= 0;i-- ) {
-                        if ( ( rvcsinrvc ( i ).length() > 0 ) || ( rvcsinrv ( i ).length() > 0 ) ) {
-                                condi = zeros ( rvcsinrv ( i ).length() + rvcsinrvc ( i ).length() );
+                        if ( ( irvcs_rvc ( i ).length() > 0 ) || ( irvcs_rv ( i ).length() > 0 ) ) {
+                                condi = zeros ( irvcs_rv ( i ).length() + irvcs_rvc ( i ).length() );
                                 // copy data in condition
-                                set_subvector ( condi,rvcsinrvc ( i ), cond );
+                                set_subvector ( condi,irvcs_rvc ( i ), get_vec ( cond, irvc_rvcs ( i ) ) );
                                 // copy data in already generated sample
-                                set_subvector ( condi,rvinrvcs ( i ), get_vec ( smp,rvcsinrv ( i ) ) );
+                                set_subvector ( condi,irvcs_rv ( i ), get_vec ( smp,irv_rvcs ( i ) ) );
 
                                 mpdfs ( i )->condition ( condi );
@@ -220 +277 @@
         };
 };
+
 #endif //MX_H

bdm/stat/libBM.cpp

@@ -151 +151 @@
                 }
         }
+        it_assert_debug(res.length()==tsize,"this rv is not fully present in crv!");
         return res;
 
+}
+
+void RV::dataind (const RV &rv2, ivec &selfi, ivec &rv2i) const {
+        //clean results
+        selfi.set_size(0);
+        rv2i.set_size(0);
+        
+        // just in case any rv is empty
+        if ((len==0)||(rv2.length()==0)){return;} 
+        
+        //find comon rv
+        ivec cids=itpp::find(this->findself(rv2)>=0);
+        
+        // index of 
+        if ( cids.length()>0 ) {
+                str str1 = tostr();
+                str str2 = rv2.tostr(); 
+                
+                ivec part1;
+                ivec part2;
+                int i,j;
+                // find common rv in strs
+                for ( j=0; j < cids.length();j++ ) {
+                        i = cids(j);
+                        part1 = itpp::find ( ( str1.ids == ids ( i ) ) & ( str1.times == times ( i ) ) );
+                        part2 = itpp::find ( ( str2.ids == ids ( i ) ) & ( str2.times == times ( i ) ) );
+                        selfi = concat ( selfi, part1 );
+                        rv2i = concat ( rv2i, part2 );
+                }
+        }
+        it_assert_debug(selfi.length() == rv2i.length(),"this should not happen!");
 }
 
@@ -206 +238 @@
                 // find ith rv in common rv
                 rvsinrv ( i ) = mpdfs ( i )->_rv().dataind ( rv );
-                rvcsinrv ( i ) = mpdfs ( i )->_rvc().dataind ( rv );
-                rvinrvcs ( i ) = rv.dataind(mpdfs ( i )->_rvc());
-                //rvcsinrv and rvinrvcs must be 1-to-1
-                it_assert_debug(rvcsinrv ( i ).length()==rvinrvcs(i).length(),"" );
+                // establish link between rv and rvc
+                rv.dataind(mpdfs(i)->_rvc(), irv_rvcs(i), irvcs_rv(i));
         }
 }

bdm/stat/libBM.h

@@ -92 +92 @@
         //! generate \c str from rv, by expanding sizes
         str tostr() const;
-        //! generate indeces into \param crv data vector that form data vector of self.
+        //! when this rv is a part of bigger rv, this function returns indeces of self in the data vector of the bigger crv.  
+        //! Then, data can be copied via: data_of_this = cdata(ind);
         ivec dataind(const RV &crv) const;
+        //! generate mutual indeces when copying data betwenn self and crv. 
+        //! Data are copied via: data_of_this(selfi) = data_of_rv2(rv2i)
+        void dataind(const RV &rv2, ivec &selfi, ivec &rv2i) const;
 
         //!access function
@@ -178 +182 @@
         }
         //! Return conditional density on the given RV, the remaining rvs will be in conditioning
-        mpdf* condition(const RV &rv){it_warning("Not implemented"); return NULL;}
+        virtual mpdf* condition(const RV &rv) const {it_warning("Not implemented"); return NULL;}
         //! Return marginal density on the given RV, the remainig rvs are intergrated out
-        epdf* marginal(const RV &rv){it_warning("Not implemented"); return NULL;}
+        epdf* marginal(const RV &rv) const {it_warning("Not implemented"); return NULL;}
 
         //! return expected value
@@ -248 +252 @@
         //!Default constructor
         mepdf (epdf &em ) :mpdf ( em._rv(),RV() ) {ep=&em;};
+        void condition(const vec &cond){}
 };
 
@@ -260 +265 @@
                 Array<ivec> rvsinrv;
                 //! Indeces of rvc in common rv
-                Array<ivec> rvcsinrv;
+                Array<ivec> irvcs_rv;
                 //! Indeces of common rv in rvc
-                Array<ivec> rvinrvcs;
+                Array<ivec> irv_rvcs;
         public:
-                compositepdf(Array<mpdf*> A0): n(A0.length()), mpdfs(A0), rvsinrv(n), rvcsinrv(n),rvinrvcs(n){};
+                compositepdf(Array<mpdf*> A0): n(A0.length()), mpdfs(A0), rvsinrv(n), irvcs_rv(n),irv_rvcs(n){};
                 //! find common rv, flag \param checkoverlap modifies whether overlaps are acceptable
                 RV getrv(bool checkoverlap=false);

bdm/stat/libEF.cpp

@@ -95 +95 @@
         
         // set mean value
-        M= iLsub*L ( xdim,end,0,xdim-1 );
+        M= L ( xdim,end,0,xdim-1 ).T()*iLsub;
         R= ldR.to_mat()  ;
 }
@@ -207 +207 @@
                 N_babies ( j ) ++;
         }
-
         // We have assigned new babies for each Particle
         // Now, we fill the resulting index such that:

bdm/stat/libEF.h

@@ -130 +130 @@
         double lognc () const;
         vec mean() const {return mu;}
-        mlnorm<sq_T>* condition ( const RV &rvn );
-        enorm<sq_T>* marginal ( const RV &rv );
+        mlnorm<sq_T>* condition ( const RV &rvn ) const ;
+        enorm<sq_T>* marginal ( const RV &rv ) const;
 //Access methods
         //! returns a pointer to the internal mean value. Use with Care!
@@ -384 +384 @@
         //! Set value of \c rvc . Result of this operation is stored in \c epdf use function \c _ep to access it.
         void condition (const vec &cond );
+        
+/*      template<class sq_M>
+        friend std::ostream &operator<< ( std::ostream &os, const mlnorm<sq_M> &ml );*/
 };
 
@@ -594 +597 @@
 
 template<class sq_T>
-enorm<sq_T>* enorm<sq_T>::marginal ( const RV &rvn ) {
+enorm<sq_T>* enorm<sq_T>::marginal ( const RV &rvn ) const {
         ivec irvn = rvn.dataind ( rv );
 
@@ -604 +607 @@
 
 template<class sq_T>
-mlnorm<sq_T>* enorm<sq_T>::condition ( const RV &rvn ) {
+mlnorm<sq_T>* enorm<sq_T>::condition ( const RV &rvn ) const {
 
         RV rvc = rv.subt ( rvn );
@@ -620 +623 @@
         int end=R.rows()-1;
         mat S11 = S.get ( 0,n, 0, n );
-        mat S12 = S.get ( rvn.count(), end, 0, n );
+        mat S12 = S.get (0, n , rvn.count(), end );
         mat S22 = S.get ( rvn.count(), end, rvn.count(), end );
 
@@ -636 +639 @@
 ///////////
 
+// template<class sq_T>
+// std::ostream &operator<< ( std::ostream &os, const mlnorm<sq_T> &ml ){
+//      os << "A:"<< ml.A<<endl;
+//      os << "mu:"<< ml.mu_const<<endl;
+//      os << "R:" << ml.est._R().to_mat() <<endl;
+//      return os;
+// }
+// ;
 
 #endif //EF_H