testregressunit.pas

maths lib with source

            XY[3*I+0,0] := I;
            XY[3*I+1,0] := I;
            XY[3*I+2,0] := I;
            XY[3*I+0,1] := TaskLevel-NoiseLevel;
            XY[3*I+1,1] := TaskLevel;
            XY[3*I+2,1] := TaskLevel+NoiseLevel;
            Inc(I);
        end;
        LRBuild(XY, 3*N, 1, Info, WT, AR);
        if Info=1 then
        begin
            LRUnpack(WT, TmpWeights, TmpI);
            V := LRRMSError(WT, XY, 3*N);
            GROtherErrors := GROtherErrors or (AbsReal(V-NoiseLevel*Sqrt(2/3))>1000*MachineEpsilon);
            V := LRAvgError(WT, XY, 3*N);
            GROtherErrors := GROtherErrors or (AbsReal(V-NoiseLevel*(2/3))>1000*MachineEpsilon);
            V := LRAvgRelError(WT, XY, 3*N);
            GROtherErrors := GROtherErrors or (AbsReal(V-(AbsReal(NoiseLevel/(TaskLevel-NoiseLevel))+AbsReal(NoiseLevel/(TaskLevel+NoiseLevel)))/3)>1000*MachineEpsilon);
        end
        else
        begin
            GROtherErrors := True;
        end;
        I:=0;
        while I<=N-1 do
        begin
            XY[3*I+0,0] := I;
            XY[3*I+1,0] := I;
            XY[3*I+2,0] := I;
            XY[3*I+0,1] := -NoiseLevel;
            XY[3*I+1,1] := 0;
            XY[3*I+2,1] := +NoiseLevel;
            Inc(I);
        end;
        LRBuild(XY, 3*N, 1, Info, WT, AR);
        if Info=1 then
        begin
            LRUnpack(WT, TmpWeights, TmpI);
            V := LRAvgRelError(WT, XY, 3*N);
            GROtherErrors := GROtherErrors or (AbsReal(V-1)>1000*MachineEpsilon);
        end
        else
        begin
            GROtherErrors := True;
        end;
        Inc(Pass);
    end;
    Pass:=1;
    while Pass<=10 do
    begin
        M := 1+RandomInteger(5);
        N := 10+RandomInteger(10);
        SetLength(XY, N-1+1, M+1);
        I:=0;
        while I<=N-1 do
        begin
            J:=0;
            while J<=M do
            begin
                XY[I,J] := 2*RandomReal-1;
                Inc(J);
            end;
            Inc(I);
        end;
        LRBuild(XY, N, M, Info, W, AR);
        if Info<0 then
        begin
            GROtherErrors := True;
            Break;
        end;
        SetLength(X1, M-1+1);
        SetLength(X2, M-1+1);
        
        //
        // Same inputs on original leads to same outputs
        // on copy created using LRCopy
        //
        UnsetLR(WT);
        LRCopy(W, WT);
        I:=0;
        while I<=M-1 do
        begin
            X1[I] := 2*RandomReal-1;
            X2[I] := X1[I];
            Inc(I);
        end;
        Y1 := LRProcess(W, X1);
        Y2 := LRProcess(W, X2);
        AllSame := Y1=Y2;
        GROtherErrors := GROtherErrors or  not AllSame;
        
        //
        // Same inputs on original leads to same outputs
        // on copy created using LRSerialize
        //
        UnsetLR(WT);
        SetLength(RA, 0+1);
        RA[0] := 0;
        RLen := 0;
        LRSerialize(W, RA, RLen);
        SetLength(RA2, RLen-1+1);
        I:=0;
        while I<=RLen-1 do
        begin
            RA2[I] := RA[I];
            Inc(I);
        end;
        LRUnserialize(RA2, WT);
        I:=0;
        while I<=M-1 do
        begin
            X1[I] := 2*RandomReal-1;
            X2[I] := X1[I];
            Inc(I);
        end;
        Y1 := LRProcess(W, X1);
        Y2 := LRProcess(WT, X2);
        AllSame := Y1=Y2;
        GROtherErrors := GROtherErrors or  not AllSame;
        Inc(Pass);
    end;
    
    //
    // TODO: Degenerate tests (when design matrix and right part are zero)
    //
    
    //
    // Final report
    //
    WasErrors := SLErrors or SLCErrors or GROptErrors or GRCovErrors or GREstErrors or GROtherErrors or GRConvErrors;
    if  not Silent then
    begin
        Write(Format('REGRESSION TEST'#13#10'',[]));
        Write(Format('STRAIGHT LINE REGRESSION:                ',[]));
        if  not SLErrors then
        begin
            Write(Format('OK'#13#10'',[]));
        end
        else
        begin
            Write(Format('FAILED'#13#10'',[]));
        end;
        Write(Format('STRAIGHT LINE REGRESSION CONVERGENCE:    ',[]));
        if  not SLCErrors then
        begin
            Write(Format('OK'#13#10'',[]));
        end
        else
        begin
            Write(Format('FAILED'#13#10'',[]));
        end;
        Write(Format('GENERAL LINEAR REGRESSION:               ',[]));
        if  not (GROptErrors or GRCovErrors or GREstErrors or GROtherErrors or GRConvErrors) then
        begin
            Write(Format('OK'#13#10'',[]));
        end
        else
        begin
            Write(Format('FAILED'#13#10'',[]));
        end;
        Write(Format('* OPTIMALITY:                            ',[]));
        if  not GROptErrors then
        begin
            Write(Format('OK'#13#10'',[]));
        end
        else
        begin
            Write(Format('FAILED'#13#10'',[]));
        end;
        Write(Format('* COV. MATRIX:                           ',[]));
        if  not GRCovErrors then
        begin
            Write(Format('OK'#13#10'',[]));
        end
        else
        begin
            Write(Format('FAILED'#13#10'',[]));
        end;
        Write(Format('* ERROR ESTIMATES:                       ',[]));
        if  not GREstErrors then
        begin
            Write(Format('OK'#13#10'',[]));
        end
        else
        begin
            Write(Format('FAILED'#13#10'',[]));
        end;
        Write(Format('* CONVERGENCE:                           ',[]));
        if  not GRConvErrors then
        begin
            Write(Format('OK'#13#10'',[]));
        end
        else
        begin
            Write(Format('FAILED'#13#10'',[]));
        end;
        Write(Format('* OTHER SUBROUTINES:                     ',[]));
        if  not GROtherErrors 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;


(*************************************************************************
Task generation. Meaningless task, just random numbers.
*************************************************************************)
procedure GenerateRandomTask(XL : Double;
     XR : Double;
     RandomX : Boolean;
     YMin : Double;
     YMax : Double;
     SMin : Double;
     SMax : Double;
     N : Integer;
     var XY : TReal2DArray;
     var S : TReal1DArray);
var
    I : Integer;
begin
    SetLength(XY, N-1+1, 1+1);
    SetLength(S, N-1+1);
    I:=0;
    while I<=N-1 do
    begin
        if RandomX then
        begin
            XY[I,0] := XL+(XR-XL)*RandomReal;
        end
        else
        begin
            XY[I,0] := XL+(XR-XL)*I/(N-1);
        end;
        XY[I,1] := YMin+(YMax-YMin)*RandomReal;
        S[I] := SMin+(SMax-SMin)*RandomReal;
        Inc(I);
    end;
end;


(*************************************************************************
Task generation.
*************************************************************************)
procedure GenerateTask(A : Double;
     B : Double;
     XL : Double;
     XR : Double;
     RandomX : Boolean;
     SMin : Double;
     SMax : Double;
     N : Integer;
     var XY : TReal2DArray;
     var S : TReal1DArray);
var
    I : Integer;
begin
    SetLength(XY, N-1+1, 1+1);
    SetLength(S, N-1+1);
    I:=0;
    while I<=N-1 do
    begin
        if RandomX then
        begin
            XY[I,0] := XL+(XR-XL)*RandomReal;
        end
        else
        begin
            XY[I,0] := XL+(XR-XL)*I/(N-1);
        end;
        S[I] := SMin+(SMax-SMin)*RandomReal;
        XY[I,1] := A+B*XY[I,0]+GenerateNormal(0, S[I]);
        Inc(I);
    end;
end;


(*************************************************************************
Task generation.
y[i] are filled based on A, B, X[I], S[I]
*************************************************************************)
procedure FillTaskWithY(A : Double;
     B : Double;
     N : Integer;
     var XY : TReal2DArray;
     var S : TReal1DArray);
var
    I : Integer;
begin
    I:=0;
    while I<=N-1 do
    begin
        XY[I,1] := A+B*XY[I,0]+GenerateNormal(0, S[I]);
        Inc(I);
    end;
end;


(*************************************************************************
Normal random numbers
*************************************************************************)
function GenerateNormal(Mean : Double; Sigma : Double):Double;
var
    U : Double;
    V : Double;
    S : Double;
    Sum : Double;
begin
    Result := Mean;
    while True 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);
            Result := Sigma*u*sum+Mean;
            Exit;
        end;
    end;
end;


(*************************************************************************
Moments estimates and their errors
*************************************************************************)
procedure CalculateMV(const X : TReal1DArray;
     N : Integer;
     var Mean : Double;
     var MeanS : Double;
     var StdDev : Double;
     var StdDevS : Double);
var
    I : Integer;
    V : Double;
    V1 : Double;
    V2 : Double;
    Variance : Double;
begin
    Mean := 0;
    MeanS := 1;
    StdDev := 0;
    StdDevS := 1;
    Variance := 0;
    if N<=1 then
    begin
        Exit;
    end;
    
    //
    // Mean
    //
    I:=0;
    while I<=N-1 do
    begin
        Mean := Mean+X[I];
        Inc(I);
    end;
    Mean := Mean/N;
    
    //
    // Variance (using corrected two-pass algorithm)
    //
    if N<>1 then
    begin
        V1 := 0;
        I:=0;
        while I<=N-1 do
        begin
            V1 := V1+Sqr(X[I]-Mean);
            Inc(I);
        end;
        V2 := 0;
        I:=0;
        while I<=N-1 do
        begin
            V2 := V2+(X[I]-Mean);
            Inc(I);
        end;
        V2 := Sqr(V2)/N;
        Variance := (V1-V2)/(N-1);
        if Variance<0 then
        begin
            Variance := 0;
        end;
        StdDev := Sqrt(Variance);
    end;
    
    //
    // Errors
    //
    MeanS := StdDev/Sqrt(N);
    StdDevS := StdDev*Sqrt(2)/Sqrt(N-1);
end;


(*************************************************************************
Unsets LR
*************************************************************************)
procedure UnsetLR(var LR : LinearModel);
var
    XY : TReal2DArray;
    Info : Integer;
    Rep : LRReport;
    I : Integer;
begin
    SetLength(XY, 5+1, 1+1);
    I:=0;
    while I<=5 do
    begin
        XY[I,0] := 0;
        XY[I,1] := 0;
        Inc(I);
    end;
    LRBuild(XY, 6, 1, Info, LR, Rep);
    Assert(Info>0);
end;


(*************************************************************************
Silent unit test
*************************************************************************)
function testregressunit_test_silent():Boolean;
begin
    Result := TestLinRegression(True);
end;


(*************************************************************************
Unit test
*************************************************************************)
function testregressunit_test():Boolean;
begin
    Result := TestLinRegression(False);
end;


end.