testcorrhtunit.pas

maths lib with source

unit testcorrhtunit;
interface
uses Math, Ap, Sysutils, gammaf, normaldistr, ibetaf, studenttdistr, correlation, correlationtests;

function TestCorrTest(Silent : Boolean):Boolean;
function testcorrhtunit_test_silent():Boolean;
function testcorrhtunit_test():Boolean;

implementation

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


function TestCorrTest(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;
    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;
    PTestOK : Boolean;
    STestOK : Boolean;
    PTestPowOK : Boolean;
    STestPowOK : Boolean;
    WasErrors : Boolean;
    I : Integer;
    J : Integer;
    K : Integer;
    Med : Double;
    Width : Double;
    V : Double;
    R : Double;
begin
    WasErrors := False;
    DetailedReport := True;
    PassCount := 20000;
    SigmaThreshold := 5;
    N := 40;
    
    //
    // Prepare place
    //
    SetLength(X, N-1+1);
    SetLength(Y, N-1+1);
    
    //
    // Prepare quantiles tables for both-tail test and left/right test.
    // Note that since Spearman test statistic has discrete distribution
    // we should carefully select quantiles. Equation P(tail<=alpha)=alpha
    // holds only for some specific alpha.
    //
    QCnt := 4;
    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;
    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 CORRELATION TEST'#13#10'',[]));
    end;
    
    //
    // Pearson correlation 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;
        GenerateNormal(N, Med, Width, X);
        GenerateNormal(N, Med, Width, Y);
        R := PearsonCorrelation(X, Y, N);
        PearsonCorrelationSignificance(R, N, 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;
    PTestOK := TmpTestOK;
    
    //
    // Report
    //
    if  not Silent and DetailedReport then
    begin
        Write(Format(''#13#10'PEARSON 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;
    
    //
    // Pearson 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 := 10;
        GenerateNormal(N, 0.0, 1.0, X);
        GenerateNormal(N, 0.0, 1.0, Y);
        APVSub(@Y[0], 0, N-1, @X[0], 0, N-1, V);
        R := PearsonCorrelation(X, Y, N);
        PearsonCorrelationSignificance(R, N, BT, LT, RT);
        if LT<0.05 then
        begin
            LPower := LPower+1/PassCount;
        end;
        GenerateNormal(N, 0.0, 1.0, Y);
        APVAdd(@Y[0], 0, N-1, @X[0], 0, N-1, V);
        R := PearsonCorrelation(X, Y, N);
        PearsonCorrelationSignificance(R, N, BT, LT, RT);
        if RT<0.05 then
        begin
            RPower := RPower+1/PassCount;
        end;
        GenerateNormal(N, 0.0, 1.0, Y);
        if RandomReal>0.5 then
        begin
            APVAdd(@Y[0], 0, N-1, @X[0], 0, N-1, V);
        end
        else
        begin
            APVSub(@Y[0], 0, N-1, @X[0], 0, N-1, V);
        end;
        R := PearsonCorrelation(X, Y, N);
        PearsonCorrelationSignificance(R, N, BT, LT, RT);
        if BT<0.05 then
        begin
            BPower := BPower+1/PassCount;
        end;
        Inc(Pass);
    end;
    
    //
    // Check
    //
    PTestPowOK := (LPower>0.95) and (RPower>0.95) and (BPower>0.95);
    
    //
    // Report
    //
    if  not Silent and DetailedReport then
    begin
        Write(Format(''#13#10'PEARSON 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 PTestPowOK then
        begin
            Write(Format('TEST PASSED'#13#10''#13#10'',[]));
        end
        else
        begin
            Write(Format('TEST FAILED'#13#10''#13#10'',[]));
        end;
    end;
    
    //
    // Spearman correlation test for validity.
    //
    
    //
    // Prepare p-table
    //
    I:=0;
    while I<=QCnt-1 do
    begin
        BTTbl[I] := 0;