root/library/doc/html/user_guide2.html @ 661

Revision 661, 16.8 kB (checked in by smidl, 15 years ago)

doc

Line 
1<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Transitional//EN" "http://www.w3.org/TR/xhtml1/DTD/xhtml1-transitional.dtd">
2<html xmlns="http://www.w3.org/1999/xhtml">
3<head>
4<meta http-equiv="Content-Type" content="text/xhtml;charset=UTF-8"/>
5<title>mixpp: BDM Use - Estimation and Bayes Rule</title>
6<link href="tabs.css" rel="stylesheet" type="text/css"/>
7<link href="doxygen.css" rel="stylesheet" type="text/css"/>
8</head>
9<body>
10<!-- Generated by Doxygen 1.6.1 -->
11<script type="text/javascript">
12<!--
13function changeDisplayState (e){
14  var num=this.id.replace(/[^[0-9]/g,'');
15  var button=this.firstChild;
16  var sectionDiv=document.getElementById('dynsection'+num);
17  if (sectionDiv.style.display=='none'||sectionDiv.style.display==''){
18    sectionDiv.style.display='block';
19    button.src='open.gif';
20  }else{
21    sectionDiv.style.display='none';
22    button.src='closed.gif';
23  }
24}
25function initDynSections(){
26  var divs=document.getElementsByTagName('div');
27  var sectionCounter=1;
28  for(var i=0;i<divs.length-1;i++){
29    if(divs[i].className=='dynheader'&&divs[i+1].className=='dynsection'){
30      var header=divs[i];
31      var section=divs[i+1];
32      var button=header.firstChild;
33      if (button!='IMG'){
34        divs[i].insertBefore(document.createTextNode(' '),divs[i].firstChild);
35        button=document.createElement('img');
36        divs[i].insertBefore(button,divs[i].firstChild);
37      }
38      header.style.cursor='pointer';
39      header.onclick=changeDisplayState;
40      header.id='dynheader'+sectionCounter;
41      button.src='closed.gif';
42      section.id='dynsection'+sectionCounter;
43      section.style.display='none';
44      section.style.marginLeft='14px';
45      sectionCounter++;
46    }
47  }
48}
49window.onload = initDynSections;
50-->
51</script>
52<div class="navigation" id="top">
53  <div class="tabs">
54    <ul>
55      <li><a href="main.html"><span>Main&nbsp;Page</span></a></li>
56      <li class="current"><a href="pages.html"><span>Related&nbsp;Pages</span></a></li>
57      <li><a href="annotated.html"><span>Classes</span></a></li>
58      <li><a href="files.html"><span>Files</span></a></li>
59    </ul>
60  </div>
61  <div class="navpath"><a class="el" href="dev_guide.html">BDM Use - in C++</a>&nbsp;&raquo&nbsp;<a class="el" href="kalman.html">Examples of (extended) Kalman filtering</a>
62  </div>
63</div>
64<div class="contents">
65
66
67<h1><a class="anchor" id="user_guide2">BDM Use - Estimation and Bayes Rule </a></h1><p>Baysian theory is predominantly used in system identification, or estimation problems. This section is concerned with recursive estimation, as implemneted in prepared scenario <code>estimator</code>.</p>
68<p>The function of the <code>estimator</code> is graphically illustrated: </p>
69<div align="center">
70<img src="inline_dotgraph_2.dot.png" alt="inline_dotgraph_2.dot" border="0" usemap="#inline_dotgraph_2.dot.map">
71<map name="inline_dotgraph_2.dot.map" id="inline_dotgraph_2.dot.map"></map>
72</div>
73<p>Here, </p>
74<ul>
75<li>Data Source is an object (class DS) providing sequential data, <img class="formulaInl" alt="$ [d_1, d_2, \ldots d_t] $" src="form_154.png"/>. </li>
76<li>Bayesian Model is an object (class BM) performing Bayesian filtering, </li>
77<li>Result Logger is an object (class logger) dedicated to storing important data from the experiment.</li>
78</ul>
79<p>Since objects datasource and the logger has already been introduced in section <a class="el" href="user_guide.html">BDM Use - System, Data, Simulation</a>, it remains to introduce object <code>Bayesian</code> <code>Model</code> (<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>).</p>
80<h2><a class="anchor" id="ug2_theory">
81Bayes rule and estimation</a></h2>
82<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>
83<p class="formulaDsp">
84<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_155.png"/>
85</p>
86<p>Since this operation can not be defined universally, the object is defined as abstract class with methods for:</p>
87<ul>
88<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_156.png"/></li>
89<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_157.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>
90<li><b> prediction </b> the object has enough information to create the one-step ahead predictor, i.e. <p class="formulaDsp">
91<img class="formulaDsp" alt="\[ f(d_{t+1}| d_1 \ldots d_{t}), \]" src="form_158.png"/>
92</p>
93 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>
94</ul>
95<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>
96<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>
97<h2><a class="anchor" id="ug2_arx_basic">
98Estimation of ARX models</a></h2>
99<p>Autoregressive models has already been introduced in <a class="el" href="user_guide.html#ug_arx_sim">Simulating autoregressive model</a> where their simulator has been presented. We will use results of simulation of the ARX datasource defined there to provide data for estimation using MemDS.</p>
100<p>The following code is from bdmtoolbox/tutorial/userguide/arx_basic_example.m </p>
101<div class="fragment"><pre class="fragment">A1.class = <span class="stringliteral">&apos;ARX&apos;</span>;
102A1.rv = y;
103A1.rgr = RVtimes([y,y],[-3,-1]) ;
104A1.options = <span class="stringliteral">&apos;logbounds,logll&apos;</span>;
105</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>
106<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>
107<ul>
108<li>logbounds - store also lower and upper bounds on estimates (obtained by calling BM::posterior().qbounds()),</li>
109<li>logll - store also loglikelihood of each step of the Bayes rule. These values are stored in given logger (<a class="el" href="user_guide.html#ug_loggers">Loggers for flexible handling of results</a>). By default, only mean values of the estimate are stored.</li>
110</ul>
111<p>Storing of the log-likelihood is useful, e.g. in model selection task when two models are compared.</p>
112<p>The bounds are useful e.g. for visualization of the results. Run of the example should provide result like the following: </p>
113<div align="center">
114<img src="arx_basic_example_small.png" alt="arx_basic_example_small.png"/>
115</div>
116 <h2><a class="anchor" id="ug2_model_sel">
117Model selection</a></h2>
118<p>In Bayesian framework, model selection is done via comparison of marginal likelihood of the recorded data. See [some theory].</p>
119<p>A trivial exammple how this can be done is presented in file bdmtoolbox/tutorial/userguide/arx_selection_example.m. The code extends the basic A1 object as follows: </p>
120<div class="fragment"><pre class="fragment">A2=A1;
121A2.constant = 0;
122
123A3=A2;
124A3.frg = 0.95;
125</pre></div><p> That is, two other ARX estimators are created,</p>
126<ul>
127<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>
128<li>A3 which is the same as A2, but assumes time-variant parameters with forgetting factor 0.95.</li>
129</ul>
130<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>
131<div class="fragment"><pre class="fragment">Model_probabilities =
132
133    0.0002    0.7318    0.2680
134</pre></div><p> Hence, the true model A2 was correctly identified as the most likely to produce this data.</p>
135<p>For this task, additional technical adjustments were needed: </p>
136<div class="fragment"><pre class="fragment">A1.name=<span class="stringliteral">&apos;A1&apos;</span>;
137A2.name=<span class="stringliteral">&apos;A2&apos;</span>;
138A2.rv_param = RV({<span class="stringliteral">&apos;a2th&apos;</span>, <span class="charliteral">&apos;r&apos;</span>},[2,1],[0,0]);
139A3.name=<span class="stringliteral">&apos;A3&apos;</span>;
140A3.rv_param = RV({<span class="stringliteral">&apos;a3th&apos;</span>, <span class="charliteral">&apos;r&apos;</span>},[2,1],[0,0]);
141</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>
142<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>
143<h2><a class="anchor" id="ug2_bm_composition">
144Composition of estimators</a></h2>
145<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>
146<p>One possible theoretical base of composition is the Marginalized particle filter, which splits the prior and the posterior in two parts: </p>
147<p class="formulaDsp">
148<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_159.png"/>
149</p>
150<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>
151<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_160.png"/>, each of them conditioned on value of a sample of <img class="formulaInl" alt="$ x_{2,t}$" src="form_161.png"/>.</p>
152<p>For example, the forgetting factor, <img class="formulaInl" alt="$ \phi $" src="form_162.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_163.png"/> decomposed as follows: </p>
153<p class="formulaDsp">
154<img class="formulaDsp" alt="\[ f(\theta_t, r_t, \phi_t) = f(\theta_t, r_t| \phi_t) f(\phi_t) \]" src="form_164.png"/>
155</p>
156<p> Note that for known trajectory of <img class="formulaInl" alt="$ \phi_t $" src="form_165.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_165.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>
157<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>
158<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>
159<div class="fragment"><pre class="fragment">%%%%%% ARX estimator conditioned on frg
160
161A1.class = <span class="stringliteral">&apos;ARXfrg&apos;</span>;
162A1.rv = y;
163A1.rgr = RVtimes([y,u],[-3,-1]) ;
164A1.options =<span class="stringliteral">&apos;logbounds,logll&apos;</span>;
165A1.frg = 0.9;
166A1.name = <span class="stringliteral">&apos;A1&apos;</span>;
167
168%%%%%% Random walk on frg - Dirichlet
169phi_pdf.class = <span class="stringliteral">&apos;mDirich&apos;</span>;         % random walk on coefficient phi
170phi_pdf.rv    = RV(<span class="stringliteral">&apos;phi&apos;</span>,2);       % 2D random walk - frg is the first element
171phi_pdf.k     = 0.01;              % width of the random walk
172phi_pdf.betac = [0.01 0.01];       % stabilizing elememnt of random walk
173
174%%%%%% Combining estimators in Marginalized particle filter
175E.class = <span class="stringliteral">&apos;MPF&apos;</span>;
176E.BM = A1;                         % ARX is the analytical part
177E.parameter_pdf = phi_pdf;         % Random walk is the parameter evolution model
178E.n = 20;                          % number of particles
179E.prior.class = <span class="stringliteral">&apos;eDirich&apos;</span>;         % prior on non-linear part
180E.prior.beta  = [1 1]; %
181E.options =<span class="stringliteral">&apos;logbounds,logll&apos;</span>;
182E.name = <span class="stringliteral">&apos;MPF&apos;</span>;
183
184M=estimator(DS,{E});
185</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>
186<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_166.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>
187<p class="formulaDsp">
188<img class="formulaDsp" alt="\[ f(\phi_t|\phi_{t-1}) = Di (\phi_{t-1}/k + \beta_c) \]" src="form_167.png"/>
189</p>
190<p> where <img class="formulaInl" alt="$ k $" src="form_92.png"/> influences the spread of the walk and <img class="formulaInl" alt="$ \beta_c $" src="form_81.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>
191<p>This example is implemented as bdmtoolbox/tutorial/userguide/frg_example.m Its typical run should look like the following: </p>
192<div align="center">
193<img src="frg_example_small.png" alt="frg_example_small.png"/>
194</div>
195 <p>Note: error bars in this case are not directly comparable with those of previous examples. The MPF class implements the qbounds function as minimum and maximum of bounds in the condidered set (even if its weight is extreemly small). Hence, the bounds of the MPF are probably larger than it should be. Nevertheless, they provide great help when designing and tuning algorithms. </p>
196</div>
197<hr size="1"/><address style="text-align: right;"><small>Generated on Wed Oct 14 23:00:33 2009 for mixpp by&nbsp;
198<a href="http://www.doxygen.org/index.html">
199<img class="footer" src="doxygen.png" alt="doxygen"/></a> 1.6.1 </small></address>
200</body>
201</html>
Note: See TracBrowser for help on using the browser.