From 2a77e690b486f4e026aa76da716a2e83fc5f00e7 Mon Sep 17 00:00:00 2001 From: marco Date: Sun, 23 Apr 2017 16:17:52 +0000 Subject: [PATCH] * 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 - --- packages/numlib/src/roo.pas | 22 ++- packages/numlib/src/spe.pas | 374 ++++++++++++++++++++++++++++++++++-- packages/numlib/src/typ.pas | 4 + 3 files changed, 381 insertions(+), 19 deletions(-) diff --git a/packages/numlib/src/roo.pas b/packages/numlib/src/roo.pas index 39da1be41f..fa6ca1c4f9 100644 --- a/packages/numlib/src/roo.pas +++ b/packages/numlib/src/roo.pas @@ -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; diff --git a/packages/numlib/src/spe.pas b/packages/numlib/src/spe.pas index 3dda7c19b9..b4201f46cb 100644 --- a/packages/numlib/src/spe.pas +++ b/packages/numlib/src/spe.pas @@ -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 diff --git a/packages/numlib/src/typ.pas b/packages/numlib/src/typ.pas index e7501cf4ff..68fea7bc6d 100644 --- a/packages/numlib/src/typ.pas +++ b/packages/numlib/src/typ.pas @@ -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}