| 33 | | protected: |
| 34 | | //! weights of the components |
| 35 | | vec w; |
| 36 | | //! Component (epdfs) |
| 37 | | Array<epdf*> Coms; |
| 38 | | public: |
| 39 | | //!Default constructor |
| 40 | | emix(RV &rv) : epdf(rv) {}; |
| 41 | | //! Set weights \c w and components \c R |
| 42 | | void set_parameters(const vec &w, const Array<epdf*> &Coms); |
| | 33 | protected: |
| | 34 | //! weights of the components |
| | 35 | vec w; |
| | 36 | //! Component (epdfs) |
| | 37 | Array<epdf*> Coms; |
| | 38 | public: |
| | 39 | //!Default constructor |
| | 40 | emix ( RV &rv ) : epdf ( rv ) {}; |
| | 41 | //! Set weights \c w and components \c R |
| | 42 | void set_parameters ( const vec &w, const Array<epdf*> &Coms ); |
| 44 | | vec sample() const; |
| 45 | | vec mean() const { |
| 46 | | int i; vec mu = zeros(rv.count()); |
| 47 | | for (i = 0;i < w.length();i++) {mu += w(i) * Coms(i)->mean(); } |
| 48 | | return mu; |
| 49 | | } |
| 50 | | double evalpdflog(const vec &val) const { |
| 51 | | int i; |
| 52 | | double sum = 0.0; |
| 53 | | for (i = 0;i < w.length();i++) {sum += w(i) * Coms(i)->evalpdflog(val);} |
| 54 | | return log(sum); |
| 55 | | }; |
| | 44 | vec sample() const; |
| | 45 | vec mean() const { |
| | 46 | int i; vec mu = zeros ( rv.count() ); |
| | 47 | for ( i = 0;i < w.length();i++ ) {mu += w ( i ) * Coms ( i )->mean(); } |
| | 48 | return mu; |
| | 49 | } |
| | 50 | double evalpdflog ( const vec &val ) const { |
| | 51 | int i; |
| | 52 | double sum = 0.0; |
| | 53 | for ( i = 0;i < w.length();i++ ) {sum += w ( i ) * Coms ( i )->evalpdflog ( val );} |
| | 54 | return log ( sum ); |
| | 55 | }; |
| 70 | | class eprod: public epdf { |
| 71 | | protected: |
| 72 | | // |
| 73 | | int n; |
| 74 | | // pointers to epdfs |
| 75 | | Array<epdf*> epdfs; |
| 76 | | Array<mpdf*> mpdfs; |
| 77 | | // |
| 78 | | Array<ivec> rvinds; |
| 79 | | Array<ivec> rvcinds; |
| 80 | | //! Indicate independence of its factors |
| 81 | | bool independent; |
| 82 | | //! Indicate internal creation of mpdfs which must be destroyed |
| 83 | | bool intermpdfs; |
| 84 | | public: |
| 85 | | //!Constructor from list of eFacs or list of mFacs |
| 86 | | eprod(Array<mpdf*> mFacs): epdf(RV()), n(mFacs.length()), epdfs(n), mpdfs(mFacs), rvinds(n), rvcinds(n) { |
| 87 | | int i; |
| 88 | | intermpdfs = false; |
| 89 | | for (i = 0;i < n;i++) { |
| 90 | | rv.add(mpdfs(i)->_rv()); //add rv to common rvs. |
| 91 | | epdfs(i) = &(mpdfs(i)->_epdf()); // add pointer to epdf |
| 92 | | }; |
| 93 | | independent = true; |
| 94 | | //test rvc of mpdfs and fill rvinds |
| 95 | | for (i = 0;i < n;i++) { |
| 96 | | // find ith rv in common rv |
| 97 | | rvinds(i) = mpdfs(i)->_rv().dataind(rv); |
| 98 | | // find ith rvc in common rv |
| 99 | | rvcinds(i) = mpdfs(i)->_rvc().dataind(rv); |
| 100 | | if (rvcinds(i).length()>0) {independent = false;} |
| 101 | | } |
| 102 | | |
| | 70 | class mprod: public mpdf { |
| | 71 | protected: |
| | 72 | // |
| | 73 | int n; |
| | 74 | // pointers to epdfs |
| | 75 | Array<epdf*> epdfs; |
| | 76 | Array<mpdf*> mpdfs; |
| | 77 | // |
| | 78 | //! Indeces of rvs in common rv |
| | 79 | Array<ivec> rvinds; |
| | 80 | //! Indeces of rvc in common rv |
| | 81 | Array<ivec> rvcinrv; |
| | 82 | //! Indeces of rvc in common rvc |
| | 83 | Array<ivec> rvcinds; |
| | 84 | //! Indicate independence of its factors |
| | 85 | bool independent; |
| | 86 | // //! Indicate internal creation of mpdfs which must be destroyed |
| | 87 | // bool intermpdfs; |
| | 88 | public: |
| | 89 | //!Constructor from list of eFacs or list of mFacs |
| | 90 | mprod ( Array<mpdf*> mFacs ) : mpdf ( RV(), RV() ), n ( mFacs.length() ), epdfs ( n ), mpdfs ( mFacs ), rvinds ( n ), rvcinrv ( n ), rvcinds ( n ) { |
| | 91 | int i; |
| | 92 | // Create rv |
| | 93 | for ( i = 0;i < n;i++ ) { |
| | 94 | rv.add ( mpdfs ( i )->_rv() ); //add rv to common rvs. |
| | 95 | epdfs ( i ) = & ( mpdfs ( i )->_epdf() ); // add pointer to epdf |
| 104 | | eprod(Array<epdf*> eFacs): epdf(RV()), n(eFacs.length()), epdfs(eFacs), mpdfs(n), rvinds(n), rvcinds(n) { |
| 105 | | int i; |
| 106 | | intermpdfs = true; |
| 107 | | for (i = 0;i < n;i++) { |
| 108 | | if (rv.add(eFacs(i)->_rv())) {it_error("Incompatible eFacs.rv!");} //add rv to common rvs. |
| 109 | | mpdfs(i) = new mepdf(*(epdfs(i))); |
| 110 | | epdfs(i) = &(mpdfs(i)->_epdf()); // add pointer to epdf |
| 111 | | }; |
| 112 | | independent = true; |
| 113 | | //test rvc of mpdfs and fill rvinds |
| 114 | | for (i = 0;i < n;i++) { |
| 115 | | // find ith rv in common rv |
| 116 | | rvinds(i) = mpdfs(i)->_rv().dataind(rv); |
| 117 | | } |
| | 97 | // Create rvc |
| | 98 | for ( i = 0;i < n;i++ ) { |
| | 99 | rvc.add ( mpdfs ( i )->_rv().subt ( rv ) ); //add rv to common rvs. |
| 120 | | double evalpdflog(const vec &val) const { |
| 121 | | int i; |
| 122 | | double res = 0.0; |
| 123 | | for (i = n - 1;i > 0;i++) { |
| 124 | | if (rvcinds(i).length() > 0) |
| 125 | | {mpdfs(i)->condition(val(rvcinds(i)));} |
| 126 | | // add logarithms |
| 127 | | res += epdfs(i)->evalpdflog(val(rvinds(i))); |
| | 102 | independent = true; |
| | 103 | //test rvc of mpdfs and fill rvinds |
| | 104 | for ( i = 0;i < n;i++ ) { |
| | 105 | // find ith rv in common rv |
| | 106 | rvinds ( i ) = mpdfs ( i )->_rv().dataind ( rv ); |
| | 107 | // find ith rvc in common rv |
| | 108 | rvcinrv ( i ) = mpdfs ( i )->_rvc().dataind ( rv ); |
| | 109 | // find ith rvc in common rv |
| | 110 | rvcinds ( i ) = mpdfs ( i )->_rvc().dataind ( rvc ); |
| | 111 | // |
| | 112 | if ( rvcinds ( i ).length() >0 ) {independent = false;} |
| | 113 | if ( rvcinds ( i ).length() >0 ) {independent = false;} |
| | 114 | } |
| | 115 | }; |
| | 116 | |
| | 117 | double evalpdflog ( const vec &val ) const { |
| | 118 | int i; |
| | 119 | double res = 0.0; |
| | 120 | for ( i = n - 1;i > 0;i++ ) { |
| | 121 | if ( rvcinds ( i ).length() > 0 ) |
| | 122 | {mpdfs ( i )->condition ( val ( rvcinds ( i ) ) );} |
| | 123 | // add logarithms |
| | 124 | res += epdfs ( i )->evalpdflog ( val ( rvinds ( i ) ) ); |
| | 125 | } |
| | 126 | return res; |
| | 127 | } |
| | 128 | vec samplecond ( const vec &cond, double &ll ) { |
| | 129 | vec smp=zeros ( rv.count() ); |
| | 130 | vec condi; |
| | 131 | for ( int i = ( n - 1 );i >= 0;i-- ) { |
| | 132 | if ( rvcinds ( i ).length() > 0 ) { |
| | 133 | condi = zeros ( rvcinrv.length() + rvcinds.length() ); |
| | 134 | // copy data in condition |
| | 135 | set_subvector ( condi,rvcinds ( i ), cond ); |
| | 136 | // copy data in already generated sample |
| | 137 | set_subvector ( condi,rvcinrv ( i ), smp ); |
| | 138 | |
| | 139 | mpdfs ( i )->condition ( condi ); |
| 131 | | vec sample () const{ |
| 132 | | vec smp=zeros(rv.count()); |
| 133 | | for (int i = (n - 1);i >= 0;i--) { |
| 134 | | if (rvcinds(i).length() > 0) |
| 135 | | {mpdfs(i)->condition(smp(rvcinds(i)));} |
| 136 | | set_subvector(smp,rvinds(i), epdfs(i)->sample()); |
| 137 | | } |
| 138 | | return smp; |
| 139 | | } |
| 140 | | vec mean() const{ |
| 141 | | vec tmp(rv.count()); |
| 142 | | if (independent) { |
| 143 | | for (int i=0;i<n;i++) { |
| 144 | | vec pom = epdfs(i)->mean(); |
| 145 | | set_subvector(tmp,rvinds(i), pom); |
| 146 | | } |
| 147 | | } |
| 148 | | else { |
| 149 | | int N=50*rv.count(); |
| 150 | | it_warning("eprod.mean() computed by sampling"); |
| 151 | | tmp = zeros(rv.count()); |
| 152 | | for (int i=0;i<N;i++){ tmp += sample();} |
| 153 | | tmp /=N; |
| 154 | | } |
| 155 | | return tmp; |
| 156 | | } |
| 157 | | ~eprod(){if (intermpdfs) {for (int i=0;i<n;i++){delete mpdfs(i);}}}; |
| | 144 | return smp; |
| | 145 | } |
| | 146 | mat samplecond ( const vec &cond, vec &ll, int N ) { |
| | 147 | mat Smp(rv.count(),N); |
| | 148 | for(int i=0;i<N;i++){Smp.set_col(i,samplecond(cond,ll(i)));} |
| | 149 | return Smp; |
| | 150 | } |
| | 151 | // vec mean() const { |
| | 152 | // vec tmp ( rv.count() ); |
| | 153 | // if ( independent ) { |
| | 154 | // for ( int i=0;i<n;i++ ) { |
| | 155 | // vec pom = epdfs ( i )->mean(); |
| | 156 | // set_subvector ( tmp,rvinds ( i ), pom ); |
| | 157 | // } |
| | 158 | // } |
| | 159 | // else { |
| | 160 | // int N=50*rv.count(); |
| | 161 | // it_warning ( "eprod.mean() computed by sampling" ); |
| | 162 | // tmp = zeros ( rv.count() ); |
| | 163 | // for ( int i=0;i<N;i++ ) { tmp += sample();} |
| | 164 | // tmp /=N; |
| | 165 | // } |
| | 166 | // return tmp; |
| | 167 | // } |
| | 168 | |
| | 169 | ~mprod() {}; |
| 164 | | protected: |
| 165 | | //! Component (epdfs) |
| 166 | | Array<mpdf*> Coms; |
| 167 | | //!Internal epdf |
| 168 | | emix Epdf; |
| 169 | | public: |
| 170 | | //!Default constructor |
| 171 | | mmix(RV &rv, RV &rvc) : mpdf(rv, rvc), Epdf(rv) {ep = &Epdf;}; |
| 172 | | //! Set weights \c w and components \c R |
| 173 | | void set_parameters(const vec &w, const Array<mpdf*> &Coms) { |
| 174 | | Array<epdf*> Eps(Coms.length()); |
| | 176 | protected: |
| | 177 | //! Component (epdfs) |
| | 178 | Array<mpdf*> Coms; |
| | 179 | //!Internal epdf |
| | 180 | emix Epdf; |
| | 181 | public: |
| | 182 | //!Default constructor |
| | 183 | mmix ( RV &rv, RV &rvc ) : mpdf ( rv, rvc ), Epdf ( rv ) {ep = &Epdf;}; |
| | 184 | //! Set weights \c w and components \c R |
| | 185 | void set_parameters ( const vec &w, const Array<mpdf*> &Coms ) { |
| | 186 | Array<epdf*> Eps ( Coms.length() ); |
