regress.pas

Delphi 的数学控件

{ **********************************************************************
  *                          Unit REGRESS.PAS                          *
  *                            Version 2.5d                            *
  *                      (c) J. Debord, July 2004                      *
  **********************************************************************
                           Regression routines
  ********************************************************************** }

unit regress;

interface

uses
  fmath, fspec, matrices, randnum, eigen, optim, simopt, mcmc, stat;

{ **********************************************************************
  Type definitions
  ********************************************************************** }

{ Algorithm for linear regression }
type
  TRegAlgo = (
    GAUSS_JORDAN,  { Gauss-Jordan solution of normal equations }
    SVD);          { Singular value decomposition }

{ Optimization algorithm for nonlinear regression }
type
  TOptAlgo = (
    NL_MARQ,       { Marquardt algorithm }
    NL_SIMP,       { Simplex algorithm }
    NL_BFGS,       { BFGS algorithm }
    NL_SA,         { Simulated annealing }
    NL_MH);        { Metropolis-Hastings }

{ Regression modes }
type
  TRegMode = (UNWEIGHTED, WEIGHTED);

{ Regression function }
type
  TRegFunc = function(X : Float; B : TVector) : Float;

{ Procedure to compute the derivatives of the regression function
  with respect to the regression parameters }
type
  TDerivProc = procedure(RegFunc : TRegFunc;
                         X, Y    : Float;
                         B, D    : TVector);

{ Test of regression }
type
  TRegTest = record
    Vr,              { Residual variance }
    R2,              { Coefficient of determination }
    R2a,             { Adjusted coeff. of determination }
    F,               { Variance ratio (explained/residual) }
    Prob : Float;    { Probability of F }
  end;

var
  MHFile : String;   { File for saving Metropolis-Hastings results }

{ **********************************************************************
  Procedures to modify the regression settings
  ********************************************************************** }

procedure SetRegAlgo(Algo : TRegAlgo);
{ ----------------------------------------------------------------------
  Sets the linear regression algorithm according to Algo, which must be
  GAUSS_JORDAN or SVD. The default algorithm is SVD.
  ---------------------------------------------------------------------- }

procedure SetOptAlgo(Algo : TOptAlgo);
{ ----------------------------------------------------------------------
  Sets the optimization algorithm according to Algo, which must be
  NL_MARQ, NL_SIMP, NL_BFGS, NL_SA or NL_MH. The default is NL_MARQ
  ---------------------------------------------------------------------- }

procedure SetFirstPoint(Index : Integer);
{ ----------------------------------------------------------------------
  Sets the index of the first data point (usually 0 or 1). The default
  value is 1.
  ---------------------------------------------------------------------- }

function GetRegAlgo : TRegAlgo;
{ ----------------------------------------------------------------------
  Returns the linear regression algorithm
  ---------------------------------------------------------------------- }

function GetOptAlgo : TOptAlgo;
{ ----------------------------------------------------------------------
  Returns the optimization algorithm
  ---------------------------------------------------------------------- }

function GetFirstPoint : Integer;
{ ----------------------------------------------------------------------
  Returns the index of the first data point
  ---------------------------------------------------------------------- }

{ **********************************************************************
  Unweighted regression routines
  **********************************************************************
  These routines fit equations to data by minimizing the sum of squared
  residuals :
                      SS = Sum [y(k) - ycalc(k)]^2

  where y(k) and ycalc(k) are respectively the observed and calculated
  value of the dependent variable for observation k. ycalc(k) is a
  function of the regression parameters b(0), b(1) ...

  The following regression types are implemented :

  * Simple linear regression :

                        y(k) = b(0) + b(1) * x(k)

  * Multiple linear regression :

    y(k) = b(0) + b(1) * x(1,k) + b(2) * x(2,k) + ... + b(Nvar) * x(Nvar,k)

  * Polynomial regression :

      y(k) = b(0) + b(1) * x(k) + b(2) * x(k)^2 + ... + b(Deg) * x(k)^Deg

  * Nonlinear regression :

                   y(k) = f[x(k), b(0), b(1), ... ]

    where f is a user-specified function.

  The following parameters are common to all routines :

  Input  : X = Vector or matrix of independent variables
           Y = Vector of dependent variable
           N = Index of the last observation
  Output : B = Regression parameters
           V = Inverse matrix of normal equations
  ********************************************************************** }

function LinFit(X, Y : TVector;
                N    : Integer;
                B    : TVector;
                V    : TMatrix) : Integer;
{ ----------------------------------------------------------------------
  Simple linear regression
  ---------------------------------------------------------------------- }

function MulFit(X        : TMatrix;
                Y        : TVector;
                N, Nvar  : Integer;
                ConsTerm : Boolean;
                B        : TVector;
                V        : TMatrix) : Integer;
{ ----------------------------------------------------------------------
  Multiple linear regression
  ----------------------------------------------------------------------
  Additional input parameters :
  Nvar     = Index of the last independent variable
  ConsTerm = Flags the presence of a constant term b(0)
  ---------------------------------------------------------------------- }

function PolFit(X, Y   : TVector;
                N, Deg : Integer;
                B      : TVector;
                V      : TMatrix) : Integer;
{ ----------------------------------------------------------------------
  Polynomial regression
  ----------------------------------------------------------------------
  Additional input parameter :
  Deg = Degree of polynomial
  ---------------------------------------------------------------------- }

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;
{ ----------------------------------------------------------------------
  Nonlinear regression
  ----------------------------------------------------------------------
  Additional input parameters :
  RegFunc        = Regression function
  DerivProc      = Procedure to compute the derivatives of RegFunc
  Lbound, Ubound = Indices of first and last function parameters
  MaxIter        = Maximum number of iterations
  Tol            = Required parameter precision
  B              = Initial parameter values
  B_min, B_max   = Lower and upper parameter bounds
  ---------------------------------------------------------------------- }

{ **********************************************************************
  Weighted regression routines
  **********************************************************************
  These routines fit equations to data by minimizing the sum of weighted
  squared residuals :

                      SWS = Sum w(k)*[y(k) - ycalc(k)]^2

  where the "weight" w(k) is inversely proportional to the variance v(k)
  of the observation y(k). v(k) is usually computed as :

                       v(k) = Vr * g[y(k)] = Vr / w(k)

  where Vr is the residual variance and g is a user-specified function
  (e.g. g[y(k)] = y(k)^2 for a constant coefficient of variation).

  Function syntax and results are the same than for unweighted regression
  except that the vector of weights (W) is passed as an additional input
  parameter.
  ********************************************************************** }

function WLinFit(X, Y, W : TVector;
                 N       : Integer;
                 B       : TVector;
                 V       : TMatrix) : Integer;

function WMulFit(X        : TMatrix;
                 Y, W     : TVector;
                 N, Nvar  : Integer;
                 ConsTerm : Boolean;
                 B        : TVector;
                 V        : TMatrix) : Integer;

function WPolFit(X, Y, W : TVector;
                 N, Deg  : Integer;
                 B       : TVector;
                 V       : TMatrix) : Integer;

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;

{ **********************************************************************
  Procedure to compute the derivatives of the regression function by
  numerical differentiation.
  ********************************************************************** }

procedure NumDeriv(RegFunc : TRegFunc;
                   X, Y    : Float;
                   B, D    : TVector);
{ ----------------------------------------------------------------------
  Input parameters  : RegFunc = Regression function
                      X, Y    = Coordinates of point
                      B       = Regression parameters

  Output parameter  : D       = Derivatives (D[I] contains the
                                derivative w.r.t. parameter B[I])
  ---------------------------------------------------------------------- }

{ **********************************************************************
  Routines to test the quality of the regression
  **********************************************************************
  These routines compute the variance-covariance matrix of the fitted
  parameters and the different statistics used to test the quality of
  the fit.

  Input parameters  : Y      = Vector of dependent variable
                      Ycalc  = Computed Y values
                      W      = Vector of weights (if any)
                      N      = Index of the last observation
                      Lbound,
                      Ubound = Indices of first & last fitted parameters
                      V      = Inverse normal equations matrix

  Output parameters : V      = Variance-covariance matrix
                      Test   = Test statistics (Vr, R2, R2a, F, Prob)
  ********************************************************************** }

procedure RegTest(Y, Ycalc          : TVector;
                  N, Lbound, Ubound : Integer;
                  V                 : TMatrix;
                  var Test          : TRegTest);
{ ----------------------------------------------------------------------
  Test of unweighted regression
  ---------------------------------------------------------------------- }

procedure WRegTest(Y, Ycalc, W       : TVector;
                   N, Lbound, Ubound : Integer;
                   V                 : TMatrix;
                   var Test          : TRegTest);
{ ----------------------------------------------------------------------
  Test of weighted regression
  ---------------------------------------------------------------------- }

{ **********************************************************************
  Test of regression parameters
  ********************************************************************** }

procedure ParamTest(B                 : TVector;
                    V                 : TMatrix;
                    N, Lbound, Ubound : Integer;
                    S, T, Prob        : TVector);
{ ----------------------------------------------------------------------
  This routine tests the significance of the parameters. It must be
  called AFTER RegTest or WRegTest since it uses the variance-covariance
  matrix.
  ----------------------------------------------------------------------
  Input parameters  : B      = Regression parameters
                      V      = Variance-covariance matrix
                      N      = Index of the last observation
                      Lbound,
                      Ubound = Indices of first & last fitted parameters
  ----------------------------------------------------------------------
  Output parameters : S      = Standard deviations of parameters
                      T      = Student's t
                      Prob   = Probabilities
  ---------------------------------------------------------------------- }

{ **********************************************************************
  Correlation and principal component analysis

  Common parameters:

  X              = matrix of variables (X[I] contains the I-th variable)
  N              = Index of the last observation
  Lbound, Ubound = Indices of first & last variables
  M              = Mean vector (M[I] = mean of X[I])
  S              = Vector of standard deviations
  V              = Variance-covariance matrix
  R              = Correlation matrix
  ********************************************************************** }

procedure VecMean(X                 : TMatrix;
                  N, Lbound, Ubound : Integer;
                  M                 : TVector);
{ ----------------------------------------------------------------------
  Computes the mean vector (M) from matrix X

  Input  : X, Lbound, Ubound
  Output : M
  ---------------------------------------------------------------------- }

procedure VecSD(X                 : TMatrix;
                N, Lbound, Ubound : Integer;
                M, S              : TVector);
{ ----------------------------------------------------------------------
  Computes the vector of standard deviations (S) from matrix X

  Input  : X, Lbound, Ubound, M
  Output : S
  ---------------------------------------------------------------------- }

procedure MatVarCov(X                 : TMatrix;
                    N, Lbound, Ubound : Integer;
                    M                 : TVector;
                    V                 : TMatrix);
{ ----------------------------------------------------------------------
  Computes the variance-covariance matrix (V) from matrix X

  Input  : X, Lbound, Ubound, M
  Output : V
  ---------------------------------------------------------------------- }

procedure MatCorrel(V              : TMatrix;
                    Lbound, Ubound : Integer;
                    R              : TMatrix);
{ ----------------------------------------------------------------------
  Computes the correlation matrix (R) from the variance-covariance
  matrix (V)

  Input  : V, Lbound, Ubound
  Output : R
  ---------------------------------------------------------------------- }

function PCA(R              : TMatrix;
             Lbound, Ubound : Integer;
             Lambda         : TVector;
             C, Rc          : TMatrix) : Integer;