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

Revision 651, 16.0 kB (checked in by mido, 15 years ago)

\doc directory cleaned a bit

Line 
1<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
2<html><head><meta http-equiv="Content-Type" content="text/html;charset=UTF-8">
3<title>mixpp: BDM Use - Estimation and Bayes Rule</title>
4<link href="tabs.css" rel="stylesheet" type="text/css">
5<link href="doxygen.css" rel="stylesheet" type="text/css">
6</head><body>
7<!-- Generated by Doxygen 1.5.9 -->
8<script type="text/javascript">
9<!--
10function changeDisplayState (e){
11  var num=this.id.replace(/[^[0-9]/g,'');
12  var button=this.firstChild;
13  var sectionDiv=document.getElementById('dynsection'+num);
14  if (sectionDiv.style.display=='none'||sectionDiv.style.display==''){
15    sectionDiv.style.display='block';
16    button.src='open.gif';
17  }else{
18    sectionDiv.style.display='none';
19    button.src='closed.gif';
20  }
21}
22function initDynSections(){
23  var divs=document.getElementsByTagName('div');
24  var sectionCounter=1;
25  for(var i=0;i<divs.length-1;i++){
26    if(divs[i].className=='dynheader'&&divs[i+1].className=='dynsection'){
27      var header=divs[i];
28      var section=divs[i+1];
29      var button=header.firstChild;
30      if (button!='IMG'){
31        divs[i].insertBefore(document.createTextNode(' '),divs[i].firstChild);
32        button=document.createElement('img');
33        divs[i].insertBefore(button,divs[i].firstChild);
34      }
35      header.style.cursor='pointer';
36      header.onclick=changeDisplayState;
37      header.id='dynheader'+sectionCounter;
38      button.src='closed.gif';
39      section.id='dynsection'+sectionCounter;
40      section.style.display='none';
41      section.style.marginLeft='14px';
42      sectionCounter++;
43    }
44  }
45}
46window.onload = initDynSections;
47-->
48</script>
49<div class="navigation" id="top">
50  <div class="tabs">
51    <ul>
52      <li><a href="main.html"><span>Main&nbsp;Page</span></a></li>
53      <li class="current"><a href="pages.html"><span>Related&nbsp;Pages</span></a></li>
54      <li><a href="annotated.html"><span>Classes</span></a></li>
55      <li><a href="files.html"><span>Files</span></a></li>
56    </ul>
57  </div>
58</div>
59<div class="contents">
60<h1><a class="anchor" name="user_guide2">BDM Use - Estimation and Bayes Rule </a></h1>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>
61The function of the <code>estimator</code> is graphically illustrated: <div align="center">
62<img src="inline_dotgraph_2.dot.png" alt="inline_dotgraph_2.dot" border="0" usemap="#inline_dotgraph_2.dot.map">
63<map name="inline_dotgraph_2.dot.map"></map>
64</div>
65<p>
66Here, <ul>
67<li>Data Source is an object (class DS) providing sequential data, <img class="formulaInl" alt="$ [d_1, d_2, \ldots d_t] $" src="form_173.png">. </li>
68<li>Bayesian Model is an object (class BM) performing Bayesian filtering, </li>
69<li>Result Logger is an object (class logger) dedicated to storing important data from the experiment.</li>
70</ul>
71Since 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>).<h2><a class="anchor" name="ug2_theory">
72Bayes rule and estimation</a></h2>
73The 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>
77Since 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>
83These 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>
84Variants 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">
85Estimation of ARX models</a></h2>
86Autoregressive 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>
87The following code is from bdmtoolbox/tutorial/userguide/arx_basic_example.m <div class="fragment"><pre class="fragment">A1.class = <span class="stringliteral">'ARX'</span>;
88A1.rv = y;
89A1.rgr = RVtimes([y,y],[-3,-1]) ;
90A1.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>
92The 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>
95Storing of the log-likelihood is useful, e.g. in model selection task when two models are compared.<p>
96The 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">
98</div>
99 <h2><a class="anchor" name="ug2_model_sel">
100Model selection</a></h2>
101In Bayesian framework, model selection is done via comparison of marginal likelihood of the recorded data. See [some theory].<p>
102A 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: <div class="fragment"><pre class="fragment">A2=A1;
103A2.constant = 0;
104
105A3=A2;
106A3.frg = 0.95;
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>
110Since 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 =
111
112    0.0002    0.7318    0.2680
113</pre></div> Hence, the true model A2 was correctly identified as the most likely to produce this data.<p>
114For this task, additional technical adjustments were needed: <div class="fragment"><pre class="fragment">A1.name=<span class="stringliteral">'A1'</span>;
115A2.name=<span class="stringliteral">'A2'</span>;
116A2.rv_param = RV({<span class="stringliteral">'a2th'</span>, <span class="charliteral">'r'</span>},[2,1],[0,0]);
117A3.name=<span class="stringliteral">'A3'</span>;
118A3.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>
120Second, 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">
121Composition of estimators</a></h2>
122Similarly 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>
123One 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>
127The 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>
128For 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>
132Extension 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>
133The 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
134
135A1.class = <span class="stringliteral">'ARXfrg'</span>;
136A1.rv = y;
137A1.rgr = RVtimes([y,u],[-3,-1]) ;
138A1.options =<span class="stringliteral">'logbounds,logll'</span>;
139A1.frg = 0.9;
140A1.name = <span class="stringliteral">'A1'</span>;
141
142%%%%%% Random walk on frg - Dirichlet
143phi_pdf.class = <span class="stringliteral">'mDirich'</span>;         % random walk on coefficient phi
144phi_pdf.rv    = RV(<span class="stringliteral">'phi'</span>,2);       % 2D random walk - frg is the first element
145phi_pdf.k     = 0.01;              % width of the random walk
146phi_pdf.betac = [0.01 0.01];       % stabilizing elememnt of random walk
147
148%%%%%% Combining estimators in Marginalized particle filter
149E.class = <span class="stringliteral">'MPF'</span>;
150E.BM = A1;                         % ARX is the analytical part
151E.parameter_pdf = phi_pdf;         % Random walk is the parameter evolution model
152E.n = 20;                          % number of particles
153E.prior.class = <span class="stringliteral">'eDirich'</span>;         % prior on non-linear part
154E.prior.beta  = [1 1]; %
155E.options =<span class="stringliteral">'logbounds,logll'</span>;
156E.name = <span class="stringliteral">'MPF'</span>;
157
158M=estimator(DS,{E});
159</pre></div><p>
160Here, 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>
161The 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>
165This 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">
167</div>
168 <p>
169Note: 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. </div>
170<hr size="1"><address style="text-align: right;"><small>Generated on Wed Oct 7 17:34:44 2009 for mixpp by&nbsp;
171<a href="http://www.doxygen.org/index.html">
172<img src="doxygen.png" alt="doxygen" align="middle" border="0"></a> 1.5.9 </small></address>
173</body>
174</html>
Note: See TracBrowser for help on using the browser.