Changeset 1025 for library/bdm/estim/arx.cpp

Timestamp:

05/31/10 15:55:10 (14 years ago)

Author:

dedecius

Message:

Class implementing partial forgetting in ARX in a way similar to ARXfrg. The ARXpartforg::bayes() method (sh|c)ould be divided into smaller functions in the future(?)

Files:

• library/bdm/estim/arx.cpp (modified) (1 diff)

library/bdm/estim/arx.cpp

@@ -291 +291 @@
         }
 }
-}
-
+
+void ARXpartialforg::bayes ( const vec &val, const vec &cond ) {
+    #define LOG2  0.69314718055995
+    vec frg = cond.right(cond.length() - rgrlen);
+    vec cond_rgr = cond.left(rgrlen);       // regression vector
+
+    int dimV = est._V().cols();
+    int nparams = dimV - 1;                 // number of parameters
+    int nalternatives = pow(2, nparams);    // number of alternatives
+    
+    // Permutation matrix
+    mat perm_matrix = ones(nalternatives, nparams);
+    int i, j, period, idx_from, idx_to, start, end;
+    for(i = 0; i < nparams; i++) {
+        idx_from = pow(2, i);
+        idx_to = 2 * pow(2, i) - 1;
+        period = pow(2, i+1);
+        j = 0;
+        start = idx_from;
+        end = idx_to;
+        while(start < pow(2, nparams)) {
+            perm_matrix.set_submatrix(start, end, i, i, 0);
+            j++;
+            start = idx_from + period * j;
+            end = idx_to + period * j;
+        }
+    }
+
+    // Array of egiws for approximation
+    Array<egiw*> egiw_array(nalternatives + 1);
+    // No. of conditioning rows in LD
+    int nalternatives_cond, position;
+
+    for(int i = 0; i < nalternatives; i++) {
+        // vector defining alternatives
+        vec vec_alt = perm_matrix.get_row(i);
+
+        // full alternative or filtered
+        if( sum(vec_alt) == vec_alt.length() )  {
+            egiw_array(i) = &alter_est;
+            continue;
+        } else if( sum(vec_alt) == 0 ) {
+            egiw_array(i) = &est;
+            continue;
+        }
+
+        nalternatives_cond = sum(vec_alt) + 1;
+        ivec vec_perm(0);                   // permutation vector
+
+        for(int j = 0; j < nparams; j++) {
+            position = dimV - j - 2;
+            if ( vec_alt(position) == 0 ) {
+                vec_perm.ins(j, position + 1);
+            }
+            else {
+                vec_perm.ins(0, position + 1);
+            }
+        }
+        vec_perm.ins(0, 0);
+
+        ldmat filt (est._V(), vec_perm);
+        ldmat alt (alter_est._V(), vec_perm);
+
+        mat tmpL(dimV, dimV);
+        tmpL.set_rows( 0, alt._L().get_rows(0, nalternatives_cond - 1) );
+        tmpL.set_rows( nalternatives_cond, filt._L().get_rows(nalternatives_cond, nparams) );
+
+        vec tmpD(dimV);
+        tmpD.set_subvector( 0, alt._D()(0, nalternatives_cond - 1) );
+        tmpD.set_subvector( nalternatives_cond, filt._D()(nalternatives_cond, nparams) );
+
+        ldmat tmpLD (tmpL, tmpD);
+
+        vec_perm = sort_index(vec_perm);
+        ldmat newLD (tmpLD, vec_perm);
+
+        egiw_array(i) = new egiw(1, newLD, alter_est._nu());
+    }
+
+    // Approximation
+    double sumVecCommon;                // frequently used term
+    vec vecNu(nalternatives);           // vector of nus of components
+    vec vecD(nalternatives);            // vector of LS reminders 
+    vec vecCommon(nalternatives);       // vector of common parts
+    mat matVecsTheta;                           // matrix whose rows are theta vects.
+
+    for (i = 0; i < nalternatives; i++) {
+        vecNu.shift_left( egiw_array(i)->_nu() );
+        vecD.shift_left( egiw_array(i)->_V()._D()(0) );
+        matVecsTheta.append_row( egiw_array(i)->est_theta()  );
+    }
+
+    vecCommon = elem_mult ( frg, elem_div(vecNu, vecD) );
+    sumVecCommon = sum(vecCommon);
+
+    // approximation of est. regr. coefficients
+    vec aprEstTheta(nparams);  aprEstTheta.zeros();
+    for (i = 0; i < nalternatives; i++) {
+        aprEstTheta +=  matVecsTheta.get_row(i) * vecCommon(i);
+    }
+    aprEstTheta /= sumVecCommon;
+    
+    // approximation of degr. of freedom
+    double aprNu;
+    double A = log( sumVecCommon );             // Term 'A' in equation
+    for ( int i = 0; i < nalternatives; i++ ) {
+        A += frg(i) * ( log( vecD(i) ) - psi( 0.5 * vecNu(i) ) );
+    }
+    aprNu = ( 1 + sqrt(1 + 4 * (A - LOG2)/3 ) ) / ( 2 * (A - LOG2) );
+
+    // approximation of LS reminder D(0,0)
+    double aprD = aprNu / sumVecCommon;
+
+    // Aproximation of covariance of LS est.
+    mat aprC = zeros(nparams, nparams);
+    for ( int i = 0; i < nalternatives; i++ ) {
+        aprC += egiw_array(i)->est_theta_cov().to_mat() * frg(i); 
+        vec tmp = ( matVecsTheta.get_row(i) - aprEstTheta );
+        aprC += vecCommon(i) * outer_product( tmp, tmp);
+    }
+
+    // Construct GiW pdf
+    ldmat aprCinv ( inv(aprC) );
+    vec D = concat( aprD, aprCinv._D() );
+    mat L = eye(dimV);
+    L.set_submatrix(1, 0, aprCinv._L() * aprEstTheta);
+    L.set_submatrix(1, 1, aprCinv._L());
+    ldmat aprLD (L, D);
+    est = egiw(1, aprLD, aprNu);
+
+    // update
+    ARX::bayes ( val, cond_rgr );
+}
+
+}
+