root/library/tests/test_util.cpp @ 706

#include "test_util.h"
#include "stat/exp_family.h"
#include <fstream>

namespace bdm {

using itpp::vec;
using itpp::mat;

std::string load_test_file ( const char *fname ) {
        char buffer[8192];
        memset ( buffer, 0, sizeof ( buffer ) );
        std::ifstream src ( fname );
        src.read ( buffer, sizeof ( buffer ) - 1 );
        return std::string ( buffer );
}

double normcoef ( const epdf *ep, const vec &xb, const vec &yb,
                  int xn, int yn ) {
        mat Pdf ( xn + 1, yn + 1 );
        vec rgr ( 2 );

        double xstep = ( xb ( 1 ) - xb ( 0 ) ) / xn;
        double ystep = ( yb ( 1 ) - yb ( 0 ) ) / yn;

        double x = xb ( 0 );
        for ( int i = 0; i <= xn; x += xstep, i++ ) {
                rgr ( 0 ) = x;
                double y = yb ( 0 );
                for ( int j = 0; j <= yn; y += ystep, j++ ) {
                        rgr ( 1 ) = y;
                        Pdf ( i, j ) = exp ( ep->evallog ( rgr ) );
                }
        }

        return sumsum ( Pdf ) * xstep * ystep;
}

vec num_mean2 ( const epdf *ep, const vec &xb, const vec &yb,
                int xn, int yn ) {
        mat Pdf ( xn + 1, yn + 1 );
        vec rgr ( 2 );

        double xstep = ( xb ( 1 ) - xb ( 0 ) ) / xn;
        double ystep = ( yb ( 1 ) - yb ( 0 ) ) / yn;

        vec Mu ( xn + 1 );
        vec Si ( yn + 1 );

        double x = xb ( 0 );
        for ( int i = 0; i <= xn; x += xstep, i++ ) {
                Mu ( i ) = x;
                rgr ( 0 ) = x;
                double y = yb ( 0 );
                for ( int j = 0; j <= yn; y += ystep, j++ ) {
                        Si ( j ) = y;
                        rgr ( 1 ) = y;
                        Pdf ( i, j ) = exp ( ep->evallog ( rgr ) );
                }
        }

        vec fm = sum ( Pdf, 2 );
        double sfm = sum ( fm );
        vec fs = sum ( Pdf, 1 );
        double sfs = sum ( fs );
        double vi0 = Mu * fm / sfm;
        double vi1 = Si * fs / sfs;
        return vec_2 ( vi0, vi1 );
}

}