Changeset 182 for bdm/stat/emix.h

Timestamp:

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

Author:

smidl

Message:

compilation error!

Files:

• bdm/stat/emix.h (modified) (8 diffs)

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