| 22 | | <a name="l00017"></a>00017 <span class="preprocessor">#include "<a class="code" href="libBM_8h.html" title="Bayesian Models (bm) that use Bayes rule to learn from observations.">libBM.h</a>"</span> |
| 23 | | <a name="l00018"></a>00018 <span class="preprocessor">#include "<a class="code" href="libEF_8h.html" title="Probability distributions for Exponential Family models.">libEF.h</a>"</span> |
| 24 | | <a name="l00019"></a>00019 <span class="preprocessor">#include "<a class="code" href="emix_8h.html" title="Probability distributions for Mixtures of pdfs.">emix.h</a>"</span> |
| 25 | | <a name="l00020"></a>00020 |
| 26 | | <a name="l00021"></a>00021 <span class="keyword">using namespace </span>itpp; |
| 27 | | <a name="l00022"></a>00022 <span class="keyword">using</span> std::string; |
| 28 | | <a name="l00023"></a>00023 |
| 29 | | <a name="l00030"></a><a class="code" href="classmerger.html">00030</a> <span class="keyword">class </span><a class="code" href="classmerger.html" title="Function for general combination of pdfs.">merger</a> : <span class="keyword">public</span> <a class="code" href="classmprod.html" title="Chain rule decomposition of epdf.">mprod</a>{ |
| 30 | | <a name="l00031"></a>00031 <span class="keyword">protected</span>: |
| 31 | | <a name="l00033"></a><a class="code" href="classmerger.html#1712341746811f2844de5499817abb46">00033</a> <a class="code" href="classenorm.html" title="Gaussian density with positive definite (decomposed) covariance matrix.">enorm<fsqmat></a> <a class="code" href="classmerger.html#1712341746811f2844de5499817abb46" title="Additional pdf on the part in condition (if undefined);.">condpdf</a>; |
| 32 | | <a name="l00035"></a><a class="code" href="classmerger.html#11eb177ef6409926b2ba3b4630bceab1">00035</a> Array<ivec> <a class="code" href="classmerger.html#11eb177ef6409926b2ba3b4630bceab1" title="Find potential overlaps in rv.">overlaps</a>; |
| 33 | | <a name="l00036"></a>00036 <span class="keyword">public</span>: |
| 34 | | <a name="l00038"></a><a class="code" href="classmerger.html#1a6a721294321c2e412ca0084c636b69">00038</a> <a class="code" href="classmerger.html#1a6a721294321c2e412ca0084c636b69" title="Default constructor.">merger</a> (<span class="keyword">const</span> Array<mpdf*> &in_sources) : <a class="code" href="classmprod.html" title="Chain rule decomposition of epdf.">mprod</a>(in_sources,true), <a class="code" href="classmerger.html#1712341746811f2844de5499817abb46" title="Additional pdf on the part in condition (if undefined);.">condpdf</a>(<a class="code" href="classmpdf.html#acb7dda792b3cd5576f39fa3129abbab" title="random variable in condition">rvc</a>), <a class="code" href="classmerger.html#11eb177ef6409926b2ba3b4630bceab1" title="Find potential overlaps in rv.">overlaps</a>(<a class="code" href="classmpdf.html#f6687c07ff07d47812dd565368ca59eb" title="modeled random variable">rv</a>.count()) { |
| 35 | | <a name="l00039"></a>00039 <span class="keywordflow">if</span> (<a class="code" href="classmpdf.html#acb7dda792b3cd5576f39fa3129abbab" title="random variable in condition">rvc</a>.<a class="code" href="classRV.html#f5c7b8bd589eef09ccdf3329a0addea0" title="Return number of scalars in the RV.">count</a>()>0) { |
| 36 | | <a name="l00040"></a>00040 vec mu = zeros(<a class="code" href="classmpdf.html#acb7dda792b3cd5576f39fa3129abbab" title="random variable in condition">rvc</a>.<a class="code" href="classRV.html#f5c7b8bd589eef09ccdf3329a0addea0" title="Return number of scalars in the RV.">count</a>()); |
| 37 | | <a name="l00041"></a>00041 mat R = 100*eye(<a class="code" href="classmpdf.html#acb7dda792b3cd5576f39fa3129abbab" title="random variable in condition">rvc</a>.<a class="code" href="classRV.html#f5c7b8bd589eef09ccdf3329a0addea0" title="Return number of scalars in the RV.">count</a>()); |
| 38 | | <a name="l00042"></a>00042 <a class="code" href="classmerger.html#1712341746811f2844de5499817abb46" title="Additional pdf on the part in condition (if undefined);.">condpdf</a>.<a class="code" href="classenorm.html#1394a65caa6e00d42e00cc99b12227af" title="Set mean value mu and covariance R.">set_parameters</a>(mu,R); |
| 39 | | <a name="l00043"></a>00043 } |
| 40 | | <a name="l00044"></a>00044 <span class="keywordflow">for</span> (<span class="keywordtype">int</span> i=0;i<<a class="code" href="classmpdf.html#f6687c07ff07d47812dd565368ca59eb" title="modeled random variable">rv</a>.<a class="code" href="classRV.html#f5c7b8bd589eef09ccdf3329a0addea0" title="Return number of scalars in the RV.">count</a>();i++){ <span class="comment">// cycle over rv</span> |
| 41 | | <a name="l00045"></a>00045 ivec tmp(0); |
| 42 | | <a name="l00046"></a>00046 <span class="keywordflow">for</span> (<span class="keywordtype">int</span> j=0;j<n;j++){<span class="comment">//cycle over rvinds</span> |
| 43 | | <a name="l00047"></a>00047 ivec <span class="keywordtype">id</span>=itpp::find(<a class="code" href="classmprod.html#c8f78fbcdb4e38ecd3197d2958082497" title="Indeces of rvs in common rv.">rvinds</a>(j)==i); |
| 44 | | <a name="l00048"></a>00048 tmp=concat(tmp,<span class="keywordtype">id</span>); <span class="comment">// add i to tmp if found</span> |
| 45 | | <a name="l00049"></a>00049 } |
| 46 | | <a name="l00050"></a>00050 <a class="code" href="classmerger.html#11eb177ef6409926b2ba3b4630bceab1" title="Find potential overlaps in rv.">overlaps</a>(i)=tmp; |
| 47 | | <a name="l00051"></a>00051 } |
| 48 | | <a name="l00052"></a>00052 } |
| 49 | | <a name="l00055"></a>00055 <span class="comment">/* vec sample(double &w, const vec &smp0){</span> |
| 50 | | <a name="l00056"></a>00056 <span class="comment"> // result</span> |
| 51 | | <a name="l00057"></a>00057 <span class="comment"> vec smp=smp0;</span> |
| 52 | | <a name="l00058"></a>00058 <span class="comment"> vec cond=condpdf.sample(); // Just like in samplecond, here it is not given!</span> |
| 53 | | <a name="l00059"></a>00059 <span class="comment"> // temporary</span> |
| 54 | | <a name="l00060"></a>00060 <span class="comment"> mat smpi=zeros(rv.count()+rvc.count(), n);</span> |
| 55 | | <a name="l00061"></a>00061 <span class="comment"> vec condi;</span> |
| 56 | | <a name="l00062"></a>00062 <span class="comment"> vec ll(n);</span> |
| 57 | | <a name="l00063"></a>00063 <span class="comment"> </span> |
| 58 | | <a name="l00064"></a>00064 <span class="comment"> // Consider arithmetic mean as a proposal density</span> |
| 59 | | <a name="l00065"></a>00065 <span class="comment"> ll = 0;</span> |
| 60 | | <a name="l00066"></a>00066 <span class="comment"> for ( int i = ( n - 1 );i >= 0;i-- ) {</span> |
| 61 | | <a name="l00067"></a>00067 <span class="comment"> if ( rvcinds ( i ).length() > 0 ) {</span> |
| 62 | | <a name="l00068"></a>00068 <span class="comment"> condi = zeros ( rvcinrv.length() + rvcinds.length() );</span> |
| 63 | | <a name="l00069"></a>00069 <span class="comment"> // copy data in condition</span> |
| 64 | | <a name="l00070"></a>00070 <span class="comment"> set_subvector ( condi,rvcinds ( i ), cond );</span> |
| 65 | | <a name="l00071"></a>00071 <span class="comment"> // copy data in already generated sample</span> |
| 66 | | <a name="l00072"></a>00072 <span class="comment"> set_subvector ( condi,rvcinrv ( i ), smp );</span> |
| 67 | | <a name="l00073"></a>00073 <span class="comment"></span> |
| 68 | | <a name="l00074"></a>00074 <span class="comment"> mpdfs ( i )->condition ( condi );</span> |
| 69 | | <a name="l00075"></a>00075 <span class="comment"> }</span> |
| 70 | | <a name="l00076"></a>00076 <span class="comment"> smpi.set_col(i, epdfs ( i )->sample());</span> |
| 71 | | <a name="l00077"></a>00077 <span class="comment"> // add ith likelihood contribution</span> |
| 72 | | <a name="l00078"></a>00078 <span class="comment"> ll+=epdfs(i)->evalpdflog(smpi);</span> |
| 73 | | <a name="l00079"></a>00079 <span class="comment"> }</span> |
| 74 | | <a name="l00080"></a>00080 <span class="comment"> // Now lets analyze samples smpi</span> |
| 75 | | <a name="l00081"></a>00081 <span class="comment"> for (int i=0;i<rv.count();i++){</span> |
| 76 | | <a name="l00082"></a>00082 <span class="comment"> // number of components</span> |
| 77 | | <a name="l00083"></a>00083 <span class="comment"> int noc = overlaps(i).length();</span> |
| 78 | | <a name="l00084"></a>00084 <span class="comment"> if (noc<2) { // only one sample in this dimension, take it</span> |
| 79 | | <a name="l00085"></a>00085 <span class="comment"> smp(i) = smpi(i,overlaps(i)(0));</span> |
| 80 | | <a name="l00086"></a>00086 <span class="comment"> }</span> |
| 81 | | <a name="l00087"></a>00087 <span class="comment"> else {</span> |
| 82 | | <a name="l00088"></a>00088 <span class="comment"> //pick dimension</span> |
| 83 | | <a name="l00089"></a>00089 <span class="comment"> double unis = UniRNG.sample();</span> |
| 84 | | <a name="l00090"></a>00090 <span class="comment"> dim = 0;</span> |
| 85 | | <a name="l00091"></a>00091 <span class="comment"> while(unis<(dim+1)/noc) {dim++;} // </span> |
| 86 | | <a name="l00092"></a>00092 <span class="comment"> //in dim we now have randomly picked dimension</span> |
| 87 | | <a name="l00093"></a>00093 <span class="comment"> smp(i) = smpi(i,overlaps(i)(dim));</span> |
| 88 | | <a name="l00094"></a>00094 <span class="comment"> }</span> |
| 89 | | <a name="l00095"></a>00095 <span class="comment"> }</span> |
| 90 | | <a name="l00096"></a>00096 <span class="comment"> </span> |
| 91 | | <a name="l00097"></a>00097 <span class="comment"> // now we have sample how to evaluate weight?</span> |
| 92 | | <a name="l00098"></a>00098 <span class="comment"> return smp;</span> |
| 93 | | <a name="l00099"></a>00099 <span class="comment"> </span> |
| 94 | | <a name="l00100"></a>00100 <span class="comment"> // copied from mprod.sample</span> |
| 95 | | <a name="l00101"></a>00101 <span class="comment"> };*/</span> |
| 96 | | <a name="l00102"></a>00102 <span class="comment">// project </span> |
| 97 | | <a name="l00104"></a><a class="code" href="classmerger.html#51e425ee1dd46d73353526e6a1c24ca3">00104</a> <span class="comment"></span> <span class="keyword">virtual</span> <a class="code" href="classmerger.html#51e425ee1dd46d73353526e6a1c24ca3" title="for future use">~merger</a>() {}; |
| 98 | | <a name="l00105"></a>00105 }; |
| 99 | | <a name="l00106"></a>00106 |
| 100 | | <a name="l00107"></a>00107 |
| 101 | | <a name="l00108"></a>00108 <span class="preprocessor">#endif // MER_H</span> |
| | 22 | <a name="l00017"></a>00017 <span class="preprocessor">#include "<a class="code" href="mixef_8h.html" title="Bayesian Filtering for mixtures of exponential family (EF) members.">mixef.h</a>"</span> |
| | 23 | <a name="l00018"></a>00018 |
| | 24 | <a name="l00019"></a>00019 <span class="keyword">using namespace </span>itpp; |
| | 25 | <a name="l00020"></a>00020 <span class="keyword">using</span> std::string; |
| | 26 | <a name="l00021"></a>00021 |
| | 27 | <a name="l00028"></a><a class="code" href="classmerger.html">00028</a> <span class="keyword">class </span><a class="code" href="classmerger.html" title="Function for general combination of pdfs.">merger</a> : <span class="keyword">public</span> <a class="code" href="classcompositepdf.html" title="Abstract composition of pdfs, a base for specific classes.">compositepdf</a>, <span class="keyword">public</span> <a class="code" href="classepdf.html" title="Probability density function with numerical statistics, e.g. posterior density.">epdf</a> { |
| | 28 | <a name="l00029"></a>00029 <span class="keyword">protected</span>: |
| | 29 | <a name="l00031"></a><a class="code" href="classmerger.html#021d1522cea70db8545d708afd32ec2f">00031</a> <a class="code" href="classMixEF.html" title="Mixture of Exponential Family Densities.">MixEF</a> <a class="code" href="classmerger.html#021d1522cea70db8545d708afd32ec2f" title="Internal mixture of EF models.">Mix</a>; |
| | 30 | <a name="l00033"></a><a class="code" href="classmerger.html#ce84de1824875b4262eefb640e8c6d69">00033</a> <span class="keywordtype">int</span> <a class="code" href="classmerger.html#ce84de1824875b4262eefb640e8c6d69" title="Number of samples used in approximation.">Ns</a>; |
| | 31 | <a name="l00035"></a><a class="code" href="classmerger.html#6acf33e1c64aaaee20689725e8110f65">00035</a> <span class="keywordtype">int</span> <a class="code" href="classmerger.html#6acf33e1c64aaaee20689725e8110f65" title="Number of components in a mixture.">Nc</a>; |
| | 32 | <a name="l00037"></a><a class="code" href="classmerger.html#e68a72e2740007ae68545e378c7701ad">00037</a> <span class="keywordtype">double</span> <a class="code" href="classmerger.html#e68a72e2740007ae68545e378c7701ad" title="Prior on the log-normal merging model.">beta</a>; |
| | 33 | <a name="l00038"></a>00038 <span class="keyword">public</span>: |
| | 34 | <a name="l00040"></a><a class="code" href="classmerger.html#3f9a955329eab65026e8322e9027d469">00040</a> <a class="code" href="classmerger.html#3f9a955329eab65026e8322e9027d469" title="Default constructor.">merger</a> ( <span class="keyword">const</span> Array<mpdf*> &S ) : |
| | 35 | <a name="l00041"></a>00041 <a class="code" href="classcompositepdf.html" title="Abstract composition of pdfs, a base for specific classes.">compositepdf</a> ( S ), <a class="code" href="classepdf.html" title="Probability density function with numerical statistics, e.g. posterior density.">epdf</a> ( getrv ( false ) ), |
| | 36 | <a name="l00042"></a>00042 <a class="code" href="classmerger.html#021d1522cea70db8545d708afd32ec2f" title="Internal mixture of EF models.">Mix</a> ( Array<<a class="code" href="classBMEF.html" title="Estimator for Exponential family.">BMEF</a>*> ( 0 ),vec ( 0 ) ) |
| | 37 | <a name="l00043"></a>00043 { <a class="code" href="classmerger.html#e68a72e2740007ae68545e378c7701ad" title="Prior on the log-normal merging model.">beta</a>=2.0; <a class="code" href="classmerger.html#ce84de1824875b4262eefb640e8c6d69" title="Number of samples used in approximation.">Ns</a>=100; <a class="code" href="classmerger.html#6acf33e1c64aaaee20689725e8110f65" title="Number of components in a mixture.">Nc</a>=10;} |
| | 38 | <a name="l00045"></a><a class="code" href="classmerger.html#79420e251308df4502c834171d518ad9">00045</a> <span class="keywordtype">void</span> <a class="code" href="classmerger.html#79420e251308df4502c834171d518ad9" title="Set internal parameters used in approximation.">set_parameters</a> ( <span class="keywordtype">double</span> beta0, <span class="keywordtype">int</span> Ns0, <span class="keywordtype">int</span> Nc0 ) { <a class="code" href="classmerger.html#e68a72e2740007ae68545e378c7701ad" title="Prior on the log-normal merging model.">beta</a>=beta0;<a class="code" href="classmerger.html#ce84de1824875b4262eefb640e8c6d69" title="Number of samples used in approximation.">Ns</a>=Ns0;<a class="code" href="classmerger.html#6acf33e1c64aaaee20689725e8110f65" title="Number of components in a mixture.">Nc</a>=Nc0;} |
| | 39 | <a name="l00047"></a><a class="code" href="classmerger.html#02d2a445e44a87f37ee3a2b9b40cffd1">00047</a> <span class="keywordtype">void</span> <a class="code" href="classmerger.html#02d2a445e44a87f37ee3a2b9b40cffd1" title="Initialize the proposal density. This function must be called before merge()!">init</a>() { |
| | 40 | <a name="l00048"></a>00048 Array<vec> Smps ( <a class="code" href="classcompositepdf.html#c1646fe41ff5934d010e54ccaa46133f" title="Number of mpdfs in the composite.">n</a> ); |
| | 41 | <a name="l00049"></a>00049 <span class="comment">//Gibbs sampling</span> |
| | 42 | <a name="l00050"></a>00050 <span class="keywordflow">for</span> ( <span class="keywordtype">int</span> i=0;i<<a class="code" href="classcompositepdf.html#c1646fe41ff5934d010e54ccaa46133f" title="Number of mpdfs in the composite.">n</a>;i++ ) {Smps ( i ) =zeros ( 0 );} |
| | 43 | <a name="l00051"></a>00051 } |
| | 44 | <a name="l00053"></a>00053 <span class="keywordtype">void</span> <a class="code" href="classmerger.html#bc23accc77ec766bce7202934e2c180c" title="Create a mixture density, make sure to call init() before the first call.">merge</a> ( <span class="keyword">const</span> <a class="code" href="classepdf.html" title="Probability density function with numerical statistics, e.g. posterior density.">epdf</a>* g0 ); |
| | 45 | <a name="l00055"></a><a class="code" href="classmerger.html#bc23accc77ec766bce7202934e2c180c">00055</a> <span class="keywordtype">void</span> <a class="code" href="classmerger.html#bc23accc77ec766bce7202934e2c180c" title="Create a mixture density, make sure to call init() before the first call.">merge</a> () {<a class="code" href="classmerger.html#bc23accc77ec766bce7202934e2c180c" title="Create a mixture density, make sure to call init() before the first call.">merge</a> ( & ( <a class="code" href="classmerger.html#021d1522cea70db8545d708afd32ec2f" title="Internal mixture of EF models.">Mix</a>.<a class="code" href="classMixEF.html#efb3e20c2151d91c4fc080b7722a2069" title="Returns a pointer to the epdf representing posterior density on parameters. Use with...">_epdf</a>() ) );}; |
| | 46 | <a name="l00056"></a>00056 |
| | 47 | <a name="l00058"></a>00058 vec <a class="code" href="classmerger.html#217c65f335c603033fa09a2ecc64bc05" title="Merge log-likelihood values.">lognorm_merge</a> ( mat &lW ); |
| | 48 | <a name="l00061"></a><a class="code" href="classmerger.html#379198c3d2063bfa63f5d1245a2511ba">00061</a> vec <a class="code" href="classmerger.html#379198c3d2063bfa63f5d1245a2511ba">sample</a> ( )<span class="keyword">const </span>{ <span class="keywordflow">return</span> <a class="code" href="classmerger.html#021d1522cea70db8545d708afd32ec2f" title="Internal mixture of EF models.">Mix</a>.<a class="code" href="classMixEF.html#efb3e20c2151d91c4fc080b7722a2069" title="Returns a pointer to the epdf representing posterior density on parameters. Use with...">_epdf</a>().<a class="code" href="classepdf.html#8019654e494bf5e458f6fb947e11b262" title="Returns a sample, from density .">sample</a>();} |
| | 49 | <a name="l00062"></a><a class="code" href="classmerger.html#8c37688902b1a1e9fa32edc5709e5a00">00062</a> <span class="keywordtype">double</span> <a class="code" href="classmerger.html#8c37688902b1a1e9fa32edc5709e5a00" title="Compute log-probability of argument val.">evalpdflog</a> ( <span class="keyword">const</span> vec &dt )<span class="keyword"> const</span>{ <span class="keywordflow">return</span> <a class="code" href="classmerger.html#021d1522cea70db8545d708afd32ec2f" title="Internal mixture of EF models.">Mix</a>.<a class="code" href="classMixEF.html#efb3e20c2151d91c4fc080b7722a2069" title="Returns a pointer to the epdf representing posterior density on parameters. Use with...">_epdf</a>().<a class="code" href="classepdf.html#6aef3eca74899692503769c18add1a4c" title="Compute log-probability of argument val.">evalpdflog</a> ( dt );} |
| | 50 | <a name="l00063"></a><a class="code" href="classmerger.html#3a46999cd77047d5d4ac6c1cb735a83e">00063</a> vec <a class="code" href="classmerger.html#3a46999cd77047d5d4ac6c1cb735a83e" title="return expected value">mean</a>()<span class="keyword">const </span>{<span class="keywordflow">return</span> <a class="code" href="classmerger.html#021d1522cea70db8545d708afd32ec2f" title="Internal mixture of EF models.">Mix</a>.<a class="code" href="classMixEF.html#efb3e20c2151d91c4fc080b7722a2069" title="Returns a pointer to the epdf representing posterior density on parameters. Use with...">_epdf</a>().<a class="code" href="classepdf.html#bf0a070a8f3e67a93604ca724638b870" title="return expected value">mean</a>();} |
| | 51 | <a name="l00065"></a><a class="code" href="classmerger.html#51e425ee1dd46d73353526e6a1c24ca3">00065</a> <span class="keyword">virtual</span> <a class="code" href="classmerger.html#51e425ee1dd46d73353526e6a1c24ca3" title="for future use">~merger</a>() {}; |
| | 52 | <a name="l00066"></a>00066 }; |
| | 53 | <a name="l00067"></a>00067 |
| | 54 | <a name="l00068"></a>00068 |
| | 55 | <a name="l00069"></a>00069 |
| | 56 | <a name="l00070"></a>00070 <span class="preprocessor">#endif // MER_H</span> |
