lazarus-ccr/applications/lazstats/source/units/mathunit.pas

unit MathUnit;

{ extract some math functions from functionslib for easier testing }

{$mode objfpc}{$H+}

interface

uses
  Classes, SysUtils;

function Beta(a, b: Double): Extended;
function BetaI(a,b,x: Double): Extended;

function GammaLn(x: double): Extended;

function tDist(x: Double; N: Integer; OneSided: Boolean): Double;
function tDensity(x: Double; N: Integer): Double;

function FDensity(x: Double; DF1, DF2: Integer): Double;

function Chi2Density(x: Double; N: Integer): Double;

implementation

uses
  Math;

function Beta(a, b: Double): Extended;
begin
  if (a > 0) and (b > 0) then
    Result := exp(GammaLn(a) + GammaLn(b) - GammaLn(a+b))
  else
    raise Exception.Create('Invalid argument for beta function.');
end;

function BetaCF(a, b, x: double): Extended;
const
  itmax = 100;
  eps = 3.0e-7;
var
  tem,qap,qam,qab,em,d: extended;
  bz,bpp,bp,bm,az,app: extended;
  am,aold,ap: extended;
  term1, term2, term3, term4, term5, term6: extended;
  m: integer;
BEGIN
  am := 1.0;
  bm := 1.0;
  az := 1.0;
  qab := a+b;
  qap := a+1.0;
  qam := a-1.0;
  bz := 1.0 - qab * x / qap;
  for m := 1 to itmax do begin
    em := m;
    tem := em+em;
    d := em * (b - m) * x / ((qam + tem) * (a + tem));
    ap := az + d * am;
    bp := bz + d * bm;
    term1 := -(a + em);
    term2 := qab + em;
    term3 := term1 * term2 * x;
    term4 := a + tem;
    term5 := qap + tem;
    term6 := term4 * term5;
    d := term3 / term6;
    app := ap + d * az;
    bpp := bp + d * bz;
    aold := az;
    am := ap/bpp;
    bm := bp/bpp;
    az := app/bpp;
    bz := 1.0;
    if ((abs(az-aold)) < (eps*abs(az))) then
      Break;
  end;

  { ShowMessage('WARNING! a or b too big, or itmax too small in betacf');}
  Result := az
end;

function BetaI(a,b,x: Double): extended;
var
  bt: extended;
  term1, term2, term3, term4, term5: extended;
begin
  if ((x < 0.0) or (x > 1.0)) then begin
     { ShowMessage('ERROR! Problem in routine BETAI');}
    Result := 0.5;
    exit;
  end;

  if ((x <= 0.0) or (x >= 1.0)) then
    bt := 0.0
  else
  begin
    term1 := GammaLn(a + b) - GammaLn(a) - GammaLn(b);
    term2 := a * ln(x);
    term3 := b * ln(1.0 - x);
    term4 := term1 + term2 + term3;
    bt    := exp(term4);
  end;
  term5 := (a + 1.0) / (a + b + 2.0);
  if x < term5 then
    Result := bt * betacf(a,b,x) / a
  else
    Result := 1.0 - bt * betacf(b,a,1.0-x) / b
end;

function GammaLn(x: double): Extended;
var
  tmp, ser: double;
  cof: array[0..5] of double;
  j: integer;
begin
   cof[0] := 76.18009173;
   cof[1] := -86.50532033;
   cof[2] := 24.01409822;
   cof[3] := -1.231739516;
   cof[4] := 0.00120858003;
   cof[5] := -0.00000536382;

   x := x - 1.0;
   tmp := x + 5.5;
   tmp := tmp - (x + 0.5) * ln(tmp);
   ser := 1.0;
   for j := 0 to 5 do
   begin
     x := x + 1.0;
     ser := ser + cof[j] / x;
   end;
   Result := -tmp + ln(2.50662827465 * ser);
end;

// Calculates the (cumulative) t distribution function for N degrees of freedom
function tDist(x: Double; N: Integer; OneSided: Boolean): Double;
begin
  Result :=  betai(0.5*N, 0.5, N/(N + sqr(x)));
  if OneSided then Result := Result * 0.5;
end;

// Returns the density curve for the t statistic with N degrees of freedom
function tDensity(x: Double; N: Integer): Double;
var
  factor: Double;
begin
  factor := exp(gammaLn((N+1)/2) - gammaLn(N/2)) / sqrt(N * pi);
  Result := factor * Power(1 + sqr(x)/N, (1-n)/2);
end;

// Returns the density curve for the F statistic for DF1 and DF2 degrees of freedom.
function FDensity(x: Double; DF1, DF2: Integer): Double;
var
  ratio1, ratio2, ratio3, ratio4: double;
  part1, part2, part3, part4, part5, part6, part7, part8, part9: double;
begin
  ratio1 := (DF1 + DF2) / 2.0;
  ratio2 := (DF1 - 2.0) / 2.0;
  ratio3 := DF1 / 2.0;
  ratio4 := DF2 / 2.0;
  part1 := exp(gammaln(ratio1));
  part2 := power(DF1, ratio3);
  part3 := power(DF2, ratio4);
  part4 := exp(gammaln(ratio3));
  part5 := exp(gammaln(ratio4));
  part6 := power(x, ratio2);
  part7 := power((x*DF1 + DF2), ratio1);
  part8 := (part1 * part2 * part3) / (part4 * part5);
  if (part7 = 0.0) then
    part9 := 0.0
  else
    part9 := part6 / part7;
  Result := part8 * part9;
end;


// Returns the density curve of the chi2 statistic for N degrees of freedom
function Chi2Density(x: Double; N: Integer): Double;
var
  factor: Double;
begin
  factor := Power(2.0, N * 0.5) * exp(gammaLN(N * 0.5));
  Result := power(x, (N-2.0) * 0.5) / (factor * exp(x * 0.5));
end;

end.