testregressunit.pas

maths lib with source

    begin
        GRConvErrors := True;
    end
    else
    begin
        GRCovErrors := GRCovErrors or (AbsReal(Ln(AR.C[0,0]/VarB))>Ln(1.2));
        GRCovErrors := GRCovErrors or (AbsReal(Ln(AR.C[1,1]/VarA))>Ln(1.2));
        GRCovErrors := GRCovErrors or (AbsReal(Ln(AR.C[0,1]/CovAB))>Ln(1.2));
        GRCovErrors := GRCovErrors or (AbsReal(Ln(AR.C[1,0]/CovAB))>Ln(1.2));
    end;
    
    //
    // General tests:
    // * basis functions - up to cubic
    // * task types:
    // * data set is noisy sine half-period with random shift
    // * tests:
    //   unpacking/packing
    //   optimality
    //   error estimates
    // * tasks:
    //   0 = noised sine
    //   1 = degenerate task with 1-of-n encoded categorical variables
    //   2 = random task with large variation (for 1-type models)
    //   3 = random task with small variation (for 1-type models)
    //
    //   Additional tasks TODO
    //   specially designed task with defective vectors which leads to
    //   the failure of the fast CV formula.
    //
    //
    ModelType:=0;
    while ModelType<=1 do
    begin
        TaskType:=0;
        while TaskType<=3 do
        begin
            if TaskType=0 then
            begin
                M1 := 1;
                M2 := 3;
            end;
            if TaskType=1 then
            begin
                M1 := 9;
                M2 := 9;
            end;
            if (TaskType=2) or (TaskType=3) then
            begin
                M1 := 9;
                M2 := 9;
            end;
            M:=M1;
            while M<=M2 do
            begin
                if TaskType=0 then
                begin
                    N1 := M+3;
                    N2 := M+20;
                end;
                if TaskType=1 then
                begin
                    N1 := 70+RandomInteger(70);
                    N2 := N1;
                end;
                if (TaskType=2) or (TaskType=3) then
                begin
                    N1 := 100;
                    N2 := N1;
                end;
                N:=N1;
                while N<=N2 do
                begin
                    SetLength(XY, N-1+1, M+1);
                    SetLength(XY0, N-1+1);
                    SetLength(S, N-1+1);
                    HStep := 0.001;
                    NoiseLevel := 0.2;
                    
                    //
                    // Prepare task
                    //
                    if TaskType=0 then
                    begin
                        I:=0;
                        while I<=N-1 do
                        begin
                            XY[I,0] := 2*RandomReal-1;
                            Inc(I);
                        end;
                        I:=0;
                        while I<=N-1 do
                        begin
                            J:=1;
                            while J<=M-1 do
                            begin
                                XY[I,J] := XY[I,0]*XY[I,J-1];
                                Inc(J);
                            end;
                            Inc(I);
                        end;
                        SinShift := RandomReal*Pi;
                        I:=0;
                        while I<=N-1 do
                        begin
                            XY0[I] := Sin(SinShift+Pi*0.5*(XY[I,0]+1));
                            XY[I,M] := XY0[I]+NoiseLevel*GenerateNormal(0, 1);
                            Inc(I);
                        end;
                    end;
                    if TaskType=1 then
                    begin
                        Assert(M=9);
                        SetLength(TA, 8+1);
                        TA[0] := 1;
                        TA[1] := 2;
                        TA[2] := 3;
                        TA[3] := 0.25;
                        TA[4] := 0.5;
                        TA[5] := 0.75;
                        TA[6] := 0.06;
                        TA[7] := 0.12;
                        TA[8] := 0.18;
                        I:=0;
                        while I<=N-1 do
                        begin
                            J:=0;
                            while J<=M-1 do
                            begin
                                XY[I,J] := 0;
                                Inc(J);
                            end;
                            XY[I,I mod 3] := 1;
                            XY[I,I div 3 mod 3] := 1;
                            XY[I,I div 9 mod 3] := 1;
                            V := APVDotProduct(@XY[I][0], 0, 8, @TA[0], 0, 8);
                            XY0[I] := V;
                            XY[I,M] := V+NoiseLevel*GenerateNormal(0, 1);
                            Inc(I);
                        end;
                    end;
                    if (TaskType=2) or (TaskType=3) then
                    begin
                        Assert(M=9);
                        SetLength(TA, 8+1);
                        TA[0] := 1;
                        TA[1] := -2;
                        TA[2] := 3;
                        TA[3] := 0.25;
                        TA[4] := -0.5;
                        TA[5] := 0.75;
                        TA[6] := -0.06;
                        TA[7] := 0.12;
                        TA[8] := -0.18;
                        I:=0;
                        while I<=N-1 do
                        begin
                            J:=0;
                            while J<=M-1 do
                            begin
                                if TaskType=2 then
                                begin
                                    XY[I,J] := 1+GenerateNormal(0, 3);
                                end
                                else
                                begin
                                    XY[I,J] := 1+GenerateNormal(0, 0.05);
                                end;
                                Inc(J);
                            end;
                            V := APVDotProduct(@XY[I][0], 0, 8, @TA[0], 0, 8);
                            XY0[I] := V;
                            XY[I,M] := V+NoiseLevel*GenerateNormal(0, 1);
                            Inc(I);
                        end;
                    end;
                    I:=0;
                    while I<=N-1 do
                    begin
                        S[I] := 1+RandomReal;
                        Inc(I);
                    end;
                    
                    //
                    // Solve (using S-variant, non-S-variant is not tested)
                    //
                    if ModelType=0 then
                    begin
                        LRBuildS(XY, S, N, M, Info, WT, AR);
                    end
                    else
                    begin
                        LRBuildZS(XY, S, N, M, Info, WT, AR);
                    end;
                    if Info<>1 then
                    begin
                        GRConvErrors := True;
                        Inc(N);
                        Continue;
                    end;
                    LRUnpack(WT, TmpWeights, TmpI);
                    
                    //
                    // LRProcess test
                    //
                    SetLength(X, M-1+1);
                    V := TmpWeights[M];
                    I:=0;
                    while I<=M-1 do
                    begin
                        X[I] := 2*RandomReal-1;
                        V := V+TmpWeights[I]*X[I];
                        Inc(I);
                    end;
                    GROtherErrors := GROtherErrors or (AbsReal(V-LRProcess(WT, X))/AbsReal(V)>1000*MachineEpsilon);
                    
                    //
                    // LRPack test
                    //
                    LRPack(TmpWeights, M, WT2);
                    SetLength(X, M-1+1);
                    I:=0;
                    while I<=M-1 do
                    begin
                        X[I] := 2*RandomReal-1;
                        Inc(I);
                    end;
                    V := LRProcess(WT, X);
                    GROtherErrors := GROtherErrors or (AbsReal(V-LRProcess(WT2, X))/AbsReal(V)>1000*MachineEpsilon);
                    
                    //
                    // Optimality test
                    //
                    K:=0;
                    while K<=M do
                    begin
                        if (ModelType=1) and (K=M) then
                        begin
                            
                            //
                            // 0-type models (with non-zero constant term)
                            // are tested for optimality of all coefficients.
                            //
                            // 1-type models (with zero constant term)
                            // are tested for optimality of non-constant terms only.
                            //
                            Inc(K);
                            Continue;
                        end;
                        F := 0;
                        FP := 0;
                        FM := 0;
                        I:=0;
                        while I<=N-1 do
                        begin
                            V := TmpWeights[M];
                            J:=0;
                            while J<=M-1 do
                            begin
                                V := V+XY[I,J]*TmpWeights[J];
                                Inc(J);
                            end;
                            F := F+Sqr((V-XY[I,M])/S[I]);
                            if K<M then
                            begin
                                VV := XY[I,K];
                            end
                            else
                            begin
                                VV := 1;
                            end;
                            FP := FP+Sqr((V+VV*HStep-XY[I,M])/S[I]);
                            FM := FM+Sqr((V-VV*HStep-XY[I,M])/S[I]);
                            Inc(I);
                        end;
                        GROptErrors := GROptErrors or (F>FP) or (F>FM);
                        Inc(K);
                    end;
                    
                    //
                    // Covariance matrix test:
                    // generate random vector, project coefficients on it,
                    // compare variance of projection with estimate provided
                    // by cov.matrix
                    //
                    SetLength(TA, EstPassCount-1+1);
                    SetLength(TB, M+1);
                    SetLength(TC, M+1);
                    SetLength(XY2, N-1+1, M+1);
                    I:=0;
                    while I<=M do
                    begin
                        TB[I] := GenerateNormal(0, 1);
                        Inc(I);
                    end;
                    EPass:=0;
                    while EPass<=EstPassCount-1 do
                    begin
                        I:=0;
                        while I<=N-1 do
                        begin
                            APVMove(@XY2[I][0], 0, M-1, @XY[I][0], 0, M-1);
                            XY2[I,M] := XY0[I]+S[I]*GenerateNormal(0, 1);
                            Inc(I);
                        end;
                        if ModelType=0 then
                        begin
                            LRBuildS(XY2, S, N, M, Info, WT, AR2);
                        end
                        else
                        begin
                            LRBuildZS(XY2, S, N, M, Info, WT, AR2);
                        end;
                        if Info<>1 then
                        begin
                            TA[EPass] := 0;
                            GRConvErrors := True;
                            Exit;
                        end;
                        LRUnpack(WT, W2, TmpI);
                        V := APVDotProduct(@TB[0], 0, M, @W2[0], 0, M);
                        TA[EPass] := V;
                        Inc(EPass);
                    end;
                    CalculateMV(TA, EstPassCount, Mean, MeanS, StdDev, StdDevS);
                    I:=0;
                    while I<=M do
                    begin
                        V := 0.0;
                        for i_ := 0 to M do
                        begin
                            V := V + TB[i_]*AR.C[i_,I];
                        end;
                        TC[I] := V;
                        Inc(I);
                    end;
                    V := APVDotProduct(@TC[0], 0, M, @TB[0], 0, M);
                    GRCovErrors := GRCovErrors or (AbsReal((Sqrt(V)-StdDev)/StdDevS)>=SigmaThreshold);
                    
                    //
                    // Test for the fast CV error:
                    // calculate CV error by definition (leaving out N
                    // points and recalculating solution).
                    //
                    // Test for the training set error
                    //
                    CVRMSError := 0;
                    CVAvgError := 0;
                    CVAvgRelError := 0;
                    RMSError := 0;
                    AvgError := 0;
                    AvgRelError := 0;
                    SetLength(XY2, N-2+1, M+1);
                    SetLength(S2, N-2+1);
                    I:=0;
                    while I<=N-2 do
                    begin
                        APVMove(@XY2[I][0], 0, M, @XY[I+1][0], 0, M);
                        S2[I] := S[I+1];
                        Inc(I);
                    end;
                    I:=0;
                    while I<=N-1 do
                    begin
                        
                        //
                        // Trn
                        //
                        V := APVDotProduct(@XY[I][0], 0, M-1, @TmpWeights[0], 0, M-1);
                        V := V+TmpWeights[M];
                        RMSError := RMSError+Sqr(V-XY[I,M]);
                        AvgError := AvgError+AbsReal(V-XY[I,M]);
                        AvgRelError := AvgRelError+AbsReal((V-XY[I,M])/XY[I,M]);
                        
                        //
                        // CV: non-defect vectors only
                        //
                        NonDefect := True;
                        K:=0;
                        while K<=AR.NCVDefects-1 do
                        begin
                            if AR.CVDefects[K]=I then
                            begin
                                NonDefect := False;
                            end;
                            Inc(K);
                        end;
                        if NonDefect then
                        begin
                            if ModelType=0 then
                            begin
                                LRBuildS(XY2, S2, N-1, M, Info2, WT, AR2);
                            end
                            else
                            begin
                                LRBuildZS(XY2, S2, N-1, M, Info2, WT, AR2);
                            end;
                            if Info2<>1 then
                            begin
                                GRConvErrors := True;
                                Inc(I);
                                Continue;
                            end;
                            LRUnpack(WT, W2, TmpI);
                            V := APVDotProduct(@XY[I][0], 0, M-1, @W2[0], 0, M-1);
                            V := V+W2[M];
                            CVRMSError := CVRMSError+Sqr(V-XY[I,M]);
                            CVAvgError := CVAvgError+AbsReal(V-XY[I,M]);
                            CVAvgRelError := CVAvgRelError+AbsReal((V-XY[I,M])/XY[I,M]);
                        end;
                        
                        //
                        // Next set
                        //
                        if I<>N-1 then
                        begin
                            APVMove(@XY2[I][0], 0, M, @XY[I][0], 0, M);
                            S2[I] := S[I];
                        end;
                        Inc(I);
                    end;
                    CVRMSError := Sqrt(CVRMSError/(N-AR.NCVDefects));
                    CVAvgError := CVAvgError/(N-AR.NCVDefects);
                    CVAvgRelError := CVAvgRelError/(N-AR.NCVDefects);
                    RMSError := Sqrt(RMSError/N);
                    AvgError := AvgError/N;
                    AvgRelError := AvgRelError/N;
                    GREstErrors := GREstErrors or (AbsReal(Ln(AR.CVRMSError/CVRMSError))>Ln(1+1.0E-5));
                    GREstErrors := GREstErrors or (AbsReal(Ln(AR.CVAvgError/CVAvgError))>Ln(1+1.0E-5));
                    GREstErrors := GREstErrors or (AbsReal(Ln(AR.CVAvgRelError/CVAvgRelError))>Ln(1+1.0E-5));
                    GREstErrors := GREstErrors or (AbsReal(Ln(AR.RMSError/RMSError))>Ln(1+1.0E-5));
                    GREstErrors := GREstErrors or (AbsReal(Ln(AR.AvgError/AvgError))>Ln(1+1.0E-5));
                    GREstErrors := GREstErrors or (AbsReal(Ln(AR.AvgRelError/AvgRelError))>Ln(1+1.0E-5));
                    Inc(N);
                end;
                Inc(M);
            end;
            Inc(TaskType);
        end;
        Inc(ModelType);
    end;
    
    //
    // Additional subroutines
    //
    Pass:=1;
    while Pass<=50 do
    begin
        N := 2;
        NoiseLevel := RandomReal+0.1;
        TaskLevel := 2*RandomReal-1;
        SetLength(XY, 3*N-1+1, 1+1);
        I:=0;
        while I<=N-1 do
        begin