regress.pas

delphi矩阵运算、回归分析等数学计算

  var
    B_sim         : TMatrix;
    I, J, ErrCode : Integer;
  begin
    DimMatrix(B_sim, Ubound, MH_SavedSim);
    for I := Lbound to Ubound do
      B[I] := B_min[I] + RanMar * (B_max[I] - B_min[I]);

    { Initialize variance-covariance matrix }
    for I := Lbound to Ubound do
      for J := Lbound to Ubound do
        if I = J then
          { The parameter range is assumed to cover 6 SD's }
          V[I,J] := Sqr((B_max[I] - B_min[I]) / 6.0)
        else
          V[I,J] := 0.0;

    ErrCode := Hastings(ObjFunc, 2.0, B, V, Lbound, Ubound, B_sim, B, F_min);

    if (ErrCode = 0) and (MHFile <> '') then
      SaveSim(B_sim, Lbound, Ubound);

    Simulate := ErrCode;
  end;

  function GenNLFit(RegFunc                    : TRegFunc;
                    DerivProc                  : TDerivProc;
                    Mode                       : TRegMode;
                    X, Y, W                    : TVector;
                    N, Lbound, Ubound, MaxIter : Integer;
                    Tol                        : Float;
                    B, B_min, B_max            : TVector;
                    V                          : TMatrix) : Integer;
  { --------------------------------------------------------------------
    General nonlinear regression routine
    -------------------------------------------------------------------- }
  var
    F_min    : Float;      { Value of objective function at minimum }
    ErrCode  : Integer;    { Error code }
    G        : TVector;    { Gradient vector }
    H        : TMatrix;    { Hessian matrix }
    ObjFunc  : TFuncNVar;  { Objective function }
    GradProc : TGradient;  { Procedure to compute gradient }
    HessProc : THessGrad;  { Procedure to compute gradient and hessian }
  begin
    SetGlobalVar(RegFunc, DerivProc, Mode, X, Y, W,
                 N, Lbound, Ubound, B_min, B_max);

    ObjFunc  := OLS_ObjFunc;
    GradProc := OLS_Gradient;
    HessProc := OLS_HessGrad;

    if Mode = WEIGHTED then
      begin
        ObjFunc  := WLS_ObjFunc;
        GradProc := WLS_Gradient;
        HessProc := WLS_HessGrad;
      end;

    ErrCode := 0;
    case OptAlgo of
      NL_MARQ : ErrCode := Marquardt(ObjFunc, HessProc, B, Lbound, Ubound,
                                     MaxIter, Tol, F_min, V);
      NL_SIMP : ErrCode := Simplex(ObjFunc, B, Lbound, Ubound,
                                   MaxIter, Tol, F_min);
      NL_BFGS : ErrCode := BFGS(ObjFunc, GradProc, B, Lbound, Ubound,
                                MaxIter, Tol, F_min, V);
      NL_SA   : ErrCode := SimAnn(ObjFunc, B, B_min, B_max, Lbound, Ubound,
                                  MaxIter, Tol, F_min);
      NL_MH   : ErrCode := Simulate(ObjFunc, B, B_min, B_max, Lbound, Ubound,
                                    F_min, V);
    end;

    if (OptAlgo in [NL_SIMP, NL_SA]) and (ErrCode = OPT_OK) then
      begin
        { Compute the Hessian matrix and its inverse }
        DimVector(G, Ubound);
        DimMatrix(H, Ubound, Ubound);
        case Mode of
          UNWEIGHTED : OLS_HessGrad(ObjFunc, B, Lbound, Ubound, G, H);
          WEIGHTED   : WLS_HessGrad(ObjFunc, B, Lbound, Ubound, G, H);
        end;
        if InvMat(H, Lbound, Ubound, V) = 0 then
          ErrCode := OPT_OK
        else
          ErrCode := OPT_SING;
      end;

    GenNLFit := ErrCode;
  end;

  function NLFit(RegFunc                    : TRegFunc;
                 DerivProc                  : TDerivProc;
                 X, Y                       : TVector;
                 N, Lbound, Ubound, MaxIter : Integer;
                 Tol                        : Float;
                 B, B_min, B_max            : TVector;
                 V                          : TMatrix) : Integer;
  var
    W : TVector;
  begin
    W := nil;
    NLFit := GenNLFit(RegFunc, DerivProc, UNWEIGHTED, X, Y, W, N,
                      Lbound, Ubound, MaxIter, Tol, B, B_min, B_max, V);
  end;

  function WNLFit(RegFunc : TRegFunc; DerivProc : TDerivProc;
                  X, Y, W : TVector; N, Lbound, Ubound, MaxIter : Integer;
                  Tol : Float; B, B_min, B_max : TVector; V : TMatrix) : Integer;
  begin
    WNLFit := GenNLFit(RegFunc, DerivProc, WEIGHTED, X, Y, W, N,
                       Lbound, Ubound, MaxIter, Tol, B, B_min, B_max, V);
  end;

  procedure GenRegTest(Mode              : TRegMode;
                       Y, Ycalc, W       : TVector;
                       N, Lbound, Ubound : Integer;
                       V                 : TMatrix;
                       var Test          : TRegTest);
  var
    Ybar     : Float;    { Average Y value }
    SSt      : Float;    { Total sum of squares }
    SSe      : Float;    { Explained sum of squares }
    SSr      : Float;    { Residual sum of squares }
    Nobs     : Integer;  { Number of observations }
    Npar     : Integer;  { Number of fitted parameters }
    Nu1, Nu2 : Integer;  { Degrees of freedom }
    I, J     : Integer;  { Loop variables }
  begin
    Nobs := N - FirstPoint + 1;
    Npar := Ubound - Lbound + 1;
    with Test do
      if Nobs > Npar then
        begin
          Ybar := Average(Y, FirstPoint, N);
          if Mode = UNWEIGHTED then
            begin
              SSt := SumSqrDif(Y, FirstPoint, N, Ybar);
              SSe := SumSqrDif(Ycalc, FirstPoint, N, Ybar);
              SSr := SumSqrDifVect(Y, Ycalc, FirstPoint, N);
            end
          else
            begin
              SSt := SumWSqrDif(Y, W, FirstPoint, N, Ybar);
              SSe := SumWSqrDif(Ycalc, W, FirstPoint, N, Ybar);
              SSr := SumWSqrDifVect(Y, Ycalc, W, FirstPoint, N);
            end;
          Nu1 := Npar - 1;
          Nu2 := Nobs - Npar;
          R2 := SSe / SSt;
          R2a := 1.0 - (1.0 - R2) * (Nobs - 1) / Nu2;
          Vr := SSr / Nu2;
          if Vr > 0.0 then
            begin
              F := (SSe / Nu1) / Vr;
              Prob := PSnedecor(Nu1, Nu2, F);
            end
          else
            begin
              F := MAXNUM;
              Prob := 0.0;
            end;
        end
      else
        begin
          Vr := 0.0;
          R2 := 1.0;
          R2a := 0.0;
          F := 0.0;
          Prob := 1.0;
        end;

    { Compute variance-covariance matrix }
    for I := Lbound to Ubound do
      for J := I to Ubound do
        V[I,J] := V[I,J] * Test.Vr;
    for I := Succ(Lbound) to Ubound do
      for J := Lbound to Pred(I) do
        V[I,J] := V[J,I];
  end;

  procedure RegTest(Y, Ycalc          : TVector;
                    N, Lbound, Ubound : Integer;
                    V                 : TMatrix;
                    var Test          : TRegTest);
  var
    W : TVector;
  begin
    W := nil;
    GenRegTest(UNWEIGHTED, Y, Ycalc, W, N, Lbound, Ubound, V, Test);
  end;

  procedure WRegTest(Y, Ycalc, W       : TVector;
                     N, Lbound, Ubound : Integer;
                     V                 : TMatrix;
                     var Test          : TRegTest);
  begin
    GenRegTest(WEIGHTED, Y, Ycalc, W, N, Lbound, Ubound, V, Test);
  end;

  procedure ParamTest(B                 : TVector;
                      V                 : TMatrix;
                      N, Lbound, Ubound : Integer;
                      S, T, Prob        : TVector);
  var
    I    : Integer;
    Nu   : Integer;  { Degrees of freedom }
    Nobs : Integer;  { Number of observations }
    Nvar : Integer;  { Number of indep. variables }
  begin
    Nobs := N - FirstPoint + 1;
    Nvar := Ubound - Lbound + 1;
    Nu := Nobs - Nvar;            { DoF = Nb points - Nb parameters }
    for I := Lbound to Ubound do
      if V[I,I] > 0.0 then
        begin
          S[I] := Sqrt(V[I,I]);
          T[I] := B[I] / S[I];
          Prob[I] := PStudent(Nu, T[I]);
        end
      else
        begin
          S[I] := 0.0;
          T[I] := 0.0;
          Prob[I] := 1.0;
        end;
  end;

  procedure VecMean(X                 : TMatrix;
                    N, Lbound, Ubound : Integer;
                    M                 : TVector);
  var
    I, K, Nobs : Integer;
    Sum        : Float;
  begin
    Nobs := N - FirstPoint + 1;
    for I := Lbound to Ubound do
      begin
        Sum := 0.0;
        for K := FirstPoint to N do
          Sum := Sum + X[I,K];
        M[I] := Sum / Nobs;
      end;
  end;

  procedure VecSD(X                 : TMatrix;
                  N, Lbound, Ubound : Integer;
                  M, S              : TVector);
  var
    I, K, Nobs : Integer;
    Sum        : Float;
  begin
    Nobs := N - FirstPoint + 1;
    for I := Lbound to Ubound do
      begin
        Sum := 0.0;
        for K := FirstPoint to N do
          Sum := Sum + Sqr(X[I,K] - M[I]);
        S[I] := Sqrt(Sum / Nobs);
      end;
  end;

  procedure MatVarCov(X                 : TMatrix;
                      N, Lbound, Ubound : Integer;
                      M                 : TVector;
                      V                 : TMatrix);
  var
    I, J, K, Nobs : Integer;
    Sum           : Float;
  begin
    Nobs := N - FirstPoint + 1;
    for I := Lbound to Ubound do
      for J := I to Ubound do
      begin
        Sum := 0.0;
        for K := FirstPoint to N do
          Sum := Sum + (X[I,K] - M[I]) * (X[J,K] - M[J]);
        V[I,J] := Sum / Nobs;
      end;
    for I := Succ(Lbound) to Ubound do
      for J := Lbound to Pred(I) do
        V[I,J] := V[J,I];
  end;

  procedure MatCorrel(V              : TMatrix;
                      Lbound, Ubound : Integer;
                      R              : TMatrix);
  var
    I, J : Integer;
  begin
    for I := Lbound to Ubound do
      begin
        R[I,I] := 1.0;
        for J := Succ(I) to Ubound do
          begin
            R[I,J] := V[I,J] / Sqrt(V[I,I] * V[J,J]);
            R[J,I] := R[I,J];
          end;
      end;
  end;

  function PCA(R              : TMatrix;
               Lbound, Ubound : Integer;
               Lambda         : TVector;
               C, Rc          : TMatrix) : Integer;
  var
    I, J, ErrCode : Integer;
    Rac           : Float;
  begin
    { Compute eigenvalues and eigenvectors of correlation matrix }
    ErrCode := Jacobi(R, Lbound, Ubound, PCA_MAXITER, PCA_TOL, C, Lambda);

    if ErrCode <> 0 then
      begin
        PCA := ErrCode;
        Exit;
      end;

    { Compute correlations between principal factors and reduced variables }
    for I := Lbound to Ubound do
      begin
        Rac := Sqrt(Lambda[I]);
        for J := Lbound to Ubound do
          Rc[I,J] := C[I,J] * Rac;
      end;

    PCA := ErrCode;
  end;

  procedure ScaleVar(X                 : TMatrix;
                     N, Lbound, Ubound : Integer;
                     M, S              : TVector;
                     Z                 : TMatrix);
  var
    I, K : Integer;
  begin
    for I := Lbound to Ubound do
      for K := FirstPoint to N do
        Z[I,K] := (X[I,K] - M[I]) / S[I];
  end;

  procedure PrinFac(Z                 : TMatrix;
                    N, Lbound, Ubound : Integer;
                    C, F              : TMatrix);
  var
    I, J, K : Integer;
  begin
    for I := Lbound to Ubound do
      for K := FirstPoint to N do
        begin
          F[I,K] := 0.0;
          for J := Lbound to Ubound do
            F[I,K] := F[I,K] + C[I,J] * Z[J,K];
        end;
  end;

begin
  MHFile     := '';
  RegAlgo    := SVD;
  OptAlgo    := NL_MARQ;
  FirstPoint := 1;
  NN         := 1;
  XX         := nil;
  YY         := nil;
  WW         := nil;
  YYcalc     := nil;
  FirstParam := 0;
  LastParam  := 1;
  ParamMin   := nil;
  ParamMax   := nil;
end.