testmwuunit.pas

maths lib with source

unit testmwuunit;
interface
uses Math, Ap, Sysutils, normaldistr, mannwhitneyu;

function TestMWUTest(Silent : Boolean):Boolean;
function testmwuunit_test_silent():Boolean;
function testmwuunit_test():Boolean;

implementation

procedure Generate(N : Integer;
     Med : Double;
     Width : Double;
     var R : TReal1DArray);forward;


function TestMWUTest(Silent : Boolean):Boolean;
var
    PassCount : Integer;
    DetailedReport : Boolean;
    SigmaThreshold : Double;
    N : Integer;
    M : Integer;
    Pass : Integer;
    DTask : Integer;
    NPos : Integer;
    NCnt : Integer;
    NTbl : TInteger1DArray;
    QCnt : Integer;
    BTQTbl : TReal1DArray;
    LRQTbl : TReal1DArray;
    BTSigmaTbl : TReal1DArray;
    LRSigmaTbl : TReal1DArray;
    BTTbl : TReal1DArray;
    LTTbl : TReal1DArray;
    RTTbl : TReal1DArray;
    BT : Double;
    LT : Double;
    RT : Double;
    LPower : Double;
    RPower : Double;
    BPower : Double;
    X : TReal1DArray;
    Y : TReal1DArray;
    TmpTestOK : Boolean;
    CTestOK : Boolean;
    PowerTestOK : Boolean;
    WasErrors : Boolean;
    I : Integer;
    J : Integer;
    K : Integer;
    Med : Double;
    Width : Double;
    V : Double;
begin
    WasErrors := False;
    DetailedReport := True;
    PassCount := 20000;
    SigmaThreshold := 5;
    N := 25+RandomInteger(15);
    M := 25+RandomInteger(15);
    
    //
    // Prepare place
    //
    SetLength(X, N-1+1);
    SetLength(Y, M-1+1);
    
    //
    // Prepare quantiles tables for both-tail test and left/right test.
    //
    QCnt := 5;
    SetLength(BTTbl, QCnt-1+1);
    SetLength(LTTbl, QCnt-1+1);
    SetLength(RTTbl, QCnt-1+1);
    SetLength(BTQTbl, QCnt-1+1);
    SetLength(LRQTbl, QCnt-1+1);
    SetLength(BTSigmaTbl, QCnt-1+1);
    SetLength(LRSigmaTbl, 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
        LRQTbl[I] := BTQTbl[I];
        Inc(I);
    end;
    I:=0;
    while I<=QCnt-1 do
    begin
        BTSigmaTbl[I] := Sqrt(BTQTbl[I]*(1-BTQTbl[I])/PassCount);
        LRSigmaTbl[I] := Sqrt(LRQTbl[I]*(1-LRQTbl[I])/PassCount);
        Inc(I);
    end;
    
    //
    // Report header
    //
    if  not Silent then
    begin
        Write(Format('TESTING MANN-WHITNEY U TEST'#13#10'',[]));
    end;
    
    //
    // Continuous distribuiton test for validity.
    //
    
    //
    // Prepare p-table
    //
    I:=0;
    while I<=QCnt-1 do
    begin
        BTTbl[I] := 0;
        LTTbl[I] := 0;
        RTTbl[I] := 0;
        Inc(I);
    end;
    
    //
    // Fill p-table
    //
    Pass:=1;
    while Pass<=PassCount do
    begin
        Med := 20*RandomReal-10;
        Width := 1+5*RandomReal;
        Generate(N, Med, Width, X);
        Generate(M, Med, Width, Y);
        MannWhitneyUTest(X, N, Y, M, BT, LT, RT);
        I:=0;
        while I<=QCnt-1 do
        begin
            if BT<=BTQTbl[I] then
            begin
                BTTbl[I] := BTTbl[I]+1/PassCount;
            end;
            if LT<=LRQTbl[I] then
            begin
                LTTbl[I] := LTTbl[I]+1/PassCount;
            end;
            if RT<=LRQTbl[I] then
            begin
                RTTbl[I] := RTTbl[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;
        if AbsReal(RTTbl[I]-LRQTbl[I])/LRSigmaTbl[I]>SigmaThreshold then
        begin
            TmpTestOK := False;
        end;
        if AbsReal(LTTbl[I]-LRQTbl[I])/LRSigmaTbl[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.     BT err.   LT err.   RT err.   '#13#10'',[]));
        I:=0;
        while I<=QCnt-1 do
        begin
            Write(Format('%1d         %3.1f std.  %3.1f std.  %3.1f std.  '#13#10'',[
                I+1,
                AbsReal(BTTbl[I]-BTQTbl[I])/BTSigmaTbl[I],
                AbsReal(LTTbl[I]-LRQTbl[I])/LRSigmaTbl[I],
                AbsReal(RTTbl[I]-LRQTbl[I])/LRSigmaTbl[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;
        LTTbl[I] := 0;
        RTTbl[I] := 0;
        Inc(I);
    end;
    
    //
    // Fill p-table
    //
    LPower := 0;
    RPower := 0;
    BPower := 0;
    Pass:=1;
    while Pass<=PassCount do
    begin
        V := 1.5;
        Generate(N, 0, 1.0, X);
        Generate(M, V, 1.0, Y);
        MannWhitneyUTest(X, N, Y, M, BT, LT, RT);
        if LT<0.05 then
        begin
            LPower := LPower+1/PassCount;
        end;
        Generate(M, -V, 1.0, Y);
        MannWhitneyUTest(X, N, Y, M, BT, LT, RT);
        if RT<0.05 then
        begin
            RPower := RPower+1/PassCount;
        end;
        Generate(M, V*(2*RandomInteger(2)-1), 1.0, Y);
        MannWhitneyUTest(X, N, Y, M, BT, LT, RT);
        if BT<0.05 then
        begin
            BPower := BPower+1/PassCount;
        end;
        Inc(Pass);
    end;
    
    //
    // Check
    //
    PowerTestOK := (LPower>0.95) and (RPower>0.95) and (BPower>0.95);
    
    //
    // Report
    //
    if  not Silent and DetailedReport then
    begin
        Write(Format(''#13#10'POWER TEST TABLE'#13#10'',[]));
        Write(Format('LEFT-TAIL TEST POWER:  %4.2f'#13#10'',[
            LPower]));
        Write(Format('RIGHT-TAIL TEST POWER: %4.2f'#13#10'',[
            RPower]));
        Write(Format('BOTH-TAIL 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 Generate(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+3*Width*V;
        end;
        Inc(i);
    end;
end;


(*************************************************************************
Silent unit test
*************************************************************************)
function testmwuunit_test_silent():Boolean;
begin
    Result := TestMWUTest(True);
end;


(*************************************************************************
Unit test
*************************************************************************)
function testmwuunit_test():Boolean;
begin
    Result := TestMWUTest(False);
end;


end.