Show
Ignore:
Timestamp:
03/03/08 13:00:32 (17 years ago)
Author:
smidl
Message:

test KF : estimation of R in KF is not possible! Likelihood of y_t is growing when R -> 0

Files:
1 modified

Legend:

Unmodified
Added
Removed
  • doc/html/libPF_8h-source.html

    r28 r32  
    1919<a name="l00015"></a>00015 <span class="preprocessor"></span> 
    2020<a name="l00016"></a>00016 <span class="preprocessor">#include &lt;itpp/itbase.h&gt;</span> 
    21 <a name="l00017"></a>00017 <span class="preprocessor">#include "../stat/libBM.h"</span> 
     21<a name="l00017"></a>00017 <span class="preprocessor">#include "../stat/libEF.h"</span> 
    2222<a name="l00018"></a>00018 <span class="preprocessor">#include "../math/libDC.h"</span> 
    2323<a name="l00019"></a>00019  
    2424<a name="l00020"></a>00020 <span class="keyword">using namespace </span>itpp; 
    2525<a name="l00021"></a>00021  
    26 <a name="l00022"></a>00022 <span class="keyword">enum</span> RESAMPLING_METHOD { MULTINOMIAL = 0, STRATIFIED = 1, SYSTEMATIC = 3 }; 
    27 <a name="l00023"></a>00023  
    28 <a name="l00029"></a><a class="code" href="classPF.html">00029</a> <span class="keyword">class </span><a class="code" href="classPF.html" title="A Particle Filter prototype.">PF</a> : <span class="keyword">public</span> <a class="code" href="classBM.html" title="Bayesian Model of the world, i.e. all uncertainty is modeled by probabilities.">BM</a> {  
    29 <a name="l00030"></a>00030 <span class="keyword">protected</span>: 
    30 <a name="l00031"></a>00031         <span class="keywordtype">int</span> n; <span class="comment">//number of particles</span> 
    31 <a name="l00032"></a>00032         vec w; <span class="comment">//particle weights</span> 
    32 <a name="l00033"></a>00033         Uniform_RNG URNG; <span class="comment">//used for resampling</span> 
    33 <a name="l00034"></a>00034          
    34 <a name="l00035"></a>00035 <span class="keyword">public</span>: 
    35 <a name="l00037"></a>00037         ivec <a class="code" href="classPF.html#a0e26b2f6a5884aca49122f3e4f0cf19" title="Returns indexes of particles that should be resampled. The ordering MUST guarantee...">resample</a>(RESAMPLING_METHOD method = SYSTEMATIC); 
    36 <a name="l00038"></a>00038         <a class="code" href="classPF.html" title="A Particle Filter prototype.">PF</a> (vec w); 
    37 <a name="l00039"></a>00039         <span class="comment">//TODO remove or implement bayes()!</span> 
    38 <a name="l00040"></a>00040         <span class="keywordtype">void</span> bayes(<span class="keyword">const</span> vec &amp;dt, <span class="keywordtype">bool</span> evell){}; 
    39 <a name="l00041"></a>00041 }; 
    40 <a name="l00042"></a>00042  
    41 <a name="l00049"></a><a class="code" href="classTrivialPF.html">00049</a> <span class="keyword">class </span><a class="code" href="classTrivialPF.html" title="Trivial particle filter with proposal density that is not conditioned on the data...">TrivialPF</a> : <span class="keyword">public</span> <a class="code" href="classPF.html" title="A Particle Filter prototype.">PF</a> { 
    42 <a name="l00050"></a>00050         Array&lt;vec&gt; ptcls; 
    43 <a name="l00051"></a>00051          
    44 <a name="l00052"></a>00052         <span class="keywordtype">bool</span> is_proposal; 
    45 <a name="l00053"></a>00053         <a class="code" href="classBM.html" title="Bayesian Model of the world, i.e. all uncertainty is modeled by probabilities.">BM</a> *prop; 
    46 <a name="l00054"></a>00054         <a class="code" href="classmpdf.html" title="Conditional probability density, e.g. modeling some dependencies.">mpdf</a> *par; 
    47 <a name="l00055"></a>00055         <a class="code" href="classmpdf.html" title="Conditional probability density, e.g. modeling some dependencies.">mpdf</a> *obs; 
    48 <a name="l00056"></a>00056          
    49 <a name="l00057"></a>00057         <span class="keyword">public</span>: 
    50 <a name="l00058"></a>00058         <a class="code" href="classTrivialPF.html" title="Trivial particle filter with proposal density that is not conditioned on the data...">TrivialPF</a>(<a class="code" href="classmpdf.html" title="Conditional probability density, e.g. modeling some dependencies.">mpdf</a> &amp;par, <a class="code" href="classmpdf.html" title="Conditional probability density, e.g. modeling some dependencies.">mpdf</a> &amp;obs, <a class="code" href="classBM.html" title="Bayesian Model of the world, i.e. all uncertainty is modeled by probabilities.">BM</a> &amp;prop, <span class="keywordtype">int</span> n0); 
    51 <a name="l00059"></a>00059         <a class="code" href="classTrivialPF.html" title="Trivial particle filter with proposal density that is not conditioned on the data...">TrivialPF</a>(<a class="code" href="classmpdf.html" title="Conditional probability density, e.g. modeling some dependencies.">mpdf</a> &amp;par, <a class="code" href="classmpdf.html" title="Conditional probability density, e.g. modeling some dependencies.">mpdf</a> &amp;obs, <span class="keywordtype">int</span> n0); 
    52 <a name="l00060"></a>00060         <span class="keywordtype">void</span> bayes(<span class="keyword">const</span> vec &amp;dt, <span class="keywordtype">bool</span> <a class="code" href="classBM.html#bf6fb59b30141074f8ee1e2f43d03129" title="If true, the filter will compute likelihood of the data record and store it in ll...">evalll</a>); 
    53 <a name="l00061"></a>00061 }; 
    54 <a name="l00062"></a>00062  
    55 <a name="l00063"></a>00063 <span class="keyword">class </span>MPF : <span class="keyword">public</span> <a class="code" href="classTrivialPF.html" title="Trivial particle filter with proposal density that is not conditioned on the data...">TrivialPF</a> { 
    56 <a name="l00064"></a>00064         Array&lt;BM&gt; Bms; 
    57 <a name="l00065"></a>00065         <span class="keyword">public</span>: 
    58 <a name="l00066"></a>00066         MPF(<a class="code" href="classBM.html" title="Bayesian Model of the world, i.e. all uncertainty is modeled by probabilities.">BM</a> &amp;B, <a class="code" href="classmpdf.html" title="Conditional probability density, e.g. modeling some dependencies.">mpdf</a> &amp;prop, <a class="code" href="classmpdf.html" title="Conditional probability density, e.g. modeling some dependencies.">mpdf</a> &amp;obs, <a class="code" href="classmpdf.html" title="Conditional probability density, e.g. modeling some dependencies.">mpdf</a> &amp;par); 
    59 <a name="l00067"></a>00067         <span class="keywordtype">void</span> bayes(vec &amp;dt);     
    60 <a name="l00068"></a>00068 }; 
    61 <a name="l00069"></a>00069  
    62 <a name="l00070"></a>00070 <span class="preprocessor">#endif // KF_H</span> 
    63 <a name="l00071"></a>00071 <span class="preprocessor"></span> 
    64 <a name="l00072"></a>00072  
    65 </pre></div><hr size="1"><address style="text-align: right;"><small>Generated on Mon Feb 18 21:48:39 2008 for mixpp by&nbsp; 
     26<a name="l00027"></a><a class="code" href="classPF.html">00027</a> <span class="keyword">class </span><a class="code" href="classPF.html" title="Trivial particle filter with proposal density equal to parameter evolution model...">PF</a> : <span class="keyword">public</span> <a class="code" href="classBM.html" title="Bayesian Model of the world, i.e. all uncertainty is modeled by probabilities.">BM</a> { 
     27<a name="l00028"></a>00028 <span class="keyword">protected</span>: 
     28<a name="l00030"></a><a class="code" href="classPF.html#2c2f44ed7a4eaa42e07bdb58d503f280">00030</a>         <span class="keywordtype">int</span> <a class="code" href="classPF.html#2c2f44ed7a4eaa42e07bdb58d503f280" title="number of particles;">n</a>; 
     29<a name="l00032"></a><a class="code" href="classPF.html#a2ac56d1e3ffbb4ff0b3f02e6399deb0">00032</a>         <a class="code" href="classeEmp.html" title="Weighted empirical density.">eEmp</a> <a class="code" href="classPF.html#a2ac56d1e3ffbb4ff0b3f02e6399deb0" title="posterior density">ePdf</a>; 
     30<a name="l00034"></a><a class="code" href="classPF.html#a97d12da4d1832c0b0c6ec5877f921f0">00034</a>         vec &amp;<a class="code" href="classPF.html#a97d12da4d1832c0b0c6ec5877f921f0" title="pointer into eEmp ">w</a>; 
     31<a name="l00036"></a><a class="code" href="classPF.html#361743a0b5b89de1a29e91d1343b2565">00036</a>         Array&lt;vec&gt; &amp;<a class="code" href="classPF.html#361743a0b5b89de1a29e91d1343b2565" title="pointer into eEmp ">samples</a>; 
     32<a name="l00038"></a><a class="code" href="classPF.html#d92ac103f88f8c21e197e90af5695a09">00038</a>         <a class="code" href="classmpdf.html" title="Conditional probability density, e.g. modeling some dependencies.">mpdf</a> &amp;<a class="code" href="classPF.html#d92ac103f88f8c21e197e90af5695a09" title="Parameter evolution model.">par</a>; 
     33<a name="l00040"></a><a class="code" href="classPF.html#dd0a687a4515333d6809147335854e77">00040</a>         <a class="code" href="classmpdf.html" title="Conditional probability density, e.g. modeling some dependencies.">mpdf</a> &amp;<a class="code" href="classPF.html#dd0a687a4515333d6809147335854e77" title="Observation model.">obs</a>; 
     34<a name="l00041"></a>00041 <span class="keyword">public</span>: 
     35<a name="l00042"></a>00042         <a class="code" href="classPF.html" title="Trivial particle filter with proposal density equal to parameter evolution model...">PF</a> ( <span class="keyword">const</span> <a class="code" href="classRV.html" title="Class representing variables, most often random variables.">RV</a> &amp;rv0, <a class="code" href="classmpdf.html" title="Conditional probability density, e.g. modeling some dependencies.">mpdf</a> &amp;par0,  <a class="code" href="classmpdf.html" title="Conditional probability density, e.g. modeling some dependencies.">mpdf</a> &amp;obs0, <span class="keywordtype">int</span> <a class="code" href="classPF.html#2c2f44ed7a4eaa42e07bdb58d503f280" title="number of particles;">n</a> ) :<a class="code" href="classBM.html" title="Bayesian Model of the world, i.e. all uncertainty is modeled by probabilities.">BM</a>(rv0), 
     36<a name="l00043"></a>00043                         n ( 1 ),<a class="code" href="classPF.html#a2ac56d1e3ffbb4ff0b3f02e6399deb0" title="posterior density">ePdf</a> ( rv0 ),<a class="code" href="classPF.html#a97d12da4d1832c0b0c6ec5877f921f0" title="pointer into eEmp ">w</a> ( <a class="code" href="classPF.html#a2ac56d1e3ffbb4ff0b3f02e6399deb0" title="posterior density">ePdf</a>._w() ),<a class="code" href="classPF.html#361743a0b5b89de1a29e91d1343b2565" title="pointer into eEmp ">samples</a> ( <a class="code" href="classPF.html#a2ac56d1e3ffbb4ff0b3f02e6399deb0" title="posterior density">ePdf</a>._samples() ), 
     37<a name="l00044"></a>00044                         <a class="code" href="classPF.html#d92ac103f88f8c21e197e90af5695a09" title="Parameter evolution model.">par</a> ( par0 ), <a class="code" href="classPF.html#dd0a687a4515333d6809147335854e77" title="Observation model.">obs</a> ( obs0 ) {}; 
     38<a name="l00045"></a>00045 <span class="comment">//      void set_parametres(mpdf &amp;par0, mpdf &amp;obs0) {par=&amp;par0;obs=&amp;obs0;};</span> 
     39<a name="l00046"></a>00046         <span class="keywordtype">void</span> set_est ( <span class="keyword">const</span> <a class="code" href="classepdf.html" title="Probability density function with numerical statistics, e.g. posterior density.">epdf</a>* &amp;epdf0 ); 
     40<a name="l00047"></a>00047         <span class="comment">//TODO remove or implement bayes()!</span> 
     41<a name="l00048"></a>00048         <span class="keywordtype">void</span> <a class="code" href="classPF.html#64f636bbd63bea9efd778214e6b631d3" title="Incremental Bayes rule.">bayes</a> ( <span class="keyword">const</span> vec &amp;dt ); 
     42<a name="l00049"></a>00049 }; 
     43<a name="l00050"></a>00050  
     44<a name="l00057"></a>00057 <span class="keyword">template</span>&lt;<span class="keyword">class</span> BM_T&gt; 
     45<a name="l00058"></a><a class="code" href="classMPF.html">00058</a> <span class="keyword">class </span><a class="code" href="classMPF.html" title="Marginalized Particle filter.">MPF</a> : <span class="keyword">public</span> <a class="code" href="classPF.html" title="Trivial particle filter with proposal density equal to parameter evolution model...">PF</a> { 
     46<a name="l00059"></a>00059         BM_T* Bms[1000]; 
     47<a name="l00060"></a>00060         <a class="code" href="classRV.html" title="Class representing variables, most often random variables.">RV</a> rvc; 
     48<a name="l00061"></a>00061         <a class="code" href="classemix.html" title="Weighted mixture of epdfs with external owned components.">emix</a> est; 
     49<a name="l00062"></a>00062 <span class="keyword">public</span>: 
     50<a name="l00064"></a><a class="code" href="classMPF.html#827a66609cf69a832535d52233f76fa0">00064</a>         <a class="code" href="classMPF.html#827a66609cf69a832535d52233f76fa0" title="Default constructor.">MPF</a> ( <span class="keyword">const</span> <a class="code" href="classRV.html" title="Class representing variables, most often random variables.">RV</a> &amp;rv0, <a class="code" href="classmpdf.html" title="Conditional probability density, e.g. modeling some dependencies.">mpdf</a> &amp;par0, <a class="code" href="classmpdf.html" title="Conditional probability density, e.g. modeling some dependencies.">mpdf</a> &amp;obs0, <span class="keywordtype">int</span> <a class="code" href="classPF.html#2c2f44ed7a4eaa42e07bdb58d503f280" title="number of particles;">n</a>, <span class="keyword">const</span> BM_T &amp;BMcond0 ) : 
     51<a name="l00065"></a>00065                         <a class="code" href="classPF.html" title="Trivial particle filter with proposal density equal to parameter evolution model...">PF</a> ( rv0,par0,obs0,n ), est(rv0,<a class="code" href="classPF.html#a97d12da4d1832c0b0c6ec5877f921f0" title="pointer into eEmp ">w</a>) { 
     52<a name="l00066"></a>00066                 <span class="comment">//</span> 
     53<a name="l00067"></a>00067                 <span class="keywordflow">if</span> (n&gt;1000) it_error(<span class="stringliteral">"increase 1000 here!"</span>); 
     54<a name="l00068"></a>00068                 <span class="keywordflow">for</span> ( <span class="keywordtype">int</span> i=0;i&lt;n;i++ ) { 
     55<a name="l00069"></a>00069                         Bms[i] = <span class="keyword">new</span> BM_T(BMcond0); <span class="comment">//copy constructor</span> 
     56<a name="l00070"></a>00070                         <a class="code" href="classepdf.html" title="Probability density function with numerical statistics, e.g. posterior density.">epdf</a>&amp; pom=Bms[i]-&gt;_epdf(); 
     57<a name="l00071"></a>00071                         est.<a class="code" href="classemix.html#599366d678152cfb0703da697b6b1de1">set_parameters</a>(i,&amp;pom); 
     58<a name="l00072"></a>00072                 } 
     59<a name="l00073"></a>00073         }; 
     60<a name="l00074"></a>00074         ~<a class="code" href="classMPF.html" title="Marginalized Particle filter.">MPF</a>() { 
     61<a name="l00075"></a>00075         } 
     62<a name="l00076"></a>00076         <span class="keywordtype">void</span> <a class="code" href="classMPF.html#55daf8e4b6553dd9f47c692de7931623" title="Incremental Bayes rule.">bayes</a> (<span class="keyword">const</span> vec &amp;dt ); 
     63<a name="l00077"></a><a class="code" href="classMPF.html#549e08268a46a250f21a33d06f19276a">00077</a>         <a class="code" href="classepdf.html" title="Probability density function with numerical statistics, e.g. posterior density.">epdf</a>&amp; <a class="code" href="classMPF.html#549e08268a46a250f21a33d06f19276a" title="Returns a pointer to the epdf representing posterior density on parameters. Use with...">_epdf</a>(){<span class="keywordflow">return</span> est;} 
     64<a name="l00078"></a>00078 }; 
     65<a name="l00079"></a>00079  
     66<a name="l00080"></a>00080 <span class="keyword">template</span>&lt;<span class="keyword">class</span> BM_T&gt; 
     67<a name="l00081"></a><a class="code" href="classMPF.html#55daf8e4b6553dd9f47c692de7931623">00081</a> <span class="keywordtype">void</span> <a class="code" href="classMPF.html#55daf8e4b6553dd9f47c692de7931623" title="Incremental Bayes rule.">MPF&lt;BM_T&gt;::bayes</a>( <span class="keyword">const</span> vec &amp;dt) { 
     68<a name="l00082"></a>00082         <span class="keywordtype">int</span> i; 
     69<a name="l00083"></a>00083         vec lls(<a class="code" href="classPF.html#2c2f44ed7a4eaa42e07bdb58d503f280" title="number of particles;">n</a>); 
     70<a name="l00084"></a>00084         ivec ind; 
     71<a name="l00085"></a>00085         <span class="keywordtype">double</span> mlls=-std::numeric_limits&lt;double&gt;::infinity(), sum=0.0; 
     72<a name="l00086"></a>00086          
     73<a name="l00087"></a>00087         <span class="keywordflow">for</span> ( i=0;i&lt;<a class="code" href="classPF.html#2c2f44ed7a4eaa42e07bdb58d503f280" title="number of particles;">n</a>;i++ ) { 
     74<a name="l00088"></a>00088         <span class="comment">//generate new samples from paramater evolution model;</span> 
     75<a name="l00089"></a>00089                 <a class="code" href="classPF.html#361743a0b5b89de1a29e91d1343b2565" title="pointer into eEmp ">samples</a>(i) = <a class="code" href="classPF.html#d92ac103f88f8c21e197e90af5695a09" title="Parameter evolution model.">par</a>.<a class="code" href="classmpdf.html#c20c796f8d0a201f0897299150e45a41" title="Returns the required moment of the epdf.">samplecond</a>( <a class="code" href="classPF.html#361743a0b5b89de1a29e91d1343b2565" title="pointer into eEmp ">samples</a>( i ), lls(i) ); 
     76<a name="l00090"></a>00090                 <span class="keywordflow">if</span> (lls(i)&gt;mlls) mlls=lls(i); <span class="comment">//find maximum</span> 
     77<a name="l00091"></a>00091         } 
     78<a name="l00092"></a>00092         <span class="comment">// compute weights </span> 
     79<a name="l00093"></a>00093         <span class="keywordflow">for</span> ( i=0;i&lt;n;i++ ) { 
     80<a name="l00094"></a>00094                 <a class="code" href="classPF.html#a97d12da4d1832c0b0c6ec5877f921f0" title="pointer into eEmp ">w</a>(i) *= exp( lls(i) - mlls); <span class="comment">// multiply w by likelihood</span> 
     81<a name="l00095"></a>00095         } 
     82<a name="l00096"></a>00096         <span class="comment">//renormalize</span> 
     83<a name="l00097"></a>00097         <span class="keywordflow">for</span> ( i=0;i&lt;n;i++ ){sum+=<a class="code" href="classPF.html#a97d12da4d1832c0b0c6ec5877f921f0" title="pointer into eEmp ">w</a>( i );}; 
     84<a name="l00098"></a>00098                 <a class="code" href="classPF.html#a97d12da4d1832c0b0c6ec5877f921f0" title="pointer into eEmp ">w</a>( i ) /=sum; <span class="comment">//?</span> 
     85<a name="l00099"></a>00099          
     86<a name="l00100"></a>00100         ind = <a class="code" href="classPF.html#a2ac56d1e3ffbb4ff0b3f02e6399deb0" title="posterior density">ePdf</a>.<a class="code" href="classeEmp.html#77268292fc4465cb73ddbfb1f2932a59" title="Function performs resampling, i.e. removal of low-weight samples and duplication...">resample</a>(); 
     87<a name="l00101"></a>00101          
     88<a name="l00102"></a>00102 } 
     89<a name="l00103"></a>00103  
     90<a name="l00104"></a>00104 <span class="preprocessor">#endif // KF_H</span> 
     91<a name="l00105"></a>00105 <span class="preprocessor"></span> 
     92<a name="l00106"></a>00106  
     93</pre></div><hr size="1"><address style="text-align: right;"><small>Generated on Thu Feb 28 16:54:40 2008 for mixpp by&nbsp; 
    6694<a href="http://www.doxygen.org/index.html"> 
    6795<img src="doxygen.png" alt="doxygen" align="middle" border="0"></a> 1.5.3 </small></address>