mirror of
				https://gitlab.com/freepascal.org/fpc/source.git
				synced 2025-11-04 02:19:22 +01:00 
			
		
		
		
	* paths from mantis #31652 by WP
- deprecated functions with duplicate in math.
    - overloaded some functions with "is nested" procvars.
    - added new functions:
        - (complete) beta function, beta()
        - incomplete beta function (betai) and its inverse (invbetai)
        - cumulative normal distribution (normaldist) and its inverse (invnormaldist)
        - chi-square distribution (chi2dist), and its inverse (invchidist)
        - Student's t distribution (tdist) and its inverse (invtdist)
        - F distribution (Fdist) and its inverse (invFdist)
        - inverse of the incomplete gamma function gammaq (invgammaq)
git-svn-id: trunk@35922 -
			
			
This commit is contained in:
		
							parent
							
								
									5165497498
								
							
						
					
					
						commit
						2a77e690b4
					
				@ -18,6 +18,9 @@
 | 
			
		||||
 | 
			
		||||
 **********************************************************************}
 | 
			
		||||
 | 
			
		||||
{$mode objfpc}{$H+}
 | 
			
		||||
{$modeswitch nestedprocvars}
 | 
			
		||||
 | 
			
		||||
Unit roo;
 | 
			
		||||
{$i direct.inc}
 | 
			
		||||
 | 
			
		||||
@ -34,6 +37,8 @@ Procedure roobin(n: ArbInt; a: complex; Var z: complex; Var term: ArbInt);
 | 
			
		||||
 | 
			
		||||
Procedure roof1r(f: rfunc1r; a, b, ae, re: ArbFloat; Var x: ArbFloat;
 | 
			
		||||
                 Var term: ArbInt);
 | 
			
		||||
Procedure roof1rn(f: rfunc1rn; a, b, ae, re: ArbFloat; Var x: ArbFloat;
 | 
			
		||||
                 Var term: ArbInt);
 | 
			
		||||
 | 
			
		||||
{Determine all zeropoints for a given n'th degree polynomal with real
 | 
			
		||||
coefficients}
 | 
			
		||||
@ -45,7 +50,7 @@ Procedure roopol(Var a: ArbFloat; n: ArbInt; Var z: complex;
 | 
			
		||||
 | 
			
		||||
Procedure rooqua(p, q: ArbFloat; Var z1, z2: complex);
 | 
			
		||||
 | 
			
		||||
{Roofnr is undocumented, but verry big}
 | 
			
		||||
{Solve a system of non-linear equations}
 | 
			
		||||
 | 
			
		||||
Procedure roofnr(f: roofnrfunc; n: ArbInt; Var x, residu: ArbFloat; re: ArbFloat;
 | 
			
		||||
                 Var term: ArbInt);
 | 
			
		||||
@ -141,13 +146,24 @@ End {roobin};
 | 
			
		||||
Procedure roof1r(f: rfunc1r; a, b, ae, re: ArbFloat; Var x: ArbFloat;
 | 
			
		||||
                 Var term: ArbInt);
 | 
			
		||||
 | 
			
		||||
  function nested_f(x: ArbFloat): ArbFloat;
 | 
			
		||||
  begin
 | 
			
		||||
    Result := f(x);
 | 
			
		||||
  end;
 | 
			
		||||
 | 
			
		||||
begin
 | 
			
		||||
  roof1rn(@nested_f, a, b, ae, re, x, term);
 | 
			
		||||
end;  
 | 
			
		||||
 | 
			
		||||
Procedure roof1rn(f: rfunc1rn; a, b, ae, re: ArbFloat; Var x: ArbFloat;
 | 
			
		||||
                 Var term: ArbInt);
 | 
			
		||||
Var fa, fb, c, fc, m, tol, w1, w2 : ArbFloat;
 | 
			
		||||
                                k : ArbInt;
 | 
			
		||||
                             stop : boolean;
 | 
			
		||||
 | 
			
		||||
Begin
 | 
			
		||||
  fa := f(a);
 | 
			
		||||
 fb := f(b);
 | 
			
		||||
  fb := f(b);
 | 
			
		||||
  If (spesgn(fa)*spesgn(fb)=1) Or (ae<0) Or (re<0)
 | 
			
		||||
   Then  {wrong input}
 | 
			
		||||
    Begin
 | 
			
		||||
@ -173,7 +189,7 @@ Begin
 | 
			
		||||
  k := 0;
 | 
			
		||||
  tol := ae+re*spemax(abs(a), abs(b));
 | 
			
		||||
  w1 := abs(b-a);
 | 
			
		||||
 stop := false;
 | 
			
		||||
  stop := false;
 | 
			
		||||
  while (abs(b-a)>tol) and (fb<>0) and (Not stop) Do
 | 
			
		||||
    Begin
 | 
			
		||||
      m := (a+b)/2;
 | 
			
		||||
 | 
			
		||||
@ -60,9 +60,17 @@ function speent(x: ArbFloat): longint;
 | 
			
		||||
{  Errorfunction ( 2/sqrt(pi)* Int(t,0,pi,exp(sqr(t)) )}
 | 
			
		||||
function speerf(x: ArbFloat): ArbFloat;
 | 
			
		||||
 | 
			
		||||
{  Errorfunction's complement ( 2/sqrt(pi)* Int(t,pi,inf,exp(sqr(t)) )}
 | 
			
		||||
{  Errorfunction's complement ( 2/sqrt(pi)* Int(t,pi,inf,exp(sqr(t))) )}
 | 
			
		||||
function speefc(x: ArbFloat): ArbFloat;
 | 
			
		||||
 | 
			
		||||
{  Calculates the cumulative normal distribution
 | 
			
		||||
   N(x) = 1/sqrt(2*pi) * Int(t, -INF, x, exp(t^2/2) ) }
 | 
			
		||||
function normaldist(x: ArbFloat): ArbFloat;
 | 
			
		||||
 | 
			
		||||
{  Inverse of cumulative normal distribution:
 | 
			
		||||
   Returns x such that y = normaldist(x) }
 | 
			
		||||
function invnormaldist(y: ArbFloat): ArbFloat;
 | 
			
		||||
 | 
			
		||||
{  Function to calculate the Gamma function ( int(t,0,inf,t^(x-1)*exp(-t)) }
 | 
			
		||||
function spegam(x: ArbFloat): ArbFloat;
 | 
			
		||||
 | 
			
		||||
@ -77,43 +85,73 @@ function gammap(s, x: ArbFloat): ArbFloat;
 | 
			
		||||
     int(t,x,inf,exp(-t)t^(s-1)) / spegam(s)  (s > 0)
 | 
			
		||||
     gammaq(s,x) = 1 - gammap(s,x) }
 | 
			
		||||
function gammaq(s, x: ArbFloat): ArbFloat;
 | 
			
		||||
function invgammaq(s, y: ArbFloat): ArbFloat;
 | 
			
		||||
 | 
			
		||||
{  Function to calculate the (complete) beta function
 | 
			
		||||
      beta(a, b) = int(t, 0, 1, t^(a-1) * (1-t)^(b-1) with a > 0, b > 0
 | 
			
		||||
      beta(a, b) = spegam(a) * spegam(b) / spegam(a + b) }
 | 
			
		||||
function beta(a, b: ArbFloat): ArbFloat;
 | 
			
		||||
 | 
			
		||||
{  Function to calculate the (regularized) incomplete beta function
 | 
			
		||||
      betai(a, b, x) = int(t, 0, x, t^(x-1) * (1-t)^(y-1) ) / beta(a,b) }
 | 
			
		||||
function betai(a, b, x: ArbFloat): ArbFloat;
 | 
			
		||||
function invbetai(a, b, y: ArbFloat; eps: ArbFloat = 0.0): ArbFloat;
 | 
			
		||||
 | 
			
		||||
{ Function to calculate the cumulative chi2 distribution with n degrees of
 | 
			
		||||
  freedom (upper tail) }
 | 
			
		||||
function chi2dist(x: ArbFloat; n: ArbInt): ArbFloat;
 | 
			
		||||
function invchi2dist(y: Arbfloat; n: ArbInt): ArbFloat;
 | 
			
		||||
 | 
			
		||||
{  Function to calculate Student's t distribution with n degrees of freedom
 | 
			
		||||
   (cumulative, upper tail if Tails = 1, else both tails }
 | 
			
		||||
type
 | 
			
		||||
  TNumTails = 1..2;
 | 
			
		||||
 | 
			
		||||
function tdist(t: ArbFloat; n: ArbInt; Tails: TNumTails): ArbFloat;
 | 
			
		||||
function invtdist(y: ArbFloat; n: ArbInt; Tails: TNumTails; eps: ArbFloat = 0.0): ArbFloat;
 | 
			
		||||
 | 
			
		||||
{  Function to calculate the cumulative F distribution function of value F
 | 
			
		||||
   with n1 and n2 degrees of freedom }
 | 
			
		||||
function Fdist(F: ArbFloat; n1, n2: ArbInt): ArbFloat;
 | 
			
		||||
function invFdist(p: ArbFloat; n1, n2: ArbInt; eps: ArbFloat = 0.0): ArbFloat;
 | 
			
		||||
 | 
			
		||||
{  "Calculates" the maximum of two ArbFloat values     }
 | 
			
		||||
function spemax(a, b: ArbFloat): ArbFloat;
 | 
			
		||||
function spemax(a, b: ArbFloat): ArbFloat; deprecated 'Use max(a,b) in unit math.';
 | 
			
		||||
 | 
			
		||||
{  Calculates the functionvalue of a polynomalfunction with n coefficients in a
 | 
			
		||||
for variable X }
 | 
			
		||||
{  Calculates the function value of a polynomial of degree n for variable x.
 | 
			
		||||
   The polynomial coefficients a are ordered from lowest to highest degree term.
 | 
			
		||||
     y = a0 + a1 x + a2 x^2 + ... + an x^n }
 | 
			
		||||
function spepol(x: ArbFloat; var a: ArbFloat; n: ArbInt): ArbFloat;
 | 
			
		||||
 | 
			
		||||
{ Calc a^b with a and b real numbers}
 | 
			
		||||
function spepow(a, b: ArbFloat): ArbFloat;
 | 
			
		||||
function spepow(a, b: ArbFloat): ArbFloat; deprecated 'Use power(a,b) in unit math.';
 | 
			
		||||
 | 
			
		||||
{ Returns sign of x (-1 for x<0, 0 for x=0 and 1 for x>0)  }
 | 
			
		||||
function spesgn(x: ArbFloat): ArbInt;
 | 
			
		||||
function spesgn(x: ArbFloat): ArbInt;  deprecated 'Use sign(x) in unit math.';
 | 
			
		||||
 | 
			
		||||
{  ArcSin(x) }
 | 
			
		||||
function spears(x: ArbFloat): ArbFloat;
 | 
			
		||||
function spears(x: ArbFloat): ArbFloat; deprecated 'Use arcsin(x) in unit math.';
 | 
			
		||||
 | 
			
		||||
{  ArcCos(x) }
 | 
			
		||||
function spearc(x: ArbFloat): ArbFloat;
 | 
			
		||||
function spearc(x: ArbFloat): ArbFloat; deprecated 'Use arccos(x) in unit math.';
 | 
			
		||||
 | 
			
		||||
{  Sinh(x) }
 | 
			
		||||
function spesih(x: ArbFloat): ArbFloat;
 | 
			
		||||
function spesih(x: ArbFloat): ArbFloat; deprecated 'Use sinh(x) in unit math.';
 | 
			
		||||
 | 
			
		||||
{  Cosh(x) }
 | 
			
		||||
function specoh(x: ArbFloat): ArbFloat;
 | 
			
		||||
function specoh(x: ArbFloat): ArbFloat; deprecated 'Use cosh(x) in unit math.';
 | 
			
		||||
 | 
			
		||||
{  Tanh(x) }
 | 
			
		||||
function spetah(x: ArbFloat): ArbFloat;
 | 
			
		||||
function spetah(x: ArbFloat): ArbFloat; deprecated 'Use tanh(x) in unit math.';
 | 
			
		||||
 | 
			
		||||
{  ArcSinH(x) }
 | 
			
		||||
function speash(x: ArbFloat): ArbFloat;
 | 
			
		||||
function speash(x: ArbFloat): ArbFloat; deprecated 'Use arcsinh(x) in unit math.';
 | 
			
		||||
 | 
			
		||||
{  ArcCosh(x) }
 | 
			
		||||
function speach(x: ArbFloat): ArbFloat;
 | 
			
		||||
function speach(x: ArbFloat): ArbFloat; deprecated 'Use arccosh(x) in unit math';
 | 
			
		||||
 | 
			
		||||
{  ArcTanH(x) }
 | 
			
		||||
function speath(x: ArbFloat): ArbFloat;
 | 
			
		||||
function speath(x: ArbFloat): ArbFloat; deprecated 'Use arctanh(x) in unit math';
 | 
			
		||||
 | 
			
		||||
{ Error numbers used within this unit:
 | 
			
		||||
 | 
			
		||||
@ -128,10 +166,22 @@ function speath(x: ArbFloat): ArbFloat;
 | 
			
		||||
  409 - function spears(s, x) is not defined for x < -1 or x > 1
 | 
			
		||||
  410 - function gammap(s, x) is not defined for s <= 0 or x < 0
 | 
			
		||||
  411 - function gammaq(s, x) is not defined for s <= 0 or x < 0
 | 
			
		||||
  412 - function beta(a, b) is not defined for a <= 0 or b <= 0
 | 
			
		||||
  413 - function betai(a, b, x) is not defined for a <= 0 or b <= 0
 | 
			
		||||
  414 - function betai(a, b, x) is not defined for x < 0 or x > 0
 | 
			
		||||
  415 - function invtdist(t, n) is not defined for t <= 0 or t >= 1 or n <= 0.
 | 
			
		||||
}
 | 
			
		||||
 | 
			
		||||
implementation
 | 
			
		||||
 | 
			
		||||
uses
 | 
			
		||||
  math, roo;
 | 
			
		||||
 | 
			
		||||
const
 | 
			
		||||
  SQRT2 = 1.4142135623730950488016887242097;   // sqrt(2)
 | 
			
		||||
  SQRT2PI = 2.506628274631000502415765284811;  // sqrt(2*pi)
 | 
			
		||||
  EXP_2 = 0.13533528323661269189399949497248;  // exp(-2)
 | 
			
		||||
 | 
			
		||||
function spebi0(x: ArbFloat): ArbFloat;
 | 
			
		||||
 | 
			
		||||
const
 | 
			
		||||
@ -896,6 +946,120 @@ begin
 | 
			
		||||
    end
 | 
			
		||||
end {speefc};
 | 
			
		||||
 | 
			
		||||
{ N(x) = 1/sqrt(2 pi) int(-INF, x, exp(t^2/2) = (1 + erf(x/sqrt(2))) / 2 }
 | 
			
		||||
function normaldist(x: ArbFloat): ArbFloat;
 | 
			
		||||
begin
 | 
			
		||||
  Result := 0.5 * (1.0 + speerf(x / SQRT2));
 | 
			
		||||
end;
 | 
			
		||||
 | 
			
		||||
function invnormaldist(y: ArbFloat): ArbFloat;
 | 
			
		||||
{ Ref.: Moshier, "Methods and programs for mathematical function" }
 | 
			
		||||
const
 | 
			
		||||
  P0: array[0..4] of ArbFloat = (
 | 
			
		||||
    -1.23916583867381258016,
 | 
			
		||||
    13.9312609387279679503,
 | 
			
		||||
    -56.6762857469070293439,
 | 
			
		||||
    98.0010754185999661536,
 | 
			
		||||
    -59.9633501014107895267);
 | 
			
		||||
  Q0: array[0..8] of ArbFloat = (
 | 
			
		||||
    -1.18331621121330003142,
 | 
			
		||||
    15.9056225126211695515,
 | 
			
		||||
    -82.0372256168333339912,
 | 
			
		||||
    200.260212380060660359,
 | 
			
		||||
    -225.462687854119370527,
 | 
			
		||||
    86.3602421390890590575,
 | 
			
		||||
    4.67627912898881538453,
 | 
			
		||||
    1.95448858338141759834,
 | 
			
		||||
    1.0);
 | 
			
		||||
  P1: array[0..8] of ArbFloat = (
 | 
			
		||||
    -8.57456785154685413611E-4,
 | 
			
		||||
    -3.50424626827848203418E-2,
 | 
			
		||||
    -1.40256079171354495875E-1,
 | 
			
		||||
    2.18663306850790267539,
 | 
			
		||||
    14.6849561928858024014,
 | 
			
		||||
    44.0805073893200834700,
 | 
			
		||||
    57.1628192246421288162,
 | 
			
		||||
    31.5251094599893866154,
 | 
			
		||||
    4.05544892305962419923);
 | 
			
		||||
  Q1: array[0..8] of Arbfloat = (
 | 
			
		||||
    -9.33259480895457427372E-4,
 | 
			
		||||
    -3.80806407691578277194E-2,
 | 
			
		||||
    -1.42182922854787788574E-1,
 | 
			
		||||
    2.50464946208309415979,
 | 
			
		||||
    15.0425385692907503408,
 | 
			
		||||
    41.3172038254672030440,
 | 
			
		||||
    45.3907635128879210584,
 | 
			
		||||
    15.7799883256466749731,
 | 
			
		||||
    1.0);
 | 
			
		||||
  P2: array[0..8] of ArbFloat = (
 | 
			
		||||
    6.23974539184983293730E-9,
 | 
			
		||||
    2.65806974686737550832E-6,
 | 
			
		||||
    3.01581553508235416007E-4,
 | 
			
		||||
    1.23716634817820021358E-2,
 | 
			
		||||
    2.01485389549179081538E-1,
 | 
			
		||||
    1.33303460815807542389,
 | 
			
		||||
    3.93881025292474443415,
 | 
			
		||||
    6.91522889068984211695,
 | 
			
		||||
    3.23774891776946035970);
 | 
			
		||||
  Q2: array[0..8] of ArbFloat = (
 | 
			
		||||
    6.79019408009981274425E-9,
 | 
			
		||||
    2.89247864745380683936E-6,
 | 
			
		||||
    3.28014464682127739104E-4,
 | 
			
		||||
    1.34204006088543189037E-2,
 | 
			
		||||
    2.16236993594496635890E-1,
 | 
			
		||||
    1.37702099489081330271,
 | 
			
		||||
    3.67983563856160859403,
 | 
			
		||||
    6.02427039364742014255,
 | 
			
		||||
    1.0);
 | 
			
		||||
var
 | 
			
		||||
  x, x0, x1: ArbFloat;
 | 
			
		||||
  yy, y2: ArbFloat;
 | 
			
		||||
  z: ArbFloat;
 | 
			
		||||
  code: Integer;
 | 
			
		||||
begin
 | 
			
		||||
  if y <= 0.0 then
 | 
			
		||||
  begin
 | 
			
		||||
    Result := -giant;
 | 
			
		||||
    exit;
 | 
			
		||||
  end;
 | 
			
		||||
 | 
			
		||||
  if y >= 1.0 then
 | 
			
		||||
  begin
 | 
			
		||||
    Result := +giant;
 | 
			
		||||
    exit;
 | 
			
		||||
  end;
 | 
			
		||||
 | 
			
		||||
  code := 1;
 | 
			
		||||
  yy := y;
 | 
			
		||||
  if yy > 1.0 - EXP_2 then begin    // EXP_2 = exp(-2)
 | 
			
		||||
    yy := 1.0 - yy;
 | 
			
		||||
    code := 0;
 | 
			
		||||
  end;
 | 
			
		||||
 | 
			
		||||
  if yy > EXP_2 then begin
 | 
			
		||||
    yy := yy - 0.5;
 | 
			
		||||
    y2 := yy * yy;
 | 
			
		||||
    x := y2 * spepol(y2, P0[0], 4) / spepol(y2, Q0[0], 8);
 | 
			
		||||
    x := (yy + yy * x) * SQRT2PI;   // SQRT2PI = sqrt(2*pi);
 | 
			
		||||
    Result := x;
 | 
			
		||||
    exit;
 | 
			
		||||
  end;
 | 
			
		||||
 | 
			
		||||
  x := sqrt(-2.0 * ln(yy));
 | 
			
		||||
  x0 := x - ln(x) / x;
 | 
			
		||||
  z := 1.0 / x;
 | 
			
		||||
 | 
			
		||||
  if x < 8.0 then
 | 
			
		||||
    x1 := z * spepol(z, P1[0], 8) / spepol(z, Q1[0], 8)
 | 
			
		||||
  else
 | 
			
		||||
    x1 := z * spepol(z, P2[0], 8) / spepol(z, Q2[0], 8);
 | 
			
		||||
 | 
			
		||||
  x := x0 - x1;
 | 
			
		||||
  if code <> 0 then
 | 
			
		||||
    x := -x;
 | 
			
		||||
  Result := x;
 | 
			
		||||
end;
 | 
			
		||||
 | 
			
		||||
function spegam(x: ArbFloat): ArbFloat;
 | 
			
		||||
const
 | 
			
		||||
 | 
			
		||||
@ -1154,6 +1318,183 @@ begin
 | 
			
		||||
    Result := gammacf(s, x);        // Use continued fraction
 | 
			
		||||
end;
 | 
			
		||||
 | 
			
		||||
{ Ref.: Moshier, "Methods and programs for mathematical functions" }
 | 
			
		||||
function invgammaq(s, y: ArbFloat): ArbFloat;
 | 
			
		||||
const
 | 
			
		||||
  NUM_IT = 30;
 | 
			
		||||
var
 | 
			
		||||
  d, y0, x0, xinit, lgm: ArbFloat;
 | 
			
		||||
  it: Integer;
 | 
			
		||||
  eps: ArbFloat;
 | 
			
		||||
begin
 | 
			
		||||
  d := 1.0 / (9 * s);
 | 
			
		||||
  y0 := invnormaldist(y);
 | 
			
		||||
  if y0 = giant then
 | 
			
		||||
    exit(0.0);
 | 
			
		||||
 | 
			
		||||
  y0 := 1.0 - d - y0 * sqrt(d);
 | 
			
		||||
  x0 := s * y0 * y0 * y0;
 | 
			
		||||
  xinit := x0;
 | 
			
		||||
  lgm := spelga(s);
 | 
			
		||||
  eps := 2.0 * MachEps;
 | 
			
		||||
 | 
			
		||||
  for it := 1 to NUM_IT do
 | 
			
		||||
  begin
 | 
			
		||||
    if (x0 <= 0.0) then   // underflow
 | 
			
		||||
      exit(0.0);
 | 
			
		||||
    y0 := gammaq(s, x0);
 | 
			
		||||
    d := (s - 1.0) * ln(x0) - x0 - lgm;
 | 
			
		||||
    if d < -lnGiant then  // underflow
 | 
			
		||||
      break;
 | 
			
		||||
    d := -exp(d);
 | 
			
		||||
    if d = 0.0 then
 | 
			
		||||
      break;
 | 
			
		||||
    d := (y0 - y) / d;
 | 
			
		||||
    x0 := x0 - d;
 | 
			
		||||
    if it <= 3 then
 | 
			
		||||
      continue;
 | 
			
		||||
    if abs(d / x0) < eps then
 | 
			
		||||
      break;
 | 
			
		||||
  end;
 | 
			
		||||
  Result := x0;
 | 
			
		||||
end;
 | 
			
		||||
 | 
			
		||||
{ Calculates the complete beta function based on its property that
 | 
			
		||||
    beta(a, b) = gamma(a) * gamma(b) / gamma(a+b)
 | 
			
		||||
  https://en.wikipedia.org/wiki/Beta_function }
 | 
			
		||||
function beta(a, b: ArbFloat): ArbFloat;
 | 
			
		||||
begin
 | 
			
		||||
  if (a <= 0) or (b <= 0) then
 | 
			
		||||
    RunError(412);
 | 
			
		||||
  Result := exp(spelga(a) + spelga(b) - spelga(a + b));
 | 
			
		||||
end;
 | 
			
		||||
 | 
			
		||||
{ Calculates the continued fraction of the incomplete beta function.
 | 
			
		||||
  Ref: https://www.encyclopediaofmath.org/index.php/Incomplete_beta-function }
 | 
			
		||||
function betaicf(a, b, x: ArbFloat): Arbfloat;
 | 
			
		||||
 | 
			
		||||
  function funca(i: Integer): ArbFloat;
 | 
			
		||||
  begin
 | 
			
		||||
    if i = 0 then Result := 0.0 else Result := 1.0;
 | 
			
		||||
  end;
 | 
			
		||||
 | 
			
		||||
  function funcb(i: Integer): ArbFloat;
 | 
			
		||||
  var
 | 
			
		||||
    am: ArbFloat;
 | 
			
		||||
    amm: ArbFloat;
 | 
			
		||||
    m: Integer;
 | 
			
		||||
  begin
 | 
			
		||||
    if i = 1 then
 | 
			
		||||
      Result := 1.0
 | 
			
		||||
    else begin
 | 
			
		||||
      m := (i-1) div 2;
 | 
			
		||||
      am := a + m;
 | 
			
		||||
      amm := am + m;
 | 
			
		||||
      if odd(i) then
 | 
			
		||||
        Result := m * (b - m) * x / ((amm - 1) * amm)
 | 
			
		||||
      else
 | 
			
		||||
        Result := -am * (am + b) * x / (amm * (amm + 1));
 | 
			
		||||
    end;
 | 
			
		||||
  end;
 | 
			
		||||
 | 
			
		||||
const
 | 
			
		||||
  MAX_IT = 100;
 | 
			
		||||
  EPS = 1E-7;
 | 
			
		||||
begin
 | 
			
		||||
  Result := CalcCF(@funca, @funcb, MAX_IT, EPS);
 | 
			
		||||
end;
 | 
			
		||||
 | 
			
		||||
function betai(a, b, x: ArbFloat): ArbFloat;
 | 
			
		||||
var
 | 
			
		||||
  factor: ArbFloat;
 | 
			
		||||
begin
 | 
			
		||||
  // Check for invalid arguments
 | 
			
		||||
  if (a <= 0) or (b <= 0) then
 | 
			
		||||
    RunError(413);
 | 
			
		||||
  if (x < 0) or (x > 1) then
 | 
			
		||||
    RunError(414);
 | 
			
		||||
 | 
			
		||||
  if (x = 0) or (x = 1) then
 | 
			
		||||
    factor := 0
 | 
			
		||||
  else
 | 
			
		||||
    factor := exp(a * ln(x) + b * ln(1.0 - x) + spelga(a + b) - spelga(a) - spelga(b));
 | 
			
		||||
 | 
			
		||||
  // The continued fraction expansion converges quickly only for
 | 
			
		||||
  //    x < (a + 1) / (a + b + 2)
 | 
			
		||||
  // For the other case, we apply the relation
 | 
			
		||||
  //    beta(a, b, x) = 1 - beta(b, a, 1-x)
 | 
			
		||||
  if x < (a + 1) / (a + b + 2) then
 | 
			
		||||
    Result := factor * betaicf(a, b, x) / a
 | 
			
		||||
  else
 | 
			
		||||
    Result := 1.0 - factor * betaicf(b, a, 1.0 - x) / b;
 | 
			
		||||
end;
 | 
			
		||||
 | 
			
		||||
{ Inverse of the incomplete beta function }
 | 
			
		||||
function invbetai(a, b, y: ArbFloat; eps: ArbFloat = 0.0): ArbFloat;
 | 
			
		||||
 | 
			
		||||
  function _betai(x: ArbFloat): ArbFloat;
 | 
			
		||||
  begin
 | 
			
		||||
    Result := betai(a, b, x) - y;
 | 
			
		||||
  end;
 | 
			
		||||
 | 
			
		||||
var
 | 
			
		||||
  term: ArbInt = 0;
 | 
			
		||||
begin
 | 
			
		||||
  if eps = 0.0 then
 | 
			
		||||
    eps := MachEps;
 | 
			
		||||
  roof1rn(@_betai, 0, 1, eps, eps, Result, term);
 | 
			
		||||
  if term = 3 then
 | 
			
		||||
    Result := NaN;
 | 
			
		||||
end;
 | 
			
		||||
 | 
			
		||||
function chi2dist(x: ArbFloat; n: ArbInt): ArbFloat;
 | 
			
		||||
begin
 | 
			
		||||
  Result := gammaQ(0.5*n, 0.5*x);
 | 
			
		||||
end;
 | 
			
		||||
 | 
			
		||||
function invchi2dist(y: Arbfloat; n: ArbInt): ArbFloat;
 | 
			
		||||
begin
 | 
			
		||||
  Result := 2.0 * invgammaQ(n/2, y);
 | 
			
		||||
//  Result := 2.0 * invgammaQ_alglib(n/2, y);
 | 
			
		||||
end;
 | 
			
		||||
 | 
			
		||||
function tdist(t: ArbFloat; n: ArbInt; Tails: TNumTails): ArbFloat;
 | 
			
		||||
begin
 | 
			
		||||
  Result := betai(0.5*n, 0.5, n/(n+t*t));
 | 
			
		||||
  if Tails = 1 then Result := Result * 0.5;
 | 
			
		||||
end;
 | 
			
		||||
 | 
			
		||||
function invtdist(y: ArbFloat; n: ArbInt; Tails: TNumTails;
 | 
			
		||||
  eps: ArbFloat = 0.0): ArbFloat;
 | 
			
		||||
var
 | 
			
		||||
  w: ArbFloat;
 | 
			
		||||
begin
 | 
			
		||||
  if (n <= 0) or (y <= 0) or (y >= 1) then
 | 
			
		||||
    RunError(415);
 | 
			
		||||
 | 
			
		||||
  if Tails = 2 then y := y * 0.5;
 | 
			
		||||
  w := invbetai(0.5*n, 0.5, 2*y, eps);
 | 
			
		||||
  Result := sqrt(n/w - n);
 | 
			
		||||
end;
 | 
			
		||||
 | 
			
		||||
// Calculates the F distribution with n1 and n2 degrees of freedom in the
 | 
			
		||||
// numerator and denominator, respectively
 | 
			
		||||
function Fdist(F: ArbFloat; n1, n2: ArbInt): ArbFloat;
 | 
			
		||||
begin
 | 
			
		||||
  Result := betai(n2*0.5, n1*0.5, n2 / (n2 + n1*F));
 | 
			
		||||
end;
 | 
			
		||||
 | 
			
		||||
// Calculates the inverse of the F distribution
 | 
			
		||||
// Ref. Moshier, "Methods and programs for mathematical functions"
 | 
			
		||||
function invFdist(p: ArbFloat; n1, n2: ArbInt; eps: ArbFloat = 0.0): ArbFloat;
 | 
			
		||||
var
 | 
			
		||||
  s: ArbFloat;
 | 
			
		||||
begin
 | 
			
		||||
  if eps = 0.0 then eps := machEps;
 | 
			
		||||
  s := invbetai(n2*0.5, n1*0.5, p, eps);
 | 
			
		||||
  Result := n2 * (1-s) / (n1 * s);
 | 
			
		||||
end;
 | 
			
		||||
 | 
			
		||||
function spepol(x: ArbFloat; var a: ArbFloat; n: ArbInt): ArbFloat;
 | 
			
		||||
var   pa : ^arfloat0;
 | 
			
		||||
       i : ArbInt;
 | 
			
		||||
@ -1419,6 +1760,7 @@ end; {speath}
 | 
			
		||||
var exitsave : pointer;
 | 
			
		||||
 | 
			
		||||
procedure MyExit;
 | 
			
		||||
{
 | 
			
		||||
const ErrorS : array[400..408,1..6] of char =
 | 
			
		||||
     ('spepow',
 | 
			
		||||
      'spebk0',
 | 
			
		||||
@ -1428,7 +1770,7 @@ const ErrorS : array[400..408,1..6] of char =
 | 
			
		||||
      'speach',
 | 
			
		||||
      'speath',
 | 
			
		||||
      'spegam',
 | 
			
		||||
      'spelga');
 | 
			
		||||
      'spelga'); }
 | 
			
		||||
 | 
			
		||||
//var ErrFil : text;
 | 
			
		||||
 | 
			
		||||
@ -1436,7 +1778,7 @@ begin
 | 
			
		||||
     ExitProc := ExitSave;
 | 
			
		||||
//     Assign(ErrFil, 'CON');
 | 
			
		||||
//    ReWrite(ErrFil);
 | 
			
		||||
     if (ExitCode>=400) AND (ExitCode<=408) then
 | 
			
		||||
     if (ExitCode>=400) AND (ExitCode<=415) then
 | 
			
		||||
       begin
 | 
			
		||||
//         write(ErrFil, 'critical error in ', ErrorS[ExitCode]);
 | 
			
		||||
         ExitCode := 201
 | 
			
		||||
 | 
			
		||||
@ -35,6 +35,9 @@ Also some stuff had to be added to get ipf running (vector object and
 | 
			
		||||
complex.inp and scale methods)
 | 
			
		||||
 }
 | 
			
		||||
 | 
			
		||||
{$mode objfpc}{$H+}
 | 
			
		||||
{$modeswitch nestedprocvars}
 | 
			
		||||
 | 
			
		||||
unit typ;
 | 
			
		||||
 | 
			
		||||
{$I DIRECT.INC}                 {Contains "global" compilerswitches which
 | 
			
		||||
@ -182,6 +185,7 @@ type
 | 
			
		||||
 | 
			
		||||
     {Standard Functions used in NumLib}
 | 
			
		||||
     rfunc1r    = Function(x : ArbFloat): ArbFloat;
 | 
			
		||||
     rfunc1rn   = Function(x : ArbFloat): ArbFloat is nested;
 | 
			
		||||
     rfunc2r    = Function(x, y : ArbFloat): ArbFloat;
 | 
			
		||||
 | 
			
		||||
     {Complex version}
 | 
			
		||||
 | 
			
		||||
		Loading…
	
		Reference in New Issue
	
	Block a user