| | 8 | V.opupdt ( dt,frg ); |
| | 9 | nu+=frg; |
| | 10 | |
| | 11 | if ( evalll ) { |
| | 12 | lnc = est.lognc(); |
| | 13 | ll = lnc - last_lognc; |
| | 14 | last_lognc = lnc; |
| | 15 | } |
| | 16 | } |
| | 17 | |
| | 18 | /*! \brief Return the best structure |
| | 19 | @param Eg a copy of GiW density that is being examined |
| | 20 | @param Eg0 a copy of prior GiW density before estimation |
| | 21 | @param Egll likelihood of the current Eg |
| | 22 | @param indeces current indeces |
| | 23 | \return best likelihood in the structure below the given one |
| | 24 | */ |
| | 25 | double egiw_bestbelow ( egiw Eg, egiw Eg0, double Egll, ivec &indeces ) { //parameter Eg is a copy! |
| | 26 | ldmat Vo = Eg._V(); //copy |
| | 27 | ldmat Vo0 = Eg._V(); //copy |
| | 28 | ldmat& Vp = Eg._V(); // pointer into Eg |
| | 29 | ldmat& Vp0 = Eg._V(); // pointer into Eg |
| | 30 | int end = Vp.rows()-1; |
| | 31 | int i; |
| | 32 | mat Li; |
| | 33 | mat Li0; |
| | 34 | double maxll=Egll; |
| | 35 | double tmpll=Egll; |
| | 36 | double belll=Egll; |
| | 37 | |
| | 38 | ivec tmpindeces; |
| | 39 | ivec maxindeces=indeces; |
| | 40 | |
| | 41 | //try to remove only one rv |
| | 42 | for ( i=0;i<end;i++ ) { |
| | 43 | //copy original |
| | 44 | Li = Vo._L(); |
| | 45 | Li0 = Vo0._L(); |
| | 46 | //remove stuff |
| | 47 | Li.del_col ( i+1 ); |
| | 48 | Li0.del_col ( i+1 ); |
| | 49 | Vp.ldform ( Li,Vo._D() ); |
| | 50 | Vp0.ldform ( Li0,Vo0._D() ); |
| | 51 | tmpll = Eg.lognc()-Eg0.lognc(); // likelihood is difference of norm. coefs. |
| | 52 | // |
| | 53 | if ( tmpll > Egll ) { //increase of the likelihood |
| | 54 | tmpindeces = indeces; |
| | 55 | tmpindeces.del ( i ); |
| | 56 | //search for a better match in this substructure |
| | 57 | belll=egiw_bestbelow ( Eg, Eg0, tmpll, tmpindeces ); |
| | 58 | if ( belll>maxll ) { //better match found |
| | 59 | maxll = belll; |
| | 60 | maxindeces = tmpindeces; |
| | 61 | } |
| | 62 | } |
| | 63 | } |
| | 64 | indeces = maxindeces; |
| | 65 | return maxll; |
| | 66 | } |
| | 67 | |
| | 68 | ivec ARX::structure_est ( egiw est0 ) { |
| | 69 | ivec ind=linspace ( 1,rv.count()-1 ); |
| | 70 | double tmp = egiw_bestbelow ( est, est0, est.lognc()- est0.lognc(), ind ); |
| | 71 | return ind; |
| | 72 | } |
