#include <oxstd.h>
#include <oxdraw.h>

#include <../lib/normtest.ox>

#import <simula>

/*-------------------------- test statistics -------------------------------*/
MVNbowmanShenton(mX)                             /* compute asymptotic form */
{
    decl ct, vskew, vkurt;

    ct = rows(mX);
    mX = IN01Data(mX);

    vskew = (sumc(mX .^ 3) / ct) * sqrt(ct / 6);
    vkurt = (sumc(mX .^ 4) / ct - 3) * sqrt(ct / 24);

return double(sumr(vskew) + sumr(vkurt))
     | double(sumr(vskew .^ 2) + sumr(vkurt .^ 2));
}
MVNdoornikHansen(mX)                              /* compute corrected form */
{
    decl ct, vskew, vkurt, newskew, newkurt;

    ct = rows(mX);
    mX = IN01Data(mX);

    vskew = sumc(mX .^ 3) / ct;
    vkurt = sumc(mX .^ 4) / ct;

    newskew = SkewSU(vskew, ct);         /* transform skewness to normality */
    newkurt = KurtGamma(vskew, vkurt, ct);  /* transform kurtosis via gamma */

return double(sumr(newskew) + sumr(newkurt))
     | double(sumr(newskew .^ 2) + sumr(newkurt .^ 2));
}
/*------------------------ ENDtest statistics -----------------------------*/

/*---------------------- SimNormal : Simulation ----------------------------*/
class SimNormal : Simulation         // inherit from simulation class
{
    decl m_mTest;                              // test statistic
    decl m_cDim;                               // dimension of data
    SimNormal(const cDim);                     // constructor
    ~SimNormal();                              // destructor
    Generate(const iRep, const cT, const mxT); // generate replication
    GetPvalues();                              // return p-values of tests
    GetTestStatistics();                       // return test statistics
    Plot(const iRep);                          // plot function
};
SimNormal::SimNormal(const cDim)
{
    m_cDim = cDim;
    Simulation(<20>, 20, 2000, TRUE, -1, <0.2,0.1,0.05,0.01>, 0);
    SetTestNames({"normal BS", "chi-sq BS", "normal DH", "chi-sq DH"});
}
SimNormal::~SimNormal()
{
}
SimNormal::Plot(const iRep)
{
    if (iRep >= 100 && imod(iRep + 1, 100) == 0)
    {
        DrawDensity(0, m_mTestVal[0][:iRep - 1], "", FALSE, TRUE, TRUE);
        DrawDensity(2, m_mTestVal[2][:iRep - 1], "", FALSE, TRUE, TRUE);

        DrawQQ(1, m_mTestVal[1][:iRep - 1], "", QQ_CHI, 2, 0);
        DrawQQ(3, m_mTestVal[3][:iRep - 1], "", QQ_CHI, 2, 0);

        ShowDrawWindow();
    }
}
SimNormal::Generate(const iRep, const cT, const mxT)
{
    decl mx = rann(cT, m_cDim);                            /* generate data */

    Plot(iRep);

    m_mTest = MVNbowmanShenton(mx)                              /* do tests */
          | MVNdoornikHansen(mx);

return 1;                                 /* 1 indicates success, 0 failure */
}
SimNormal::GetPvalues()
{
    decl m = new matrix[4][1];
                                                   /* chi^2(2n) distributed */
    m[<1,3>][] = tailchi(m_mTest[<1,3>][], 2 * m_cDim);
                                                     /* N(0,2n) distributed */
    m[<0,2>][] = 2 * tailn( fabs(m_mTest[<0,2>][] / sqrt(2 * m_cDim)) );

    return m;
}
SimNormal::GetTestStatistics()
{
    return m_mTest;
}
/*-------------------- END SimNormal : Simulation --------------------------*/

main()
{
    decl experiment = new SimNormal(1);                       // create object

    experiment.Simulate();                                   // do simulations

    delete experiment;                                        // remove object
}

