Show
Ignore:
Timestamp:
02/11/09 14:15:39 (16 years ago)
Author:
smidl
Message:

doc

Files:
1 modified

Legend:

Unmodified
Added
Removed

  • doc/html/libPF_8h-source.html

    r264 r269  
    2626<a name="l00019"></a>00019 <span class="keyword">namespace </span>bdm{ 
    2727<a name="l00020"></a>00020  
    28 <a name="l00027"></a><a class="code" href="classbdm_1_1PF.html">00027</a> <span class="keyword">class </span><a class="code" href="classbdm_1_1PF.html" title="Trivial particle filter with proposal density equal to parameter evolution model...">PF</a> : <span class="keyword">public</span> <a class="code" href="classbdm_1_1BM.html" title="Bayesian Model of the world, i.e. all uncertainty is modeled by probabilities.">BM</a> { 
     28<a name="l00027"></a><a class="code" href="classbdm_1_1PF.html">00027</a> <span class="keyword">class </span><a class="code" href="classbdm_1_1PF.html" title="Trivial particle filter with proposal density equal to parameter evolution model...">PF</a> : <span class="keyword">public</span> <a class="code" href="classbdm_1_1BM.html" title="Bayesian Model of a system, i.e. all uncertainty is modeled by probabilities.">BM</a> { 
    2929<a name="l00028"></a>00028 <span class="keyword">protected</span>: 
    3030<a name="l00030"></a><a class="code" href="classbdm_1_1PF.html#eeafaf9b8ad75fe62ee9fd6369e3f7fe">00030</a>         <span class="keywordtype">int</span> <a class="code" href="classbdm_1_1PF.html#eeafaf9b8ad75fe62ee9fd6369e3f7fe" title="number of particles;">n</a>; 
    … …  
    3535<a name="l00040"></a><a class="code" href="classbdm_1_1PF.html#c58b8fa634272c3f48845a9020ba55aa">00040</a>         <a class="code" href="classbdm_1_1mpdf.html" title="Conditional probability density, e.g. modeling some dependencies.">mpdf</a> &amp;<a class="code" href="classbdm_1_1PF.html#c58b8fa634272c3f48845a9020ba55aa" title="Observation model.">obs</a>; 
    3636<a name="l00041"></a>00041 <span class="keyword">public</span>: 
    37 <a name="l00043"></a><a class="code" href="classbdm_1_1PF.html#df9e7dc54dbe1c78a4aa503c42f7f8c4">00043</a>         <a class="code" href="classbdm_1_1PF.html#df9e7dc54dbe1c78a4aa503c42f7f8c4" title="Default constructor.">PF</a> ( <span class="keyword">const</span> <a class="code" href="classbdm_1_1RV.html" title="Class representing variables, most often random variables.">RV</a> &amp;rv0, <a class="code" href="classbdm_1_1mpdf.html" title="Conditional probability density, e.g. modeling some dependencies.">mpdf</a> &amp;par0,  <a class="code" href="classbdm_1_1mpdf.html" title="Conditional probability density, e.g. modeling some dependencies.">mpdf</a> &amp;obs0, <span class="keywordtype">int</span> n0 ) :<a class="code" href="classbdm_1_1BM.html" title="Bayesian Model of the world, i.e. all uncertainty is modeled by probabilities.">BM</a> ( rv0 ), 
     37<a name="l00043"></a><a class="code" href="classbdm_1_1PF.html#df9e7dc54dbe1c78a4aa503c42f7f8c4">00043</a>         <a class="code" href="classbdm_1_1PF.html#df9e7dc54dbe1c78a4aa503c42f7f8c4" title="Default constructor.">PF</a> ( <span class="keyword">const</span> <a class="code" href="classbdm_1_1RV.html" title="Class representing variables, most often random variables.">RV</a> &amp;rv0, <a class="code" href="classbdm_1_1mpdf.html" title="Conditional probability density, e.g. modeling some dependencies.">mpdf</a> &amp;par0,  <a class="code" href="classbdm_1_1mpdf.html" title="Conditional probability density, e.g. modeling some dependencies.">mpdf</a> &amp;obs0, <span class="keywordtype">int</span> n0 ) :<a class="code" href="classbdm_1_1BM.html" title="Bayesian Model of a system, i.e. all uncertainty is modeled by probabilities.">BM</a> ( rv0 ), 
    3838<a name="l00044"></a>00044                         <a class="code" href="classbdm_1_1PF.html#eeafaf9b8ad75fe62ee9fd6369e3f7fe" title="number of particles;">n</a> ( n0 ),<a class="code" href="classbdm_1_1PF.html#dc049265b9086cad7071f98d00a2b9af" title="posterior density">est</a> ( rv0,<a class="code" href="classbdm_1_1PF.html#eeafaf9b8ad75fe62ee9fd6369e3f7fe" title="number of particles;">n</a> ),<a class="code" href="classbdm_1_1PF.html#f5149d5522d1095d39240c4c607f61a3" title="pointer into eEmp ">_w</a> ( <a class="code" href="classbdm_1_1PF.html#dc049265b9086cad7071f98d00a2b9af" title="posterior density">est</a>.<a class="code" href="classbdm_1_1PF.html#f5149d5522d1095d39240c4c607f61a3" title="pointer into eEmp ">_w</a>() ),<a class="code" href="classbdm_1_1PF.html#914bd66025692c4018dbd482cb3c47c1" title="pointer into eEmp ">_samples</a> ( <a class="code" href="classbdm_1_1PF.html#dc049265b9086cad7071f98d00a2b9af" title="posterior density">est</a>.<a class="code" href="classbdm_1_1PF.html#914bd66025692c4018dbd482cb3c47c1" title="pointer into eEmp ">_samples</a>() ), 
    3939<a name="l00045"></a>00045                         <a class="code" href="classbdm_1_1PF.html#cf3a1b2a407012e47ac878e3aa2fbf34" title="Parameter evolution model.">par</a> ( par0 ), <a class="code" href="classbdm_1_1PF.html#c58b8fa634272c3f48845a9020ba55aa" title="Observation model.">obs</a> ( obs0 ) {}; 
    … …  
    118118<a name="l00136"></a>00136         } 
    119119<a name="l00137"></a>00137  
    120 <a name="l00139"></a><a class="code" href="classbdm_1_1MPF.html#82b5a34d9ed0e78452f98d2ecbf1e93c">00139</a>         <a class="code" href="classbdm_1_1BM.html" title="Bayesian Model of the world, i.e. all uncertainty is modeled by probabilities.">BM</a>* <a class="code" href="classbdm_1_1MPF.html#82b5a34d9ed0e78452f98d2ecbf1e93c" title="Access function.">_BM</a>(<span class="keywordtype">int</span> i){<span class="keywordflow">return</span> Bms[i];} 
     120<a name="l00139"></a><a class="code" href="classbdm_1_1MPF.html#82b5a34d9ed0e78452f98d2ecbf1e93c">00139</a>         <a class="code" href="classbdm_1_1BM.html" title="Bayesian Model of a system, i.e. all uncertainty is modeled by probabilities.">BM</a>* <a class="code" href="classbdm_1_1MPF.html#82b5a34d9ed0e78452f98d2ecbf1e93c" title="Access function.">_BM</a>(<span class="keywordtype">int</span> i){<span class="keywordflow">return</span> Bms[i];} 
    121121<a name="l00140"></a>00140 }; 
    122122<a name="l00141"></a>00141  
    … …  
    132132<a name="l00151"></a>00151 <span class="preprocessor"></span>        <span class="keywordflow">for</span> ( i=0;i&lt;<a class="code" href="classbdm_1_1PF.html#eeafaf9b8ad75fe62ee9fd6369e3f7fe" title="number of particles;">n</a>;i++ ) {    
    133133<a name="l00152"></a>00152                 <span class="comment">//generate new samples from paramater evolution model;</span> 
    134 <a name="l00153"></a>00153                 <a class="code" href="classbdm_1_1PF.html#914bd66025692c4018dbd482cb3c47c1" title="pointer into eEmp ">_samples</a> ( i ) = <a class="code" href="classbdm_1_1PF.html#cf3a1b2a407012e47ac878e3aa2fbf34" title="Parameter evolution model.">par</a>.<a class="code" href="classbdm_1_1mpdf.html#e4848a428d8ef0549c6e4a9ed386d9f2" title="Returns a sample from the density conditioned on cond, .">samplecond</a> ( <a class="code" href="classbdm_1_1PF.html#914bd66025692c4018dbd482cb3c47c1" title="pointer into eEmp ">_samples</a> ( i ), llsP ( i ) ); 
    135 <a name="l00154"></a>00154                 Bms[i]-&gt;condition ( <a class="code" href="classbdm_1_1PF.html#914bd66025692c4018dbd482cb3c47c1" title="pointer into eEmp ">_samples</a> ( i ) ); 
    136 <a name="l00155"></a>00155                 Bms[i]-&gt;bayes ( dt ); 
    137 <a name="l00156"></a>00156                 lls ( i ) = Bms[i]-&gt;_ll(); <span class="comment">// lls above is also in proposal her must be lls(i) =, not +=!!</span> 
    138 <a name="l00157"></a>00157                 <span class="keywordflow">if</span> ( lls ( i ) &gt;mlls ) mlls=lls ( i ); <span class="comment">//find maximum likelihood (for numerical stability)</span> 
    139 <a name="l00158"></a>00158         } 
    140 <a name="l00159"></a>00159  
    141 <a name="l00160"></a>00160         <span class="keywordtype">double</span> sum_w=0.0; 
    142 <a name="l00161"></a>00161         <span class="comment">// compute weights</span> 
    143 <a name="l00162"></a>00162 <span class="preprocessor">        #pragma omp parallel for</span> 
    144 <a name="l00163"></a>00163 <span class="preprocessor"></span>        <span class="keywordflow">for</span> ( i=0;i&lt;n;i++ ) { 
    145 <a name="l00164"></a>00164                 <a class="code" href="classbdm_1_1PF.html#f5149d5522d1095d39240c4c607f61a3" title="pointer into eEmp ">_w</a> ( i ) *= exp ( lls ( i ) - mlls ); <span class="comment">// multiply w by likelihood</span> 
    146 <a name="l00165"></a>00165                 sum_w+=<a class="code" href="classbdm_1_1PF.html#f5149d5522d1095d39240c4c607f61a3" title="pointer into eEmp ">_w</a>(i); 
    147 <a name="l00166"></a>00166         } 
    148 <a name="l00167"></a>00167  
    149 <a name="l00168"></a>00168         <span class="keywordflow">if</span> ( sum_w  &gt;0.0 ) { 
    150 <a name="l00169"></a>00169                 <a class="code" href="classbdm_1_1PF.html#f5149d5522d1095d39240c4c607f61a3" title="pointer into eEmp ">_w</a> /=sum_w; <span class="comment">//?</span> 
    151 <a name="l00170"></a>00170         } <span class="keywordflow">else</span> { 
    152 <a name="l00171"></a>00171                 cout&lt;&lt;<span class="stringliteral">"sum(w)==0"</span>&lt;&lt;endl; 
    153 <a name="l00172"></a>00172         } 
    154 <a name="l00173"></a>00173  
     134<a name="l00153"></a>00153                 <a class="code" href="classbdm_1_1PF.html#914bd66025692c4018dbd482cb3c47c1" title="pointer into eEmp ">_samples</a> ( i ) = <a class="code" href="classbdm_1_1PF.html#cf3a1b2a407012e47ac878e3aa2fbf34" title="Parameter evolution model.">par</a>.<a class="code" href="classbdm_1_1mpdf.html#f0c1db6fcbb3aae2dd6123884457a367" title="Returns a sample from the density conditioned on cond, .">samplecond</a> ( <a class="code" href="classbdm_1_1PF.html#914bd66025692c4018dbd482cb3c47c1" title="pointer into eEmp ">_samples</a> ( i ) ); 
     135<a name="l00154"></a>00154                 llsP ( i ) = <a class="code" href="classbdm_1_1PF.html#cf3a1b2a407012e47ac878e3aa2fbf34" title="Parameter evolution model.">par</a>.<a class="code" href="classbdm_1_1mpdf.html#05e843fd11c410a99dad2b88c55aca80" title="access function">_e</a>()-&gt;<a class="code" href="classbdm_1_1epdf.html#a8e39e2ff5e5cf5471bc159530d3b2d3" title="Compute log-probability of argument val.">evallog</a>(<a class="code" href="classbdm_1_1PF.html#914bd66025692c4018dbd482cb3c47c1" title="pointer into eEmp ">_samples</a>(i)); 
     136<a name="l00155"></a>00155                 Bms[i]-&gt;condition ( <a class="code" href="classbdm_1_1PF.html#914bd66025692c4018dbd482cb3c47c1" title="pointer into eEmp ">_samples</a> ( i ) ); 
     137<a name="l00156"></a>00156                 Bms[i]-&gt;bayes ( dt ); 
     138<a name="l00157"></a>00157                 lls ( i ) = Bms[i]-&gt;_ll(); <span class="comment">// lls above is also in proposal her must be lls(i) =, not +=!!</span> 
     139<a name="l00158"></a>00158                 <span class="keywordflow">if</span> ( lls ( i ) &gt;mlls ) mlls=lls ( i ); <span class="comment">//find maximum likelihood (for numerical stability)</span> 
     140<a name="l00159"></a>00159         } 
     141<a name="l00160"></a>00160  
     142<a name="l00161"></a>00161         <span class="keywordtype">double</span> sum_w=0.0; 
     143<a name="l00162"></a>00162         <span class="comment">// compute weights</span> 
     144<a name="l00163"></a>00163 <span class="preprocessor">        #pragma omp parallel for</span> 
     145<a name="l00164"></a>00164 <span class="preprocessor"></span>        <span class="keywordflow">for</span> ( i=0;i&lt;n;i++ ) { 
     146<a name="l00165"></a>00165                 <a class="code" href="classbdm_1_1PF.html#f5149d5522d1095d39240c4c607f61a3" title="pointer into eEmp ">_w</a> ( i ) *= exp ( lls ( i ) - mlls ); <span class="comment">// multiply w by likelihood</span> 
     147<a name="l00166"></a>00166                 sum_w+=<a class="code" href="classbdm_1_1PF.html#f5149d5522d1095d39240c4c607f61a3" title="pointer into eEmp ">_w</a>(i); 
     148<a name="l00167"></a>00167         } 
     149<a name="l00168"></a>00168  
     150<a name="l00169"></a>00169         <span class="keywordflow">if</span> ( sum_w  &gt;0.0 ) { 
     151<a name="l00170"></a>00170                 <a class="code" href="classbdm_1_1PF.html#f5149d5522d1095d39240c4c607f61a3" title="pointer into eEmp ">_w</a> /=sum_w; <span class="comment">//?</span> 
     152<a name="l00171"></a>00171         } <span class="keywordflow">else</span> { 
     153<a name="l00172"></a>00172                 cout&lt;&lt;<span class="stringliteral">"sum(w)==0"</span>&lt;&lt;endl; 
     154<a name="l00173"></a>00173         } 
    155155<a name="l00174"></a>00174  
    156 <a name="l00175"></a>00175         <span class="keywordtype">double</span> eff = 1.0/ ( <a class="code" href="classbdm_1_1PF.html#f5149d5522d1095d39240c4c607f61a3" title="pointer into eEmp ">_w</a>*<a class="code" href="classbdm_1_1PF.html#f5149d5522d1095d39240c4c607f61a3" title="pointer into eEmp ">_w</a> ); 
    157 <a name="l00176"></a>00176         <span class="keywordflow">if</span> ( eff &lt; ( 0.3*n ) ) { 
    158 <a name="l00177"></a>00177                 ind = <a class="code" href="classbdm_1_1PF.html#dc049265b9086cad7071f98d00a2b9af" title="posterior density">est</a>.<a class="code" href="classbdm_1_1eEmp.html#f06ce255de5dbb2313f52ee51f82ba3d" title="Function performs resampling, i.e. removal of low-weight samples and duplication...">resample</a>(); 
    159 <a name="l00178"></a>00178                 <span class="comment">// Resample Bms!</span> 
    160 <a name="l00179"></a>00179  
    161 <a name="l00180"></a>00180 <span class="preprocessor">                #pragma omp parallel for</span> 
    162 <a name="l00181"></a>00181 <span class="preprocessor"></span>                <span class="keywordflow">for</span> ( i=0;i&lt;n;i++ ) { 
    163 <a name="l00182"></a>00182                         <span class="keywordflow">if</span> ( ind ( i ) !=i ) {<span class="comment">//replace the current Bm by a new one</span> 
    164 <a name="l00183"></a>00183                                 <span class="comment">//fixme this would require new assignment operator</span> 
    165 <a name="l00184"></a>00184                                 <span class="comment">// *Bms[i] = *Bms[ind ( i ) ];</span> 
    166 <a name="l00185"></a>00185  
    167 <a name="l00186"></a>00186                                 <span class="comment">// poor-man's solution: replicate constructor here</span> 
    168 <a name="l00187"></a>00187                                 <span class="comment">// copied from MPF::MPF</span> 
    169 <a name="l00188"></a>00188                                 <span class="keyword">delete</span> Bms[i]; 
    170 <a name="l00189"></a>00189                                 Bms[i] = <span class="keyword">new</span> BM_T ( *Bms[ind ( i ) ] ); <span class="comment">//copy constructor</span> 
    171 <a name="l00190"></a>00190                                 <span class="keyword">const</span> <a class="code" href="classbdm_1_1epdf.html" title="Probability density function with numerical statistics, e.g. posterior density.">epdf</a>&amp; pom=Bms[i]-&gt;_epdf(); 
    172 <a name="l00191"></a>00191                                 jest.set_elements ( i,1.0/n,&amp;pom ); 
    173 <a name="l00192"></a>00192                         } 
    174 <a name="l00193"></a>00193                 }; 
    175 <a name="l00194"></a>00194                 cout &lt;&lt; <span class="charliteral">'.'</span>; 
    176 <a name="l00195"></a>00195         } 
    177 <a name="l00196"></a>00196 } 
    178 <a name="l00197"></a>00197  
    179 <a name="l00198"></a>00198 } 
    180 <a name="l00199"></a>00199 <span class="preprocessor">#endif // KF_H</span> 
    181 <a name="l00200"></a>00200 <span class="preprocessor"></span> 
    182 <a name="l00201"></a>00201  
     156<a name="l00175"></a>00175  
     157<a name="l00176"></a>00176         <span class="keywordtype">double</span> eff = 1.0/ ( <a class="code" href="classbdm_1_1PF.html#f5149d5522d1095d39240c4c607f61a3" title="pointer into eEmp ">_w</a>*<a class="code" href="classbdm_1_1PF.html#f5149d5522d1095d39240c4c607f61a3" title="pointer into eEmp ">_w</a> ); 
     158<a name="l00177"></a>00177         <span class="keywordflow">if</span> ( eff &lt; ( 0.3*n ) ) { 
     159<a name="l00178"></a>00178                 ind = <a class="code" href="classbdm_1_1PF.html#dc049265b9086cad7071f98d00a2b9af" title="posterior density">est</a>.<a class="code" href="classbdm_1_1eEmp.html#f06ce255de5dbb2313f52ee51f82ba3d" title="Function performs resampling, i.e. removal of low-weight samples and duplication...">resample</a>(); 
     160<a name="l00179"></a>00179                 <span class="comment">// Resample Bms!</span> 
     161<a name="l00180"></a>00180  
     162<a name="l00181"></a>00181 <span class="preprocessor">                #pragma omp parallel for</span> 
     163<a name="l00182"></a>00182 <span class="preprocessor"></span>                <span class="keywordflow">for</span> ( i=0;i&lt;n;i++ ) { 
     164<a name="l00183"></a>00183                         <span class="keywordflow">if</span> ( ind ( i ) !=i ) {<span class="comment">//replace the current Bm by a new one</span> 
     165<a name="l00184"></a>00184                                 <span class="comment">//fixme this would require new assignment operator</span> 
     166<a name="l00185"></a>00185                                 <span class="comment">// *Bms[i] = *Bms[ind ( i ) ];</span> 
     167<a name="l00186"></a>00186  
     168<a name="l00187"></a>00187                                 <span class="comment">// poor-man's solution: replicate constructor here</span> 
     169<a name="l00188"></a>00188                                 <span class="comment">// copied from MPF::MPF</span> 
     170<a name="l00189"></a>00189                                 <span class="keyword">delete</span> Bms[i]; 
     171<a name="l00190"></a>00190                                 Bms[i] = <span class="keyword">new</span> BM_T ( *Bms[ind ( i ) ] ); <span class="comment">//copy constructor</span> 
     172<a name="l00191"></a>00191                                 <span class="keyword">const</span> <a class="code" href="classbdm_1_1epdf.html" title="Probability density function with numerical statistics, e.g. posterior density.">epdf</a>&amp; pom=Bms[i]-&gt;_epdf(); 
     173<a name="l00192"></a>00192                                 jest.set_elements ( i,1.0/n,&amp;pom ); 
     174<a name="l00193"></a>00193                         } 
     175<a name="l00194"></a>00194                 }; 
     176<a name="l00195"></a>00195                 cout &lt;&lt; <span class="charliteral">'.'</span>; 
     177<a name="l00196"></a>00196         } 
     178<a name="l00197"></a>00197 } 
     179<a name="l00198"></a>00198  
     180<a name="l00199"></a>00199 } 
     181<a name="l00200"></a>00200 <span class="preprocessor">#endif // KF_H</span> 
     182<a name="l00201"></a>00201 <span class="preprocessor"></span> 
     183<a name="l00202"></a>00202  
    183184</pre></div></div> 
    184 <hr size="1"><address style="text-align: right;"><small>Generated on Fri Feb 6 19:49:46 2009 for mixpp by&nbsp; 
     185<hr size="1"><address style="text-align: right;"><small>Generated on Wed Feb 11 10:20:05 2009 for mixpp by&nbsp; 
    185186<a href="http://www.doxygen.org/index.html"> 
    186187<img src="doxygen.png" alt="doxygen" align="middle" border="0"></a> 1.5.6 </small></address>