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

/*!
  \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 =     heta_1 \psi_1 +     heta_2 + \psi_2 +\ldots +     heta_n \psi_n + r e_t
\f]
where unknown parameters \c rv are \f$[    heta 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

            odo 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;
        vec samplepred (  const vec &cond ) ;
        void flatten ( const BMEF* B , double weight );
    //! Conditioned version of the predictor
    enorm<ldmat>* epredictor ( const vec &rgr ) const;
        estudent<ldmat>* epredictor_student ( 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() 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;
    }
    //!@}

    /*! Create object from the following structure

    \code
    class = 'ARX';
    rgr = RV({'names',...},[sizes,...],[times,...]);   % description of regressor variables
    --- optional fields ---    
    prior = configuration of bdm::egiw;                % any offspring of eqiw for prior density, bdm::egiw::from_setting
    alternative = configuration of bdm::egiw;          % any offspring of eqiw for alternative density in stabilized estimation of prior density    
    constant = [];                                     % 0/1 switch if the constant term is modelled or not
    --- inherited fields ---
    bdm::BMEF::from_setting
    \endcode
    If the optional fields are not given, they will be filled as follows:
    \code
    prior = posterior;                                % when prior is not given the posterior is used (TODO it is unclear)
    alternative = prior;                              % when alternative is not given the prior is used
    constant = 1;                                     % constant term is modelled on 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 );

//! \brief ARX moidel with parameters in LS form
class ARXls : public BMEF{
  public:
        egw_ls<ldmat> est;
        
        egw_ls<ldmat> alternative;
        
        const egw_ls<ldmat>& posterior() {return est;};
        
        void bayes(const vec &dt, const vec &psi){
//        ldmat &Pbeta = est.P;
//        ldmat &Palpha = alternative.P;
          
          
        }
};

/*! \brief ARX model conditined by knowledge of the forgetting factor
\f[ f(    heta| 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 model with fixed maxent forgetting on increments,
 * \f[ f(    heta| d_1 \ldots d_t , \phi_t) \f]
 * 
 * The symbol \f$ \phi \f$ is not interpreted as exponentila forgetting but foreggting of incoming data!!
 */
class ARXmaxent : public ARX {
protected:
        double maxent_frg;
public:
        ARXmaxent() : ARX() {};
        //! copy constructor
        ARXmaxent ( const ARXmaxent &A0 ) : ARX ( A0 ),maxent_frg(A0.maxent_frg) {};
        virtual ARXmaxent* _copy() const {
                ARXmaxent *A = new ARXmaxent ( *this );
                return A;
        }
        
        void bayes ( const vec &val, const vec &cond ) {
                ARX::bayes_weighted ( val, cond, maxent_frg );
        }
        void from_setting(const Setting &set){
                ARX::from_setting(set);
                maxent_frg=frg;
                frg = 1.0;
        }
};
        
UIREGISTER ( ARXmaxent );

/*! \brief ARX model with fixed-length window - old entries are removed
 * \f[ f(    heta| d_1 \ldots d_t) \f]
 * 
 */
class ARXwin : public ARX {
protected:
        mat Y;
        mat Cond;
        
        int win_length;
        ldmat V0;
        double nu0;
public:
        ARXwin() : ARX() {};
        //! copy constructor
        void set_parameters(const int win_length0){ 
                win_length=win_length0;
        }
        ARXwin ( const ARXwin &A0 ) : ARX(A0), Y( A0.Y), Cond(A0.Cond) {};
        
        virtual ARXwin* _copy() const {
                ARXwin *A = new ARXwin ( *this );
                return A;
        }
        
        void bayes ( const vec &val, const vec &cond ) {
                
                if(cond.size()>0)
                {
                        // fill window
                        Y.append_col(val);
                        Cond.append_col(cond);          
                        if (Y.cols()>win_length){
                                // shift the buffer
                                Y=Y.get_cols(1,Y.cols()-1);
                                Cond=Cond.get_cols(1,Cond.cols()-1);
                        }
                        
                        est._V()=V0;
                        est._nu()=nu0;
                        for ( int t = 0; t < Y.cols(); t++ ) {
                                ARX::bayes ( Y.get_col ( t ), Cond.get_col ( t ) );
                        }
                }
                else
                {
                        Y.append_col(val);

                        if (Y.cols()>win_length){
                                // shift the buffer
                                Y=Y.get_cols(1,Y.cols()-1);                             
                        }

                        est._V()=V0;
                        est._nu()=nu0;

                        for ( int t = 0; t < Y.cols(); t++ ) {
                                ARX::bayes (Y.get_col ( t ));
                        }
                }
                
        }

        void from_setting(const Setting &set){
                ARX::from_setting(set);
                UI::get(win_length,set,"win_length",UI::compulsory);
        }
        void validate() {
                ARX::validate();
                V0=est._V();
                nu0=est._nu();
        }
};
        
UIREGISTER ( ARXwin );
        


/*! \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