testregressunit.pas

maths lib with source

unit testregressunit;
interface
uses Math, Ap, Sysutils, descriptivestatistics, gammaf, normaldistr, igammaf, reflections, bidiagonal, qr, lq, blas, rotations, bdsvd, svd, linreg;

function TestLinRegression(Silent : Boolean):Boolean;
function testregressunit_test_silent():Boolean;
function testregressunit_test():Boolean;

implementation

procedure GenerateRandomTask(XL : Double;
     XR : Double;
     RandomX : Boolean;
     YMin : Double;
     YMax : Double;
     SMin : Double;
     SMax : Double;
     N : Integer;
     var XY : TReal2DArray;
     var S : TReal1DArray);forward;
procedure GenerateTask(A : Double;
     B : Double;
     XL : Double;
     XR : Double;
     RandomX : Boolean;
     SMin : Double;
     SMax : Double;
     N : Integer;
     var XY : TReal2DArray;
     var S : TReal1DArray);forward;
procedure FillTaskWithY(A : Double;
     B : Double;
     N : Integer;
     var XY : TReal2DArray;
     var S : TReal1DArray);forward;
function GenerateNormal(Mean : Double; Sigma : Double):Double;forward;
procedure CalculateMV(const X : TReal1DArray;
     N : Integer;
     var Mean : Double;
     var MeanS : Double;
     var StdDev : Double;
     var StdDevS : Double);forward;
procedure UnsetLR(var LR : LinearModel);forward;


function TestLinRegression(Silent : Boolean):Boolean;
var
    SigmaThreshold : Double;
    MaxN : Integer;
    MaxM : Integer;
    PassCount : Integer;
    EstPassCount : Integer;
    N : Integer;
    I : Integer;
    J : Integer;
    K : Integer;
    TmpI : Integer;
    Pass : Integer;
    EPass : Integer;
    M : Integer;
    TaskType : Integer;
    ModelType : Integer;
    M1 : Integer;
    M2 : Integer;
    N1 : Integer;
    N2 : Integer;
    Info : Integer;
    Info2 : Integer;
    XY : TReal2DArray;
    XY2 : TReal2DArray;
    S : TReal1DArray;
    S2 : TReal1DArray;
    W2 : TReal1DArray;
    X : TReal1DArray;
    TA : TReal1DArray;
    TB : TReal1DArray;
    TC : TReal1DArray;
    XY0 : TReal1DArray;
    TmpWeights : TReal1DArray;
    W : LinearModel;
    WT : LinearModel;
    WT2 : LinearModel;
    X1 : TReal1DArray;
    X2 : TReal1DArray;
    RA : TReal1DArray;
    RA2 : TReal1DArray;
    Y1 : Double;
    Y2 : Double;
    RLen : Integer;
    AllSame : Boolean;
    EA : Double;
    EB : Double;
    VarATested : Double;
    VarBTested : Double;
    A : Double;
    B : Double;
    VarA : Double;
    VarB : Double;
    A2 : Double;
    B2 : Double;
    CovAB : Double;
    CorrAB : Double;
    P : Double;
    QCnt : Integer;
    QTbl : TReal1DArray;
    QVals : TReal1DArray;
    QSigma : TReal1DArray;
    AR : LRReport;
    AR2 : LRReport;
    F : Double;
    FP : Double;
    FM : Double;
    V : Double;
    VV : Double;
    CVRMSError : Double;
    CVAvgError : Double;
    CVAvgRelError : Double;
    RMSError : Double;
    AvgError : Double;
    AvgRelError : Double;
    NonDefect : Boolean;
    SinShift : Double;
    TaskLevel : Double;
    NoiseLevel : Double;
    HStep : Double;
    Sigma : Double;
    Mean : Double;
    MeanS : Double;
    StdDev : Double;
    StdDevS : Double;
    SLCErrors : Boolean;
    SLErrors : Boolean;
    GRCovErrors : Boolean;
    GROptErrors : Boolean;
    GREstErrors : Boolean;
    GROtherErrors : Boolean;
    GRConvErrors : Boolean;
    WasErrors : Boolean;
    i_ : Integer;
begin
    
    //
    // Primary settings
    //
    MaxN := 40;
    MaxM := 5;
    PassCount := 3;
    EstPassCount := 1000;
    SigmaThreshold := 7;
    SLErrors := False;
    SLCErrors := False;
    GRCovErrors := False;
    GROptErrors := False;
    GREstErrors := False;
    GROtherErrors := False;
    GRConvErrors := False;
    WasErrors := False;
    
    //
    // Quantiles table setup
    //
    QCnt := 5;
    SetLength(QTbl, QCnt-1+1);
    SetLength(QVals, QCnt-1+1);
    SetLength(QSigma, QCnt-1+1);
    QTbl[0] := 0.5;
    QTbl[1] := 0.25;
    QTbl[2] := 0.10;
    QTbl[3] := 0.05;
    QTbl[4] := 0.025;
    I:=0;
    while I<=QCnt-1 do
    begin
        QSigma[I] := Sqrt(QTbl[I]*(1-QTbl[I])/EstPassCount);
        Inc(I);
    end;
    
    //
    // Other setup
    //
    SetLength(TA, EstPassCount-1+1);
    SetLength(TB, EstPassCount-1+1);
    
    //
    // Test straight line regression
    //
    N:=2;
    while N<=MaxN do
    begin
        
        //
        // Fail/pass test
        //
        GenerateRandomTask(-1, 1, False, -1, 1, 1, 2, N, XY, S);
        LRLineS(XY, S, N, Info, A, B, VarA, VarB, CovAB, CorrAB, P);
        SLCErrors := SLCErrors or (Info<>1);
        GenerateRandomTask(+1, 1, False, -1, 1, 1, 2, N, XY, S);
        LRLineS(XY, S, N, Info, A, B, VarA, VarB, CovAB, CorrAB, P);
        SLCErrors := SLCErrors or (Info<>-3);
        GenerateRandomTask(-1, 1, False, -1, 1, -1, -1, N, XY, S);
        LRLineS(XY, S, N, Info, A, B, VarA, VarB, CovAB, CorrAB, P);
        SLCErrors := SLCErrors or (Info<>-2);
        GenerateRandomTask(-1, 1, False, -1, 1, 2, 1, 2, XY, S);
        LRLineS(XY, S, 1, Info, A, B, VarA, VarB, CovAB, CorrAB, P);
        SLCErrors := SLCErrors or (Info<>-1);
        
        //
        // Multipass tests
        //
        Pass:=1;
        while Pass<=PassCount do
        begin
            
            //
            // Test S variant against non-S variant
            //
            EA := 2*RandomReal-1;
            EB := 2*RandomReal-1;
            GenerateTask(EA, EB, -5*RandomReal, +5*RandomReal, RandomReal>0.5, 1, 1, N, XY, S);
            LRLineS(XY, S, N, Info, A, B, VarA, VarB, CovAB, CorrAB, P);
            LRLine(XY, N, Info2, A2, B2);
            if (Info<>1) or (Info2<>1) then
            begin
                SLCErrors := True;
            end
            else
            begin
                SLErrors := SLErrors or (AbsReal(A-A2)>100*MachineEpsilon) or (AbsReal(B-B2)>100*MachineEpsilon);
            end;
            
            //
            // Test for A/B
            //
            // Generate task with exact, non-perturbed y[i],
            // then make non-zero s[i]
            //
            EA := 2*RandomReal-1;
            EB := 2*RandomReal-1;
            GenerateTask(EA, EB, -5*RandomReal, +5*RandomReal, N>4, 0.0, 0.0, N, XY, S);
            I:=0;
            while I<=N-1 do
            begin
                S[I] := 1+RandomReal;
                Inc(I);
            end;
            LRLineS(XY, S, N, Info, A, B, VarA, VarB, CovAB, CorrAB, P);
            if Info<>1 then
            begin
                SLCErrors := True;
            end
            else
            begin
                SLErrors := SLErrors or (AbsReal(A-EA)>0.001) or (AbsReal(B-EB)>0.001);
            end;
            
            //
            // Test for VarA, VarB, P (P is being tested only for N>2)
            //
            I:=0;
            while I<=QCnt-1 do
            begin
                QVals[I] := 0;
                Inc(I);
            end;
            EA := 2*RandomReal-1;
            EB := 2*RandomReal-1;
            GenerateTask(EA, EB, -5*RandomReal, +5*RandomReal, N>4, 1.0, 2.0, N, XY, S);
            LRLineS(XY, S, N, Info, A, B, VarA, VarB, CovAB, CorrAB, P);
            if Info<>1 then
            begin
                SLCErrors := True;
                Inc(Pass);
                Continue;
            end;
            VarATested := VarA;
            VarBTested := VarB;
            EPass:=0;
            while EPass<=EstPassCount-1 do
            begin
                
                //
                // Generate
                //
                FillTaskWithY(EA, EB, N, XY, S);
                LRLineS(XY, S, N, Info, A, B, VarA, VarB, CovAB, CorrAB, P);
                if Info<>1 then
                begin
                    SLCErrors := True;
                    Inc(EPass);
                    Continue;
                end;
                
                //
                // A, B, P
                // (P is being tested for uniformity, additional p-tests are below)
                //
                TA[EPass] := A;
                TB[EPass] := B;
                I:=0;
                while I<=QCnt-1 do
                begin
                    if P<=QTbl[I] then
                    begin
                        QVals[I] := QVals[I]+1/EstPassCount;
                    end;
                    Inc(I);
                end;
                Inc(EPass);
            end;
            CalculateMV(TA, EstPassCount, Mean, MeanS, StdDev, StdDevS);
            SLErrors := SLErrors or (AbsReal(Mean-EA)/MeanS>=SigmaThreshold);
            SLErrors := SLErrors or (AbsReal(StdDev-Sqrt(VarATested))/StdDevS>=SigmaThreshold);
            CalculateMV(TB, EstPassCount, Mean, MeanS, StdDev, StdDevS);
            SLErrors := SLErrors or (AbsReal(Mean-EB)/MeanS>=SigmaThreshold);
            SLErrors := SLErrors or (AbsReal(StdDev-Sqrt(VarBTested))/StdDevS>=SigmaThreshold);
            if N>2 then
            begin
                I:=0;
                while I<=QCnt-1 do
                begin
                    if AbsReal(QTbl[I]-QVals[I])/QSigma[I]>SigmaThreshold then
                    begin
                        SLErrors := True;
                    end;
                    Inc(I);
                end;
            end;
            
            //
            // Additional tests for P: correlation with fit quality
            //
            if N>2 then
            begin
                GenerateTask(EA, EB, -5*RandomReal, +5*RandomReal, False, 0.0, 0.0, N, XY, S);
                I:=0;
                while I<=N-1 do
                begin
                    S[I] := 1+RandomReal;
                    Inc(I);
                end;
                LRLineS(XY, S, N, Info, A, B, VarA, VarB, CovAB, CorrAB, P);
                if Info<>1 then
                begin
                    SLCErrors := True;
                    Inc(Pass);
                    Continue;
                end;
                SLErrors := SLErrors or (P<0.999);
                GenerateTask(0, 0, -5*RandomReal, +5*RandomReal, False, 1.0, 1.0, N, XY, S);
                I:=0;
                while I<=N-1 do
                begin
                    if I mod 2=0 then
                    begin
                        XY[I,1] := +5.0;
                    end
                    else
                    begin
                        XY[I,1] := -5.0;
                    end;
                    Inc(I);
                end;
                if N mod 2<>0 then
                begin
                    XY[N-1,1] := 0;
                end;
                LRLineS(XY, S, N, Info, A, B, VarA, VarB, CovAB, CorrAB, P);
                if Info<>1 then
                begin
                    SLCErrors := True;
                    Inc(Pass);
                    Continue;
                end;
                SLErrors := SLErrors or (P>0.001);
            end;
            Inc(Pass);
        end;
        Inc(N);
    end;
    
    //
    // General regression tests:
    //
    
    //
    // Simple linear tests (small sample, optimum point, covariance)
    //
    N:=3;
    while N<=MaxN do
    begin
        SetLength(S, N-1+1);
        
        //
        // Linear tests:
        // a. random points, sigmas
        // b. no sigmas
        //
        SetLength(XY, N-1+1, 1+1);
        I:=0;
        while I<=N-1 do
        begin
            XY[I,0] := 2*RandomReal-1;
            XY[I,1] := 2*RandomReal-1;
            S[I] := 1+RandomReal;
            Inc(I);
        end;
        LRBuildS(XY, S, N, 1, Info, WT, AR);
        if Info<>1 then
        begin
            GRConvErrors := True;
            Inc(N);
            Continue;
        end;
        LRUnpack(WT, TmpWeights, TmpI);
        LRLineS(XY, S, N, Info2, A, B, VarA, VarB, CovAB, CorrAB, P);
        GROptErrors := GROptErrors or (AbsReal(A-TmpWeights[1])>1000000*MachineEpsilon);
        GROptErrors := GROptErrors or (AbsReal(B-TmpWeights[0])>100000*MachineEpsilon);
        GRCovErrors := GRCovErrors or (AbsReal(VarA-AR.C[1,1])>100000*MachineEpsilon);
        GRCovErrors := GRCovErrors or (AbsReal(VarB-AR.C[0,0])>100000*MachineEpsilon);
        GRCovErrors := GRCovErrors or (AbsReal(CovAB-AR.C[1,0])>100000*MachineEpsilon);
        GRCovErrors := GRCovErrors or (AbsReal(CovAB-AR.C[0,1])>100000*MachineEpsilon);
        LRBuild(XY, N, 1, Info, WT, AR);
        if Info<>1 then
        begin
            GRConvErrors := True;
            Inc(N);
            Continue;
        end;
        LRUnpack(WT, TmpWeights, TmpI);
        LRLine(XY, N, Info2, A, B);
        GROptErrors := GROptErrors or (AbsReal(A-TmpWeights[1])>100000*MachineEpsilon);
        GROptErrors := GROptErrors or (AbsReal(B-TmpWeights[0])>100000*MachineEpsilon);
        Inc(N);
    end;
    
    //
    // S covariance versus S-less covariance.
    // Slightly skewed task, large sample size.
    // Will S-less subroutine estimate covariance matrix good enough?
    //
    N := 1000+RandomInteger(3000);
    Sigma := 0.1+RandomReal*1.9;
    SetLength(XY, N-1+1, 1+1);
    SetLength(S, N-1+1);
    I:=0;
    while I<=N-1 do
    begin
        XY[I,0] := 1.5*RandomReal-0.5;
        XY[I,1] := 1.2*XY[I,0]-0.3+GenerateNormal(0, Sigma);
        S[I] := Sigma;
        Inc(I);
    end;
    LRBuild(XY, N, 1, Info, WT, AR);
    LRLineS(XY, S, N, Info2, A, B, VarA, VarB, CovAB, CorrAB, P);
    if (Info<>1) or (Info2<>1) then