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