Changeset 178

Timestamp:

10/15/08 19:04:09 (16 years ago)

Author:

smidl

Message:

Changes in basic structures! new methods

Location:

bdm/stat

Files:

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

bdm/stat/emix.cpp

@@ -3 +3 @@
 using namespace itpp;
 
-void emix::set_parameters ( const vec &w0, const Array<epdf*> &Coms0 ) {
-        w = w0/sum(w0);
+void emix::set_parameters ( const vec &w0, const Array<epdf*> &Coms0, bool copy ) {
+        w = w0/sum ( w0 );
         int i;
         for ( i=0;i<w.length();i++ ) {
                 it_assert_debug ( rv.equal ( Coms0 ( i )->_rv() ),"RVs do not match!" );
         }
-        Coms = Coms0;
+        if ( copy ) {
+                Coms.set_length(Coms0.length());
+                for ( i=0;i<w.length();i++ ) {*Coms ( i ) =*Coms0 ( i );}
+                destroyComs=true;
+        }
+        else {
+                Coms = Coms0;
+                destroyComs=false;
+        }
 }
 
@@ -37 +45 @@
 //                      rvc.add ( mpdfs ( i )->_rvc().subt ( rv ) ); //add rv to common rvs.
 //              };
-//              
+//
 // //           independent = true;
 //              //test rvc of mpdfs and fill rvinds

bdm/stat/emix.h

@@ -36 +36 @@
         //! Component (epdfs)
         Array<epdf*> Coms;
+        //!Flag if owning Coms
+        bool destroyComs;
 public:
         //!Default constructor
-        emix ( RV &rv ) : epdf ( rv ) {};
-        //! Set weights \c w and components \c R
-        void set_parameters ( const vec &w, const Array<epdf*> &Coms );
+        emix (const RV &rv ) : epdf ( rv ) {};
+        //! Set weights \c w and components \c Coms , Coms are not copied!
+        void set_parameters ( const vec &w, const Array<epdf*> &Coms, bool copy=true );
 
         vec sample() const;
@@ -58 +60 @@
         //! returns a pointer to the internal mean value. Use with Care!
         vec& _w() {return w;}
+        virtual ~emix(){if (destroyComs){for(int i=0;i<Coms.length();i++){delete Coms(i);}}}
+        //! Auxiliary function for taking ownership of the Coms()
+        void ownComs(){destroyComs=true;}
 };
 
@@ -102 +107 @@
                 ll = 0;
                 for ( int i = ( n - 1 );i >= 0;i-- ) {
-                        if ( rvcinds ( i ).length() > 0 ) {
-                                condi = zeros ( rvcsinrv.length() + rvcinds.length() );
+                        if (( rvcinds ( i ).length() > 0 )||( rvcsinrv ( i ).length() > 0 )) {
+                                condi = zeros ( rvcsinrv(i).length() + rvcinds(i).length() );
                                 // copy data in condition
                                 set_subvector ( condi,rvcinds ( i ), cond );
                                 // copy data in already generated sample
-                                set_subvector ( condi,rvcsinrv ( i ), smp );
+                                set_subvector ( condi,rvinrvcs ( i ), get_vec(smp,rvcsinrv(i)) );
 
                                 mpdfs ( i )->condition ( condi );

bdm/stat/libBM.cpp

@@ -157 +157 @@
 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
+        ivec valid;
+        if (tsize>0) {valid= itpp::find ( res == -1 );} //-1 => value not found => it remains
         return ( *this ) ( valid ); //keep those that were not found in rv2
 }
@@ -196 +197 @@
 void compositepdf::setrvc(const RV &rv, RV &rvc){
         for (int i = 0;i < n;i++ ) {
-                rvc.add ( mpdfs ( i )->_rvc().subt ( rv ) ); //add rv to common rvc
+                RV rvx = mpdfs ( i )->_rvc().subt ( rv );
+                rvc.add ( rvx ); //add rv to common rvc
         };
 }
@@ -202 +204 @@
 void compositepdf::setindices(const RV &rv){
         for (int i = 0;i < n;i++ ) {
-                        // find ith rv in common rv
+                // 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(),"" );
         }
 }

bdm/stat/libBM.h

@@ -140 +140 @@
 };
 
+class mpdf;
 
 //! Probability density function with numerical statistics, e.g. posterior density.
@@ -173 +174 @@
                 return x;
         }
+        //! 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;}
+        //! Return marginal density on the given RV, the remainig rvs are intergrated out
+        epdf* marginal(const RV &rv){it_warning("Not implemented"); return NULL;}
 
         //! return expected value
@@ -183 +188 @@
         //! modifier function - useful when copying epdfs
         void _renewrv(const RV &in_rv){rv=in_rv;}
+        //!
 };
 
@@ -213 +219 @@
         };
         //! Update \c ep so that it represents this mpdf conditioned on \c rvc = cond
-        virtual void condition ( const vec &cond ) {};
+        virtual void condition ( const vec &cond ) {it_error("Not implemented");};
 
         //! Shortcut for conditioning and evaluation of the internal epdf. In some cases,  this operation can be implemented efficiently.
@@ -252 +258 @@
                 //! Indeces of rvc in common rv
                 Array<ivec> rvcsinrv;
+                //! Indeces of common rv in rvc
+                Array<ivec> rvinrvcs;
         public:
-                compositepdf(Array<mpdf*> A0): n(A0.length()), mpdfs(A0), rvsinrv(n), rvcsinrv(n){};
+                compositepdf(Array<mpdf*> A0): n(A0.length()), mpdfs(A0), rvsinrv(n), rvcsinrv(n),rvinrvcs(n){};
                 RV getrv(bool checkoverlap=false);
                 void setrvc(const RV &rv, RV &rvc);
@@ -329 +337 @@
         vec logpred_m(const mat &dt)const{vec tmp(dt.cols());for(int i=0;i<dt.cols();i++){tmp(i)=logpred(dt.get_col(i));}return tmp;}
         
+        //!Constructs a predictive density (marginal density on data)
+        virtual epdf* predictor(const RV &rv){it_error("Not implemented");return NULL;};
+        
         //! Destructor for future use;
         virtual ~BM() {};

bdm/stat/libEF.cpp

@@ -49 +49 @@
 #define logpi 1.144729885849400163877476189
 #define log2pi 1.83787706640935
+#define Inf std::numeric_limits<double>::infinity()
 
         double nkG = 0.5* xdim* ( -nPsi *log2pi + sum ( log ( D ( xdim,D.length()-1 ) ) ) );
@@ -58 +59 @@
                      - 0.5*xdim* ( m*log2 + 0.5* ( xdim-1 ) *log2pi )  - lg;
 
+        it_assert_debug(((-nkG-nkW)>-Inf) && ((-nkG-nkW)<Inf), "ARX improper");
         return -nkG-nkW;
 }
@@ -244 +246 @@
 }
 
-void eEmp::set_parameters ( const vec &w0, epdf* epdf0 ) {
+void eEmp::set_parameters ( const vec &w0, const epdf* epdf0 ) {
         //it_assert_debug(rv==epdf0->rv(),"Wrong epdf0");
         w=w0;
@@ -253 +255 @@
         for ( int i=0;i<n;i++ ) {samples ( i ) =epdf0->sample();}
 }
+
+void eEmp::set_samples (  const epdf* epdf0 ) {
+        //it_assert_debug(rv==epdf0->rv(),"Wrong epdf0");
+        w=1;
+        w/=sum ( w );//renormalize
+
+        for ( int i=0;i<n;i++ ) {samples ( i ) =epdf0->sample();}
+}
+

bdm/stat/libEF.h

@@ -44 +44 @@
         virtual double lognc() const =0;
         //!TODO decide if it is really needed
-        virtual void dupdate ( mat &v ) {it_error ( "Not implemneted" );};
+        virtual void dupdate ( mat &v ) {it_error ( "Not implemented" );};
         //!Evaluate normalized log-probability
-        virtual double evalpdflog_nn ( const vec &val ) const{it_error ( "Not implemneted" );return 0.0;};
+        virtual double evalpdflog_nn ( const vec &val ) const{it_error ( "Not implemented" );return 0.0;};
         //!Evaluate normalized log-probability
         virtual double evalpdflog ( const vec &val ) const {return evalpdflog_nn ( val )-lognc();}
@@ -90 +90 @@
         //original Bayes
         void bayes ( const vec &dt );
-        //!Flatten the posterior
-        virtual void flatten ( BMEF * B ) {it_error ( "Not implemented" );}
-};
+        //!Flatten the posterior according to the given BMEF (of the same type!)
+        virtual void flatten ( const BMEF * B ) {it_error ( "Not implemented" );}
+        //!Flatten the posterior as if to keep nu0 data
+        virtual void flatten ( double nu0 ) {it_error ( "Not implemented" );}
+};
+
+template<class sq_T>
+class mlnorm;
 
 /*!
@@ -100 +105 @@
 */
 template<class sq_T>
-
 class enorm : public eEF {
 protected:
@@ -110 +114 @@
         int dim;
 public:
-//      enorm() :eEF() {};
         //!Default constructor
-        enorm ( RV &rv );
+        enorm ( const RV &rv );
         //! Set mean value \c mu and covariance \c R
         void set_parameters ( const vec &mu,const sq_T &R );
-        //! tupdate in exponential form (not really handy)
-        void tupdate ( double phi, mat &vbar, double nubar );
+        ////! tupdate in exponential form (not really handy)
+        //void tupdate ( double phi, mat &vbar, double nubar );
         //! dupdate in exponential form (not really handy)
         void dupdate ( mat &v,double nu=1.0 );
@@ -124 +127 @@
         mat sample ( int N ) const;
         double eval ( const vec &val ) const ;
-        double evalpdflog ( const vec &val ) const;
+        double evalpdflog_nn ( const vec &val ) const;
         double lognc () const;
         vec mean() const {return mu;}
-
+        mlnorm<sq_T>* condition ( const RV &rvn );
+        enorm<sq_T>* marginal ( const RV &rv );
 //Access methods
         //! returns a pointer to the internal mean value. Use with Care!
@@ -139 +143 @@
 
         //! access method
-        mat getR () {return R.to_mat();}
+//      mat getR () {return R.to_mat();}
 };
 
@@ -215 +219 @@
         };
         //!access function
-        vec& _beta() {return beta;}
+        vec& _beta()  {return beta;}
         //!Set internal parameters
-        void set_parameters(const vec &beta0){
-                if(beta0.length()!=beta.length()){
-                        it_assert_debug(rv.length()==1,"Undefined");
-                        rv.set_size(0,beta0.length());
+        void set_parameters ( const vec &beta0 ) {
+                if ( beta0.length() !=beta.length() ) {
+                        it_assert_debug ( rv.length() ==1,"Undefined" );
+                        rv.set_size ( 0,beta0.length() );
                 }
                 beta= beta0;
@@ -258 +262 @@
                 return pred.lognc()-lll;
         }
-        void flatten ( BMEF* B ) {
-                eDirich* E=dynamic_cast<eDirich*> ( B );
+        void flatten ( const BMEF* B ) {
+                const eDirich* E=dynamic_cast<const eDirich*> ( B );
                 // sum(beta) should be equal to sum(B.beta)
-                const vec &Eb=E->_beta();
+                const vec &Eb=const_cast<eDirich*> ( E )->_beta();
                 est.pow ( sum ( beta ) /sum ( Eb ) );
                 if ( evalll ) {last_lognc=est.lognc();}
@@ -267 +271 @@
         const epdf& _epdf() const {return est;};
         void set_parameters ( const vec &beta0 ) {
-                est.set_parameters(beta0);
+                est.set_parameters ( beta0 );
                 rv = est._rv();
-                if(evalll){last_lognc=est.lognc();}
+                if ( evalll ) {last_lognc=est.lognc();}
         }
 };
@@ -359 +363 @@
  \brief Normal distributed linear function with linear function of mean value;
 
- Mean value \f$mu=A*rvc\f$.
+ Mean value \f$mu=A*rvc+mu_0\f$.
 */
 template<class sq_T>
@@ -366 +370 @@
         enorm<sq_T> epdf;
         mat A;
+        vec mu_const;
         vec& _mu; //cached epdf.mu;
 public:
         //! Constructor
-        mlnorm ( RV &rv,RV &rvc );
+        mlnorm (const RV &rv, const RV &rvc );
         //! Set \c A and \c R
-        void set_parameters ( const  mat &A, const sq_T &R );
+        void set_parameters ( const  mat &A, const vec &mu0, const sq_T &R );
         //!Generate one sample of the posterior
-        vec samplecond ( vec &cond, double &lik );
+        vec samplecond (const vec &cond, double &lik );
         //!Generate matrix of samples of the posterior
-        mat samplecond ( vec &cond, vec &lik, int n );
+        mat samplecond (const vec &cond, vec &lik, int n );
         //! Set value of \c rvc . Result of this operation is stored in \c epdf use function \c _ep to access it.
-        void condition ( vec &cond );
+        void condition (const vec &cond );
 };
 
@@ -453 +458 @@
         //! Default constructor
         eEmp ( const RV &rv0 ,int n0 ) :epdf ( rv0 ),n ( n0 ),w ( n ),samples ( n ) {};
+        //! Set samples and weights
+        void set_parameters ( const vec &w0, const epdf* pdf0 );
         //! Set sample
-        void set_parameters ( const vec &w0, epdf* pdf0 );
+        void set_samples ( const epdf* pdf0 );
         //! Potentially dangerous, use with care.
         vec& _w()  {return w;};
@@ -476 +483 @@
 
 template<class sq_T>
-enorm<sq_T>::enorm ( RV &rv ) :eEF ( rv ), mu ( rv.count() ),R ( rv.count() ),dim ( rv.count() ) {};
+enorm<sq_T>::enorm ( const RV &rv ) :eEF ( rv ), mu ( rv.count() ),R ( rv.count() ),dim ( rv.count() ) {};
 
 template<class sq_T>
@@ -490 +497 @@
 };
 
-template<class sq_T>
-void enorm<sq_T>::tupdate ( double phi, mat &vbar, double nubar ) {
-        //
-};
+// template<class sq_T>
+// void enorm<sq_T>::tupdate ( double phi, mat &vbar, double nubar ) {
+//      //
+// };
 
 template<class sq_T>
@@ -531 +538 @@
 
 template<class sq_T>
-double enorm<sq_T>::evalpdflog ( const vec &val ) const {
+double enorm<sq_T>::evalpdflog_nn ( const vec &val ) const {
         // 1.83787706640935 = log(2pi)
-        return  -0.5* ( +R.invqform ( mu-val ) ) - lognc();
+        return  -0.5* ( R.invqform ( mu-val ) );// - lognc();
 };
 
@@ -539 +546 @@
 inline double enorm<sq_T>::lognc () const {
         // 1.83787706640935 = log(2pi)
-        return -0.5* ( R.cols() * 1.83787706640935 +R.logdet() );
-};
-
-template<class sq_T>
-mlnorm<sq_T>::mlnorm ( RV &rv0,RV &rvc0 ) :mEF ( rv0,rvc0 ),epdf ( rv0 ),A ( rv0.count(),rv0.count() ),_mu ( epdf._mu() ) {
+        return 0.5* ( R.cols() * 1.83787706640935 +R.logdet() );
+};
+
+template<class sq_T>
+mlnorm<sq_T>::mlnorm ( const RV &rv0, const RV &rvc0 ) :mEF ( rv0,rvc0 ),epdf ( rv0 ),A ( rv0.count(),rv0.count() ),_mu ( epdf._mu() ) {
         ep =&epdf;
 }
 
 template<class sq_T>
-void mlnorm<sq_T>::set_parameters ( const mat &A0, const sq_T &R0 ) {
+void mlnorm<sq_T>::set_parameters ( const mat &A0, const vec &mu0, const sq_T &R0 ) {
         epdf.set_parameters ( zeros ( rv.count() ),R0 );
         A = A0;
+        mu_const = mu0;
 }
 
 template<class sq_T>
-vec mlnorm<sq_T>::samplecond ( vec &cond, double &lik ) {
+vec mlnorm<sq_T>::samplecond (const  vec &cond, double &lik ) {
         this->condition ( cond );
         vec smp = epdf.sample();
@@ -562 +570 @@
 
 template<class sq_T>
-mat mlnorm<sq_T>::samplecond ( vec &cond, vec &lik, int n ) {
+mat mlnorm<sq_T>::samplecond (const vec &cond, vec &lik, int n ) {
         int i;
         int dim = rv.count();
@@ -579 +587 @@
 
 template<class sq_T>
-void mlnorm<sq_T>::condition ( vec &cond ) {
-        _mu = A*cond;
+void mlnorm<sq_T>::condition (const vec &cond ) {
+        _mu = A*cond + mu_const;
 //R is already assigned;
 }
 
+template<class sq_T>
+enorm<sq_T>* enorm<sq_T>::marginal ( const RV &rvn ) {
+        ivec irvn = rvn.dataind ( rv );
+
+        sq_T Rn ( R,irvn );
+        enorm<sq_T>* tmp = new enorm<sq_T>( rvn );
+        tmp->set_parameters ( mu ( irvn ), Rn );
+        return tmp;
+}
+
+template<class sq_T>
+mlnorm<sq_T>* enorm<sq_T>::condition ( const RV &rvn ) {
+
+        RV rvc = rv.subt ( rvn );
+        it_assert_debug ( ( rvc.count() +rvn.count() ==rv.count() ),"wrong rvn" );
+        //Permutation vector of the new R
+        ivec irvn = rvn.dataind ( rv );
+        ivec irvc = rvc.dataind ( rv );
+        ivec perm=concat ( irvn , irvc );
+        sq_T Rn ( R,perm );
+
+        //fixme - could this be done in general for all sq_T?
+        mat S=R.to_mat();
+        //fixme
+        int n=rvn.count()-1;
+        int end=R.rows()-1;
+        mat S11 = S.get ( 0,n, 0, n );
+        mat S12 = S.get ( rvn.count(), end, 0, n );
+        mat S22 = S.get ( rvn.count(), end, rvn.count(), end );
+
+        vec mu1 = mu ( irvn );
+        vec mu2 = mu ( irvc );
+        mat A=S12*inv ( S22 );
+        sq_T R_n ( S11 - A *S12.T() );
+
+        mlnorm<sq_T>* tmp=new mlnorm<sq_T> ( rvn,rvc );
+
+        tmp->set_parameters ( A,mu1-A*mu2,R_n );
+        return tmp;
+}
+
 ///////////