| 80 | | <p>The object <a class="el" href="classbdm_1_1BM.html" title="Bayesian Model of a system, i.e. all uncertainty is modeled by probabilities.">bdm::BM</a> is basic software image of the Bayes rule: </p> |
| 81 | | <p class="formulaDsp"> |
| 82 | | <img class="formulaDsp" alt="\[ f(x_t|d_1\ldots d_t) \propto f(d_t|x_t,d_1\ldots d_{t-1}) f(x_t| d_1\ldots d_{t-1}) \]" src="form_174.png"/> |
| 83 | | </p> |
| 84 | | <p>Since this operation can not be defined universally, the object is defined as abstract class with methods for:</p> |
| 85 | | <ul> |
| 86 | | <li><b> Bayes rule </b> as defined above, operation <a class="el" href="classbdm_1_1BM.html#a60b1779a577367c369a932cabd3a6188" title="Incremental Bayes rule.">bdm::BM::bayes()</a> which expects to get the current data record <code>dt</code>, <img class="formulaInl" alt="$ d_t $" src="form_175.png"/></li> |
| 87 | | <li><b> log-likelihood </b> i.e. numerical value of <img class="formulaInl" alt="$ f(d_t|d_1\ldots d_{t-1})$" src="form_176.png"/> as a typical side-product, since it is required in denominator of the above formula. For some models, computation of this value may require extra effort, hence it computation can be suppressed by setting BM::set_evalll(false).</li> |
| 88 | | <li><b> prediction </b> the object has enough information to create the one-step ahead predictor, i.e. <p class="formulaDsp"> |
| 89 | | <img class="formulaDsp" alt="\[ f(d_{t+1}| d_1 \ldots d_{t}), \]" src="form_180.png"/> |
| 90 | | </p> |
| 91 | | this object can be either created <a class="el" href="classbdm_1_1BM.html#a598b25e3f3d96a5bc00a5faeb5b3c912" title="Constructs conditional density of 1-step ahead predictor .">bdm::BM::predictor()</a>, sometimes it is enought only a few values of prediction hence construction of the full predictor would be too expensive operation. For this situation were designed cheaper operations <a class="el" href="classbdm_1_1BM.html#a50257e0c1e5b5c73153ea6e716ad8ae0">bdm::BM::logpred()</a> or <a class="el" href="classbdm_1_1BM.html#a688d7a2aced1e06aa1c468d73a9e5eba" title="Constructs a predictive density .">bdm::BM::epredictor()</a>. These are only basic operations, see full documentation for full range of defined operations.</li> |
| 92 | | </ul> |
| 93 | | <p>These operation are abstract, i.e. not implemented for the general class. Implementation of these operations is heavily dependent on the specific class of prior pdf, or its approximations. We can identify only a few principal approaches to this problem. For example, analytical estimation which is possible within sufficient the Exponential Family, or estimation when both prior and posterior are approximated by empirical densities. These approaches are first level of descendants of class <code>BM</code>, classes <a class="el" href="classbdm_1_1BMEF.html" title="Estimator for Exponential family.">bdm::BMEF</a> and <a class="el" href="classbdm_1_1PF.html" title="Trivial particle filter with proposal density equal to parameter evolution model...">bdm::PF</a>, repectively.</p> |
| 94 | | <p>Variants of these approaches are implemented as descendats of these level-two classes. This way, each estimation method (represented a class) is fitted in its place in the tree of approximations. This is useful even from software point of view, since related approximations have common methods and data fields.</p> |
| 95 | | <h2><a class="anchor" id="ug2_arx_basic"> |
| | 73 | The object <a class="el" href="classbdm_1_1BM.html" title="Bayesian Model of a system, i.e. all uncertainty is modeled by probabilities.">bdm::BM</a> is basic software image of the Bayes rule: <p class="formulaDsp"> |
| | 74 | <img class="formulaDsp" alt="\[ f(x_t|d_1\ldots d_t) \propto f(d_t|x_t,d_1\ldots d_{t-1}) f(x_t| d_1\ldots d_{t-1}) \]" src="form_174.png"> |
| | 75 | <p> |
| | 76 | <p> |
| | 77 | Since this operation can not be defined universally, the object is defined as abstract class with methods for:<ul> |
| | 78 | <li><b> Bayes rule </b> as defined above, operation <a class="el" href="classbdm_1_1BM.html#60b1779a577367c369a932cabd3a6188" title="Incremental Bayes rule.">bdm::BM::bayes()</a> which expects to get the current data record <code>dt</code>, <img class="formulaInl" alt="$ d_t $" src="form_175.png"></li><li><b> log-likelihood </b> i.e. numerical value of <img class="formulaInl" alt="$ f(d_t|d_1\ldots d_{t-1})$" src="form_176.png"> as a typical side-product, since it is required in denominator of the above formula. For some models, computation of this value may require extra effort, hence it computation can be suppressed by setting BM::set_evalll(false).</li><li><b> prediction </b> the object has enough information to create the one-step ahead predictor, i.e. <p class="formulaDsp"> |
| | 79 | <img class="formulaDsp" alt="\[ f(d_{t+1}| d_1 \ldots d_{t}), \]" src="form_180.png"> |
| | 80 | <p> |
| | 81 | this object can be either created <a class="el" href="classbdm_1_1BM.html#598b25e3f3d96a5bc00a5faeb5b3c912" title="Constructs conditional density of 1-step ahead predictor .">bdm::BM::predictor()</a>, sometimes it is enought only a few values of prediction hence construction of the full predictor would be too expensive operation. For this situation were designed cheaper operations <a class="el" href="classbdm_1_1BM.html#50257e0c1e5b5c73153ea6e716ad8ae0">bdm::BM::logpred()</a> or <a class="el" href="classbdm_1_1BM.html#688d7a2aced1e06aa1c468d73a9e5eba" title="Constructs a predictive density .">bdm::BM::epredictor()</a>. These are only basic operations, see full documentation for full range of defined operations.</li></ul> |
| | 82 | <p> |
| | 83 | These operation are abstract, i.e. not implemented for the general class. Implementation of these operations is heavily dependent on the specific class of prior pdf, or its approximations. We can identify only a few principal approaches to this problem. For example, analytical estimation which is possible within sufficient the Exponential Family, or estimation when both prior and posterior are approximated by empirical densities. These approaches are first level of descendants of class <code>BM</code>, classes <a class="el" href="classbdm_1_1BMEF.html" title="Estimator for Exponential family.">bdm::BMEF</a> and <a class="el" href="classbdm_1_1PF.html" title="Trivial particle filter with proposal density equal to parameter evolution model...">bdm::PF</a>, repectively.<p> |
| | 84 | Variants of these approaches are implemented as descendats of these level-two classes. This way, each estimation method (represented a class) is fitted in its place in the tree of approximations. This is useful even from software point of view, since related approximations have common methods and data fields.<h2><a class="anchor" name="ug2_arx_basic"> |
| 102 | | A1.options = <span class="stringliteral">'logbounds,logll'</span>; |
| 103 | | </pre></div><p> This is the minimal configuration of an ARX estimator. Optional elements of <a class="el" href="classbdm_1_1ARX.html#a9637412df898048bafaefee9dc7e9f6c">bdm::ARX::from_setting()</a> were set using their default values:</p> |
| 104 | | <p>The first three fileds are self explanatory, they identify which data are predicted (field <code>rv</code>) and which are in regressor (field <code>rgr</code>). The field <code>options</code> is a string of options passed to the object. In particular, class <code>BM</code> understand only options related to storing results:</p> |
| 105 | | <ul> |
| 106 | | <li>logbounds - store also lower and upper bounds on estimates (obtained by calling BM::posterior().qbounds()),</li> |
| 107 | | <li>logll - store also loglikelihood of each step of the Bayes rule. These values are stored in given logger (ug_loggers). By default, only mean values of the estimate are stored.</li> |
| 108 | | </ul> |
| 109 | | <p>Storing of the log-likelihood is useful, e.g. in model selection task when too models are compared.</p> |
| 110 | | <p>The bounds are useful e.g. for visualization of the results. Run of the example should provide result like the following: </p> |
| 111 | | <div align="center"> |
| 112 | | <img src="arx_basic_example_small.png" alt="arx_basic_example_small.png"/> |
| | 90 | A1.options = <span class="stringliteral">'logbounds,logll'</span>; |
| | 91 | </pre></div> This is the minimal configuration of an ARX estimator. Optional elements of <a class="el" href="classbdm_1_1ARX.html#9637412df898048bafaefee9dc7e9f6c">bdm::ARX::from_setting()</a> were set using their default values:<p> |
| | 92 | The first three fileds are self explanatory, they identify which data are predicted (field <code>rv</code>) and which are in regressor (field <code>rgr</code>). The field <code>options</code> is a string of options passed to the object. In particular, class <code>BM</code> understand only options related to storing results:<ul> |
| | 93 | <li>logbounds - store also lower and upper bounds on estimates (obtained by calling BM::posterior().qbounds()),</li><li>logll - store also loglikelihood of each step of the Bayes rule. These values are stored in given logger (ug_loggers). By default, only mean values of the estimate are stored.</li></ul> |
| | 94 | <p> |
| | 95 | Storing of the log-likelihood is useful, e.g. in model selection task when two models are compared.<p> |
| | 96 | The bounds are useful e.g. for visualization of the results. Run of the example should provide result like the following: <div align="center"> |
| | 97 | <img src="arx_basic_example_small.png" alt="arx_basic_example_small.png"> |
| 123 | | </pre></div><p> That is, two other ARX estimators are created,</p> |
| 124 | | <ul> |
| 125 | | <li>A2 which is the same as A1 except it does not model constant term in the linear regression. Note that if the constant was set to zero, then this is the correct model.</li> |
| 126 | | <li>A3 which is the same as A2, but assumes time-variant parameters with forgetting factor 0.95.</li> |
| 127 | | </ul> |
| 128 | | <p>Since all estimator were configured to store values of marginal log-likelihood, we can easily compare them by computint total log-likelihood for each of them and converting them to probabilities. Typically, the results should look like: </p> |
| 129 | | <div class="fragment"><pre class="fragment">Model_probabilities = |
| | 107 | </pre></div> That is, two other ARX estimators are created,<ul> |
| | 108 | <li>A2 which is the same as A1 except it does not model constant term in the linear regression. Note that if the constant was set to zero, then this is the correct model.</li><li>A3 which is the same as A2, but assumes time-variant parameters with forgetting factor 0.95.</li></ul> |
| | 109 | <p> |
| | 110 | Since all estimator were configured to store values of marginal log-likelihood, we can easily compare them by computint total log-likelihood for each of them and converting them to probabilities. Typically, the results should look like: <div class="fragment"><pre class="fragment">Model_probabilities = |
| 132 | | </pre></div><p> Hence, the true model A2 was correctly identified as the most likely to produce this data.</p> |
| 133 | | <p>For this task, additional technical adjustments were needed: </p> |
| 134 | | <div class="fragment"><pre class="fragment">A1.name=<span class="stringliteral">'A1'</span>; |
| 135 | | A2.name=<span class="stringliteral">'A2'</span>; |
| 136 | | A2.rv_param = RV({<span class="stringliteral">'a2th'</span>, <span class="charliteral">'r'</span>},[2,1],[0,0]); |
| 137 | | A3.name=<span class="stringliteral">'A3'</span>; |
| 138 | | A3.rv_param = RV({<span class="stringliteral">'a3th'</span>, <span class="charliteral">'r'</span>},[2,1],[0,0]); |
| 139 | | </pre></div><p> First, in order to distinguish the estimators from each other, the estimators were given names. Hence, the results will be logged with prefix given by the name, such as M.A1ll for field <code>ll</code>.</p> |
| 140 | | <p>Second, if the parameters of a ARX model are not specified, they are automatically named <code>theta</code> and <code>r</code>. However, in this case, <code>A1</code> and <code>A2</code> differ in size, hence their random variables differ and can not use the same name. Therefore, we have explicitly used another names (RVs) of the parameters.</p> |
| 141 | | <h2><a class="anchor" id="ug2_bm_composition"> |
| | 113 | </pre></div> Hence, the true model A2 was correctly identified as the most likely to produce this data.<p> |
| | 114 | For this task, additional technical adjustments were needed: <div class="fragment"><pre class="fragment">A1.name=<span class="stringliteral">'A1'</span>; |
| | 115 | A2.name=<span class="stringliteral">'A2'</span>; |
| | 116 | A2.rv_param = RV({<span class="stringliteral">'a2th'</span>, <span class="charliteral">'r'</span>},[2,1],[0,0]); |
| | 117 | A3.name=<span class="stringliteral">'A3'</span>; |
| | 118 | A3.rv_param = RV({<span class="stringliteral">'a3th'</span>, <span class="charliteral">'r'</span>},[2,1],[0,0]); |
| | 119 | </pre></div> First, in order to distinguish the estimators from each other, the estimators were given names. Hence, the results will be logged with prefix given by the name, such as M.A1ll for field <code>ll</code>.<p> |
| | 120 | Second, if the parameters of a ARX model are not specified, they are automatically named <code>theta</code> and <code>r</code>. However, in this case, <code>A1</code> and <code>A2</code> differ in size, hence their random variables differ and can not use the same name. Therefore, we have explicitly used another names (RVs) of the parameters.<h2><a class="anchor" name="ug2_bm_composition"> |
| 143 | | <p>Similarly to mpdfs which could be composed via <code>mprod</code>, the Bayesian models can be composed. However, justification of this step is less clear than in the case of epdfs.</p> |
| 144 | | <p>One possible theoretical base of composition is the Marginalized particle filter, which splits the prior and the posterior in two parts: </p> |
| 145 | | <p class="formulaDsp"> |
| 146 | | <img class="formulaDsp" alt="\[ f(x_t|d_1\ldots d_t)=f(x_{1,t}|x_{2,t},d_1\ldots d_t)f(x_{2,t}|d_1\ldots d_t) \]" src="form_181.png"/> |
| 147 | | </p> |
| 148 | | <p> each of these parts is estimated using different approach. The first part is assumed to be analytically tractable, while the second is approximated using empirical approximation.</p> |
| 149 | | <p>The whole algorithm runs by parallel evaluation of many <code>BMs</code> for estimation of <img class="formulaInl" alt="$ x_{1,t}$" src="form_182.png"/>, each of them conditioned on value of a sample of <img class="formulaInl" alt="$ x_{2,t}$" src="form_191.png"/>.</p> |
| 150 | | <p>For example, the forgetting factor, <img class="formulaInl" alt="$ \phi $" src="form_135.png"/> of an ARX model can be considered to be unknown. Then, the whole parameter space is <img class="formulaInl" alt="$ [\theta_t, r_t, \phi_t]$" src="form_184.png"/> decomposed as follows: </p> |
| 151 | | <p class="formulaDsp"> |
| 152 | | <img class="formulaDsp" alt="\[ f(\theta_t, r_t, \phi_t) = f(\theta_t, r_t| \phi_t) f(\phi_t) \]" src="form_185.png"/> |
| 153 | | </p> |
| 154 | | <p> Note that for known trajectory of <img class="formulaInl" alt="$ \phi_t $" src="form_186.png"/> the standard ARX estimator can be used if we find a way how to feed the changing <img class="formulaInl" alt="$ \phi_t $" src="form_186.png"/> into it. This is achieved by a trivial extension using inheritance method <a class="el" href="classbdm_1_1BM.html#a6799f4b16a6a59ed58b1d0d6e17116f4" title="Substitute val for rvc.">bdm::BM::condition()</a>.</p> |
| 155 | | <p>Extension of standard ARX estimator to conditional estimator is implemented as class <a class="el" href="classbdm_1_1ARXfrg.html">bdm::ARXfrg</a>. The only difference from standard ARX is that this object will change its forgetting factor via method ARXfrg::condition(). Existence of this function is assumed by the MPF estimator. Informally, the name 'ARXfrg' means: "if anybody calls your condition(0.9), it tells you new value of forgetting factor".</p> |
| 156 | | <p>The MPF estimator is implemented by class <a class="el" href="classbdm_1_1MPF.html" title="Marginalized Particle filter.">bdm::MPF</a>. In the toolbox, it can be constructed as follows: </p> |
| 157 | | <div class="fragment"><pre class="fragment">%%%%%% ARX estimator conditioned on frg |
| | 122 | Similarly to mpdfs which could be composed via <code>mprod</code>, the Bayesian models can be composed. However, justification of this step is less clear than in the case of epdfs.<p> |
| | 123 | One possible theoretical base of composition is the Marginalized particle filter, which splits the prior and the posterior in two parts: <p class="formulaDsp"> |
| | 124 | <img class="formulaDsp" alt="\[ f(x_t|d_1\ldots d_t)=f(x_{1,t}|x_{2,t},d_1\ldots d_t)f(x_{2,t}|d_1\ldots d_t) \]" src="form_181.png"> |
| | 125 | <p> |
| | 126 | each of these parts is estimated using different approach. The first part is assumed to be analytically tractable, while the second is approximated using empirical approximation.<p> |
| | 127 | The whole algorithm runs by parallel evaluation of many <code>BMs</code> for estimation of <img class="formulaInl" alt="$ x_{1,t}$" src="form_182.png">, each of them conditioned on value of a sample of <img class="formulaInl" alt="$ x_{2,t}$" src="form_191.png">.<p> |
| | 128 | For example, the forgetting factor, <img class="formulaInl" alt="$ \phi $" src="form_135.png"> of an ARX model can be considered to be unknown. Then, the whole parameter space is <img class="formulaInl" alt="$ [\theta_t, r_t, \phi_t]$" src="form_184.png"> decomposed as follows: <p class="formulaDsp"> |
| | 129 | <img class="formulaDsp" alt="\[ f(\theta_t, r_t, \phi_t) = f(\theta_t, r_t| \phi_t) f(\phi_t) \]" src="form_185.png"> |
| | 130 | <p> |
| | 131 | Note that for known trajectory of <img class="formulaInl" alt="$ \phi_t $" src="form_186.png"> the standard ARX estimator can be used if we find a way how to feed the changing <img class="formulaInl" alt="$ \phi_t $" src="form_186.png"> into it. This is achieved by a trivial extension using inheritance method <a class="el" href="classbdm_1_1BM.html#6799f4b16a6a59ed58b1d0d6e17116f4" title="Substitute val for rvc.">bdm::BM::condition()</a>.<p> |
| | 132 | Extension of standard ARX estimator to conditional estimator is implemented as class <a class="el" href="classbdm_1_1ARXfrg.html">bdm::ARXfrg</a>. The only difference from standard ARX is that this object will change its forgetting factor via method ARXfrg::condition(). Existence of this function is assumed by the MPF estimator. Informally, the name 'ARXfrg' means: "if anybody calls your condition(0.9), it tells you new value of forgetting factor".<p> |
| | 133 | The MPF estimator is implemented by class <a class="el" href="classbdm_1_1MPF.html" title="Marginalized Particle filter.">bdm::MPF</a>. In the toolbox, it can be constructed as follows: <div class="fragment"><pre class="fragment">%%%%%% ARX estimator conditioned on frg |
| 183 | | </pre></div><p>Here, the configuration structure <code>A1</code> is a description of an ARX model, as used in previous examples, the only difference is in its name 'ARXfrg'.</p> |
| 184 | | <p>The configuration structure <code>phi_pdf</code> defines random walk on the forgetting factor. It was chosen as Dirichlet, hence it will produce 2-dimensional vector of <img class="formulaInl" alt="$[\phi, 1-\phi]$" src="form_192.png"/>. The class <code>ARXfrg</code> was designed to read only the first element of its condition. The random walk of type mDirich is: </p> |
| 185 | | <p class="formulaDsp"> |
| 186 | | <img class="formulaDsp" alt="\[ f(\phi_t|\phi_{t-1}) = Di (\phi_{t-1}/k + \beta_c) \]" src="form_193.png"/> |
| 187 | | </p> |
| 188 | | <p> where <img class="formulaInl" alt="$ k $" src="form_89.png"/> influences the spread of the walk and <img class="formulaInl" alt="$ \beta_c $" src="form_189.png"/> has the role of stabilizing, to avoid traps of corner cases such as [0,1] and [1,0]. Its influence on the results is quite important.</p> |
| 189 | | <p>This example is implemented as bdmtoolbox/tutorial/userguide/frg_example.m Its typical run should look like the following: </p> |
| 190 | | <div align="center"> |
| 191 | | <img src="frg_example_small.png" alt="frg_example_small.png"/> |
| | 159 | </pre></div><p> |
| | 160 | Here, the configuration structure <code>A1</code> is a description of an ARX model, as used in previous examples, the only difference is in its name 'ARXfrg'.<p> |
| | 161 | The configuration structure <code>phi_pdf</code> defines random walk on the forgetting factor. It was chosen as Dirichlet, hence it will produce 2-dimensional vector of <img class="formulaInl" alt="$[\phi, 1-\phi]$" src="form_192.png">. The class <code>ARXfrg</code> was designed to read only the first element of its condition. The random walk of type mDirich is: <p class="formulaDsp"> |
| | 162 | <img class="formulaDsp" alt="\[ f(\phi_t|\phi_{t-1}) = Di (\phi_{t-1}/k + \beta_c) \]" src="form_193.png"> |
| | 163 | <p> |
| | 164 | where <img class="formulaInl" alt="$ k $" src="form_89.png"> influences the spread of the walk and <img class="formulaInl" alt="$ \beta_c $" src="form_189.png"> has the role of stabilizing, to avoid traps of corner cases such as [0,1] and [1,0]. Its influence on the results is quite important.<p> |
| | 165 | This example is implemented as bdmtoolbox/tutorial/userguide/frg_example.m Its typical run should look like the following: <div align="center"> |
| | 166 | <img src="frg_example_small.png" alt="frg_example_small.png"> |
