IMPLEMENTATION MODULE SpezFunkt4;
(*------------------------------------------------------------------------*)
(* Verschiedene Dichte- und Wahrscheinlichkeitsfunktionen *)
(* Various density and probability functions *)
(*------------------------------------------------------------------------*)
(* Letzte Veraenderung: *)
(* *)
(* 27.08.16, MRi: Erstellen der ersten Version *)
(* 03.09.15, MRi: NormCDF,DisLandau & DenLandau eingefuegt *)
(*------------------------------------------------------------------------*)
(* Offene Punkte: *)
(* *)
(* - Ausfuerhlicheres testen der Routinen *)
(*------------------------------------------------------------------------*)
(* Based on work found in dtmath086 Pascal library with is in parts *)
(* based on Cephes library (http://www.moshier.net) *)
(* Translated and modified by M.Riedl, August 2016 *)
(*------------------------------------------------------------------------*)
(* Testroutine vorhanden (TestSpF4 und TestGam.p) *)
(*------------------------------------------------------------------------*)
(* Implementation : Michael Riedl *)
(* Licence : GNU Lesser General Public License (LGPL) *)
(*------------------------------------------------------------------------*)
(* $Id: SpezFunkt4.mod,v 1.1 2016/09/03 22:08:48 mriedl Exp mriedl $ *)
FROM TestReal IMPORT Real8NaNsignaled,Real8INF;
FROM LongMath IMPORT pi,sqrt,ln,exp,power;
FROM SpezFunkt2 IMPORT dLnGamma,ErrorF,CompErrorF,IGamma,IBeta;
FROM LMathLib IMPORT MinLog,sqrt2;
VAR NAN,INF : LONGREAL;
PROCEDURE DExpo(A,X : LONGREAL) : LONGREAL;
VAR Y : LONGREAL;
BEGIN
IF (A <= 0.0) OR (X < 0.0) THEN
RETURN NAN;
END;
Y := - A*X;
IF (Y < MinLog) THEN
RETURN 0.0;
END;
RETURN A*exp(Y);
END DExpo;
PROCEDURE FExpo(A,X : LONGREAL) : LONGREAL;
VAR Y : LONGREAL;
BEGIN
IF (A <= 0.0) OR (X < 0.0) THEN
RETURN NAN;;
END;
Y := - A * X;
IF (Y < MinLog) THEN
RETURN 1.0;
END;
RETURN 1.0 - exp(Y);
END FExpo;
PROCEDURE DBeta(A,B : LONGREAL;
X : LONGREAL) : LONGREAL;
VAR L : LONGREAL;
BEGIN
IF (A <= 0.0) OR (B <= 0.0) OR (X < 0.0) OR (X > 1.0) THEN
RETURN NAN;
END;
IF (X = 0.0) THEN
IF (A < 1.0) THEN
RETURN INF;
ELSE
RETURN 0.0;
END;
END;
IF (X = 1.0) THEN
IF (B < 1.0) THEN
RETURN INF;
ELSE
RETURN 0.0;
END;
END;
L := dLnGamma(A + B) - dLnGamma(A) - dLnGamma(B) +
(A - 1.0)*ln(X) + (B - 1.0)*ln(1.0 - X);
IF (L < MinLog) THEN
RETURN 0.0;
ELSE
RETURN exp(L);
END;
END DBeta;
PROCEDURE FBeta(A,B : LONGREAL;
X : LONGREAL) : LONGREAL;
BEGIN
RETURN IBeta(A,B,X);
END FBeta;
PROCEDURE FBinom(N : INTEGER;
P : LONGREAL;
K : INTEGER) : LONGREAL;
BEGIN
IF (P < 0.0) OR (P > 1.0) OR (N <= 0) OR (N < K) THEN
RETURN NAN;
ELSIF (K = 0) THEN
RETURN power(1.0 - P,VAL(LONGREAL,N));
ELSIF (K = N) THEN
RETURN 1.0;
ELSE
RETURN 1.0 - IBeta(VAL(LONGREAL,K + 1),VAL(LONGREAL,N - K),P);
END;
END FBinom;
PROCEDURE DStudent(Nu : INTEGER;
X : LONGREAL) : LONGREAL;
VAR L,P,Q : LONGREAL;
BEGIN
IF (Nu < 1) THEN
RETURN NAN;
END;
P := 0.5*(VAL(LONGREAL,Nu + 1));
Q := 0.5* VAL(LONGREAL,Nu);
L := dLnGamma(P) - dLnGamma(Q) - 0.5*ln(VAL(LONGREAL,Nu)*pi) -
P*ln(1.0 + X*X / VAL(LONGREAL,Nu));
IF (L < MinLog) THEN
RETURN 0.0;
ELSE
RETURN exp(L);
END;
END DStudent;
PROCEDURE FStudent(Nu : INTEGER;
X : LONGREAL) : LONGREAL;
VAR F : LONGREAL;
BEGIN
IF (Nu < 1) THEN
RETURN NAN;
ELSIF (X = 0.0) THEN
RETURN 0.5;
ELSE
F := 0.5*IBeta(0.5*VAL(LONGREAL,Nu),0.5,VAL(LONGREAL,Nu) /
(VAL(LONGREAL,Nu) + X*X));
IF (X < 0.0) THEN RETURN F ELSE RETURN 1.0 - F; END;
END;
END FStudent;
PROCEDURE PStudent(Nu : INTEGER;
X : LONGREAL) : LONGREAL;
BEGIN
IF (Nu < 1) THEN
RETURN NAN;
ELSE
RETURN IBeta(0.5*VAL(LONGREAL,Nu),0.5,VAL(LONGREAL,Nu) /
(VAL(LONGREAL,Nu) + X*X));
END;
END PStudent;
PROCEDURE DSnedecor(Nu1,Nu2 : INTEGER;
X : LONGREAL) : LONGREAL;
VAR L,P1,P2,R,S : LONGREAL;
BEGIN
IF (Nu1 < 1) OR (Nu2 < 1) OR (X <= 0.0) THEN
RETURN NAN;;
END;
R := VAL(LONGREAL,Nu1) / VAL(LONGREAL,Nu2);
P1 := 0.5*VAL(LONGREAL,Nu1);
P2 := 0.5*VAL(LONGREAL,Nu2);
S := P1 + P2;
L := dLnGamma(S) - dLnGamma(P1) - dLnGamma(P2) +
P1 * ln(R) + (P1 - 1.0) * ln(X) - S * ln(1.0 + R * X);
IF (L < MinLog) THEN
RETURN 0.0;
ELSE
RETURN exp(L);
END;
END DSnedecor;
PROCEDURE FSnedecor(Nu1,Nu2 : INTEGER;
X : LONGREAL) : LONGREAL;
VAR nu1,nu2 : LONGREAL;
BEGIN
IF (Nu1 < 1) OR (Nu2 < 1) OR (X <= 0.0) THEN
RETURN NAN;
ELSE
nu1:=VAL(LONGREAL,Nu1);
nu2:=VAL(LONGREAL,Nu2);
RETURN 1.0 - IBeta(0.5*nu2,0.5*nu1,nu2 / (nu2 + nu1*X));
END;
END FSnedecor;
PROCEDURE PSnedecor(Nu1,Nu2 : INTEGER;
X : LONGREAL) : LONGREAL;
VAR nu1,nu2 : LONGREAL;
BEGIN
IF (Nu1 < 1) OR (Nu2 < 1) OR (X <= 0.0) THEN
RETURN NAN;
ELSE
nu1:=VAL(LONGREAL,Nu1);
nu2:=VAL(LONGREAL,Nu2);
RETURN IBeta(0.5*nu2,0.5*nu1,nu2 / (nu2 + nu1*X));
END;
END PSnedecor;
PROCEDURE DGamma(A,B : LONGREAL;
X : LONGREAL) : LONGREAL;
VAR L : LONGREAL;
BEGIN
IF (A <= 0.0) OR (B <= 0.0) OR (X < 0.0) THEN
RETURN NAN;
END;
IF (X = 0.0) THEN
IF (A < 1.0) THEN
RETURN INF;
ELSIF (A = 1.0) THEN
RETURN B;
ELSE
RETURN 0.0;
END;
END;
L := A * ln(B) - dLnGamma(A) + (A - 1.0)*ln(X) - B*X;
IF (L < MinLog) THEN
RETURN 0.0;
ELSE
RETURN exp(L);
END;
END DGamma;
PROCEDURE FGamma(A,B : LONGREAL;
X : LONGREAL) : LONGREAL;
BEGIN
RETURN IGamma(A, B*X);
END FGamma;
PROCEDURE PPoisson(Mu : LONGREAL;
K : INTEGER) : LONGREAL;
VAR P : LONGREAL;
i : INTEGER;
BEGIN
IF (Mu <= 0.0) OR (K < 0) THEN
RETURN NAN;
ELSIF (VAL(LONGREAL,- Mu) < MinLog) THEN
RETURN 0.0;
ELSIF (K = 0) THEN
RETURN exp(- Mu);
ELSE
P := Mu;
FOR i:=2 TO K DO (* P = Mu^K / K! *)
P:=P * (Mu / VAL(LONGREAL,i));
END;
RETURN exp(- Mu)*P;
END;
END PPoisson;
PROCEDURE FPoisson(Mu : LONGREAL;
K : INTEGER) : LONGREAL;
BEGIN
IF (Mu <= 0.0) OR (K < 0) THEN
RETURN NAN;
ELSIF (K = 0) THEN
IF ((- Mu) < MinLog) THEN
RETURN 0.0;
ELSE
RETURN exp(- Mu);
END;
ELSE
RETURN 1.0 - IGamma(VAL(LONGREAL,K + 1), Mu);
END;
END FPoisson;
PROCEDURE FNorm(X : LONGREAL) : LONGREAL;
BEGIN
RETURN 0.5*(1.0 + ErrorF(0.5*sqrt2*X));
END FNorm;
PROCEDURE PNorm(X : LONGREAL) : LONGREAL;
VAR A : LONGREAL;
BEGIN
A := ABS(X);
IF (A = 0.0) THEN
RETURN 1.0;
ELSIF (A < 1.0) THEN
RETURN 1.0 - ErrorF(0.5*sqrt2*A)
ELSE
RETURN CompErrorF(0.5*sqrt2*A);
END;
END PNorm;
PROCEDURE NormCDF(Z : LONGREAL) : LONGREAL;
(*----------------------------------------------------------------*)
(* 30.03.86, AMi : Last revision by Alan Miller *)
(* 01.09.16, MRi : Translation to Modula-2 *)
(*----------------------------------------------------------------*)
CONST P0 = 2.202068679123761E+02; P1 = 2.212135961699311E+02;
P2 = 1.120792914978709E+02; P3 = 3.391286607838300E+01;
P4 = 6.373962203531650E+00; P5 = 7.003830644436881E-01;
P6 = 3.526249659989109E-02;
Q0 = 4.404137358247522E+02; Q1 = 7.938265125199484E+02;
Q2 = 6.373336333788311E+02; Q3 = 2.965642487796737E+02;
Q4 = 8.678073220294608E+01; Q5 = 1.606417757920695E+01;
Q6 = 1.755667163182642E+00; Q7 = 8.838834764831844E-02;
Cutoff = 7.071067811865475;
RootPi = 2.506628274631001;
VAR p,expntl,zabs : LONGREAL;
BEGIN
zabs := ABS(Z);
IF (zabs > 37.0) THEN
p := 0.0
ELSE
expntl := exp(-0.5*zabs*zabs);
IF (zabs < Cutoff) THEN (* |z| < cutoff = 10/sqrt(2) *)
p:=expntl*((((((
P6*zabs + P5)*zabs + P4)*zabs + P3)*
zabs + P2)*zabs + P1)*zabs + P0) /
(((((((Q7*zabs + Q6)*zabs + Q5)*zabs + Q4)*
zabs + Q3)*zabs + Q2)*zabs + Q1)*zabs + Q0);
ELSE (* |z| >= cutoff. *)
p := expntl / (zabs + 1.0 / (zabs + 2.00 / (zabs + 3.0 /
(zabs + 4.0 / (zabs + 0.65))))) / RootPi;
END;
END;
IF (Z > 0.0) THEN
p := 1.0 - p
END;
RETURN p;
END NormCDF;
PROCEDURE DKhi2(Nu : INTEGER;
X : LONGREAL) : LONGREAL;
BEGIN
RETURN DGamma(0.5*VAL(LONGREAL,Nu),0.5,X)
END DKhi2;
PROCEDURE FKhi2(Nu : INTEGER;
X : LONGREAL) : LONGREAL;
BEGIN
IF (Nu < 1) OR (X <= 0.0) THEN
RETURN NAN;
ELSE
RETURN IGamma(0.5*VAL(LONGREAL,Nu), 0.5*X);
END;
END FKhi2;
PROCEDURE PKhi2(Nu : INTEGER;
X : LONGREAL) : LONGREAL;
BEGIN
IF (Nu < 1) OR (X <= 0.0) THEN
RETURN NAN;
ELSE
RETURN 1.0 - IGamma(0.5*VAL(LONGREAL,Nu), 0.5*X);
END;
END PKhi2;
PROCEDURE DisLandau(X : LONGREAL) : LONGREAL;
TYPE V0t4 = ARRAY [0..4] OF LONGREAL;
V0t3 = ARRAY [0..3] OF LONGREAL;
V1t3 = ARRAY [1..3] OF LONGREAL;
V1t2 = ARRAY [1..2] OF LONGREAL;
CONST P1 = V0t4{ 2.514091491E-01, -6.250580444E-02,
1.458381230E-02, -2.108817737E-03,
7.411247290E-04 };
P2 = V0t3{ 2.868328584E-01, 3.564363231E-01,
1.523518695E-01, 2.251304883E-02 };
P3 = V0t3{ 2.868329066E-01, 3.003828436E-01,
9.950951941E-02, 8.733827185E-03 };
P4 = V0t3{ 1.000351630E+00, 4.503592498E+00,
1.085883880E+01, 7.536052269E+00 };
P5 = V0t3{ 1.000006517E+00, 4.909414111E+01,
8.505544753E+01, 1.532153455E+02 };
P6 = V0t3{ 1.000000983E+00, 1.329868456E+02,
9.162149244E+02, -9.605054274E+02 };
Q1 = V0t4{ 1.000000000E+00, -5.571175625E-03,
6.225310236E-02, -3.137378427E-03,
1.931496439E-03 };
Q2 = V0t3{ 1.000000000E+00, 6.191136137E-01,
1.720721448E-01, 2.278594771E-02 };
Q3 = V0t3{ 1.000000000E+00, 4.237190502E-01,
1.095631512E-01, 8.693851567E-03 };
Q4 = V0t3{ 1.000000000E+00, 5.539969678E+00,
1.933581111E+01, 2.721321508E+01 };
Q5 = V0t3{ 1.000000000E+00, 5.009928881E+01,
1.399819104E+02, 4.200002909E+02 };
Q6 = V0t3{ 1.000000000E+00, 1.339887843E+02,
1.055990413E+03, 5.532224619E+02 };
A1 = V1t3{-4.583333333E-01, 6.675347222E-01,
-1.641741416E+00 };
A2 = V1t3{ 1.000000000E+00, -4.227843351E-01,
-2.043403138E+00 };
VAR u,v,Res : LONGREAL;
BEGIN
v:=X;
IF (v < -5.5) THEN
u := exp(v + 1.0);
Res := 0.3989422803*exp( -1.0 / u)*sqrt(u)*
(1.0 + (A1[1] + (A1[2] + A1[3]*u)*u)*u);
ELSIF (v < -1.0) THEN
u := exp( -(v + 1.0));
Res := (exp( - u)/sqrt(u))*
(P1[0] + (P1[1] + (P1[2] + (P1[3] + P1[4]*v)*v)*v)*v) /
(Q1[0] + (Q1[1] + (Q1[2] + (Q1[3] + Q1[4]*v)*v)*v)*v);
ELSIF (v < 1.0) THEN
Res := (P2[0] + (P2[1] + (P2[2] + P2[3]*v)*v)*v) /
(Q2[0] + (Q2[1] + (Q2[2] + Q2[3]*v)*v)*v);
ELSIF (v < 4.0) THEN
Res := (P3[0] + (P3[1] + (P3[2] + P3[3]*v)*v)*v) /
(Q3[0] + (Q3[1] + (Q3[2] + Q3[3]*v)*v)*v);
ELSIF (v < 12.0) THEN
u := 1.0 / v;
Res := (P4[0] + (P4[1] + (P4[2] + P4[3]*u)*u)*u) /
(Q4[0] + (Q4[1] + (Q4[2] + Q4[3]*u)*u)*u);
ELSIF (v < 50.0) THEN
u := 1.0 / v;
Res := (P5[0] + (P5[1] + (P5[2] + P5[3]*u)*u)*u) /
(Q5[0] + (Q5[1] + (Q5[2] + Q5[3]*u)*u)*u);
ELSIF (v < 300.0) THEN
u := 1.0 / v;
Res := (P6[0] + (P6[1] + (P6[2] + P6[3]*u)*u)*u) /
(Q6[0] + (Q6[1] + (Q6[2] + Q6[3]*u)*u)*u);
ELSE
u := 1.0 / (v - v*ln(v) / (v + 1.0));
Res := 1.0 - (A2[1] + (A2[2] + A2[3]*u)*u)*u;
END;
RETURN Res;
END DisLandau;
PROCEDURE DenLandau(X : LONGREAL) : LONGREAL;
TYPE V0t4 = ARRAY [0..4] OF LONGREAL;
V1t3 = ARRAY [1..3] OF LONGREAL;
V1t2 = ARRAY [1..2] OF LONGREAL;
CONST A1 = V1t3{ 4.166666667E-02, -1.996527778E-02,
2.709538966E-02 };
A2 = V1t2{-1.845568670E+00, -4.284640743E+00 };
P1 = V0t4{ 4.259894875E-01, -1.249762550E-01,
3.984243700E-02, -6.298287635E-03,
1.511162253E-03 };
P2 = V0t4{ 1.788541609E-01, 1.173957403E-01,
1.488850518E-02, -1.394989411E-03,
1.283617211E-04 };
P3 = V0t4{ 1.788544503E-01, 9.359161662E-02,
6.325387654E-03, 6.611667319E-05,
-2.031049101E-06 };
P4 = V0t4{ 9.874054407E-01, 1.186723273E+02,
8.492794360E+02, -7.437792444E+02,
4.270262186E+02 };
P5 = V0t4{ 1.003675074E+00, 1.675702434E+02,
4.789711289E+03, 2.121786767E+04,
-2.232494910E+04 };
P6 = V0t4{ 1.000827619E+00, 6.649143136E+02,
6.297292665E+04, 4.755546998E+05,
-5.743609109E+06 };
Q1 = V0t4{ 1.000000000E+00, -3.388260629E-01,
9.594393323E-02, -1.608042283E-02,
3.778942063E-03 };
Q2 = V0t4{ 1.000000000E+00, 7.428795082E-01,
3.153932961E-01, 6.694219548E-02,
8.790609714E-03 };
Q3 = V0t4{ 1.000000000E+00, 6.097809921E-01,
2.560616665E-01, 4.746722384E-02,
6.957301675E-03 };
Q4 = V0t4{ 1.000000000E+00, 1.068615961E+02,
3.376496214E+02, 2.016712389E+03,
1.597063511E+03 };
Q5 = V0t4{ 1.000000000E+00, 1.569424537E+02,
3.745310488E+03, 9.834698876E+03,
6.692428357E+04 };
Q6 = V0t4{ 1.000000000E+00, 6.514101098E+02,
5.697473333E+04, 1.659174725E+05,
-2.815759939E+06 };
VAR v,u,Res : LONGREAL;
BEGIN
v:=X;
IF (v < -5.5) THEN
u := exp(v + 1.0);
Res := 0.3989422803*(exp( -1.0 / u) / sqrt(u))*
(1.0 + (A1[1] + (A1[2] + A1[3]*u)*u)*u)
ELSIF (v < -1.0) THEN
u := exp(-v - 1.0);
Res := exp( - u)*sqrt(u)*
(P1[0] + (P1[1] + (P1[2] + (P1[3] + P1[4]*v)*v)*v)*v) /
(Q1[0] + (Q1[1] + (Q1[2] + (Q1[3] + Q1[4]*v)*v)*v)*v);
ELSIF (v < 1.0) THEN
Res := (P2[0] + (P2[1] + (P2[2] + (P2[3] + P2[4]*v)*v)*v)*v) /
(Q2[0] + (Q2[1] + (Q2[2] + (Q2[3] + Q2[4]*v)*v)*v)*v);
ELSIF (v < 5.0) THEN
Res := (P3[0] + (P3[1] + (P3[2] + (P3[3] + P3[4]*v)*v)*v)*v) /
(Q3[0] + (Q3[1] + (Q3[2] + (Q3[3] + Q3[4]*v)*v)*v)*v);
ELSIF (v < 12.0) THEN
u := 1.0 / v;
Res := u*u*(P4[0] + (P4[1] + (P4[2] + (P4[3] + P4[4]*u)*u)*u)*u) /
(Q4[0] + (Q4[1] + (Q4[2] + (Q4[3] + Q4[4]*u)*u)*u)*u)
ELSIF (v < 50.0) THEN
u := 1.0 / v;
Res := u*u*(P5[0] + (P5[1] + (P5[2] + (P5[3] + P5[4]*u)*u)*u)*u) /
(Q5[0] + (Q5[1] + (Q5[2] + (Q5[3] + Q5[4]*u)*u)*u)*u);
ELSIF (v < 300.0) THEN
u := 1.0 / v;
Res := u*u*(P6[0] + (P6[1] + (P6[2] + (P6[3] + P6[4]*u)*u)*u)*u) /
(Q6[0] + (Q6[1] + (Q6[2] + (Q6[3] + Q6[4]*u)*u)*u)*u);
ELSE
u := 1.0 / (v - v*ln(v) / (v + 1.0));
Res := u*u*(1.0 + (A2[1] + A2[2]*u)*u);
END;
RETURN Res;
END DenLandau;
BEGIN
NAN := Real8NaNsignaled();
INF := Real8INF();
END SpezFunkt4.