root/matlab/mex/linefit.cpp @ 248

#include <itpp/itmex.h>

#include "../../bdm/estim/arx.h"

using namespace itpp;

void mexFunction ( int n_output, mxArray *output[], int n_input, const mxArray *input[] ) {
        // Check the number of inputs and output arguments
        if ( n_output<1 ) mexErrMsgTxt ( "Wrong number of output 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()+1; //add a constant
        RV thr ( "1","{ theta_r }",vec_1 ( npsi ),vec_1 ( 0 ) );

        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 );
        // 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 -----------------------------

        // Create output vectors
        output[0] = mxCreateDoubleMatrix ( 1,th.length(), mxREAL );
        // Convert the IT++ format to Matlab format for output
        vec2mxArray ( th, output[0] );

        if ( n_output>1 ) {
                output[1] = mxCreateDoubleMatrix ( 1,1, mxREAL );
                double2mxArray ( ilow, output[1] );
        }
        if ( n_output>1 ) {
                output[2] = mxCreateDoubleMatrix ( 1,1, mxREAL );
                double2mxArray ( ihi, output[2] );
        }
        
}