nonlinear regression routines.txt

VB写的非线性回归源程序

  General nonlinear regression routine
    -------------------------------------------------------------------- }
  var
    F_min : Float;         { Value of objective function at minimum }
    ErrCode : Integer;     { Error code }
    G : PVector;           { Gradient vector }
    H : PMatrix;           { 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;
        DelVector(G, Ubound);
        DelMatrix(H, Ubound, Ubound);
      end;

    GenNLFit := ErrCode;
  end;

  function NLFit(RegFunc : TRegFunc; DerivProc : TDerivProc;
                 X, Y : PVector; N, Lbound, Ubound, MaxIter : Integer;
                 Tol : Float; B, B_min, B_max : PVector; V : PMatrix) : Integer;
  var
    W : PVector;
  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 : PVector; N, Lbound, Ubound, MaxIter : Integer;
                  Tol : Float; B, B_min, B_max : PVector; V : PMatrix) : 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 : PVector;
                       N, Lbound, Ubound : Integer; V : PMatrix;
                       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 : PVector; N, Lbound, Ubound : Integer;
                    V : PMatrix; var Test : TRegTest);
  var
    W : PVector;
  begin
    W := nil;
    GenRegTest(UNWEIGHTED, Y, Ycalc, W, N, Lbound, Ubound, V, Test);
  end;

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

  procedure ParamTest(B : PVector; V : PMatrix; N, Lbound, Ubound : Integer;
                      S, T, Prob : PVector);
  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 : PMatrix; N, Lbound, Ubound : Integer; M : PVector);
  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 : PMatrix; N, Lbound, Ubound : Integer; M, S : PVector);
  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 : PMatrix; N, Lbound, Ubound : Integer; M : PVector; V : PMatrix);
  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 : PMatrix; Lbound, Ubound : Integer; R : PMatrix);
  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 : PMatrix; Lbound, Ubound : Integer;
               Lambda : PVector; C, Rc : PMatrix) : 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 : PMatrix; N, Lbound, Ubound : Integer;
                     M, S : PVector; Z : PMatrix);
  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 : PMatrix; N, Lbound, Ubound : Integer;
                    C, F : PMatrix);
  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.