Changeset 802

Timestamp:

02/04/10 20:45:38 (15 years ago)

Author:

smidl

Message:

new class estudent + tests

Location:

library

Files:

• bdm/base/bdmbase.cpp (modified) (1 diff)
• bdm/base/bdmbase.h (modified) (1 diff)
• bdm/estim/arx.cpp (modified) (1 diff)
• bdm/stat/exp_family.h (modified) (4 diffs)
• tests/testsuite/epdf_test.cpp (modified) (1 diff)
• tests/testsuite/estudent.cfg (added)

library/bdm/base/bdmbase.cpp

@@ -479 +479 @@
 void epdf::log_write() const {
         if ( log_level == 10 ) {
-                to_setting ( logrec->L.log_to_setting ( logrec->ids ( 0 ) ) );
+                UI::save(this,  logrec->L.log_to_setting ( logrec->ids ( 0 ) ) );
         } else {
                 if ( log_level > 0 ) {

library/bdm/base/bdmbase.h

@@ -538 +538 @@
         virtual vec variance() const = 0;
 
+        //! return expected covariance -- default is diag(variance)!!
+        virtual mat covariance() const {return diag(variance());};
+        
         //! Lower and upper bounds of \c percentage % quantile, returns mean-2*sigma as default
         virtual void qbounds ( vec &lb, vec &ub, double percentage = 0.95 ) const {

library/bdm/estim/arx.cpp

@@ -215 +215 @@
 
 void ARX::from_setting ( const Setting &set ) {
+        BMEF::from_setting(set);
         shared_ptr<RV> yrv_ = UI::build<RV> ( set, "rv", UI::compulsory );
         shared_ptr<RV> rrv = UI::build<RV> ( set, "rgr", UI::compulsory );

library/bdm/stat/exp_family.h

@@ -209 +209 @@
                 return diag ( R.to_mat() );
         }
-//      mlnorm<sq_T>* condition ( const RV &rvn ) const ; <=========== fails to cmpile. Why?
+        mat covariance() const {
+                return R.to_mat();
+        }
+        //      mlnorm<sq_T>* condition ( const RV &rvn ) const ; <=========== fails to cmpile. Why?
         shared_ptr<pdf> condition ( const RV &rvn ) const;
 
@@ -250 +253 @@
 SHAREDPTR2 ( enorm, fsqmat );
 
+//forward declaration
+class mstudent;
+
+/*! distribution of multivariate Student t density
+
+Based on article by Genest and Zidek, 
+*/
+template<class sq_T>
+class estudent : public eEF{
+        protected:
+                //! mena value
+                vec mu;
+                //! matrix H
+                sq_T H;
+                //! degrees of freedom
+                double delta;
+        public:
+                double evallog_nn(const vec &val) const{
+                        double tmp = -0.5*H.logdet() - 0.5*(delta + dim) * log(1+ H.invqform(val - mu)/delta);
+                        return tmp;
+                }
+                double lognc() const {
+                        //log(pi) = 1.14472988584940
+                        double tmp = -lgamma(0.5*(delta+dim))+lgamma(0.5*delta) + 0.5*dim*(log(delta) + 1.14472988584940);
+                        return tmp;
+                }
+                void marginal (const RV &rvm, estudent<sq_T> &marg) const {
+                        ivec ind = rvm.findself_ids(rv); // indeces of rvm in rv
+                        marg._mu() = mu(ind);
+                        marg._H() = sq_T(H,ind);
+                        marg._delta() = delta;
+                        marg.validate();
+                }
+                shared_ptr<epdf> marginal(const RV &rvm) const {
+                        shared_ptr<estudent<sq_T> > tmp = new estudent<sq_T>;
+                        marginal(rvm, *tmp);
+                        return tmp;
+                }
+                vec sample() const NOT_IMPLEMENTED(vec(0))
+                
+                vec mean() const {return mu;}
+                mat covariance() const {return delta/(delta-2)*H.to_mat();}
+                vec variance() const {return diag(covariance());}
+                        //! \name access
+                //! @{
+                        //! access function
+                        vec& _mu() {return mu;}
+                        //! access function
+                        sq_T& _H() {return H;}
+                        //! access function
+                        double& _delta() {return delta;}
+                        //!@}
+                        //! todo
+                        void from_setting(const Setting &set){
+                                epdf::from_setting(set);
+                                mat H0;
+                                UI::get(H0,set, "H");
+                                H= H0; // conversion!!
+                                UI::get(delta,set,"delta");
+                                UI::get(mu,set,"mu");
+                        }
+                        void to_setting(Setting &set) const{
+                                epdf::to_setting(set);
+                                UI::save(H.to_mat(), set, "H");
+                                UI::save(delta, set, "delta");
+                                UI::save(mu, set, "mu");
+                        }
+                        void validate() {
+                                dim = H.rows();
+                        }
+};
+UIREGISTER2(estudent,fsqmat);
+UIREGISTER2(estudent,ldmat);
+UIREGISTER2(estudent,chmat);
 
 /*!
@@ -301 +378 @@
         };
 
+        //! marginal density (only student for now)
+        shared_ptr<epdf> marginal(const RV &rvm) const {
+                bdm_assert(dimx==1, "Not supported");
+                //TODO - this is too trivial!!!
+                ivec ind = rvm.findself_ids(rv);
+                if (min(ind)==0) { //assume it si 
+                        shared_ptr<estudent<ldmat> > tmp = new estudent<ldmat>;
+                        mat M;
+                        ldmat Vz;
+                        ldmat Lam;
+                        factorize(M,Vz,Lam);
+                        
+                        tmp->_mu() = M.get_col(0);
+                        ldmat H;
+                        Vz.inv(H);
+                        H *=Lam._D()(0)/nu;
+                        tmp->_H() = H;
+                        tmp->_delta() = nu;
+                        tmp->validate();
+                        return tmp;
+                }
+                return NULL;
+        }
         //! \name Access attributes
         //!@{
@@ -361 +461 @@
         void validate() {
                 // check sizes, rvs etc.
+                // also check if RV are meaningful!!!
+                // meaningful =  rv for theta  and rv for r are split!
         }
         void log_register ( bdm::logger& L, const string& prefix );

library/tests/testsuite/epdf_test.cpp

@@ -14 +14 @@
 TEST ( enorm_test ) {
         epdf_harness::test_config ( "enorm.cfg" );
+}
+
+TEST ( estudent_test ) {
+        epdf_harness::test_config ( "estudent.cfg" );
 }