Changeset 33 for bdm/stat

Show
Ignore:
Timestamp:
03/05/08 16:01:56 (16 years ago)
Author:
smidl
Message:

Oprava PF a MPF + jejich implementace pro pmsm system

Location:
bdm/stat
Files:
5 modified

Legend:

Unmodified
Added
Removed
  • bdm/stat/libBM.h

    r32 r33  
    3131        //! len = number of individual rvs 
    3232        int len; 
     33        //! Vector of unique IDs  
    3334        ivec ids; 
     35        //! Vector of sizes 
    3436        ivec sizes; 
     37        //! Vector of shifts from current time 
    3538        ivec times; 
    36         ivec obs; 
     39        //! Array of names 
    3740        Array<std::string> names; 
    3841 
    3942private: 
     43        //! auxiliary function used in constructor 
    4044        void init ( ivec in_ids, Array<std::string> in_names, ivec in_sizes, ivec in_times); 
    4145public: 
     
    7882class fnc { 
    7983protected: 
     84        //! Length of the output vector 
    8085        int dimy; 
    8186public: 
     
    8388        virtual vec eval ( const vec &cond ) { 
    8489                return vec ( 0 ); 
    85         }; //Fixme: virtual? 
     90        };  
    8691 
    8792        //! access function 
     
    97102class epdf { 
    98103protected: 
     104        //! Identified of the random variable 
    99105        RV rv; 
    100106public: 
     
    140146        //! Returns the required moment of the epdf 
    141147//      virtual fnc moment ( const int order = 1 ); 
    142         //! Returns a sample from the density conditioned on \c cond, $x \sim epdf(rv|cond)$. \param ll is a return value of log-likelihood of the sample. 
     148        //! Returns a sample 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. 
    143149        virtual vec samplecond ( vec &cond, double &ll ) {this->condition(cond);vec temp= ep->sample();ll=ep->evalpdflog(temp);return temp;}; 
    144  
     150        //! Update \c ep so that it represents this mpdf conditioned on \c rvc = cond 
    145151        virtual void condition ( const vec &cond ) {}; 
    146152         
     153        //! Shortcut for conditioning and evaluation of the internal epdf. In some cases,  this operation can be implemented efficiently. 
    147154        virtual double evalcond (const vec &dt, const vec &cond ) {this->condition(cond);return ep->eval(dt);}; 
    148155 
     
    238245class BMcond { 
    239246protected: 
     247        //! Identificator of the conditioning variable 
    240248        RV rvc; 
    241249public: 
    242250        //! Substitute \c val for \c rvc. 
    243251        virtual void condition ( const vec &val ) =0; 
     252        //! Default constructor 
    244253        BMcond(RV &rv0):rvc(rv0){}; 
     254        //! Destructor for future use 
    245255        virtual ~BMcond(){}; 
     256        //! access function 
    246257        const RV& _rvc() const {return rvc;} 
    247258}; 
  • bdm/stat/libDS.h

    r19 r33  
    3939        void write(vec &ut,ivec &indexes){it_error("MemDS::write is not supported");} 
    4040        void step(); 
     41        //!Default constructor 
    4142        MemDS(mat &Dat, ivec &rowid, ivec &delays); 
    4243}; 
  • bdm/stat/libEF.h

    r32 r33  
    2525//! Global Uniform_RNG 
    2626extern Uniform_RNG UniRNG; 
     27//! Global Normal_RNG 
    2728extern Normal_RNG NorRNG; 
     29//! Global Gamma_RNG 
    2830extern Gamma_RNG GamRNG; 
    29  
    3031 
    3132/*! 
     
    4142        //! default constructor 
    4243        eEF ( const RV &rv ) :epdf ( rv ) {}; 
    43  
     44        //!TODO decide if it is really needed 
    4445        virtual void tupdate ( double phi, mat &vbar, double nubar ) {}; 
    45  
     46        //!TODO decide if it is really needed 
    4647        virtual void dupdate ( mat &v,double nu=1.0 ) {}; 
    4748}; 
    4849 
     50/*! 
     51* \brief Exponential family model. 
     52 
     53* More?... 
     54*/ 
     55 
    4956class mEF : public mpdf { 
    5057 
    5158public: 
     59        //! Default constructor 
    5260        mEF ( const RV &rv0, const RV &rvc0 ) :mpdf ( rv0,rvc0 ) {}; 
    5361}; 
     
    7482public: 
    7583//      enorm() :eEF() {}; 
    76  
     84        //!Default constructor 
    7785        enorm ( RV &rv ); 
     86        //! Set mean value \c mu and covariance \c R 
    7887        void set_parameters ( const vec &mu,const sq_T &R ); 
    7988        //! tupdate in exponential form (not really handy) 
    8089        void tupdate ( double phi, mat &vbar, double nubar ); 
     90        //! dupdate in exponential form (not really handy) 
    8191        void dupdate ( mat &v,double nu=1.0 ); 
    8292 
    8393        vec sample() const; 
     94        //! TODO is it used? 
    8495        mat sample ( int N ) const; 
    8596        double eval ( const vec &val ) const ; 
     
    106117 Multvariate Gamma density as product of independent univariate densities. 
    107118 \f[ 
    108  f(x|a,b) = \prod f(x_i|a_i,b_i) 
     119 f(x|\alpha,\beta) = \prod f(x_i|\alpha_i,\beta_i) 
    109120 \f] 
    110121*/ 
     
    112123class egamma : public eEF { 
    113124protected: 
     125        //! Vector \f$\alpha\f$ 
    114126        vec alpha; 
     127        //! Vector \f$\beta\f$ 
    115128        vec beta; 
    116129public : 
     
    120133        void set_parameters ( const vec &a, const vec &b ) {alpha=a,beta=b;}; 
    121134        vec sample() const; 
     135        //! TODO: is it used anywhere? 
    122136        mat sample ( int N ) const; 
    123137        double evalpdflog ( const vec &val ) const; 
     
    126140        vec mean()const {vec pom(alpha); pom/=beta; return pom;} 
    127141}; 
    128  
     142/* 
    129143//! Weighted mixture of epdfs with external owned components. 
    130144class emix : public epdf { 
     
    140154        vec sample() {it_error ( "Not implemented" );return 0;} 
    141155}; 
     156*/ 
    142157 
    143158//!  Uniform distributed density on a rectangular support 
     
    152167        vec distance; 
    153168//! normalizing coefficients 
    154         double nk,lnk; 
    155 public: 
     169        double nk; 
     170//! cache of log( \c nk ) 
     171        double lnk; 
     172public: 
     173        //! Defualt constructor 
    156174        euni ( const RV rv ) :epdf ( rv ) {} 
    157175        double eval ( const vec &val ) const  {return nk;} 
     
    161179                return low+distance*smp; 
    162180        } 
     181        //! set values of \c low and \c high 
    163182        void set_parameters ( const vec &low0, const vec &high0 ) { 
    164183                distance = high0-low0; 
     
    180199template<class sq_T> 
    181200class mlnorm : public mEF { 
     201        //! Internal epdf that arise by conditioning on \c rvc 
    182202        enorm<sq_T> epdf; 
    183203        vec* _mu; //cached epdf.mu; 
     
    186206        //! Constructor 
    187207        mlnorm ( RV &rv,RV &rvc ); 
     208        //! Set \c A and \c R 
    188209        void set_parameters ( const  mat &A, const sq_T &R ); 
    189210        //!Generate one sample of the posterior 
     
    191212        //!Generate matrix of samples of the posterior 
    192213        mat samplecond ( vec &cond, vec &lik, int n ); 
     214        //! Set value of \c rvc . Result of this operation is stored in \c epdf use function \c _ep to access it. 
    193215        void condition ( vec &cond ); 
    194216}; 
     
    197219\brief  Gamma random walk 
    198220 
    199 Mean value, $\mu$, of this density is given by \c rvc . 
     221Mean value, \f$\mu\f$, of this density is given by \c rvc . 
    200222Standard deviation of the random walk is proportional to one $k$-th the mean. 
    201 This is achieved by setting $\alpha=k$ and $\beta=k/\mu$. 
    202  
    203 The standard deviation of the walk is then: $\mu/\sqrt(k)$. 
     223This is achieved by setting \f$\alpha=k\f$ and \f$\beta=k/\mu\f$. 
     224 
     225The standard deviation of the walk is then: \f$\mu/\sqrt(k)\f$. 
    204226*/ 
    205227class mgamma : public mEF { 
     228        //! Internal epdf that arise by conditioning on \c rvc 
    206229        egamma epdf; 
     230        //! Constant $k$ 
    207231        double k; 
     232        //! cache of epdf.beta 
    208233        vec* _beta; 
    209234 
     
    211236        //! Constructor 
    212237        mgamma ( const RV &rv,const RV &rvc ); 
     238        //! Set value of \c k 
    213239        void set_parameters ( double k ); 
    214240        //!Generate one sample of the posterior 
     
    230256        //! Number of particles 
    231257        int n; 
     258        //! Sample weights $w$ 
    232259        vec w; 
     260        //! Samples \f$x^{(i)}, i=1..n\f$ 
    233261        Array<vec> samples; 
    234262public: 
     263        //! Default constructor 
    235264        eEmp ( const RV &rv0 ,int n0) :epdf ( rv0 ),n(n0),w(n),samples(n) {}; 
     265        //! Set sample 
    236266        void set_parameters ( const vec &w0, epdf* pdf0 ); 
    237267        //! Potentially dangerous, use with care. 
    238268        vec& _w()  {return w;}; 
     269        //! access function 
    239270        Array<vec>& _samples() {return samples;}; 
    240271        //! Function performs resampling, i.e. removal of low-weight samples and duplication of high-weight samples such that the new samples represent the same density. 
    241272        ivec resample ( RESAMPLING_METHOD method = SYSTEMATIC ); 
     273        //! inherited operation : NOT implemneted 
    242274        vec sample() const {it_error ( "Not implemented" );return 0;} 
     275        //! inherited operation : NOT implemneted 
    243276        double evalpdflog(const vec &val) const {it_error ( "Not implemented" );return 0.0;} 
    244277        vec mean()const {vec pom=zeros(rv.count());  
  • bdm/stat/libFN.cpp

    r22 r33  
    33 
    44using std::endl; 
    5  
    6 linfn linfn::evalsome ( ivec &rvind ) 
    7 { 
    8         return *this; 
    9 } 
    105 
    116bilinfn::bilinfn ( const RV &rvx0, const RV &rvu0, const mat &A0, const mat &B0 ) : diffbifn ( rvx0,rvu0 ) 
  • bdm/stat/libFN.h

    r22 r33  
    1010// 
    1111// 
     12#ifndef FN_H 
     13#define FN_H 
     14 
    1215#include <itpp/itbase.h> 
    1316#include "libBM.h" 
     
    1821class constfn : public fnc 
    1922{ 
    20                 RV rv; 
     23                //! value of the function 
    2124                vec val; 
    2225 
    23  
    2426        public: 
    25                 vec eval() {return val;}; 
    26                 vec eval ( vec &cond ) {return val;}; 
     27                //vec eval() {return val;}; 
     28                //! inherited 
     29                vec eval ( const vec &cond ) {return val;}; 
    2730                //!Default constructor 
    28                 constfn ( const vec &val0 ) :rv(),val ( val0 ) {}; 
     31                constfn ( const vec &val0 ) :val ( val0 ) {}; 
    2932}; 
    3033 
     
    3235class linfn: public fnc 
    3336{ 
     37                //! Identification of $x$ 
    3438                RV rv; 
    35                 ivec indexlist; // needed by evalsome 
     39                //! Matrix A 
    3640                mat A; 
     41                //! Matrix B 
    3742                vec B; 
    3843        public : 
    39                 vec eval ( vec &cond ) {it_assert_debug ( cond.length() ==rv.count(), "linfn::eval Wrong cond." );return A*cond+B;}; 
     44                vec eval (const vec &cond ) {it_assert_debug ( cond.length() ==rv.count(), "linfn::eval Wrong cond." );return A*cond+B;}; 
    4045 
    41                 linfn evalsome ( ivec &rvind ); 
    42                 linfn ( const RV &rv0 ) :rv ( rv0 ),A ( eye ( rv0.count() ) ),B ( zeros ( rv0.count() ) ) { indexlist=rv.indexlist();}; 
    43                 linfn ( const RV &rv0, const mat &A0 ) : rv ( rv0 ), A ( A0 ), B ( zeros ( rv0.count() ) ) { indexlist=rv.indexlist();}; 
    44                 linfn ( const RV &rv0, const mat &A0, const vec &B0 ) :rv ( rv0 ), A ( A0 ), B ( B0 ) { indexlist=rv.indexlist();}; 
     46//              linfn evalsome ( ivec &rvind ); 
     47                //!default constructor 
     48                linfn ( const RV &rv0 ) :rv ( rv0 ),A ( eye ( rv0.count() ) ),B ( zeros ( rv0.count() ) ) { }; 
     49                //! Set values of \c A and \c B 
     50                void set_parameters ( const mat &A0 , const vec &B0 ) {A=A0; B=B0;}; 
    4551}; 
    4652 
     
    5864{ 
    5965        protected: 
    60                 RV rvx,rvu; 
     66                //! Indentifier of the first rv. 
     67                RV rvx; 
     68                //! Indentifier of the second rv. 
     69                RV rvu; 
     70                //! cache for rvx.count() 
    6171                int dimx; 
     72                //! cache for rvu.count() 
    6273                int dimu; 
    6374        public: 
     
    7182                //! Evaluates $f(x0,u0)$ 
    7283                virtual vec eval ( const vec &x0, const vec &u0 ) {return zeros ( dimy );}; 
    73                 //! Evaluates \f$A=\frac{d}{dx}f(x,u)|_{x0,u0}\f$ and writes result into \c A . @param full denotes that even unchanged entries are to be rewritten. When, false only the changed elements are computed. 
     84                //! Evaluates \f$A=\frac{d}{dx}f(x,u)|_{x0,u0}\f$ and writes result into \c A . @param full denotes that even unchanged entries are to be rewritten. When, false only the changed elements are computed. @param x0 numeric value of $x$, @param u0 numeric value of $u$ @param A a place where the result will be stored. 
    7485                virtual void dfdx_cond ( const vec &x0, const vec &u0, mat &A , bool full=true ) {}; 
    75                 //! Evaluates \f$A=\frac{d}{du}f(x,u)|_{x0,u0}\f$ and writes result into \c A . @param full denotes that even unchanged entries are to be rewritten. When, false only the changed elements are computed. 
     86                //! Evaluates \f$A=\frac{d}{du}f(x,u)|_{x0,u0}\f$ and writes result into \c A . @param full denotes that even unchanged entries are to be rewritten. When, false only the changed elements are computed.        @param x0 numeric value of $x$, @param u0 numeric value of $u$ @param A a place where the result will be stored. 
    7687                virtual void dfdu_cond ( const vec &x0, const vec &u0, mat &A, bool full=true ) {}; 
    7788                //!Default constructor (dimy is not set!) 
     
    94105                //! Default constructor 
    95106                bilinfn ( const RV &rvx0, const RV &rvu0 ) : diffbifn ( rvx0,rvu0 ) ,A ( eye ( dimx ) ),B ( zeros ( dimx,dimu ) )       {}; 
    96                 // 
     107                //! Alternative constructor 
    97108                bilinfn ( const RV &rvx0, const RV &rvu0, const mat &A0, const mat &B0 ); 
    98                 // 
     109                //! 
    99110                void dfdx_cond ( const vec &x0, const vec &u0, mat &F, bool full ) 
    100111                { 
     
    102113                        if ( full ) F=A;        //else : nothing has changed no need to regenerate 
    103114                } 
    104                 // 
     115                //! 
    105116                void dfdu_cond ( const vec &x0, const vec &u0, mat &F,  bool full=true ) 
    106117                { 
     
    109120                } 
    110121}; 
     122 
     123#endif // FN_H