testjbunit.pas

maths lib with source

unit testjbunit;
interface
uses Math, Ap, Sysutils, jarquebera;

function TestJBTest(Silent : Boolean):Boolean;
function testjbunit_test_silent():Boolean;
function testjbunit_test():Boolean;

implementation

procedure GenerateNonNormal1(N : Integer;
     Med : Double;
     Width : Double;
     var R : TReal1DArray);forward;
procedure GenerateNonNormal2(N : Integer;
     Med : Double;
     Width : Double;
     var R : TReal1DArray);forward;
procedure GenerateNormal(N : Integer;
     Mean : Double;
     Sigma : Double;
     var R : TReal1DArray);forward;


function TestJBTest(Silent : Boolean):Boolean;
var
    PassCount : Integer;
    DetailedReport : Boolean;
    SigmaThreshold : Double;
    N : Integer;
    Pass : Integer;
    DTask : Integer;
    NPos : Integer;
    NCnt : Integer;
    NTbl : TInteger1DArray;
    QCnt : Integer;
    BTQTbl : TReal1DArray;
    BTSigmaTbl : TReal1DArray;
    BTTbl : TReal1DArray;
    BT : Double;
    BPower : Double;
    X : TReal1DArray;
    TmpTestOK : Boolean;
    CTestOK : Boolean;
    PowerTestOK : Boolean;
    WasErrors : Boolean;
    I : Integer;
    J : Integer;
    K : Integer;
    Mean : Double;
    Sigma : Double;
    V : Double;
begin
    WasErrors := False;
    DetailedReport := True;
    PassCount := 20000;
    SigmaThreshold := 5;
    N := 100+RandomInteger(50);
    
    //
    // Prepare place
    //
    SetLength(X, N-1+1);
    
    //
    // Prepare quantiles tables for both-tail test and left/right test.
    // Note that since sign test statistic has discrete distribution
    // we should carefully select quantiles. Equation P(tail<=alpha)=alpha
    // holds only for some specific alpha.
    //
    QCnt := 5;
    SetLength(BTTbl, QCnt-1+1);
    SetLength(BTQTbl, QCnt-1+1);
    SetLength(BTSigmaTbl, QCnt-1+1);
    BTQTbl[0] := 0.25;
    BTQTbl[1] := 0.15;
    BTQTbl[2] := 0.10;
    BTQTbl[3] := 0.05;
    BTQTbl[4] := 0.01;
    I:=0;
    while I<=QCnt-1 do
    begin
        BTSigmaTbl[I] := Sqrt(BTQTbl[I]*(1-BTQTbl[I])/PassCount);
        Inc(I);
    end;
    
    //
    // Report header
    //
    if  not Silent then
    begin
        Write(Format('TESTING JARQUE BERA TEST'#13#10'',[]));
    end;
    
    //
    // Continuous distribuiton test for validity.
    //
    
    //
    // Prepare p-table
    //
    I:=0;
    while I<=QCnt-1 do
    begin
        BTTbl[I] := 0;
        Inc(I);
    end;
    
    //
    // Fill p-table
    //
    Pass:=1;
    while Pass<=PassCount do
    begin
        Mean := 20*RandomReal-10;
        Sigma := 1+2*RandomReal;
        GenerateNormal(N, Mean, Sigma, X);
        JarqueBeraTest(X, N, BT);
        I:=0;
        while I<=QCnt-1 do
        begin
            if BT<=BTQTbl[I] then
            begin
                BTTbl[I] := BTTbl[I]+1/PassCount;
            end;
            Inc(I);
        end;
        Inc(Pass);
    end;
    
    //
    // Check
    //
    TmpTestOK := True;
    I:=0;
    while I<=QCnt-1 do
    begin
        if AbsReal(BTTbl[I]-BTQTbl[I])/BTSigmaTbl[I]>SigmaThreshold then
        begin
            TmpTestOK := False;
        end;
        Inc(I);
    end;
    CTestOK := TmpTestOK;
    
    //
    // Report
    //
    if  not Silent and DetailedReport then
    begin
        Write(Format(''#13#10'CONTINUOUS TEST TABLE'#13#10'',[]));
        Write(Format('Q no.     err.'#13#10'',[]));
        I:=0;
        while I<=QCnt-1 do
        begin
            Write(Format('%1d         %3.1f std.'#13#10'',[
                I+1,
                AbsReal(BTTbl[I]-BTQTbl[I])/BTSigmaTbl[I]]));
            Inc(I);
        end;
        if TmpTestOK then
        begin
            Write(Format('TEST PASSED'#13#10''#13#10'',[]));
        end
        else
        begin
            Write(Format('TEST FAILED'#13#10''#13#10'',[]));
        end;
    end;
    
    //
    // Power test
    //
    
    //
    // Prepare p-table
    //
    I:=0;
    while I<=QCnt-1 do
    begin
        BTTbl[I] := 0;
        Inc(I);
    end;
    
    //
    // Fill p-table
    //
    BPower := 0;
    Pass:=1;
    while Pass<=PassCount do
    begin
        if Pass mod 2=0 then
        begin
            GenerateNonNormal1(N, 20*RandomReal-10, 5, X);
        end
        else
        begin
            GenerateNonNormal2(N, 20*RandomReal-10, 5, X);
        end;
        JarqueBeraTest(X, N, BT);
        if BT<0.05 then
        begin
            BPower := BPower+1/PassCount;
        end;
        Inc(Pass);
    end;
    
    //
    // Check
    //
    PowerTestOK := BPower>0.95;
    
    //
    // Report
    //
    if  not Silent and DetailedReport then
    begin
        Write(Format(''#13#10'POWER TEST TABLE'#13#10'',[]));
        Write(Format('TEST POWER:  %4.2f'#13#10'',[
            BPower]));
        if PowerTestOK then
        begin
            Write(Format('TEST PASSED'#13#10''#13#10'',[]));
        end
        else
        begin
            Write(Format('TEST FAILED'#13#10''#13#10'',[]));
        end;
    end;
    
    //
    // Final report
    //
    WasErrors :=  not CTestOK or  not PowerTestOK;
    if  not Silent then
    begin
        Write(Format('CONTINUOUS VALIDITY TEST:                ',[]));
        if CTestOK then
        begin
            Write(Format('OK'#13#10'',[]));
        end
        else
        begin
            Write(Format('FAILED'#13#10'',[]));
        end;
        Write(Format('POWER TEST:                              ',[]));
        if PowerTestOK then
        begin
            Write(Format('OK'#13#10'',[]));
        end
        else
        begin
            Write(Format('FAILED'#13#10'',[]));
        end;
        if WasErrors then
        begin
            Write(Format('TEST SUMMARY: FAILED'#13#10'',[]));
        end
        else
        begin
            Write(Format('TEST SUMMARY: PASSED'#13#10'',[]));
        end;
        Write(Format(''#13#10''#13#10'',[]));
    end;
    Result :=  not WasErrors;
end;


procedure GenerateNonNormal1(N : Integer;
     Med : Double;
     Width : Double;
     var R : TReal1DArray);
var
    I : Integer;
    V : Double;
begin
    i:=0;
    while i<=N-1 do
    begin
        R[I] := Med+Width*(2*RandomInteger(2)-1);
        Inc(i);
    end;
end;


procedure GenerateNonNormal2(N : Integer;
     Med : Double;
     Width : Double;
     var R : TReal1DArray);
var
    I : Integer;
    V : Double;
begin
    i:=0;
    while i<=N-1 do
    begin
        V := 2*RandomReal-1;
        if V>0 then
        begin
            R[I] := Med+Width*V;
        end
        else
        begin
            R[I] := Med+2*Width*V;
        end;
        Inc(i);
    end;
end;


procedure GenerateNormal(N : Integer;
     Mean : Double;
     Sigma : Double;
     var R : TReal1DArray);
var
    I : Integer;
    U : Double;
    V : Double;
    S : Double;
    Sum : Double;
begin
    i := 0;
    while i<n do
    begin
        u := (2*RandomInteger(2)-1)*RandomReal;
        v := (2*RandomInteger(2)-1)*RandomReal;
        sum := u*u+v*v;
        if (sum<1) and (sum>0) then
        begin
            sum := sqrt(-2*ln(sum)/sum);
            if i<n then
            begin
                r[i] := Sigma*u*sum+Mean;
            end;
            if i+1<n then
            begin
                r[i+1] := Sigma*v*sum+Mean;
            end;
            i := i+2;
        end;
    end;
end;


(*************************************************************************
Silent unit test
*************************************************************************)
function testjbunit_test_silent():Boolean;
begin
    Result := TestJBTest(True);
end;


(*************************************************************************
Unit test
*************************************************************************)
function testjbunit_test():Boolean;
begin
    Result := TestJBTest(False);
end;


end.