root/matlab/mex/linefit2.cpp @ 133

Revision 133, 1.9 kB (checked in by smidl, 16 years ago)

linefit pro Martina

  • Property svn:eol-style set to native
RevLine 
[133]1#include <itpp/itmex.h>
2
3#include "../../bdm/estim/arx.h"
4
5using namespace itpp;
6
7void mexFunction ( int n_output, mxArray *output[], int n_input, const mxArray *input[] ) {
8        // Check the number of inputs and output arguments
9        if ( n_output<1 ) mexErrMsgTxt ( "Wrong number of output variables!" );
10        if ( n_input!=3 ) mexErrMsgTxt ( "Wrong number of input variables!" );
11
12        // Convert input variables to IT++ format
13        mat Data = mxArray2mat ( input[0] );
14        if (Data.rows()!=2)
15                { if (Data.cols()==2) Data=Data.T(); else mexErrMsgTxt("Data are not 2D!");}
16        int ilow = mxArray2int ( input[1] )-1;
17        int ihi = mxArray2int ( input[2] )-1; //coreection for different indeces
18
19        // ------------------ Start of routine ---------------------------
20        // ------------------ Start of routine ---------------------------
21        int ndat=Data.cols();
22        int npsi=Data.rows()+1; //add a constant
23        RV thr ( "1","{ theta_r }",vec_1 ( npsi ),vec_1 ( 0 ) );
24
25        mat V0=1e-8*eye ( npsi ); V0 ( 0,0 ) =10.0;
26        vec rgr(npsi);
27        double nu0=2;
28        // fitting a linear part => third coef is "1"
29        rgr(2) = 1.0;
30
31        // RESULTS
32        mat Tlls=zeros(ihi,ihi);
33        // AR model
34        ARX Ar ( thr,V0,nu0 );
35
36        //Initialize
37        for (int i=ilow; i<ihi; i++) {
38                rgr.set_subvector(0,Data.get_col ( i ));
39                Ar.bayes(rgr);
40        }
41        Ar.get_parameters(V0,nu0);
42        V0 /=V0(1,1);
43        nu0=3;
44        Ar.set_parameters(V0,nu0);
45               
46        // FLATTEN
47        for (int i=ilow; i<ihi; i++) {
48                Ar.set_parameters(V0,nu0);
49                for (int t=i+1; t<ihi; t++) {
50                        rgr.set_subvector(0,Data.get_col ( t ));
51                        Ar.bayes(rgr);
52                        Tlls(i,t)=Ar._tll();
53                }
54        }
55
56        max_index(Tlls,ilow,ihi);
57
58        // ------------------ End of routine -----------------------------
59
60        output[0] = mxCreateDoubleMatrix ( 2,1, mxREAL );
61        ivec2mxArray (vec_2(ilow+1,ihi+1) , output[0] );
62
63        if ( n_output>1 ) {
64        // Create output vectors
65        output[1] = mxCreateDoubleMatrix ( Tlls.rows(),Tlls.cols(), mxREAL );
66        // Convert the IT++ format to Matlab format for output
67        mat2mxArray ( Tlls, output[1] );
68        }
69       
70}
Note: See TracBrowser for help on using the browser.