root/library/tests/egiw_test.cpp @ 436

#define BDMLIB // not an ideal way to prevent double registration of UI factories...
#include "base/bdmbase.h"
#include "base/user_info.h"
#include "stat/exp_family.h"
#include "itpp_ext.h"
#include "epdf_harness.h"
#include "mat_checks.h"
#include "UnitTest++.h"

const double epsilon = 0.00001;

namespace UnitTest
{

inline void CheckClose(TestResults &results, const itpp::vec &expected,
                       const itpp::vec &actual, double tolerance,
                       TestDetails const &details) {
    if (!AreClose(expected, actual, tolerance)) { 
        MemoryOutStream stream;
        stream << "Expected " << expected << " +/- " << tolerance << " but was " << actual;

        results.OnTestFailure(details, stream.GetText());
    }
}

}

using namespace bdm;

template<>
const ParticularUI<egiw> &ParticularUI<egiw>::factory(
    ParticularUI<egiw>("egiw"));

TEST(test_egiw) {
    RV::clear_all();
    UIFile in("egiw.cfg");
    Array<epdf_harness *> input;
    UI::get(input, in, "data");
    int sz = input.size();
    CHECK(sz > 0);
    for (int i = 0; i < sz; ++i) {
        input(i)->test();
    }
}

TEST(test_egiw_1_2) {
    // Setup model
    double mu = 1.1; //unit step parametr
    double b = 3.0; // sequence of <1 -1 1 -1...>
    double s = 0.1;


    // TEST 1x1 EGIW
    mat V(3, 3);
    V(0, 0) = pow(mu, 2) + pow(b, 2) + s;
    V(1, 0) = mu;
    V(2, 0) = b;

    V(0, 1) = V(1, 0);
    V(1, 1) = 1.0;
    V(2, 1) = 0.0;

    V(0, 2) = V(2, 0);
    V(1, 2) = V(2, 1);
    V(2, 2) = 1.0;

    double nu = 20;

    egiw E(1, nu * V, nu);    
    CHECK_CLOSE(vec("1.1 3.0 0.142857"), E.mean(), epsilon);
    CHECK_CLOSE(7.36731, E.lognc(), epsilon);

    int n = 100;
    vec rgr(3);

    mat Tmp(2 * n, n);

    double summ = 0.0;
    for (int k = 0; k < n; k++) { // ALL b
        rgr(1) = 1 + k * (1.0 / n) * 4.0;
        for (int i = 0; i < 2*n; i++) { //ALL mu
            rgr(0) = -2 + i * (1.0 / n) * 3.0;
            for (int j = 0; j < n; j++) { // All sigma
                rgr(2) = (j + 1) * (1.0 / n) * 2.0;

                Tmp(i, j) = E.evallog(rgr);
            }
        }
        summ += sumsum(exp(Tmp)) / n / n / n * 3.0 * 2.0 * 4.0;
    }

    CHECK_CLOSE(1.0, summ, 0.1);
}