root/bdm/estim/arx.cpp @ 204

#include "arx.h"

using namespace itpp;

void ARX::bayes ( const vec &dt, const double w ) {
        double lnc;

        if ( frg<1.0 ) {
                est.pow ( frg );
                if ( evalll ) {
                        last_lognc = est.lognc();
                }
        }
        V.opupdt ( dt,w );
        nu+=w;

        // log(sqrt(2*pi)) = 0.91893853320467
        if ( evalll ) {
                lnc = est.lognc();
                ll = lnc - last_lognc - 0.91893853320467;
                last_lognc = lnc;
        }
}

double ARX::logpred ( const vec &dt ) const {
        egiw pred ( est );
        ldmat &V=pred._V();
        double &nu=pred._nu();

        double lll;

        if ( frg<1.0 ) {
                pred.pow ( frg );
                lll = pred.lognc();
        }
        else//should be save: last_lognc is changed only by bayes;
                if ( evalll ) {lll=last_lognc;}
                else{lll=pred.lognc();}

        V.opupdt ( dt,1.0 );
        nu+=1.0;
        // log(sqrt(2*pi)) = 0.91893853320467
        return pred.lognc()-lll - 0.91893853320467;
}

ARX* ARX::_copy_ ( bool changerv ) {
        ARX* Tmp=new ARX ( *this );
        if ( changerv ) {Tmp->rv.newids(); Tmp->est._renewrv ( Tmp->rv );}
        return Tmp;
}

void ARX::set_statistics ( const BMEF* B0 ) {
        const ARX* A0=dynamic_cast<const ARX*> ( B0 );

        it_assert_debug ( V.rows() ==A0->V.rows(),"ARX::set_statistics Statistics  differ" );
        set_parameters ( A0->V,A0->nu );
}

enorm<ldmat>* ARX::predictor ( const RV &rv, const vec &rgr ) const {
        mat mu ( rv.count(), V.rows()-rv.count() );
        mat R ( rv.count(),rv.count() );
        enorm<ldmat>* tmp;
        tmp=new enorm<ldmat> ( rv );

        est.mean_mat ( mu,R ); //mu =
        //correction for student-t  -- TODO check if correct!!
        //R*=nu/(nu-2);
        mat p_mu=mu.T() *rgr;   //the result is one column
        tmp->set_parameters ( p_mu.get_col ( 0 ),ldmat ( R ) );
        return tmp;
}

mlnorm<ldmat>* ARX::predictor ( const RV &rv, const RV &rvc ) const {
        int dif=V.rows() - rv.count() - rvc.count();
        it_assert_debug ( ( dif==0 ) || ( dif==1 ), "Give RVs do not match" );

        mat mu ( rv.count(), V.rows()-rv.count() );
        mat R ( rv.count(),rv.count() );
        mlnorm<ldmat>* tmp;
        tmp=new mlnorm<ldmat> ( rv,rvc );

        est.mean_mat ( mu,R ); //mu =
        mu = mu.T();
        //correction for student-t  -- TODO check if correct!!

        if ( dif==0 ) { // no constant term
                tmp->set_parameters ( mu, zeros ( rv.count() ), ldmat ( R ) );
        }
        else {
                //Assume the constant term is the last one:
                tmp->set_parameters ( mu.get_cols (0,mu.cols()-2 ), mu.get_col ( mu.cols()-1 ), ldmat ( R ) );
        }
        return tmp;
}

mlstudent* ARX::predictor_student ( const RV &rv, const RV &rvc ) const {
        int dif=V.rows() - rv.count() - rvc.count();
        it_assert_debug ( ( dif==0 ) || ( dif==1 ), "Give RVs do not match" );

        mat mu ( rv.count(), V.rows()-rv.count() );
        mat R ( rv.count(),rv.count() );
        mlstudent* tmp;
        tmp=new mlstudent ( rv,rvc );

        est.mean_mat ( mu,R ); //
        mu = mu.T();
        
        int xdim = rv.count();
        int end = V._L().rows()-1;
        ldmat Lam ( V._L() ( xdim,end,xdim,end ), V._D() ( xdim,end ) );  //exp val of R


        if ( dif==0 ) { // no constant term
                tmp->set_parameters ( mu, zeros ( rv.count() ), ldmat ( R ), Lam);
        }
        else {
                //Assume the constant term is the last one:
                tmp->set_parameters ( mu.get_cols (0,mu.cols()-2 ), mu.get_col ( mu.cols()-1 ), ldmat ( R ), Lam);
        }
        return tmp;
}

/*! \brief Return the best structure
@param Eg a copy of GiW density that is being examined
@param Eg0 a copy of prior GiW density before estimation
@param Egll likelihood of the current Eg
@param indeces current indeces
\return best likelihood in the structure below the given one
*/
double egiw_bestbelow ( egiw Eg, egiw Eg0, double Egll, ivec &indeces ) { //parameter Eg is a copy!
        ldmat Vo = Eg._V(); //copy
        ldmat Vo0 = Eg._V(); //copy
        ldmat& Vp = Eg._V(); // pointer into Eg
        ldmat& Vp0 = Eg._V(); // pointer into Eg
        int end = Vp.rows()-1;
        int i;
        mat Li;
        mat Li0;
        double maxll=Egll;
        double tmpll=Egll;
        double belll=Egll;

        ivec tmpindeces;
        ivec maxindeces=indeces;


        cout << "bb:(" << indeces <<") ll=" << Egll <<endl;

        //try to remove only one rv
        for ( i=0;i<end;i++ ) {
                //copy original
                Li = Vo._L();
                Li0 = Vo0._L();
                //remove stuff
                Li.del_col ( i+1 );
                Li0.del_col ( i+1 );
                Vp.ldform ( Li,Vo._D() );
                Vp0.ldform ( Li0,Vo0._D() );
                tmpll = Eg.lognc()-Eg0.lognc(); // likelihood is difference of norm. coefs.

                cout << "i=(" << i <<") ll=" << tmpll <<endl;

                //
                if ( tmpll > Egll ) { //increase of the likelihood
                        tmpindeces = indeces;
                        tmpindeces.del ( i );
                        //search for a better match in this substructure
                        belll=egiw_bestbelow ( Eg, Eg0, tmpll, tmpindeces );
                        if ( belll>maxll ) { //better match found
                                maxll = belll;
                                maxindeces = tmpindeces;
                        }
                }
        }
        indeces = maxindeces;
        return maxll;
}

ivec ARX::structure_est ( egiw est0 ) {
        ivec ind=linspace ( 1,rv.count()-1 );
        egiw_bestbelow ( est, est0, est.lognc()- est0.lognc(), ind );
        return ind;
}