teststudentttestsunit.pas

maths lib with source

unit teststudentttestsunit;
interface
uses Math, Ap, Sysutils, gammaf, normaldistr, ibetaf, studenttdistr, studentttests;

function TestStudentT(Silent : Boolean):Boolean;
function teststudentttestsunit_test_silent():Boolean;
function teststudentttestsunit_test():Boolean;

implementation

procedure TestStudentTGenNormRnd(N : Integer;
     Mean : Double;
     Sigma : Double;
     var R : TReal1DArray);forward;


function TestStudentT(Silent : Boolean):Boolean;
var
    Pass : Integer;
    PassCount : Integer;
    MaxN : Integer;
    N : Integer;
    M : Integer;
    I : Integer;
    J : Integer;
    QCnt : Integer;
    X : TReal1DArray;
    Y : TReal1DArray;
    QTbl : TReal1DArray;
    PTbl : TReal1DArray;
    LPTbl : TReal1DArray;
    RPTbl : TReal1DArray;
    SigmaTbl : TReal1DArray;
    Bt : Double;
    LT : Double;
    RT : Double;
    V : Double;
    WasErrors : Boolean;
    SigmaThreshold : Double;
begin
    WasErrors := False;
    MaxN := 1000;
    PassCount := 20000;
    SigmaThreshold := 5;
    SetLength(X, MaxN-1+1);
    SetLength(Y, MaxN-1+1);
    QCnt := 8;
    SetLength(PTbl, QCnt-1+1);
    SetLength(LPTbl, QCnt-1+1);
    SetLength(RPTbl, QCnt-1+1);
    SetLength(QTbl, QCnt-1+1);
    SetLength(SigmaTbl, QCnt-1+1);
    QTbl[0] := 0.25;
    QTbl[1] := 0.15;
    QTbl[2] := 0.10;
    QTbl[3] := 0.05;
    QTbl[4] := 0.04;
    QTbl[5] := 0.03;
    QTbl[6] := 0.02;
    QTbl[7] := 0.01;
    I:=0;
    while I<=QCnt-1 do
    begin
        SigmaTbl[I] := Sqrt(QTbl[I]*(1-QTbl[I])/PassCount);
        Inc(I);
    end;
    if  not Silent then
    begin
        Write(Format('TESTING STUDENT T'#13#10'',[]));
    end;
    
    //
    // 1-sample test
    //
    if  not Silent then
    begin
        Write(Format('Testing 1-sample test for 1-type errors'#13#10'',[]));
    end;
    N := 15;
    I:=0;
    while I<=QCnt-1 do
    begin
        PTbl[I] := 0;
        LPTbl[I] := 0;
        RPTbl[I] := 0;
        Inc(I);
    end;
    Pass:=1;
    while Pass<=PassCount do
    begin
        
        //
        // Both tails
        //
        TestStudentTGenNormRnd(N, 0.0, 1+4*RandomReal, X);
        StudentTTest1(X, N, 0, BT, LT, RT);
        I:=0;
        while I<=QCnt-1 do
        begin
            if BT<=QTbl[I] then
            begin
                PTbl[I] := PTbl[I]+1/PassCount;
            end;
            Inc(I);
        end;
        
        //
        // Left tail
        //
        TestStudentTGenNormRnd(N, 0.5*RandomReal, 1+4*RandomReal, X);
        StudentTTest1(X, N, 0, BT, LT, RT);
        I:=0;
        while I<=QCnt-1 do
        begin
            if LT<=QTbl[I] then
            begin
                LPTbl[I] := LPTbl[I]+1/PassCount;
            end;
            Inc(I);
        end;
        
        //
        // Right tail
        //
        TestStudentTGenNormRnd(N, -0.5*RandomReal, 1+4*RandomReal, X);
        StudentTTest1(X, N, 0, BT, LT, RT);
        I:=0;
        while I<=QCnt-1 do
        begin
            if RT<=QTbl[I] then
            begin
                RPTbl[I] := RPTbl[I]+1/PassCount;
            end;
            Inc(I);
        end;
        Inc(Pass);
    end;
    if  not Silent then
    begin
        Write(Format('Expect. Both    Left    Right'#13#10'',[]));
        I:=0;
        while I<=QCnt-1 do
        begin
            Write(Format('%4.1f%%   %4.1f%%   %4.1f%%   %4.1f%%   '#13#10'',[
                QTbl[I]*100,
                PTbl[I]*100,
                LPTbl[I]*100,
                RPTbl[I]*100]));
            Inc(I);
        end;
    end;
    I:=0;
    while I<=QCnt-1 do
    begin
        WasErrors := WasErrors or (AbsReal(PTbl[I]-QTbl[I])/SigmaTbl[I]>SigmaThreshold);
        WasErrors := WasErrors or ((LPTbl[I]-QTbl[I])/SigmaTbl[I]>SigmaThreshold);
        WasErrors := WasErrors or ((RPTbl[I]-QTbl[I])/SigmaTbl[I]>SigmaThreshold);
        Inc(I);
    end;
    
    //
    // 2-sample test
    //
    if  not Silent then
    begin
        Write(Format(''#13#10''#13#10'Testing 2-sample test for 1-type errors'#13#10'',[]));
    end;
    I:=0;
    while I<=QCnt-1 do
    begin
        PTbl[I] := 0;
        LPTbl[I] := 0;
        RPTbl[I] := 0;
        Inc(I);
    end;
    Pass:=1;
    while Pass<=PassCount do
    begin
        N := 5+RandomInteger(20);
        M := 5+RandomInteger(20);
        V := 1+4*RandomReal;
        
        //
        // Both tails
        //
        TestStudentTGenNormRnd(N, 0.0, V, X);
        TestStudentTGenNormRnd(M, 0.0, V, Y);
        StudentTTest2(X, N, Y, M, BT, LT, RT);
        I:=0;
        while I<=QCnt-1 do
        begin
            if BT<=QTbl[I] then
            begin
                PTbl[I] := PTbl[I]+1/PassCount;
            end;
            Inc(I);
        end;
        
        //
        // Left tail
        //
        TestStudentTGenNormRnd(N, 0.5*RandomReal, V, X);
        TestStudentTGenNormRnd(M, 0.0, V, Y);
        StudentTTest2(X, N, Y, M, BT, LT, RT);
        I:=0;
        while I<=QCnt-1 do
        begin
            if LT<=QTbl[I] then
            begin
                LPTbl[I] := LPTbl[I]+1/PassCount;
            end;
            Inc(I);
        end;
        
        //
        // Right tail
        //
        TestStudentTGenNormRnd(N, -0.5*RandomReal, V, X);
        TestStudentTGenNormRnd(M, 0.0, V, Y);
        StudentTTest2(X, N, Y, M, BT, LT, RT);
        I:=0;
        while I<=QCnt-1 do
        begin
            if RT<=QTbl[I] then
            begin
                RPTbl[I] := RPTbl[I]+1/PassCount;
            end;
            Inc(I);
        end;
        Inc(Pass);
    end;
    if  not Silent then
    begin
        Write(Format('Expect. Both    Left    Right'#13#10'',[]));
        I:=0;
        while I<=QCnt-1 do
        begin
            Write(Format('%4.1f%%   %4.1f%%   %4.1f%%   %4.1f%%   '#13#10'',[
                QTbl[I]*100,
                PTbl[I]*100,
                LPTbl[I]*100,
                RPTbl[I]*100]));
            Inc(I);
        end;
    end;
    I:=0;
    while I<=QCnt-1 do
    begin
        WasErrors := WasErrors or (AbsReal(PTbl[I]-QTbl[I])/SigmaTbl[I]>SigmaThreshold);
        WasErrors := WasErrors or ((LPTbl[I]-QTbl[I])/SigmaTbl[I]>SigmaThreshold);
        WasErrors := WasErrors or ((RPTbl[I]-QTbl[I])/SigmaTbl[I]>SigmaThreshold);
        Inc(I);
    end;
    
    //
    // Unequal variance test
    //
    if  not Silent then
    begin
        Write(Format(''#13#10''#13#10'Testing unequal variance test for 1-type errors'#13#10'',[]));
    end;
    I:=0;
    while I<=QCnt-1 do
    begin
        PTbl[I] := 0;
        LPTbl[I] := 0;
        RPTbl[I] := 0;
        Inc(I);
    end;
    Pass:=1;
    while Pass<=PassCount do
    begin
        N := 15+RandomInteger(20);
        M := 15+RandomInteger(20);
        
        //
        // Both tails
        //
        TestStudentTGenNormRnd(N, 0.0, 1+4*RandomReal, X);
        TestStudentTGenNormRnd(M, 0.0, 1+4*RandomReal, Y);
        UnequalVarianceTTest(X, N, Y, M, BT, LT, RT);
        I:=0;
        while I<=QCnt-1 do
        begin
            if BT<=QTbl[I] then
            begin
                PTbl[I] := PTbl[I]+1/PassCount;
            end;
            Inc(I);
        end;
        
        //
        // Left tail
        //
        TestStudentTGenNormRnd(N, 0.5*RandomReal, 1+4*RandomReal, X);
        TestStudentTGenNormRnd(M, 0.0, 1+4*RandomReal, Y);
        UnequalVarianceTTest(X, N, Y, M, BT, LT, RT);
        I:=0;
        while I<=QCnt-1 do
        begin
            if LT<=QTbl[I] then
            begin
                LPTbl[I] := LPTbl[I]+1/PassCount;
            end;
            Inc(I);
        end;
        
        //
        // Right tail
        //
        TestStudentTGenNormRnd(N, -0.5*RandomReal, 1+4*RandomReal, X);
        TestStudentTGenNormRnd(M, 0.0, 1+4*RandomReal, Y);
        UnequalVarianceTTest(X, N, Y, M, BT, LT, RT);
        I:=0;
        while I<=QCnt-1 do
        begin
            if RT<=QTbl[I] then
            begin
                RPTbl[I] := RPTbl[I]+1/PassCount;
            end;
            Inc(I);
        end;
        Inc(Pass);
    end;
    if  not Silent then
    begin
        Write(Format('Expect. Both    Left    Right'#13#10'',[]));
        I:=0;
        while I<=QCnt-1 do
        begin
            Write(Format('%4.1f%%   %4.1f%%   %4.1f%%   %4.1f%%   '#13#10'',[
                QTbl[I]*100,
                PTbl[I]*100,
                LPTbl[I]*100,
                RPTbl[I]*100]));
            Inc(I);
        end;
    end;
    I:=0;
    while I<=QCnt-1 do
    begin
        WasErrors := WasErrors or (AbsReal(PTbl[I]-QTbl[I])/SigmaTbl[I]>SigmaThreshold);
        WasErrors := WasErrors or ((LPTbl[I]-QTbl[I])/SigmaTbl[I]>SigmaThreshold);
        WasErrors := WasErrors or ((RPTbl[I]-QTbl[I])/SigmaTbl[I]>SigmaThreshold);
        Inc(I);
    end;
    
    //
    //
    //
    if  not Silent then
    begin
        if WasErrors then
        begin
            Write(Format('TEST FAILED'#13#10'',[]));
        end
        else
        begin
            Write(Format('TEST PASSED'#13#10'',[]));
        end;
    end;
    Result :=  not WasErrors;
end;


procedure TestStudentTGenNormRnd(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 teststudentttestsunit_test_silent():Boolean;
begin
    Result := TestStudentT(True);
end;


(*************************************************************************
Unit test
*************************************************************************)
function teststudentttestsunit_test():Boolean;
begin
    Result := TestStudentT(False);
end;


end.