fspec.pas

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

{ **********************************************************************
  *                            Unit FSPEC.PAS                          *
  *                             Version 1.2d                           *
  *                       (c) J. Debord, May 2003                      *
  **********************************************************************
             Special functions and probability distributions
  **********************************************************************
  Error handling: The global variable MathError returns the error code
  from the last function evaluation. It must be checked immediately
  after a function call:

       Y := f(X);        (* f is one of the functions of the library *)
       if MathError = FN_OK then ...

  The possible error codes are the following:

     ---------------------------------------------
     Error code   Value  Significance
     ---------------------------------------------
     FN_OK          0    No error
     FN_DOMAIN     -1    Argument domain error
     FN_SING       -2    Function singularity
     FN_OVERFLOW   -3    Overflow range error
     FN_UNDERFLOW  -4    Underflow range error
     FN_TLOSS      -5    Total loss of precision
     FN_PLOSS      -6    Partial loss of precision
     ---------------------------------------------

  **********************************************************************
  Most functions are translated C code from Cephes math library
  by S. Moshier (http://www.moshier.net)
  ********************************************************************** }

unit fspec;

interface

uses
  fmath;

{ ----------------------------------------------------------------------
  Table of factorials
  ---------------------------------------------------------------------- }

const
  NFACT = 33;

var
  FactArray : array[0..NFACT] of Float =
    (1.0,
     1.0,
     2.0,
     6.0,
     24.0,
     120.0,
     720.0,
     5040.0,
     40320.0,
     362880.0,
     3628800.0,
     39916800.0,
     479001600.0,
     6227020800.0,
     87178291200.0,
     1307674368000.0,
     20922789888000.0,
     355687428096000.0,
     6402373705728000.0,
     121645100408832000.0,
     2432902008176640000.0,
     51090942171709440000.0,
     1124000727777607680000.0,
     25852016738884976640000.0,
     620448401733239439360000.0,
     15511210043330985984000000.0,
     403291461126605635584000000.0,
     10888869450418352160768000000.0,
     304888344611713860501504000000.0,
     8841761993739701954543616000000.0,
     265252859812191058636308480000000.0,
     8222838654177922817725562880000000.0,
     263130836933693530167218012160000000.0,
     8683317618811886495518194401280000000.0);

{ ----------------------------------------------------------------------
  Special functions
  ---------------------------------------------------------------------- }

function Fact(N : Integer) : Float;         { Factorial }
function Binomial(N, K : Integer) : Float;  { Binomial coef. C(N,K) }
function Gamma(X : Float) : Float;          { Gamma function }
function SgnGamma(X : Float) : Integer;     { Sign of Gamma function }
function IGamma(A, X : Float) : Float;      { Incomplete Gamma function }
function JGamma(A, X : Float) : Float;      { Complement of IGamma }
function Beta(X, Y : Float) : Float;        { Beta function }
function IBeta(A, B, X : Float) : Float;    { Incomplete Beta function }
function Erf(X : Float) : Float;            { Error function }
function Erfc(X : Float) : Float;           { Complement of Erf }

function LnGamma(X : Float)   : Float;   overload;  { Log(|Gamma(X)|) }
function LnGamma(Z : Complex) : Complex; overload;

{ ----------------------------------------------------------------------
  Binomial distribution with probability P and number of repetitions N
  ---------------------------------------------------------------------- }

function PBinom(N : Integer; P : Float; K : Integer) : Float; { Prob(X = K) }
function FBinom(N : Integer; P : Float; K : Integer) : Float; { Prob(X <= K) }

{ ----------------------------------------------------------------------
  Poisson distribution with mean Mu
  ---------------------------------------------------------------------- }

function PPoisson(Mu : Float; K : Integer) : Float;  { Prob(X = K) }
function FPoisson(Mu : Float; K : Integer) : Float;  { Prob(X <= K) }

{ ----------------------------------------------------------------------
  Standard normal distribution
  ---------------------------------------------------------------------- }

function DNorm(X : Float) : Float;    { Density of standard normal }
function FNorm(X : Float) : Float;    { Prob(U <= X) }
function PNorm(X : Float) : Float;    { Prob(|U| >= |X|) }
function InvNorm(P : Float) : Float;  { Inverse of FNorm : returns X
                                        such that Prob(U <= X) = P}

{ ----------------------------------------------------------------------
  Student distribution with Nu d.o.f.
  ---------------------------------------------------------------------- }

function DStudent(Nu : Integer; X : Float) : Float;  { Density of t }
function FStudent(Nu : Integer; X : Float) : Float;  { Prob(t <= X) }
function PStudent(Nu : Integer; X : Float) : Float;  { Prob(|t| >= |X|) }

{ ----------------------------------------------------------------------
  Khi-2 distribution with Nu d.o.f.
  ---------------------------------------------------------------------- }

function DKhi2(Nu : Integer; X : Float) : Float;  { Density of Khi2 }
function FKhi2(Nu : Integer; X : Float) : Float;  { Prob(Khi2 <= X) }
function PKhi2(Nu : Integer; X : Float) : Float;  { Prob(Khi2 >= X) }

{ ----------------------------------------------------------------------
  Fisher-Snedecor distribution with Nu1 and Nu2 d.o.f.
  ---------------------------------------------------------------------- }

function DSnedecor(Nu1, Nu2 : Integer; X : Float) : Float;  { Density of F }
function FSnedecor(Nu1, Nu2 : Integer; X : Float) : Float;  { Prob(F <= X) }
function PSnedecor(Nu1, Nu2 : Integer; X : Float) : Float;  { Prob(F >= X) }

{ ----------------------------------------------------------------------
  Exponential distribution
  ---------------------------------------------------------------------- }

function DExpo(A, X : Float) : Float;  { Density of exponential distrib. }
function FExpo(A, X : Float) : Float;  { Prob( <= X) }

{ ----------------------------------------------------------------------
  Beta distribution
  ---------------------------------------------------------------------- }

function DBeta(A, B, X : Float) : Float;   { Density of Beta distribution }
function FBeta(A, B, X : Float) : Float;   { Prob( <= X) }

{ ----------------------------------------------------------------------
  Gamma distribution
  ---------------------------------------------------------------------- }

function DGamma(A, B, X : Float) : Float;  { Density of Gamma distribution }
function FGamma(A, B, X : Float) : Float;  { Prob( <= X) }

{ ********************************************************************** }

implementation

{ ----------------------------------------------------------------------
  Error handling function
  ---------------------------------------------------------------------- }

  function DefaultVal(ErrCode : Integer; Default : Float) : Float;
  { Sets the global variable MathError and the function default value
    according to the error code }
  begin
    MathError := ErrCode;
    DefaultVal := Default;
  end;

{ ----------------------------------------------------------------------
  Special functions
  ---------------------------------------------------------------------- }

const { Used by IGamma and IBeta }
  BIG    = 9.223372036854775808E18;
  BIGINV = 1.084202172485504434007E-19;

type
  TabCoef = array[0..9] of Float;

  function PolEvl(var X : Float; Coef : TabCoef; N : Integer) : Float;
{ ----------------------------------------------------------------------
  Evaluates polynomial of degree N:

			2	   N
    y  =  C  + C x + C x  +...+ C x
	   0	1     2 	 N

  Coefficients are stored in reverse order:

  Coef[0] = C  , ..., Coef[N] = C
             N                   0

  The function P1Evl() assumes that Coef[N] = 1.0 and is
  omitted from the array. Its calling arguments are
  otherwise the same as PolEvl().
  ---------------------------------------------------------------------- }
  var
    Ans : Float;
    I : Integer;
  begin
    Ans := Coef[0];
    for I := 1 to N do
      Ans := Ans * X + Coef[I];
    PolEvl := Ans;
  end;

  function P1Evl(var X : Float; Coef : TabCoef; N : Integer) : Float;
{ ----------------------------------------------------------------------
  Evaluate polynomial when coefficient of X is 1.0.
  Otherwise same as PolEvl.
  ---------------------------------------------------------------------- }
  var
    Ans : Float;
    I : Integer;
  begin
    Ans := X + Coef[0];
    for I := 1 to N - 1 do
      Ans := Ans * X + Coef[I];
    P1Evl := Ans;
  end;

  function SgnGamma(X : Float) : Integer;
  begin
    if X > 0.0 then
      SgnGamma := 1
    else if Odd(Trunc(Abs(X))) then
      SgnGamma := 1
    else
      SgnGamma := - 1;
  end;

  function Stirf(X : Float) : Float;
  { Stirling's formula for the gamma function
    Gamma(x) = Sqrt(2*Pi) x^(x-.5) exp(-x) (1 + 1/x P(1/x))
    where P(x) is a polynomial }
  const
    STIR : TabCoef = (
        7.147391378143610789273E-4,
      - 2.363848809501759061727E-5,
      - 5.950237554056330156018E-4,
        6.989332260623193171870E-5,
        7.840334842744753003862E-4,
      - 2.294719747873185405699E-4,
      - 2.681327161876304418288E-3,
        3.472222222230075327854E-3,
        8.333333333333331800504E-2,
        0);

  var
    W, P : Float;
  begin
    W := 1.0 / X;
    if X > 1024.0 then
      begin
        P := 6.97281375836585777429E-5 * W + 7.84039221720066627474E-4;
        P := P * W - 2.29472093621399176955E-4;
        P := P * W - 2.68132716049382716049E-3;
        P := P * W + 3.47222222222222222222E-3;
        P := P * W + 8.33333333333333333333E-2;
      end
    else
      P := PolEvl(W, STIR, 8);
    Stirf := SQRT2PI * Exp((X - 0.5) * Ln(X) - X) * (1.0 + W * P);
  end;

  function GamSmall(X1, Z : Float) : Float;
  { Gamma function for small values of the argument }
  const
    S : TabCoef = (
      - 1.193945051381510095614E-3,
        7.220599478036909672331E-3,
      - 9.622023360406271645744E-3,
      - 4.219773360705915470089E-2,
        1.665386113720805206758E-1,
      - 4.200263503403344054473E-2,
      - 6.558780715202540684668E-1,
        5.772156649015328608253E-1,
        1.000000000000000000000E0,
        0);

    SN : TabCoef = (
        1.133374167243894382010E-3,
        7.220837261893170325704E-3,
        9.621911155035976733706E-3,
      - 4.219773343731191721664E-2,
      - 1.665386113944413519335E-1,
      - 4.200263503402112910504E-2,
        6.558780715202536547116E-1,
        5.772156649015328608727E-1,
      - 1.000000000000000000000E0,
        0);

  var
    P : Float;
  begin
    if X1 = 0.0 then
      begin
        GamSmall := DefaultVal(FN_SING, MAXNUM);
        Exit;
      end;
    if X1 < 0.0 then
      begin
        X1 := - X1;
        P := PolEvl(X1, SN, 8);
      end
    else
      P := PolEvl(X1, S, 8);
    GamSmall := Z / (X1 * P);
  end;

  function StirfL(X : Float) : Float;
  { Approximate Ln(Gamma) by Stirling's formula, for X >= 13 }
  const
    P : TabCoef = (
        4.885026142432270781165E-3,
      - 1.880801938119376907179E-3,
        8.412723297322498080632E-4,
      - 5.952345851765688514613E-4,
        7.936507795855070755671E-4,
      - 2.777777777750349603440E-3,
        8.333333333333331447505E-2,
        0, 0, 0);

  var
    Q, W : Float;
  begin
    Q := Ln(X) * (X - 0.5) - X;
    Q := Q + LNSQRT2PI;
    if X > 1.0E+10 then
      StirfL := Q
    else
      begin
        W := 1.0 / Sqr(X);
        StirfL := Q + PolEvl(W, P, 6) / X;
      end;
  end;

  function Gamma(X : Float) : Float;
  const
    P : TabCoef = (
      4.212760487471622013093E-5,
      4.542931960608009155600E-4,
      4.092666828394035500949E-3,
      2.385363243461108252554E-2,
      1.113062816019361559013E-1,
      3.629515436640239168939E-1,
      8.378004301573126728826E-1,
      1.000000000000000000009E0,
      0, 0);

    Q : TabCoef = (
      - 1.397148517476170440917E-5,
        2.346584059160635244282E-4,
      - 1.237799246653152231188E-3,
      - 7.955933682494738320586E-4,
        2.773706565840072979165E-2,
      - 4.633887671244534213831E-2,
      - 2.243510905670329164562E-1,
        4.150160950588455434583E-1,
        9.999999999999999999908E-1,
        0);

  var
    SgnGam, N : Integer;
    A, X1, Z : Float;
  begin
    MathError := FN_OK;
    SgnGam := SgnGamma(X);

    if (X = 0.0) or ((X < 0.0) and (Frac(X) = 0.0)) then
      begin
        Gamma := DefaultVal(FN_SING, SgnGam * MAXNUM);
        Exit;
      end;

    if X > MAXGAM then
      begin
        Gamma := DefaultVal(FN_OVERFLOW, MAXNUM);
        Exit;
      end;

    A := Abs(X);
    if A > 13.0 then
      begin
        if X < 0.0 then
          begin
            N := Trunc(A);
            Z := A - N;
            if Z > 0.5 then
              begin
                N := N + 1;
                Z := A - N;
              end;
            Z := Abs(A * Sin(PI * Z)) * Stirf(A);
            if Z <= PI / MAXNUM then
              begin
                Gamma := DefaultVal(FN_OVERFLOW, SgnGam * MAXNUM);
                Exit;
              end;
            Z := PI / Z;
          end
        else
          Z := Stirf(X);
        Gamma := SgnGam * Z;
      end
    else
      begin
        Z := 1.0;
        X1 := X;
        while X1 >= 3.0 do
          begin
            X1 := X1 - 1.0;
            Z := Z * X1;
          end;
        while X1 < - 0.03125 do
          begin
            Z := Z / X1;
            X1 := X1 + 1.0;
          end;
        if X1 <= 0.03125 then
          Gamma := GamSmall(X1, Z)
        else
          begin
            while X1 < 2.0 do
              begin
                Z := Z / X1;
                X1 := X1 + 1.0;
              end;
            if (X1 = 2.0) or (X1 = 3.0) then
              Gamma := Z
            else
              begin
                X1 := X1 - 2.0;
                Gamma := Z * PolEvl(X1, P, 7) / PolEvl(X1, Q, 8);
              end;
          end;
      end;
  end;

  function LnGamma(X : Float) : Float; overload;
  const
    P : TabCoef = (
      - 2.163690827643812857640E3,
      - 8.723871522843511459790E4,
      - 1.104326814691464261197E6,
      - 6.111225012005214299996E6,
      - 1.625568062543700591014E7,
      - 2.003937418103815175475E7,
      - 8.875666783650703802159E6,
        0, 0, 0);

    Q : TabCoef = (
      - 5.139481484435370143617E2,
      - 3.403570840534304670537E4,
      - 6.227441164066219501697E5,
      - 4.814940379411882186630E6,
      - 1.785433287045078156959E7,
      - 3.138646407656182662088E7,
      - 2.099336717757895876142E7,
        0, 0, 0);

  var
    N, SgnGam : Integer;
    A, X1, Z : Float;
  begin
    MathError := FN_OK;
    SgnGam := SgnGamma(X);

    if (X = 0.0) or ((X < 0.0) and (Frac(X) = 0.0)) then
      begin
        LnGamma := DefaultVal(FN_SING, MAXNUM);
        Exit;
      end;

    if X > MAXLGM then
      begin
        LnGamma := DefaultVal(FN_OVERFLOW, MAXNUM);
        Exit;
      end;

    A := Abs(X);
    if A > 34.0 then
      begin
        if X < 0.0 then
          begin
            N := Trunc(A);
            Z := A - N;
            if Z > 0.5 then
              begin
                N := N + 1;
                Z := N - A;
              end;