root/library/bdm/estim/arx.h @ 1036

/*!
  \file
  \brief Bayesian Filtering for generalized autoregressive (ARX) model
  \author Vaclav Smidl.

  -----------------------------------
  BDM++ - C++ library for Bayesian Decision Making under Uncertainty

  Using IT++ for numerical operations
  -----------------------------------
*/

#ifndef AR_H
#define AR_H

#include "../math/functions.h"
#include "../stat/exp_family.h"
#include "../base/user_info.h"
//#include "../estim/kalman.h"
#include "arx_straux.h"

namespace bdm {

/*!
* \brief Linear Autoregressive model with Gaussian noise

Regression of the following kind:
\f[
y_t = \theta_1 \psi_1 + \theta_2 + \psi_2 +\ldots + \theta_n \psi_n + r e_t
\f]
where unknown parameters \c rv are \f$[\theta r]\f$, regression vector \f$\psi=\psi(y_{1:t},u_{1:t})\f$ is a known function of past outputs and exogeneous variables \f$u_t\f$. Distrubances \f$e_t\f$ are supposed to be normally distributed:
\f[
e_t \sim \mathcal{N}(0,1).
\f]

See \ref tut_arx for mathematical treatment.

The easiest way how to use the class is:
\include arx_simple.cpp

        \todo sort out constant terms - bayes should accept vec without additional 1s
*/
class ARX: public BMEF {
protected:
        //! switch if constant is modelled or not
        bool have_constant;
        //! vector of dyadic update
        vec dyad;
        //! RV of regressor
        RV rgr;
        //! length of the regressor (without optional constant)
        int rgrlen;
        //! posterior density
        egiw est;
        //! Alternative estimate of parameters, used in stabilized forgetting, see [Kulhavy]
        egiw alter_est;
public:
        //! \name Constructors
        //!@{
        ARX ( const double frg0 = 1.0 ) : BMEF ( frg0 ),  have_constant ( true ), dyad(), rgrlen(),est(), alter_est() {};
        ARX ( const ARX &A0 ) : BMEF ( A0 ),  have_constant ( A0.have_constant ), dyad ( A0.dyad ),rgrlen(A0.rgrlen), est ( A0.est ), alter_est ( A0.alter_est ) { };

        ARX* _copy() const;

        void set_frg ( double frg0 ) {
                frg = frg0;
        }
        void set_constant ( bool const0 ) {
                have_constant = const0;
        }
        void set_statistics ( int dimy0, const ldmat V0, double nu0 = -1.0 ) {
                est.set_parameters ( dimy0, V0, nu0 );
                est.validate();
                last_lognc = est.lognc();
                dimy = dimy0;
        }
        //!@}

        //! Set sufficient statistics
        void set_statistics ( const BMEF* BM0 );

        //!\name Mathematical operations
        //!@{

        //! Weighted Bayes \f$ dt = [y_t psi_t] \f$.
        void bayes_weighted ( const vec &yt, const vec &cond = empty_vec, const double w = 1.0 );
        void bayes ( const vec &yt, const vec &cond = empty_vec ) {
                bayes_weighted ( yt, cond, 1.0 );
        };
        double logpred ( const vec &yt, const vec &cond ) const;
        void flatten ( const BMEF* B , double weight );
        //! Conditioned version of the predictor
        enorm<ldmat>* epredictor ( const vec &rgr ) const;
        //! conditional version of the predictor
        template<class sq_T>
        shared_ptr<mlnorm<sq_T> > ml_predictor() const;
        //! fast version of predicto
        template<class sq_T>
        void ml_predictor_update ( mlnorm<sq_T> &pred ) const;
        mlstudent* predictor_student() const;
        //! Brute force structure estimation.\return indices of accepted regressors.
        ivec structure_est ( const egiw &Eg0 );
        //! Smarter structure estimation by Ludvik Tesar.\return indices of accepted regressors.
        ivec structure_est_LT ( const egiw &Eg0 );
        //! reduce structure to the given ivec of matrix V
        void reduce_structure(ivec &inds_in_V){
                ldmat V = posterior()._V();
                if (max(inds_in_V)>=V.rows()) {bdm_error("Incompatible structure");}
                
                ldmat newV(V,inds_in_V);
                est.set_parameters(dimy,newV, posterior()._nu());
                
                if (have_constant){
                        ivec rgr_elem= find(inds_in_V<(V.rows()-1)); // < -- find non-constant 
                        rgr = rgr.subselect(rgr_elem);
                        rgrlen = rgr_elem.length();
                } else{
                        rgr = rgr.subselect(inds_in_V);
                }
                validate();
        }
        //!@}

        //!\name Access attributes
        //!@{
        //! return correctly typed posterior (covariant return)
        const egiw& posterior() const {
                return est;
        }
        //!@}

        /*! UI for ARX estimator

        \code
        class = 'ARX';
        yrv   = RV({names_of_dt} )                 // description of output variables
        rgr   = RV({names_of_regressors}, [-1,-2]} // description of regressor variables
        constant = 1;                              // 0/1 switch if the constant term is modelled or not

        --- optional ---
        prior = {class='egiw',...};                // Prior density, when given default is used instead
        alternative = {class='egiw',...};          // Alternative density in stabilized estimation, when not given prior is used

        frg = 1.0;                                 // forgetting, default frg=1.0

        rv  = RV({names_of_parameters}}            // description of parametetr names
                                                                                           // default: [""]
        \endcode
        */
        void from_setting ( const Setting &set );

        void validate() {
                BMEF::validate();       
                est.validate();
                
                // When statistics is defined, it has priority
                if(posterior()._dimx()>0) {//statistics is assigned
                        dimy = posterior()._dimx();
                        rgrlen=posterior()._V().rows() - dimy - int ( have_constant == true );
                        dimc = rgrlen;
                } else{  // statistics is not assigned - build it from dimy and rgrlen
                        bdm_assert(dimy>0,"No way to validate egiw: empty statistics and empty dimy");
                        est.set_parameters(dimy, zeros(dimy+rgrlen+int(have_constant==true)));
                        set_prior_default(est);
                }
                if (alter_est.dimension()==0) alter_est=est;

                dyad = ones ( est._V().rows() );
        }
        //! function sets prior and alternative density
        void set_prior ( const epdf *prior );
        //! build default prior and alternative when all values are set
        void set_prior_default ( egiw &prior ) {
                //assume
                vec dV0 ( prior._V().rows() );
                dV0.set_subvector ( 0, prior._dimx() - 1, 1.0 );
                if (dV0.length()>prior._dimx())
                        dV0.set_subvector ( prior._dimx(), dV0.length() - 1, 1e-5 );
                
                prior.set_parameters ( prior._dimx(), ldmat ( dV0 ) );
                prior.validate();
        }

        void to_setting ( Setting &set ) const
        {                       
                BMEF::to_setting( set ); // takes care of rv, yrv, rvc
                UI::save(rgr, set, "rgr");
                int constant = have_constant ? 1 : 0;
                UI::save(constant, set, "constant");
                UI::save(&alter_est, set, "alternative");
                UI::save(&posterior(), set, "posterior");
                
        } 
        //! access function
        RV & _rgr() {return rgr;}
        bool _have_constant() {return have_constant;}
        int _rgrlen() {return rgrlen;}
};
UIREGISTER ( ARX );
SHAREDPTR ( ARX );

/*! ARX model conditined by knowledge of the forgetting factor
\f[ f(\theta| d_1 \ldots d_t , \phi_t) \f]

The symbol \f$ \phi \f$ is assumed to be the last of the conditioning variables.
*/
class ARXfrg : public ARX {
public:
        ARXfrg() : ARX() {};
        //! copy constructor
        ARXfrg ( const ARXfrg &A0 ) : ARX ( A0 ) {};
        virtual ARXfrg* _copy() const {
                ARXfrg *A = new ARXfrg ( *this );
                return A;
        }

        void bayes ( const vec &val, const vec &cond ) {
                bdm_assert_debug(cond.size()>rgrlen, "ARXfrg: Insufficient conditioning, frg not given.");
                frg = cond ( rgrlen); // the first part after rgrlen
                ARX::bayes ( val, cond.left(rgrlen) );
        }
        void validate() {
                ARX::validate();
                rvc.add ( RV ( "{phi }", vec_1 ( 1 ) ) );
                dimc += 1;
        }
};
UIREGISTER ( ARXfrg );


/*!
* \brief ARX with partial forgetting

Implements partial forgetting when <tt>bayes(const vec &yt, const vec &cond=empty_vec)</tt> is called, where \c cond is a vector <em>(regressor', forg.factor')</em>.

Note, that the weights have the same order as hypotheses in partial forgetting, and follow this scheme:
\li 0 means that the parameter doesn't change,
\li 1 means that the parameter varies.

The combinations of parameters are described binary:
\f{bmatrix}[
  0 & 0 & 0 & \ldots \\
  1 & 0 & 0 & \ldots \\
  0 & 1 & 0 & \ldots \\
  1 & 1 & 0 & \ldots \\
  \vdots & \vdots & \vdots & \vdots 
\f{bmatrix}]
Notice, that each n-th column has altering n-tuples of 1's and 0's, n = 0,...,number of params. Hence, the first forg. factor relates to the situation when no parameter changes, the second one when the first parameter changes etc.

See ARX class for more information about ARX.
*/
class ARXpartialforg : public ARX {
public:
        ARXpartialforg() : ARX(1.0) {};
        //! copy constructor
        ARXpartialforg ( const ARXpartialforg &A0 ) : ARX ( A0 ) {};
        virtual ARXpartialforg* _copy() const {
                ARXpartialforg *A = new ARXpartialforg ( *this );
                return A;
        }

    void bayes ( const vec &val, const vec &cond ); 
        void validate() {
                ARX::validate();
                int philen = 1 << (est._V().cols() - 1);
                rvc.add ( RV ( "{phi }", vec_1(philen) ) ); // pocet 2^parametru
                dimc += philen;
        }
};
UIREGISTER ( ARXpartialforg );


////////////////////
template<class sq_T>
shared_ptr< mlnorm<sq_T> > ARX::ml_predictor ( ) const {
        shared_ptr< mlnorm<sq_T> > tmp = new mlnorm<sq_T> ( );
        tmp->set_rv ( yrv );
        tmp->set_rvc ( _rvc() );

        ml_predictor_update ( *tmp );
        tmp->validate();
        return tmp;
}

template<class sq_T>
void ARX::ml_predictor_update ( mlnorm<sq_T> &pred ) const {
        mat mu ( dimy, posterior()._V().rows() - dimy );
        mat R ( dimy, dimy );

        est.mean_mat ( mu, R ); //mu =
        mu = mu.T();
        //correction for student-t  -- TODO check if correct!!

        if ( have_constant ) { // constant term
                //Assume the constant term is the last one:
                pred.set_parameters ( mu.get_cols ( 0, mu.cols() - 2 ), mu.get_col ( mu.cols() - 1 ), sq_T ( R ) );
        } else {
                pred.set_parameters ( mu, zeros ( dimy ), sq_T ( R ) );
        }
}

};
#endif // AR_H