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 | |
---|
16 | #include "../stat/libFN.h" |
---|
17 | #include "../stat/libEF.h" |
---|
18 | #include "../user_info.h" |
---|
19 | |
---|
20 | namespace bdm { |
---|
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 | |
---|
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 | |
---|
39 | */ |
---|
40 | class ARX: public BMEF { |
---|
41 | protected: |
---|
42 | //!size of output variable (needed in regressors) |
---|
43 | int dimx; |
---|
44 | //!description of modelled data \f$ y_t \f$ in the likelihood function |
---|
45 | //! Do NOT access directly, only via \c get_yrv(). |
---|
46 | RV _yrv; |
---|
47 | //! Posterior estimate of \f$\theta,r\f$ in the form of Normal-inverse Wishart density |
---|
48 | egiw est; |
---|
49 | //! cached value of est.V |
---|
50 | ldmat &V; |
---|
51 | //! cached value of est.nu |
---|
52 | double ν |
---|
53 | public: |
---|
54 | //! \name Constructors |
---|
55 | //!@{ |
---|
56 | ARX ( const double frg0=1.0 ) : BMEF ( frg0 ),est (), V ( est._V() ), nu ( est._nu() ) {}; |
---|
57 | ARX ( const ARX &A0 ) : BMEF (),est (), V ( est._V() ), nu ( est._nu() ) { |
---|
58 | set_statistics ( A0.dimx,A0.V,A0.nu ); |
---|
59 | set_parameters(A0.frg); |
---|
60 | }; |
---|
61 | ARX* _copy_() const; |
---|
62 | void set_parameters ( double frg0 ) {frg=frg0;} |
---|
63 | void set_statistics ( int dimx0, const ldmat V0, double nu0=-1.0 ) {est.set_parameters ( dimx0,V0,nu0 );last_lognc=est.lognc();dimx=dimx0;} |
---|
64 | //!@} |
---|
65 | |
---|
66 | // //! Set parameters given by moments, \c mu (mean of theta), \c R (mean of R) and \c C (variance of theta) |
---|
67 | // void set_parameters ( const vec &mu, const mat &R, const mat &C, double dfm){}; |
---|
68 | //! Set sufficient statistics |
---|
69 | void set_statistics ( const BMEF* BM0 ); |
---|
70 | // //! Returns sufficient statistics |
---|
71 | // void get_parameters ( mat &V0, double &nu0 ) {V0=est._V().to_mat(); nu0=est._nu();} |
---|
72 | //!\name Mathematical operations |
---|
73 | //!@{ |
---|
74 | |
---|
75 | //! Weighted Bayes \f$ dt = [y_t psi_t] \f$. |
---|
76 | void bayes ( const vec &dt, const double w ); |
---|
77 | void bayes ( const vec &dt ) {bayes ( dt,1.0 );}; |
---|
78 | double logpred ( const vec &dt ) const; |
---|
79 | void flatten ( const BMEF* B ) { |
---|
80 | const ARX* A=dynamic_cast<const ARX*> ( B ); |
---|
81 | // nu should be equal to B.nu |
---|
82 | est.pow ( A->nu/nu ); |
---|
83 | if ( evalll ) {last_lognc=est.lognc();} |
---|
84 | } |
---|
85 | //! Conditioned version of the predictor |
---|
86 | enorm<ldmat>* epredictor ( const vec &rgr ) const; |
---|
87 | //! Predictor for empty regressor |
---|
88 | enorm<ldmat>* epredictor() const { |
---|
89 | it_assert_debug ( dimx==V.rows()-1,"Regressor is not only 1" ); |
---|
90 | return epredictor ( vec_1 ( 1.0 ) ); |
---|
91 | } |
---|
92 | //! conditional version of the predictor |
---|
93 | mlnorm<ldmat>* predictor() const; |
---|
94 | mlstudent* predictor_student() const; |
---|
95 | //! Brute force structure estimation.\return indeces of accepted regressors. |
---|
96 | ivec structure_est ( egiw Eg0 ); |
---|
97 | //!@} |
---|
98 | |
---|
99 | //!\name Access attributes |
---|
100 | //!@{ |
---|
101 | const egiw* _e() const {return &est ;}; |
---|
102 | const egiw& posterior() const {return est;} |
---|
103 | //!@} |
---|
104 | |
---|
105 | //!\name Connection |
---|
106 | //!@{ |
---|
107 | void set_drv ( const RV &drv0 ) {drv=drv0;} |
---|
108 | RV& get_yrv() { |
---|
109 | //if yrv is not ready create it |
---|
110 | if ( _yrv._dsize() !=dimx ) { |
---|
111 | int i=0; |
---|
112 | while ( _yrv._dsize() <dimx ) {_yrv.add ( drv ( vec_1 ( i ) ) );i++;} |
---|
113 | } |
---|
114 | //yrv should be ready by now |
---|
115 | it_assert_debug ( _yrv._dsize() ==dimx,"incompatible drv" ); |
---|
116 | return _yrv; |
---|
117 | } |
---|
118 | //!@} |
---|
119 | |
---|
120 | // TODO dokumentace - aktualizovat |
---|
121 | /*! UI for ARX estimator |
---|
122 | |
---|
123 | The ARX is constructed from a structure with fields: |
---|
124 | \code |
---|
125 | estimator = { |
---|
126 | type = "ARX"; |
---|
127 | y = {type="rv", ...} // description of output variables |
---|
128 | rgr = {type="rv", ...} // description of regressor variables |
---|
129 | constant = true; // boolean switch if the constant term is modelled or not |
---|
130 | |
---|
131 | //optional fields |
---|
132 | dV0 = [1e-3, 1e-5, 1e-5, 1e-5]; |
---|
133 | // default: 1e-3 for y, 1e-5 for rgr |
---|
134 | nu0 = 6; // default: rgrlen + 2 |
---|
135 | frg = 1.0; // forgetting, default frg=1.0 |
---|
136 | }; |
---|
137 | \endcode |
---|
138 | |
---|
139 | The estimator will assign names of the posterior in the form ["theta_i" and "r_i"] |
---|
140 | */ |
---|
141 | void from_setting( const Setting &root ); |
---|
142 | |
---|
143 | // TODO dodelat void to_setting( Setting &root ) const; |
---|
144 | }; |
---|
145 | |
---|
146 | UIREGISTER(ARX); |
---|
147 | |
---|
148 | } |
---|
149 | |
---|
150 | #endif // AR_H |
---|
151 | |
---|