Show
Ignore:
Timestamp:
02/24/09 14:12:46 (16 years ago)
Author:
smidl
Message:

progress...

Files:
1 modified

Legend:

Unmodified
Added
Removed

  • doc/html/libPF_8h-source.html

    r275 r280  
    8080<a name="l00047"></a>00047                         { <a class="code" href="classbdm_1_1PF.html#521e9621d3b5e1274275f323691afdaf" title="Parameter evolution model.">par</a> = par0; <a class="code" href="classbdm_1_1PF.html#d6e7a62fba1e0a0d73c9b87f4fb683ec" title="Observation model.">obs</a>=obs0; <a class="code" href="classbdm_1_1PF.html#dc049265b9086cad7071f98d00a2b9af" title="posterior density">est</a>.<a class="code" href="classbdm_1_1eEmp.html#82320074a9b0ad7e1bb33a6e885b65d7" title="Set samples and weights.">set_parameters</a> ( ones(n0),epdf0 ); }; 
    8181<a name="l00048"></a>00048         <span class="keywordtype">void</span> set_parameters ( mpdf *par0, mpdf *obs0, <span class="keywordtype">int</span> n0 )  
    82 <a name="l00049"></a>00049                         { <a class="code" href="classbdm_1_1PF.html#521e9621d3b5e1274275f323691afdaf" title="Parameter evolution model.">par</a> = par0; <a class="code" href="classbdm_1_1PF.html#d6e7a62fba1e0a0d73c9b87f4fb683ec" title="Observation model.">obs</a>=obs0; <a class="code" href="classbdm_1_1PF.html#eeafaf9b8ad75fe62ee9fd6369e3f7fe" title="number of particles;">n</a>=n0; }; 
     82<a name="l00049"></a>00049                         { <a class="code" href="classbdm_1_1PF.html#521e9621d3b5e1274275f323691afdaf" title="Parameter evolution model.">par</a> = par0; <a class="code" href="classbdm_1_1PF.html#d6e7a62fba1e0a0d73c9b87f4fb683ec" title="Observation model.">obs</a>=obs0; <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>.set_n(<a class="code" href="classbdm_1_1PF.html#eeafaf9b8ad75fe62ee9fd6369e3f7fe" title="number of particles;">n</a>);}; 
    8383<a name="l00050"></a>00050         <span class="keywordtype">void</span> set_statistics (<span class="keyword">const</span> vec w0, epdf *epdf0){<a class="code" href="classbdm_1_1PF.html#dc049265b9086cad7071f98d00a2b9af" title="posterior density">est</a>.set_parameters ( w0,epdf0 );}; 
    8484<a name="l00053"></a>00053         <span class="keywordtype">void</span> <a class="code" href="classbdm_1_1PF.html#6f1988db4c3f602d187a6c15ec89cb1e">set_est</a> ( <span class="keyword">const</span> epdf &amp;epdf0 ); 
    … …  
    103103<a name="l00080"></a>00080                                 Coms ( <a class="code" href="classbdm_1_1PF.html#f5149d5522d1095d39240c4c607f61a3" title="pointer into eEmp ">_w</a>.length() ) { 
    104104<a name="l00081"></a>00081                 }; 
    105 <a name="l00082"></a>00082  
    106 <a name="l00083"></a>00083                 <span class="keywordtype">void</span> set_elements ( <span class="keywordtype">int</span> &amp;i, <span class="keywordtype">double</span> wi, <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>* ep ) 
    107 <a name="l00084"></a>00084                 {<a class="code" href="classbdm_1_1PF.html#f5149d5522d1095d39240c4c607f61a3" title="pointer into eEmp ">_w</a> ( i ) =wi; Coms ( i ) =ep;}; 
    108 <a name="l00085"></a>00085  
     105<a name="l00082"></a>00082                 <span class="keywordtype">void</span> set_elements ( <span class="keywordtype">int</span> &amp;i, <span class="keywordtype">double</span> wi, <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>* ep ) 
     106<a name="l00083"></a>00083                 {<a class="code" href="classbdm_1_1PF.html#f5149d5522d1095d39240c4c607f61a3" title="pointer into eEmp ">_w</a> ( i ) =wi; Coms ( i ) =ep;}; 
     107<a name="l00084"></a>00084  
     108<a name="l00085"></a>00085                 <span class="keywordtype">void</span> set_n(<span class="keywordtype">int</span> <a class="code" href="classbdm_1_1PF.html#eeafaf9b8ad75fe62ee9fd6369e3f7fe" title="number of particles;">n</a>){E.set_n(n); Coms.set_length(n);} 
    109109<a name="l00086"></a>00086                 vec mean()<span class="keyword"> const </span>{ 
    110110<a name="l00087"></a>00087                         <span class="comment">// ugly</span> 
    … …  
    134134<a name="l00113"></a><a class="code" href="classbdm_1_1MPF.html#ca0e773df05fd70cf8ef3a7f1b3e42ef">00113</a>         <a class="code" href="classbdm_1_1MPF.html#ca0e773df05fd70cf8ef3a7f1b3e42ef" title="Default constructor.">MPF</a> ( <a class="code" href="classbdm_1_1mpdf.html" title="Conditional probability density, e.g. modeling some dependencies.">mpdf</a> *par0, <a class="code" href="classbdm_1_1mpdf.html" title="Conditional probability density, e.g. modeling some dependencies.">mpdf</a> *obs0, <span class="keywordtype">int</span> <a class="code" href="classbdm_1_1PF.html#eeafaf9b8ad75fe62ee9fd6369e3f7fe" title="number of particles;">n</a>, <span class="keyword">const</span> BM_T &amp;BMcond0 ) : <a class="code" href="classbdm_1_1PF.html" title="Trivial particle filter with proposal density equal to parameter evolution model...">PF</a> (), jest ( <a class="code" href="classbdm_1_1PF.html#dc049265b9086cad7071f98d00a2b9af" title="posterior density">est</a> ) { 
    135135<a name="l00114"></a>00114                 PF::set_parameters(par0,obs0,n); 
    136 <a name="l00115"></a>00115                 <span class="comment">//</span> 
    137 <a name="l00116"></a>00116                 <span class="comment">//TODO test if rv and BMcond.rv are compatible.</span> 
    138 <a name="l00117"></a>00117 <span class="comment">//              rv.add ( rvlin );</span> 
    139 <a name="l00118"></a>00118                 <span class="comment">//</span> 
    140 <a name="l00119"></a>00119  
    141 <a name="l00120"></a>00120                 <span class="keywordflow">if</span> ( n&gt;10000 ) {it_error ( <span class="stringliteral">"increase 10000 here!"</span> );} 
    142 <a name="l00121"></a>00121  
    143 <a name="l00122"></a>00122                 <span class="keywordflow">for</span> ( <span class="keywordtype">int</span> i=0;i&lt;n;i++ ) { 
    144 <a name="l00123"></a>00123                         Bms[i] = <span class="keyword">new</span> BM_T ( BMcond0 ); <span class="comment">//copy constructor</span> 
    145 <a name="l00124"></a>00124                         <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;posterior(); 
    146 <a name="l00125"></a>00125                         jest.set_elements ( i,1.0/n,&amp;pom ); 
    147 <a name="l00126"></a>00126                 } 
    148 <a name="l00127"></a>00127         }; 
    149 <a name="l00128"></a>00128  
    150 <a name="l00129"></a>00129         ~<a class="code" href="classbdm_1_1MPF.html" title="Marginalized Particle filter.">MPF</a>() { 
    151 <a name="l00130"></a>00130         } 
    152 <a name="l00131"></a>00131  
    153 <a name="l00132"></a>00132         <span class="keywordtype">void</span> <a class="code" href="classbdm_1_1MPF.html#286d040770d08bd7ff416cea617b1b14" title="Incremental Bayes rule.">bayes</a> ( <span class="keyword">const</span> vec &amp;dt ); 
    154 <a name="l00133"></a>00133         <span class="keyword">const</span> epdf&amp; posterior()<span class="keyword"> const </span>{<span class="keywordflow">return</span> jest;} 
    155 <a name="l00134"></a>00134         <span class="keyword">const</span> epdf* _e()<span class="keyword"> const </span>{<span class="keywordflow">return</span> &amp;jest;} <span class="comment">//Fixme: is it useful?</span> 
    156 <a name="l00136"></a><a class="code" href="classbdm_1_1MPF.html#dcecdaf2acbbee51acf3018a70989a7e">00136</a> <span class="comment"></span>        <span class="keywordtype">void</span> <a class="code" href="classbdm_1_1MPF.html#dcecdaf2acbbee51acf3018a70989a7e" title="Set postrior of rvc to samples from epdf0. Statistics of Bms are not re-computed!...">set_est</a> ( <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; epdf0 ) { 
    157 <a name="l00137"></a>00137                 <a class="code" href="classbdm_1_1PF.html#6f1988db4c3f602d187a6c15ec89cb1e">PF::set_est</a> ( epdf0 );  <span class="comment">// sample params in condition</span> 
    158 <a name="l00138"></a>00138                 <span class="comment">// copy conditions to BMs</span> 
    159 <a name="l00139"></a>00139  
    160 <a name="l00140"></a>00140                 <span class="keywordflow">for</span> ( <span class="keywordtype">int</span> i=0;i&lt;<a class="code" href="classbdm_1_1PF.html#eeafaf9b8ad75fe62ee9fd6369e3f7fe" title="number of particles;">n</a>;i++ ) {Bms[i]-&gt;condition ( <a class="code" href="classbdm_1_1PF.html#914bd66025692c4018dbd482cb3c47c1" title="pointer into eEmp ">_samples</a> ( i ) );} 
    161 <a name="l00141"></a>00141         } 
    162 <a name="l00142"></a>00142  
    163 <a name="l00144"></a><a class="code" href="classbdm_1_1MPF.html#82b5a34d9ed0e78452f98d2ecbf1e93c">00144</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];} 
    164 <a name="l00145"></a>00145 }; 
    165 <a name="l00146"></a>00146  
    166 <a name="l00147"></a>00147 <span class="keyword">template</span>&lt;<span class="keyword">class</span> BM_T&gt; 
    167 <a name="l00148"></a><a class="code" href="classbdm_1_1MPF.html#286d040770d08bd7ff416cea617b1b14">00148</a> <span class="keywordtype">void</span> <a class="code" href="classbdm_1_1MPF.html" title="Marginalized Particle filter.">MPF&lt;BM_T&gt;::bayes</a> ( <span class="keyword">const</span> vec &amp;dt ) { 
    168 <a name="l00149"></a>00149         <span class="keywordtype">int</span> i; 
    169 <a name="l00150"></a>00150         vec lls ( <a class="code" href="classbdm_1_1PF.html#eeafaf9b8ad75fe62ee9fd6369e3f7fe" title="number of particles;">n</a> ); 
    170 <a name="l00151"></a>00151         vec llsP ( <a class="code" href="classbdm_1_1PF.html#eeafaf9b8ad75fe62ee9fd6369e3f7fe" title="number of particles;">n</a> ); 
    171 <a name="l00152"></a>00152         ivec ind; 
    172 <a name="l00153"></a>00153         <span class="keywordtype">double</span> mlls=-std::numeric_limits&lt;double&gt;::infinity(); 
    173 <a name="l00154"></a>00154  
    174 <a name="l00155"></a>00155 <span class="preprocessor">#pragma omp parallel for</span> 
    175 <a name="l00156"></a>00156 <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++ ) { 
    176 <a name="l00157"></a>00157                 <span class="comment">//generate new samples from paramater evolution model;</span> 
    177 <a name="l00158"></a>00158                 <a class="code" href="classbdm_1_1PF.html#914bd66025692c4018dbd482cb3c47c1" title="pointer into eEmp ">_samples</a> ( i ) = <a class="code" href="classbdm_1_1PF.html#521e9621d3b5e1274275f323691afdaf" title="Parameter evolution model.">par</a>-&gt;<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 ) ); 
    178 <a name="l00159"></a>00159                 llsP ( i ) = <a class="code" href="classbdm_1_1PF.html#521e9621d3b5e1274275f323691afdaf" title="Parameter evolution model.">par</a>-&gt;<a class="code" href="classbdm_1_1mpdf.html#05e843fd11c410a99dad2b88c55aca80">_e</a>()-&gt;<a class="code" href="classbdm_1_1epdf.html#deab266d63c236c277538867d5c3f249" 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 ) ); 
    179 <a name="l00160"></a>00160                 Bms[i]-&gt;condition ( <a class="code" href="classbdm_1_1PF.html#914bd66025692c4018dbd482cb3c47c1" title="pointer into eEmp ">_samples</a> ( i ) ); 
    180 <a name="l00161"></a>00161                 Bms[i]-&gt;bayes ( dt ); 
    181 <a name="l00162"></a>00162                 lls ( i ) = Bms[i]-&gt;_ll(); <span class="comment">// lls above is also in proposal her must be lls(i) =, not +=!!</span> 
    182 <a name="l00163"></a>00163                 <span class="keywordflow">if</span> ( lls ( i ) &gt;mlls ) mlls=lls ( i ); <span class="comment">//find maximum likelihood (for numerical stability)</span> 
    183 <a name="l00164"></a>00164         } 
    184 <a name="l00165"></a>00165  
    185 <a name="l00166"></a>00166         <span class="keywordtype">double</span> sum_w=0.0; 
    186 <a name="l00167"></a>00167         <span class="comment">// compute weights</span> 
    187 <a name="l00168"></a>00168 <span class="preprocessor">#pragma omp parallel for</span> 
    188 <a name="l00169"></a>00169 <span class="preprocessor"></span>        <span class="keywordflow">for</span> ( i=0;i&lt;n;i++ ) { 
    189 <a name="l00170"></a>00170                 <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> 
    190 <a name="l00171"></a>00171                 sum_w+=<a class="code" href="classbdm_1_1PF.html#f5149d5522d1095d39240c4c607f61a3" title="pointer into eEmp ">_w</a> ( i ); 
    191 <a name="l00172"></a>00172         } 
    192 <a name="l00173"></a>00173  
    193 <a name="l00174"></a>00174         <span class="keywordflow">if</span> ( sum_w  &gt;0.0 ) { 
    194 <a name="l00175"></a>00175                 <a class="code" href="classbdm_1_1PF.html#f5149d5522d1095d39240c4c607f61a3" title="pointer into eEmp ">_w</a> /=sum_w; <span class="comment">//?</span> 
    195 <a name="l00176"></a>00176         } 
    196 <a name="l00177"></a>00177         <span class="keywordflow">else</span> { 
    197 <a name="l00178"></a>00178                 cout&lt;&lt;<span class="stringliteral">"sum(w)==0"</span>&lt;&lt;endl; 
    198 <a name="l00179"></a>00179         } 
    199 <a name="l00180"></a>00180  
     136<a name="l00115"></a>00115                 jest.set_n(n); 
     137<a name="l00116"></a>00116                 <span class="comment">//</span> 
     138<a name="l00117"></a>00117                 <span class="comment">//TODO test if rv and BMcond.rv are compatible.</span> 
     139<a name="l00118"></a>00118 <span class="comment">//              rv.add ( rvlin );</span> 
     140<a name="l00119"></a>00119                 <span class="comment">//</span> 
     141<a name="l00120"></a>00120  
     142<a name="l00121"></a>00121                 <span class="keywordflow">if</span> ( n&gt;10000 ) {it_error ( <span class="stringliteral">"increase 10000 here!"</span> );} 
     143<a name="l00122"></a>00122  
     144<a name="l00123"></a>00123                 <span class="keywordflow">for</span> ( <span class="keywordtype">int</span> i=0;i&lt;n;i++ ) { 
     145<a name="l00124"></a>00124                         Bms[i] = <span class="keyword">new</span> BM_T ( BMcond0 ); <span class="comment">//copy constructor</span> 
     146<a name="l00125"></a>00125                         <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;posterior(); 
     147<a name="l00126"></a>00126                         jest.set_elements ( i,1.0/n,&amp;pom ); 
     148<a name="l00127"></a>00127                 } 
     149<a name="l00128"></a>00128         }; 
     150<a name="l00129"></a>00129  
     151<a name="l00130"></a>00130         ~<a class="code" href="classbdm_1_1MPF.html" title="Marginalized Particle filter.">MPF</a>() { 
     152<a name="l00131"></a>00131         } 
     153<a name="l00132"></a>00132  
     154<a name="l00133"></a>00133         <span class="keywordtype">void</span> <a class="code" href="classbdm_1_1MPF.html#286d040770d08bd7ff416cea617b1b14" title="Incremental Bayes rule.">bayes</a> ( <span class="keyword">const</span> vec &amp;dt ); 
     155<a name="l00134"></a>00134         <span class="keyword">const</span> epdf&amp; posterior()<span class="keyword"> const </span>{<span class="keywordflow">return</span> jest;} 
     156<a name="l00135"></a>00135         <span class="keyword">const</span> epdf* _e()<span class="keyword"> const </span>{<span class="keywordflow">return</span> &amp;jest;} <span class="comment">//Fixme: is it useful?</span> 
     157<a name="l00137"></a><a class="code" href="classbdm_1_1MPF.html#dcecdaf2acbbee51acf3018a70989a7e">00137</a> <span class="comment"></span>        <span class="keywordtype">void</span> <a class="code" href="classbdm_1_1MPF.html#dcecdaf2acbbee51acf3018a70989a7e" title="Set postrior of rvc to samples from epdf0. Statistics of Bms are not re-computed!...">set_est</a> ( <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; epdf0 ) { 
     158<a name="l00138"></a>00138                 <a class="code" href="classbdm_1_1PF.html#6f1988db4c3f602d187a6c15ec89cb1e">PF::set_est</a> ( epdf0 );  <span class="comment">// sample params in condition</span> 
     159<a name="l00139"></a>00139                 <span class="comment">// copy conditions to BMs</span> 
     160<a name="l00140"></a>00140  
     161<a name="l00141"></a>00141                 <span class="keywordflow">for</span> ( <span class="keywordtype">int</span> i=0;i&lt;<a class="code" href="classbdm_1_1PF.html#eeafaf9b8ad75fe62ee9fd6369e3f7fe" title="number of particles;">n</a>;i++ ) {Bms[i]-&gt;condition ( <a class="code" href="classbdm_1_1PF.html#914bd66025692c4018dbd482cb3c47c1" title="pointer into eEmp ">_samples</a> ( i ) );} 
     162<a name="l00142"></a>00142         } 
     163<a name="l00143"></a>00143  
     164<a name="l00145"></a><a class="code" href="classbdm_1_1MPF.html#82b5a34d9ed0e78452f98d2ecbf1e93c">00145</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];} 
     165<a name="l00146"></a>00146 }; 
     166<a name="l00147"></a>00147  
     167<a name="l00148"></a>00148 <span class="keyword">template</span>&lt;<span class="keyword">class</span> BM_T&gt; 
     168<a name="l00149"></a><a class="code" href="classbdm_1_1MPF.html#286d040770d08bd7ff416cea617b1b14">00149</a> <span class="keywordtype">void</span> <a class="code" href="classbdm_1_1MPF.html" title="Marginalized Particle filter.">MPF&lt;BM_T&gt;::bayes</a> ( <span class="keyword">const</span> vec &amp;dt ) { 
     169<a name="l00150"></a>00150         <span class="keywordtype">int</span> i; 
     170<a name="l00151"></a>00151         vec lls ( <a class="code" href="classbdm_1_1PF.html#eeafaf9b8ad75fe62ee9fd6369e3f7fe" title="number of particles;">n</a> ); 
     171<a name="l00152"></a>00152         vec llsP ( <a class="code" href="classbdm_1_1PF.html#eeafaf9b8ad75fe62ee9fd6369e3f7fe" title="number of particles;">n</a> ); 
     172<a name="l00153"></a>00153         ivec ind; 
     173<a name="l00154"></a>00154         <span class="keywordtype">double</span> mlls=-std::numeric_limits&lt;double&gt;::infinity(); 
     174<a name="l00155"></a>00155  
     175<a name="l00156"></a>00156 <span class="preprocessor">#pragma omp parallel for</span> 
     176<a name="l00157"></a>00157 <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++ ) { 
     177<a name="l00158"></a>00158                 <span class="comment">//generate new samples from paramater evolution model;</span> 
     178<a name="l00159"></a>00159                 <a class="code" href="classbdm_1_1PF.html#914bd66025692c4018dbd482cb3c47c1" title="pointer into eEmp ">_samples</a> ( i ) = <a class="code" href="classbdm_1_1PF.html#521e9621d3b5e1274275f323691afdaf" title="Parameter evolution model.">par</a>-&gt;<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 ) ); 
     179<a name="l00160"></a>00160                 llsP ( i ) = <a class="code" href="classbdm_1_1PF.html#521e9621d3b5e1274275f323691afdaf" title="Parameter evolution model.">par</a>-&gt;<a class="code" href="classbdm_1_1mpdf.html#05e843fd11c410a99dad2b88c55aca80">_e</a>()-&gt;<a class="code" href="classbdm_1_1epdf.html#deab266d63c236c277538867d5c3f249" 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 ) ); 
     180<a name="l00161"></a>00161                 Bms[i]-&gt;condition ( <a class="code" href="classbdm_1_1PF.html#914bd66025692c4018dbd482cb3c47c1" title="pointer into eEmp ">_samples</a> ( i ) ); 
     181<a name="l00162"></a>00162                 Bms[i]-&gt;bayes ( dt ); 
     182<a name="l00163"></a>00163                 lls ( i ) = Bms[i]-&gt;_ll(); <span class="comment">// lls above is also in proposal her must be lls(i) =, not +=!!</span> 
     183<a name="l00164"></a>00164                 <span class="keywordflow">if</span> ( lls ( i ) &gt;mlls ) mlls=lls ( i ); <span class="comment">//find maximum likelihood (for numerical stability)</span> 
     184<a name="l00165"></a>00165         } 
     185<a name="l00166"></a>00166  
     186<a name="l00167"></a>00167         <span class="keywordtype">double</span> sum_w=0.0; 
     187<a name="l00168"></a>00168         <span class="comment">// compute weights</span> 
     188<a name="l00169"></a>00169 <span class="preprocessor">#pragma omp parallel for</span> 
     189<a name="l00170"></a>00170 <span class="preprocessor"></span>        <span class="keywordflow">for</span> ( i=0;i&lt;n;i++ ) { 
     190<a name="l00171"></a>00171                 <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> 
     191<a name="l00172"></a>00172                 sum_w+=<a class="code" href="classbdm_1_1PF.html#f5149d5522d1095d39240c4c607f61a3" title="pointer into eEmp ">_w</a> ( i ); 
     192<a name="l00173"></a>00173         } 
     193<a name="l00174"></a>00174  
     194<a name="l00175"></a>00175         <span class="keywordflow">if</span> ( sum_w  &gt;0.0 ) { 
     195<a name="l00176"></a>00176                 <a class="code" href="classbdm_1_1PF.html#f5149d5522d1095d39240c4c607f61a3" title="pointer into eEmp ">_w</a> /=sum_w; <span class="comment">//?</span> 
     196<a name="l00177"></a>00177         } 
     197<a name="l00178"></a>00178         <span class="keywordflow">else</span> { 
     198<a name="l00179"></a>00179                 cout&lt;&lt;<span class="stringliteral">"sum(w)==0"</span>&lt;&lt;endl; 
     199<a name="l00180"></a>00180         } 
    200200<a name="l00181"></a>00181  
    201 <a name="l00182"></a>00182         <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> ); 
    202 <a name="l00183"></a>00183         <span class="keywordflow">if</span> ( eff &lt; ( 0.3*n ) ) { 
    203 <a name="l00184"></a>00184                 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>(); 
    204 <a name="l00185"></a>00185                 <span class="comment">// Resample Bms!</span> 
    205 <a name="l00186"></a>00186  
    206 <a name="l00187"></a>00187 <span class="preprocessor">#pragma omp parallel for</span> 
    207 <a name="l00188"></a>00188 <span class="preprocessor"></span>                <span class="keywordflow">for</span> ( i=0;i&lt;n;i++ ) { 
    208 <a name="l00189"></a>00189                         <span class="keywordflow">if</span> ( ind ( i ) !=i ) {<span class="comment">//replace the current Bm by a new one</span> 
    209 <a name="l00190"></a>00190                                 <span class="comment">//fixme this would require new assignment operator</span> 
    210 <a name="l00191"></a>00191                                 <span class="comment">// *Bms[i] = *Bms[ind ( i ) ];</span> 
    211 <a name="l00192"></a>00192  
    212 <a name="l00193"></a>00193                                 <span class="comment">// poor-man's solution: replicate constructor here</span> 
    213 <a name="l00194"></a>00194                                 <span class="comment">// copied from MPF::MPF</span> 
    214 <a name="l00195"></a>00195                                 <span class="keyword">delete</span> Bms[i]; 
    215 <a name="l00196"></a>00196                                 Bms[i] = <span class="keyword">new</span> BM_T ( *Bms[ind ( i ) ] ); <span class="comment">//copy constructor</span> 
    216 <a name="l00197"></a>00197                                 <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;posterior(); 
    217 <a name="l00198"></a>00198                                 jest.set_elements ( i,1.0/n,&amp;pom ); 
    218 <a name="l00199"></a>00199                         } 
    219 <a name="l00200"></a>00200                 }; 
    220 <a name="l00201"></a>00201                 cout &lt;&lt; <span class="charliteral">'.'</span>; 
    221 <a name="l00202"></a>00202         } 
    222 <a name="l00203"></a>00203 } 
    223 <a name="l00204"></a>00204  
    224 <a name="l00205"></a>00205 } 
    225 <a name="l00206"></a>00206 <span class="preprocessor">#endif // KF_H</span> 
    226 <a name="l00207"></a>00207 <span class="preprocessor"></span> 
    227 <a name="l00208"></a>00208  
     201<a name="l00182"></a>00182  
     202<a name="l00183"></a>00183         <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> ); 
     203<a name="l00184"></a>00184         <span class="keywordflow">if</span> ( eff &lt; ( 0.3*n ) ) { 
     204<a name="l00185"></a>00185                 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>(); 
     205<a name="l00186"></a>00186                 <span class="comment">// Resample Bms!</span> 
     206<a name="l00187"></a>00187  
     207<a name="l00188"></a>00188 <span class="preprocessor">#pragma omp parallel for</span> 
     208<a name="l00189"></a>00189 <span class="preprocessor"></span>                <span class="keywordflow">for</span> ( i=0;i&lt;n;i++ ) { 
     209<a name="l00190"></a>00190                         <span class="keywordflow">if</span> ( ind ( i ) !=i ) {<span class="comment">//replace the current Bm by a new one</span> 
     210<a name="l00191"></a>00191                                 <span class="comment">//fixme this would require new assignment operator</span> 
     211<a name="l00192"></a>00192                                 <span class="comment">// *Bms[i] = *Bms[ind ( i ) ];</span> 
     212<a name="l00193"></a>00193  
     213<a name="l00194"></a>00194                                 <span class="comment">// poor-man's solution: replicate constructor here</span> 
     214<a name="l00195"></a>00195                                 <span class="comment">// copied from MPF::MPF</span> 
     215<a name="l00196"></a>00196                                 <span class="keyword">delete</span> Bms[i]; 
     216<a name="l00197"></a>00197                                 Bms[i] = <span class="keyword">new</span> BM_T ( *Bms[ind ( i ) ] ); <span class="comment">//copy constructor</span> 
     217<a name="l00198"></a>00198                                 <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;posterior(); 
     218<a name="l00199"></a>00199                                 jest.set_elements ( i,1.0/n,&amp;pom ); 
     219<a name="l00200"></a>00200                         } 
     220<a name="l00201"></a>00201                 }; 
     221<a name="l00202"></a>00202                 cout &lt;&lt; <span class="charliteral">'.'</span>; 
     222<a name="l00203"></a>00203         } 
     223<a name="l00204"></a>00204 } 
     224<a name="l00205"></a>00205  
     225<a name="l00206"></a>00206 } 
     226<a name="l00207"></a>00207 <span class="preprocessor">#endif // KF_H</span> 
     227<a name="l00208"></a>00208 <span class="preprocessor"></span> 
     228<a name="l00209"></a>00209  
    228229</pre></div></div> 
    229 <hr size="1"><address style="text-align: right;"><small>Generated on Mon Feb 16 10:06:24 2009 for mixpp by&nbsp; 
     230<hr size="1"><address style="text-align: right;"><small>Generated on Wed Feb 18 17:38:40 2009 for mixpp by&nbsp; 
    230231<a href="http://www.doxygen.org/index.html"> 
    231232<img src="doxygen.png" alt="doxygen" align="middle" border="0"></a> 1.5.6 </small></address>