regress.pas

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

{ ----------------------------------------------------------------------
  Performs a principal component analysis of the correlation matrix R
  ----------------------------------------------------------------------
  Input  : R, Lbound, Ubound
  Output : Lambda = Eigenvalues of the correlation matrix
                    (in descending order)
           C      = Eigenvectors of the correlation matrix
                    (C[I] is the I-th eigenvector)
           Rc     = Correlations between principal factors and variables
                    (Rc[I,J] is the correlation coefficient between
                     factor I and variable J)
  ----------------------------------------------------------------------
  Possible results : MAT_OK       : No error
                     MAT_NON_CONV : Non-convergence of eigenvalue
                                    determination
  ----------------------------------------------------------------------
  NB : This procedure destroys the original matrix R
  ---------------------------------------------------------------------- }

procedure ScaleVar(X                 : TMatrix;
                   N, Lbound, Ubound : Integer;
                   M, S              : TVector;
                   Z                 : TMatrix);
{ ----------------------------------------------------------------------
  Scales a set of variables by subtracting means and dividing by SD's
  ----------------------------------------------------------------------
  Input  : X, N, Lbound, Ubound, M, S
  Output : Z = matrix of scaled variables (Z[I] contains the I-th var.)
  ---------------------------------------------------------------------- }

procedure PrinFac(Z                 : TMatrix;
                  N, Lbound, Ubound : Integer;
                  C, F              : TMatrix);
{ ----------------------------------------------------------------------
  Computes principal factors
  ----------------------------------------------------------------------
  Input  : Z, N, Lbound, Ubound
           C = matrix of eigenvectors from PCA
  Output : F = matrix of principal factors (F[I] contains the I-th factor)
  ---------------------------------------------------------------------- }

implementation

{ Constants for eigenvalue determination in PCA }
const
  PCA_MAXITER = 100;      { Max number of iterations }
  PCA_TOL     = 1.0E-6;   { Required precision }
  MAX_FUNC    = 1.0E+30;  { Max. value for objective function
                            (used to prevent overflow) }
{ Default settings }
var
  RegAlgo    : TRegAlgo;  { Linear regression algorithm }
  OptAlgo    : TOptAlgo;  { Optimization algorithms }
  FirstPoint : Integer;   { Index of first data point }

{ Global variables used by the nonlinear regression routines }
var
  NN         : Integer;     { Number of observations }
  XX         : TVector;     { X coordinates }
  YY         : TVector;     { Y coordinates }
  WW         : TVector;     { Weights }
  YYcalc     : TVector;     { Estimated Y values }
  FirstParam : Integer;     { Index of first fitted parameter }
  LastParam  : Integer;     { Index of last fitted parameter }
  ParamMin   : TVector;     { Lower bounds on parameters }
  ParamMax   : TVector;     { Higher bounds on parameters }
  RegFunc1   : TRegFunc;    { Regression function }
  DerivProc1 : TDerivProc;  { Derivation procedure }

  function TolSVD(N : Integer) : Float;
  { This function sets the relative threshold below which a singular value
    is considered zero. N is the number of observations. }
  begin
    TolSVD := N * MACHEP;
  end;

  procedure SetRegAlgo(Algo : TRegAlgo);
  begin
    RegAlgo := Algo;
  end;

  procedure SetOptAlgo(Algo : TOptAlgo);
  begin
    OptAlgo := Algo;
  end;

  procedure SetFirstPoint(Index : Integer);
  begin
    if Index >= 0 then
      FirstPoint := Index;
  end;

  function GetRegAlgo : TRegAlgo;
  begin
    GetRegAlgo := RegAlgo;
  end;

  function GetOptAlgo : TOptAlgo;
  begin
    GetOptAlgo := OptAlgo;
  end;

  function GetFirstPoint : Integer;
  begin
    GetFirstPoint := FirstPoint;
  end;

  function GenLinFit(Mode    : TRegMode;
                     X, Y, W : TVector;
                     N       : Integer;
                     B       : TVector;
                     V       : TMatrix) : Integer;
{ ----------------------------------------------------------------------
  General linear regression routine
  ---------------------------------------------------------------------- }
  var
    WX, S, SX, SY, SX2, SXY, D : Float;
    K : Integer;
  begin
    S := 0.0;
    SX := 0.0;
    SY := 0.0;
    SX2 := 0.0;
    SXY := 0.0;
    if Mode = UNWEIGHTED then
      begin
        S := N - FirstPoint + 1;
        for K := FirstPoint to N do
          begin
            SX := SX + X[K];
            SY := SY + Y[K];
            SX2 := SX2 + Sqr(X[K]);
            SXY := SXY + X[K] * Y[K];
          end;
      end
    else
      begin
        for K := FirstPoint to N do
          begin
            WX := W[K] * X[K];
            S := S + W[K];
            SX := SX + WX;
            SY := SY + W[K] * Y[K];
            SX2 := SX2 + WX * X[K];
            SXY := SXY + WX * Y[K];
          end;
      end;
    D := S * SX2 - Sqr(SX);
    if D <= 0.0 then
      GenLinFit := MAT_SINGUL
    else
      begin
        V[0,0] := SX2 / D;
        V[0,1] := - SX / D;
        V[1,0] := V[0,1];
        V[1,1] := S / D;
        B[0] := V[0,0] * SY + V[0,1] * SXY;
        B[1] := V[1,0] * SY + V[1,1] * SXY;
        GenLinFit := MAT_OK;
      end;
  end;

  function LinFit(X, Y : TVector;
                  N    : Integer;
                  B    : TVector;
                  V    : TMatrix) : Integer;
  var
    W : TVector;
  begin
    W := nil;
    LinFit := GenLinFit(UNWEIGHTED, X, Y, W, N, B, V);
  end;

  function WLinFit(X, Y, W : TVector;
                   N       : Integer;
                   B       : TVector;
                   V       : TMatrix) : Integer;
  begin
    WLinFit := GenLinFit(WEIGHTED, X, Y, W, N, B, V);
  end;

  function Gauss_GenMulFit(Mode     : TRegMode;
                           X        : TMatrix;
                           Y, W     : TVector;
                           N, Nvar  : Integer;
                           ConsTerm : Boolean;
                           B        : TVector;
                           V        : TMatrix) : Integer;
{ ----------------------------------------------------------------------
  General multiple linear regression routine (Gauss-Jordan algorithm)
  ---------------------------------------------------------------------- }
  var
    A       : TMatrix;  { Matrix of normal equations }
    G       : TVector;  { Constant vector }
    Det     : Float;    { Determinant of A }
    I, J, K : Integer;  { Loop variables }
    WX      : Float;
  begin
    DimMatrix(A, Nvar, Nvar);
    DimVector(G, Nvar);

    { If constant term, set line 0 and column 0 of matrix A,
      and element 0 of vecteur G }
    if ConsTerm then
      begin
        if Mode = UNWEIGHTED then
          begin
            A[0,0] := Int(N - FirstPoint + 1);
            for K := FirstPoint to N do
              begin
                for J := 1 to Nvar do
                  A[0,J] := A[0,J] + X[J,K];
                G[0] := G[0] + Y[K];
              end;
          end
        else
          begin
            for K := FirstPoint to N do
              begin
                A[0,0] := A[0,0] + W[K];
                for J := 1 to Nvar do
                  A[0,J] := A[0,J] + W[K] * X[J,K];
                G[0] := G[0] + W[K] * Y[K];
              end;
          end;
        for J := 1 to Nvar do
          A[J,0] := A[0,J];
      end;

    { Set other elements of A and G }
    if Mode = UNWEIGHTED then
      for K := FirstPoint to N do
        for I := 1 to Nvar do
          begin
            for J := I to Nvar do
              A[I,J] := A[I,J] + X[I,K] * X[J,K];
            G[I] := G[I] + X[I,K] * Y[K];
          end
    else
      for K := FirstPoint to N do
        for I := 1 to Nvar do
          begin
            WX := W[K] * X[I,K];
            for J := I to Nvar do
              A[I,J] := A[I,J] + WX * X[J,K];
            G[I] := G[I] + WX * Y[K];
          end;

    { Fill in symmetric matrix }
    for I := 2 to Nvar do
      for J := 1 to Pred(I) do
        A[I,J] := A[J,I];

    { Solve normal equations }
    if ConsTerm then
      Gauss_GenMulFit := GaussJordan(A, G, 0, Nvar, V, B, Det)
    else
      Gauss_GenMulFit := GaussJordan(A, G, 1, Nvar, V, B, Det);
  end;

  function SVD_GenMulFit(Mode     : TRegMode;
                         X        : TMatrix;
                         Y, W     : TVector;
                         N, Nvar  : Integer;
                         ConsTerm : Boolean;
                         B        : TVector;
                         V        : TMatrix) : Integer;
{ ----------------------------------------------------------------------
  General multiple linear regression routine (SVD algorithm)
  ---------------------------------------------------------------------- }
  var
    U       : TMatrix;  { Matrix of independent variables for SVD }
    Z       : TVector;  { Vector of dependent variables for SVD }
    S       : TVector;  { Singular values }
    S2inv   : TVector;  { Inverses of squared singular values }
    V1      : TMatrix;  { Orthogonal matrix from SVD }
    Lbound  : Integer;  { Lower bound of U matrix in both dims. }
    Ubound  : Integer;  { Upper bound of U matrix in 1st dim. }
    I, J, K : Integer;  { Loop variables }
    Sigma   : Float;    { Square root of weight }
    Sum     : Float;    { Element of variance-covariance matrix }
    ErrCode : Integer;  { Error code }
  begin
    if ConsTerm then
      begin
        Lbound := 0;
        Ubound := N - FirstPoint;
      end
    else
      begin
        Lbound := 1;
        Ubound := N - FirstPoint + 1;
      end;

    { Dimension arrays }
    DimMatrix(U, Ubound, Nvar);
    DimVector(Z, Ubound);
    DimVector(S, Nvar);
    DimVector(S2inv, Nvar);
    DimMatrix(V1, Nvar, Nvar);

    { ----------------------------------------------------------
      Prepare arrays for SVD :
      If constant term, use U[0..(N - FirstPoint), 0..Nvar]
                        and Z[0..(N - FirstPoint)]
      Else              use U[1..(N - FirstPoint + 1), 1..Nvar]
                        and Z[1..(N - FirstPoint + 1)]
      ---------------------------------------------------------- }
    if Mode = UNWEIGHTED then
      for I := Lbound to Ubound do
        begin
          K := I - Lbound + FirstPoint;
          Z[I] := Y[K];
          if ConsTerm then
            U[I,0] := 1.0;
          for J := 1 to Nvar do
            U[I,J] := X[J,K];
        end
    else
      for I := Lbound to Ubound do
        begin
          K := I - Lbound + FirstPoint;
          Sigma := Sqrt(W[K]);
          Z[I] := Y[K] * Sigma;
          if ConsTerm then
            U[I,0] := Sigma;
          for J := 1 to Nvar do
            U[I,J] := X[J,K] * Sigma;
        end;

    { Perform singular value decomposition }
    ErrCode := SV_Decomp(U, Lbound, Ubound, Nvar, S, V1);

    if ErrCode = MAT_OK then
      begin
        { Set the lowest singular values to zero }
        SV_SetZero(S, Lbound, Nvar, TolSVD(N - FirstPoint + 1));

        { Solve the system }
        SV_Solve(U, S, V1, Z, Lbound, Ubound, Nvar, B);

        { Compute variance-covariance matrix }
        for I := Lbound to Nvar do
          if S[I] > 0.0 then
            S2inv[I] := 1.0 / Sqr(S[I])
          else
            S2inv[I] := 0.0;
        for I := Lbound to Nvar do
          for J := Lbound to I do
            begin
              Sum := 0.0;
              for K := Lbound to Nvar do
                Sum := Sum + V1[I,K] * V1[J,K] * S2inv[K];
              V[I,J] := Sum;
              V[J,I] := Sum;
            end;
      end;

    SVD_GenMulFit := ErrCode;
  end;

  function GenMulFit(Mode     : TRegMode;
                     X        : TMatrix;
                     Y, W     : TVector;
                     N, Nvar  : Integer;
                     ConsTerm : Boolean;
                     B        : TVector;
                     V        : TMatrix) : Integer;
{ ----------------------------------------------------------------------
  General multiple linear regression routine
  ---------------------------------------------------------------------- }
  var
    ErrCode : Integer;
  begin