[97] | 1 | /*! |
---|
| 2 | \file |
---|
| 3 | \brief Bayesian Filtering for generalized autoregressive (ARX) model |
---|
| 4 | \author Vaclav Smidl. |
---|
| 5 | |
---|
| 6 | ----------------------------------- |
---|
| 7 | BDM++ - C++ library for Bayesian Decision Making under Uncertainty |
---|
| 8 | |
---|
| 9 | Using IT++ for numerical operations |
---|
| 10 | ----------------------------------- |
---|
| 11 | */ |
---|
| 12 | |
---|
| 13 | #ifndef AR_H |
---|
| 14 | #define AR_H |
---|
| 15 | |
---|
[384] | 16 | #include "../math/functions.h" |
---|
| 17 | #include "../stat/exp_family.h" |
---|
| 18 | #include "../base/user_info.h" |
---|
[97] | 19 | |
---|
[270] | 20 | namespace bdm { |
---|
[97] | 21 | |
---|
| 22 | /*! |
---|
| 23 | * \brief Linear Autoregressive model with Gaussian noise |
---|
| 24 | |
---|
| 25 | Regression of the following kind: |
---|
| 26 | \f[ |
---|
| 27 | y_t = \theta_1 \psi_1 + \theta_2 + \psi_2 +\ldots + \theta_n \psi_n + r e_t |
---|
| 28 | \f] |
---|
| 29 | 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: |
---|
| 30 | \f[ |
---|
| 31 | e_t \sim \mathcal{N}(0,1). |
---|
| 32 | \f] |
---|
| 33 | |
---|
[271] | 34 | See \ref tut_arx for mathematical treatment. |
---|
| 35 | |
---|
| 36 | The easiest way how to use the class is: |
---|
| 37 | \include arx_simple.cpp |
---|
| 38 | |
---|
[384] | 39 | \todo sort out constant terms - bayes should accept vec without additional 1s |
---|
[97] | 40 | */ |
---|
[170] | 41 | class ARX: public BMEF { |
---|
[97] | 42 | protected: |
---|
[270] | 43 | //!size of output variable (needed in regressors) |
---|
| 44 | int dimx; |
---|
| 45 | //!description of modelled data \f$ y_t \f$ in the likelihood function |
---|
| 46 | //! Do NOT access directly, only via \c get_yrv(). |
---|
| 47 | RV _yrv; |
---|
[97] | 48 | //! Posterior estimate of \f$\theta,r\f$ in the form of Normal-inverse Wishart density |
---|
| 49 | egiw est; |
---|
| 50 | //! cached value of est.V |
---|
| 51 | ldmat &V; |
---|
| 52 | //! cached value of est.nu |
---|
| 53 | double ν |
---|
| 54 | public: |
---|
[270] | 55 | //! \name Constructors |
---|
| 56 | //!@{ |
---|
[477] | 57 | ARX ( const double frg0 = 1.0 ) : BMEF ( frg0 ), est (), V ( est._V() ), nu ( est._nu() ) {}; |
---|
| 58 | ARX ( const ARX &A0 ) : BMEF (), est (), V ( est._V() ), nu ( est._nu() ) { |
---|
| 59 | set_statistics ( A0.dimx, A0.V, A0.nu ); |
---|
| 60 | set_parameters ( A0.frg ); |
---|
[286] | 61 | }; |
---|
[283] | 62 | ARX* _copy_() const; |
---|
[477] | 63 | void set_parameters ( double frg0 ) { |
---|
| 64 | frg = frg0; |
---|
| 65 | } |
---|
| 66 | void set_statistics ( int dimx0, const ldmat V0, double nu0 = -1.0 ) { |
---|
| 67 | est.set_parameters ( dimx0, V0, nu0 ); |
---|
| 68 | last_lognc = est.lognc(); |
---|
| 69 | dimx = dimx0; |
---|
| 70 | } |
---|
[270] | 71 | //!@} |
---|
[170] | 72 | |
---|
| 73 | // //! Set parameters given by moments, \c mu (mean of theta), \c R (mean of R) and \c C (variance of theta) |
---|
| 74 | // void set_parameters ( const vec &mu, const mat &R, const mat &C, double dfm){}; |
---|
[145] | 75 | //! Set sufficient statistics |
---|
[170] | 76 | void set_statistics ( const BMEF* BM0 ); |
---|
[270] | 77 | // //! Returns sufficient statistics |
---|
| 78 | // void get_parameters ( mat &V0, double &nu0 ) {V0=est._V().to_mat(); nu0=est._nu();} |
---|
| 79 | //!\name Mathematical operations |
---|
| 80 | //!@{ |
---|
| 81 | |
---|
| 82 | //! Weighted Bayes \f$ dt = [y_t psi_t] \f$. |
---|
[170] | 83 | void bayes ( const vec &dt, const double w ); |
---|
[477] | 84 | void bayes ( const vec &dt ) { |
---|
| 85 | bayes ( dt, 1.0 ); |
---|
| 86 | }; |
---|
[170] | 87 | double logpred ( const vec &dt ) const; |
---|
[270] | 88 | void flatten ( const BMEF* B ) { |
---|
[477] | 89 | const ARX* A = dynamic_cast<const ARX*> ( B ); |
---|
[170] | 90 | // nu should be equal to B.nu |
---|
[477] | 91 | est.pow ( A->nu / nu ); |
---|
| 92 | if ( evalll ) { |
---|
| 93 | last_lognc = est.lognc(); |
---|
| 94 | } |
---|
[170] | 95 | } |
---|
[262] | 96 | //! Conditioned version of the predictor |
---|
[270] | 97 | enorm<ldmat>* epredictor ( const vec &rgr ) const; |
---|
| 98 | //! Predictor for empty regressor |
---|
[286] | 99 | enorm<ldmat>* epredictor() const { |
---|
[565] | 100 | bdm_assert_debug ( dimx == V.rows() - 1, "Regressor is not only 1" ); |
---|
[286] | 101 | return epredictor ( vec_1 ( 1.0 ) ); |
---|
| 102 | } |
---|
[262] | 103 | //! conditional version of the predictor |
---|
[270] | 104 | mlnorm<ldmat>* predictor() const; |
---|
| 105 | mlstudent* predictor_student() const; |
---|
[97] | 106 | //! Brute force structure estimation.\return indeces of accepted regressors. |
---|
[170] | 107 | ivec structure_est ( egiw Eg0 ); |
---|
[577] | 108 | //! Smarter structure estimation by Ludvik Tesar.\return indeces of accepted regressors. |
---|
| 109 | ivec structure_est_LT ( egiw Eg0 ); |
---|
[270] | 110 | //!@} |
---|
| 111 | |
---|
| 112 | //!\name Access attributes |
---|
| 113 | //!@{ |
---|
[477] | 114 | const egiw& posterior() const { |
---|
| 115 | return est; |
---|
| 116 | } |
---|
[270] | 117 | //!@} |
---|
| 118 | |
---|
| 119 | //!\name Connection |
---|
| 120 | //!@{ |
---|
[477] | 121 | void set_drv ( const RV &drv0 ) { |
---|
| 122 | drv = drv0; |
---|
| 123 | } |
---|
[565] | 124 | |
---|
[270] | 125 | RV& get_yrv() { |
---|
| 126 | //if yrv is not ready create it |
---|
[477] | 127 | if ( _yrv._dsize() != dimx ) { |
---|
| 128 | int i = 0; |
---|
| 129 | while ( _yrv._dsize() < dimx ) { |
---|
| 130 | _yrv.add ( drv ( vec_1 ( i ) ) ); |
---|
| 131 | i++; |
---|
| 132 | } |
---|
[270] | 133 | } |
---|
| 134 | //yrv should be ready by now |
---|
[565] | 135 | bdm_assert_debug ( _yrv._dsize() == dimx, "incompatible drv" ); |
---|
[270] | 136 | return _yrv; |
---|
| 137 | } |
---|
| 138 | //!@} |
---|
[357] | 139 | |
---|
| 140 | // TODO dokumentace - aktualizovat |
---|
| 141 | /*! UI for ARX estimator |
---|
| 142 | |
---|
| 143 | The ARX is constructed from a structure with fields: |
---|
| 144 | \code |
---|
| 145 | estimator = { |
---|
| 146 | type = "ARX"; |
---|
| 147 | y = {type="rv", ...} // description of output variables |
---|
| 148 | rgr = {type="rv", ...} // description of regressor variables |
---|
| 149 | constant = true; // boolean switch if the constant term is modelled or not |
---|
| 150 | |
---|
| 151 | //optional fields |
---|
[477] | 152 | dV0 = [1e-3, 1e-5, 1e-5, 1e-5]; |
---|
[357] | 153 | // default: 1e-3 for y, 1e-5 for rgr |
---|
| 154 | nu0 = 6; // default: rgrlen + 2 |
---|
| 155 | frg = 1.0; // forgetting, default frg=1.0 |
---|
| 156 | }; |
---|
| 157 | \endcode |
---|
| 158 | |
---|
| 159 | The estimator will assign names of the posterior in the form ["theta_i" and "r_i"] |
---|
| 160 | */ |
---|
[477] | 161 | void from_setting ( const Setting &set ); |
---|
[357] | 162 | |
---|
[97] | 163 | }; |
---|
| 164 | |
---|
[477] | 165 | UIREGISTER ( ARX ); |
---|
[529] | 166 | SHAREDPTR ( ARX ); |
---|
[357] | 167 | |
---|
[254] | 168 | } |
---|
[97] | 169 | |
---|
| 170 | #endif // AR_H |
---|
| 171 | |
---|