Changeset 133

Timestamp:

07/07/08 19:45:36 (17 years ago)

Author:

smidl

Message:

linefit pro Martina

Location:

matlab/mex

Files:

• linefit.cpp (copied) (copied from matlab/mex/arx.cpp) (3 diffs)
• linefit2.cpp (added)

matlab/mex/linefit.cpp

@@ -1 @@
-#include <itpp/itcomm.h>
 #include <itpp/itmex.h>
 
@@ -9 +8 @@
         // Check the number of inputs and output arguments
         if ( n_output<1 ) mexErrMsgTxt ( "Wrong number of output variables!" );
-        if ( n_input!=1 ) mexErrMsgTxt ( "Wrong number of input variables!" );
+        if ( n_input!=3 ) mexErrMsgTxt ( "Wrong number of input variables!" );
 
         // Convert input variables to IT++ format
         mat Data = mxArray2mat ( input[0] );
+        if (Data.rows()!=2)
+                { if (Data.cols()==2) Data=Data.T(); else mexErrMsgTxt("Data are not 2D!");}
+        int ilow = mxArray2int ( input[1] );
+        int ihi = mxArray2int ( input[2] );
 
         // ------------------ Start of routine ---------------------------
         int ndat=Data.cols();
-        int npsi=Data.rows();
+        int npsi=Data.rows()+1; //add a constant
         RV thr ( "1","{ theta_r }",vec_1 ( npsi ),vec_1 ( 0 ) );
 
-        mat V0=1e-8*eye ( npsi ); V0 ( 0,0 ) *=100.0;
-        double nu0=npsi;
+        mat V0=1e-8*eye ( npsi ); V0 ( 0,0 ) =10.0;
+        vec rgr(npsi);
+        double nu0=npsi+2;
+        // fitting a linear part => third coef is "1"
+        rgr(2) = 1.0;
+
+        // AR model
         ARX Ar ( thr,V0,nu0 );
-        // estimate
-        for (int i=0;i<ndat;i++ ) {Ar.bayes ( Data.get_col ( i ) );     }
+        // start at i0
+        bool gohi=true; bool golow=true;
+        // certain period 
+        for (int i=ilow;i<ihi;i++){
+                rgr.set_subvector(0,Data.get_col ( i ));
+                Ar.bayes ( rgr );
+        }
+        double tll=Ar._tll();
+        // go hi
+        do {
+                rgr.set_subvector(0,Data.get_col ( ++ihi ));
+                Ar.bayes (rgr);
+                if (Ar._tll()>tll) tll=Ar._tll();
+                else gohi=false;
+        } while (gohi);
+        // go low
+        do { 
+                rgr.set_subvector(0,Data.get_col ( --ilow ));
+                Ar.bayes (rgr);
+                if (Ar._tll()>tll) tll=Ar._tll();
+                else golow=false;
+        } while (golow);
 
         vec th = Ar._epdf().mean();
+
 
         // ------------------ End of routine -----------------------------
@@ -35 +64 @@
 
         if ( n_output>1 ) {
-                double ll=Ar._tll();
                 output[1] = mxCreateDoubleMatrix ( 1,1, mxREAL );
-                double2mxArray ( ll, output[1] );
+                double2mxArray ( ilow, output[1] );
+        }
+        if ( n_output>1 ) {
+                output[2] = mxCreateDoubleMatrix ( 1,1, mxREAL );
+                double2mxArray ( ihi, output[2] );
         }
         
-        if ( n_output>2 ) {
-                ivec str = Ar.structure_est ( egiw ( thr,V0,nu0 ) );
-                output[2] = mxCreateDoubleMatrix ( 1,str.length(), mxREAL );
-                ivec2mxArray ( str, output[2] );
-        }
 }