StatLib.mod.m2pp
514 lines (438 with data), 18.5 kB
IMPLEMENTATION MODULE StatLib;
(*========================================================================*)
(* HINWEIS : Bitte nur die Datei StatLib.mod.m2pp editieren *)
(*========================================================================*)
(* Es sind 2 Versionen enthalten die mit *)
(* *)
(* m2pp -D __{Parameter}__ < StatLib.mod.m2pp > StatLib.mod *)
(* *)
(* mit Parameter = {XDS|GM2} erzeugt werden koennen. *)
(* *)
(* XDS : Der Export von Objekten aus lokalen Modulen wird mit *)
(* "IMPORT" angezeigt - das scheint nicht ganz "regelkonform" *)
(* GM2 : Normales "EXPORT" wird genutzt (ELSE-Fall). *)
(* *)
(* ansonsten gibt es keine Aenderungen am Quellcode *)
(*------------------------------------------------------------------------*)
(* Modul fuer einige Statistikfunktionen. *)
(* Modul for some statistical functions. *)
(*------------------------------------------------------------------------*)
(* Letzte Bearbeitung: *)
(* *)
(* 94/95, MRi: Erstellen der ersten Versionen *)
(* 26.10.96, MRi: Durchsicht *)
(* 16.09.15, MRi: Durchsicht und Bereinigung *)
(* 30.09.15, MRi: Einfuegen von Fisher und Student *)
(* 01.10.15, MRi: Fisher mit FORTRAN Routine geprueft - ist in Ordnung *)
(* 02.10.15, MRi: Huelle fur Fisher (FisherM2) fuer Fortran eingefuegt *)
(* Abfragen in Q eingefuegt so dass p \in (0,1) liegen *)
(* kann. Prozedur Q0 geloescht und in InvProbInt umbenannt *)
(* sowie StudentT und Wilxocon entsprechend angepasst. *)
(* 05.10.15, MRi: InvNorm eingefuegt *)
(* 07.10.15, MRi: MODULE InverseNorm & WilcoxonDaten zum Kapseln der *)
(* entsprechenden Daten eingefuegt. *)
(* WilcoxonDaten ist komplett neu, Parameterliste von *)
(* Wilcoxon um Ua,Uxy,NInv,iFehl erweitert *)
(*------------------------------------------------------------------------*)
(* Offene Punkte *)
(* *)
(* - In der Prozedure Wilcoxon muessen folgende Punkte geklaert werden *)
(* = Durchsicht ob Berechung der kritischen Werte fuer grosse nx,ny *)
(* stabil ist *)
(* = Ist die Beruecksichtigung der zwei Werte fuer die Inversionen gut *)
(*------------------------------------------------------------------------*)
(* Implementation : Michael Riedl *)
(* Licence : GNU Lesser General Public License (LGPL) *)
(*------------------------------------------------------------------------*)
(* $Id: StatLib.mod.m2pp,v 1.5 2018/01/16 09:19:51 mriedl Exp mriedl $ *)
FROM Errors IMPORT ErrOut,Fehlerflag;
IMPORT Errors;
FROM LongMath IMPORT pi,sqrt,exp,ln,arctan;
FROM LMathLib IMPORT sqrt2pi,MachEps,CardPot;
FROM Storage IMPORT ALLOCATE,DEALLOCATE;
FROM Deklera IMPORT MaxPMat,SIZELONGREAL,SIZECARDINAL;
FROM Integral IMPORT Romberg;
FROM MatLib IMPORT SortVek;
IMPORT TIO;
CONST debug = FALSE;
PROCEDURE CDF(x : LONGREAL) : LONGREAL;
VAR val,sum,sa : LONGREAL;
i : CARDINAL;
BEGIN
sum:=x; val:=x; i:=1;
REPEAT
val:=val * (x*x / (2.0*VAL(LONGREAL,i) + 1.0));
sa:=sum;
sum:=sum + val;
INC(i);
UNTIL (sum = sa);
RETURN 0.5 + (sum / sqrt(2.0*pi))*exp(-0.5*(x*x));
END CDF;
PROCEDURE phi(x : LONGREAL) : LONGREAL; (* fuer InvProbInt *)
BEGIN
RETURN exp(-0.5*x*x) / sqrt2pi;
END phi;
PROCEDURE InvProbInt(p : LONGREAL;
genau : LONGREAL) : LONGREAL;
(*----------------------------------------------------------------*)
(* Die Prozedure berechnert das Integral der Wahrscheinlichkeits- *)
(* funktion \Phi_{0,1} durch numerische Integration in Subinter- *)
(* vallen und Adaption der Schrittweiten im Intervall [0..p) *)
(*----------------------------------------------------------------*)
CONST MaxIter = 46;
VAR S,Sa,dXp,Xu,Xo : LONGREAL;
iter : CARDINAL;
BEGIN
IF (genau < 1000.0*MachEps) THEN genau:=1000.0*MachEps; END;
IF (p < 0.000005) THEN
Errors.ErrOut(
'P < 0,000005 nicht sinnvoll (StatLib.InvProbInt) !');
RETURN -3.915;
ELSIF (p < 0.5) THEN (* p = 0,5 -> InvProbInt = 0 *)
RETURN -InvProbInt(1.0 - p,genau); (* \Phi(-z) = 1 - \Phi(z) *)
ELSIF (p > 0.999995) THEN
Errors.ErrOut(
'P > 0,999995 nicht sinnvoll (StatLib.InvProbInt) !');
RETURN 3.915;
END;
S := 0.5; Sa := S;
dXp := 0.025;
Xu := 0.0; Xo := dXp;
LOOP
S:=S + Romberg(Xu,Xo,1.0E-03*genau,MaxIter,iter,phi);
IF (iter > MaxIter) THEN
Errors.ErrOut("Integral unsicher (StatLib.InvProbInt) !");
END;
IF debug THEN
TIO.WrStr("p,S,Xu,Xo,dXp = ");
TIO.WrLngReal(p,9,4);
TIO.WrLngReal(S,9,4);
TIO.WrLngReal(Xu,12,7);
TIO.WrLngReal(Xo,12,7);
TIO.WrLngReal(dXp,12,7);
TIO.WrLn;
END;
IF (S > p) THEN
dXp := dXp / 2.0; (* Halbiere die Schrittweite *)
IF (dXp < 0.1*genau) THEN EXIT END;
S := Sa; (* Setze S zurueck *)
Xo := Xu + dXp; (* Neues Intervall *)
ELSE
Xu := Xo; (* Naechstes Teilintervall *)
Xo := Xo + dXp;
Sa := S;
END;
END; (* LOOP *)
RETURN Xu + 0.25*dXp;
END InvProbInt;
PROCEDURE StudentT(n : CARDINAL;
Alpha : LONGREAL) : LONGREAL;
(* Berechnet den studentschen t-Faktor *)
PROCEDURE G1(x : LONGREAL) : LONGREAL;
BEGIN
RETURN 0.25*(x*x*x + x);
END G1;
PROCEDURE G2(x : LONGREAL) : LONGREAL;
VAR x2,x3 : LONGREAL;
BEGIN
x2 := x*x; x3:=x2*x;
RETURN (5.0*x3*x2 + 16.0*x3 + 3.0*x) / 96.0;
END G2;
PROCEDURE G3(x : LONGREAL) : LONGREAL;
VAR x2,x3 : LONGREAL;
BEGIN
x2 := x*x; x3:=x2*x;
RETURN (3.0*x3*x3*x + 19.0*x3*x2 + 17.0*x3 - 15.0*x) / 384.0;
END G3;
PROCEDURE G4(x : LONGREAL) : LONGREAL;
VAR x2,x3,x6 : LONGREAL;
BEGIN
x2 := x*x; x3:=x2*x; x6:=x3*x3;
RETURN (79.0*x6*x3 + 778.0*x6*x + 1482.0*x3*x2 - 1920.0*x3 - 945.0*x)
/ 92160.0;
END G4;
VAR p,Xp,xn : LONGREAL;
BEGIN
IF (n < 6) THEN
Errors.ErrOut('Werte ungenau für n <= 5 (StudentT) !');
END;
IF (Alpha >= 1.0) OR (Alpha <= 0.0) THEN
Errors.ErrOut('Unzulässige Eingabe für Alpha (StudentT) !');
END;
p := 1.0 - Alpha;
Xp := InvProbInt(p,1.0E-04);
xn := VAL(LONGREAL,n);
RETURN Xp + G1(Xp) / xn +
G2(Xp) / (xn*xn) +
G3(Xp) / (xn*xn*xn) +
G4(Xp) / (xn*xn*xn*xn);
END StudentT;
PROCEDURE Fisher(m,n : INTEGER;
x : LONGREAL) : LONGREAL;
CONST rezpi = 1.0 / pi; (* 0.3183098862... *)
VAR a,b,i,j : INTEGER;
w,y,z,zk,d,p : LONGREAL;
BEGIN
IF (m > 300) THEN m:= 300; END;
IF (n > 5000) THEN n:=5000; END;
a := 2*(m DIV 2) - m + 2;
b := 2*(n DIV 2) - n + 2;
w := (x*VAL(LONGREAL,m)) / VAL(LONGREAL,n);
z := 1.0 / (1.0 + w);
IF (a = 1) THEN
IF (b = 1) THEN
p := sqrt(w);
d := rezpi*z / p;
p := 2.0*rezpi*arctan(p);
ELSE
p := sqrt(w*z);
d := 0.5*p*z/w;
END;
ELSE
IF (b = 1) THEN
p := sqrt(z);
d := 0.5*z*p;
p := 1.0 - p;
ELSE
d := z*z;
p := w*z;
END;
END;
y := 2.0*w / z;
IF (a = 1) THEN
FOR j:=b+2 TO n BY 2 DO
d := d*(1.0 + 1.0 / VAL(LONGREAL,j-2))*z;
p := p + d*y / VAL(LONGREAL,j-1);
END;
ELSE
zk := CardPot(z,((n-1) DIV 2));
d := d * zk*(VAL(LONGREAL,n) / VAL(LONGREAL,b));
p := p*zk + w*z*(zk - 1.0) / (z - 1.0);
END;
y := w*z;
z := 2.0 / z;
b := n - 2;
FOR i:=a+2 TO m BY 2 DO
j := i+b;
d := (d*y*VAL(LONGREAL,j)) / VAL(LONGREAL,i-2);
p :=p - z * d / VAL(LONGREAL,j);
END;
IF (p > 1.0) THEN p:=1.0; ELSIF (p < 0.0) THEN p:=0.0; END;
RETURN p;
END Fisher;
PROCEDURE FisherM2(VAR m,n : INTEGER;
VAR x : LONGREAL) : LONGREAL;
BEGIN (* FORTRAN Schnittstelle *)
RETURN Fisher(m,n,x);
END FisherM2;
PROCEDURE Student(df : INTEGER;
t : LONGREAL) : LONGREAL;
BEGIN
RETURN Fisher(1,df,t*t);
END Student;
MODULE InverseNorm; (* Kapselt Daten für InvNorm *)
IMPORT ErrOut,sqrt,ln;
<* IF (__XDS__) THEN *>
IMPORT InvNorm;
<* ELSE *>
EXPORT InvNorm;
<* END *>
TYPE Vekt4 = ARRAY [1..4] OF LONGREAL;
Vekt5 = ARRAY [1..5] OF LONGREAL;
Vekt6 = ARRAY [1..6] OF LONGREAL;
CONST A = Vekt6{
-3.96968302866538E+01, 2.20946098424521E+02,
-2.75928510446969E+02, 1.38357751867269E+02,
-3.06647980661472E+01, 2.50662827745924E+00
};
CONST B = Vekt5{
-5.44760987982241E+01, 1.61585836858041E+02,
-1.55698979859887E+02, 6.68013118877197E+01,
-1.32806815528857E+01
};
CONST C = Vekt6{
-7.78489400243029E-03, -3.22396458041136E-01,
-2.40075827716184E+00, -2.54973253934373E+00,
4.37466414146497E+00, 2.93816398269878E+00
};
CONST D = Vekt4{
7.78469570904146E-03, 3.22467129070040E-01,
2.44513413714300E+00, 3.75440866190742E+00
};
PROCEDURE InvNorm(P : LONGREAL) : LONGREAL;
(*----------------------------------------------------------------*)
(* Normal Inverse *)
(*----------------------------------------------------------------*)
CONST lowP = 0.02425;
highP = 1.0 - lowP;
VAR z,q,r : LONGREAL;
BEGIN
IF (P <= 0.0) THEN
ErrOut("Warnung: P nicht im Interval [0,1) (StatLib.InvNorm)");
z := -7.531918931
ELSIF (P < lowP) THEN
q:=sqrt(-2.0*ln(P));
z:= (((((C[1]*q + C[2])*q + C[3])*q + C[4])*q + C[5])*q + C[6]) /
(((( D[1]*q + D[2])*q + D[3])*q + D[4])*q + 1.0);
ELSIF (P < highP) THEN
q := P - 0.5;
r := q*q;
z :=(((((A[1]*r + A[2])*r + A[3])*r + A[4])*r + A[5])*r + A[6])*q /
(((((B[1]*r + B[2])*r + B[3])*r + B[4])*r + B[5])*r + 1.0);
ELSIF (P < 1.0) THEN
q:=sqrt(-2.0*ln(1.0 - P));
z:= - (((((C[1]*q + C[2])*q + C[3])*q + C[4])*q + C[5])*q + C[6]) /
(((( D[1]*q + D[2])*q + D[3])*q + D[4])*q + 1.0);
ELSE (* P >= 1.0 *)
ErrOut("Warnung: P nicht im Interval [0,1) (StatLib.InvNorm)");
z := 8.125890658;
END;
RETURN z;
END InvNorm;
END InverseNorm;
MODULE WilcoxonDaten; (* Kapselt Daten fuer Wilcoxon-Test *)
(*------------------------------------------------------------------------*)
(* Kritische Werte des Wilcoxon Vorzeichen-Rang Tests *)
(* *)
(* 06.10.15, MRi: Erstellen der 1. Version. *)
(*------------------------------------------------------------------------*)
EXPORT A010,A005,A002,A001;
TYPE Vekt50 = ARRAY [1..50] OF LONGREAL;
CONST A010 = Vekt50{ (* Alpah = 0.10 *)
-1.0, -1.0, -1.0, -1.0, 0.0,
2.0, 3.0, 5.0, 8.0, 10.0,
13.0, 17.0, 21.0, 25.0, 30.0,
35.0, 41.0, 47.0, 53.0, 60.0,
67.0, 75.0, 83.0, 91.0, 100.0,
110.0, 119.0, 130.0, 140.0, 151.0,
163.0, 175.0, 187.0, 200.0, 213.0,
227.0, 241.0, 256.0, 271.0, 286.0,
302.0, 319.0, 336.0, 353.0, 371.0,
389.0, 407.0, 426.0, 446.0, 466.0
};
CONST A005 = Vekt50{ (* Alpah = 0.05 *)
-1.0, -1.0, -1.0, -1.0, 0.0,
0.0, 2.0, 3.0, 5.0, 8.0,
10.0, 13.0, 17.0, 21.0, 25.0,
29.0, 34.0, 40.0, 46.0, 52.0,
58.0, 65.0, 73.0, 81.0, 89.0,
98.0, 107.0, 116.0, 126.0, 137.0,
147.0, 159.0, 170.0, 182.0, 195.0,
208.0, 221.0, 235.0, 249.0, 264.0,
279.0, 294.0, 310.0, 327.0, 343.0,
361.0, 378.0, 396.0, 415.0, 434.0
};
CONST A002 = Vekt50{ (* Alpah = 0.02 *)
-1.0, -1.0, -1.0, -1.0, -1.0,
-1.0, 0.0, 1.0, 3.0, 5.0,
7.0, 9.0, 12.0, 15.0, 19.0,
23.0, 27.0, 32.0, 37.0, 43.0,
49.0, 55.0, 62.0, 69.0, 76.0,
84.0, 92.0, 101.0, 110.0, 120.0,
130.0, 140.0, 151.0, 162.0, 173.0,
185.0, 198.0, 211.0, 224.0, 238.0,
252.0, 266.0, 280.0, 296.0, 312.0,
328.0, 345.0, 362.0, 379.0, 397.0
};
CONST A001 = Vekt50{ (* Alpah = 0.01 *)
-1.0, -1.0, -1.0, -1.0, -1.0,
-1.0, -1.0, 0.0, 1.0, 3.0,
5.0, 7.0, 9.0, 12.0, 15.0,
19.0, 23.0, 27.0, 32.0, 37.0,
42.0, 48.0, 54.0, 61.0, 68.0,
75.0, 83.0, 91.0, 100.0, 109.0,
118.0, 128.0, 138.0, 148.0, 159.0,
171.0, 182.0, 194.0, 207.0, 220.0,
233.0, 247.0, 261.0, 276.0, 291.0,
307.0, 322.0, 339.0, 355.0, 373.0
};
END WilcoxonDaten;
PROCEDURE Wilcoxon(VAR X,Y : ARRAY OF LONGREAL;
nx,ny : CARDINAL;
Alpha : LONGREAL;
VAR NInv : CARDINAL;
VAR Uxy,Ua : LONGREAL;
VAR iFehl : CARDINAL) : BOOLEAN;
VAR XY : POINTER TO ARRAY [1..MaxPMat] OF LONGREAL;
index : POINTER TO ARRAY [1..MaxPMat] OF CARDINAL;
i,ii,N,N2,iy,ix : CARDINAL;
ninvy,ninvx : CARDINAL;
Nx,Ny,Za,F : LONGREAL;
BEGIN
iFehl:=0;
ALLOCATE(XY ,nx*ny*SIZELONGREAL);
ALLOCATE(index,nx*ny*SIZECARDINAL);
IF (XY = NIL) OR (index = NIL) THEN
iFehl:=98;
Fehlerflag:='Kein Freispeicher vorhanden (StatLib.Wilcoxon) !';
RETURN FALSE;
END;
ii:=1; (* Setze X und Y in XY ein. *)
FOR i:=0 TO nx-1 DO XY^[ii]:=X[i]; INC(ii); END;
FOR i:=0 TO ny-1 DO XY^[ii]:=Y[i]; INC(ii); END;
N := nx + ny;
SortVek(XY^,index^,N,-1);
iy:=0; ninvy:=0;
FOR i:=1 TO N DO (* Bestimme die Anzahl der (x,y) Inversionen. *)
IF (index^[i] > nx) THEN
INC(iy);
ELSE
INC(ninvy,iy);
END;
END;
ix:=0; ninvx:=0;
FOR i:=1 TO N DO (* Bestimme die Anzahl der (y,x) Inversionen. *)
IF (index^[i] <= nx) THEN (* !!! *)
INC(ix);
ELSE
INC(ninvx,ix);
END;
END;
(* Nehme den groesseren Wert wobei gilt (nx*xy - ninvy = ninvx) => *)
(* eigentlich mueeste man also nur ninvx oder ninvy berechen, das *)
(* ist hier nur als "Uberpr"ufung gedacht. *)
IF (ninvy > ninvx) THEN NInv:=ninvy; ELSE NInv:=ninvx; END;
IF ((ninvx + ninvy) # nx*ny) THEN
(* Hmmm - watten nu - dat sollnich vorkommen ... *)
(* Fehler bei der Berechnung der Inversionen (StatLib.Wilcoxon) *)
Fehlerflag:='Interner Fehler (StatLib.Wilcoxon) !';
iFehl:=99;
END;
N2:=((nx*ny) DIV 2);
Uxy := ABS(VAL(LONGREAL,NInv - N2));
(* Das ist noch nicht perfekt. *)
IF (N <= 50) AND (Alpha = 0.10) THEN
Ua := A010[N];
ELSIF (N <= 50) AND (Alpha = 0.05) THEN
Ua := A005[N];
ELSIF (N <= 50) AND (Alpha = 0.02) THEN
Ua := A002[N];
ELSIF (N <= 50) AND (Alpha = 0.01) THEN
Ua := A001[N];
ELSE
IF (N <= 50) THEN
iFehl := 1;
Fehlerflag:=
'Kein Vorgabewert fuer kritischen Wert Ua (StatLib.Wilcoxon) !';
END;
Nx := VAL(LONGREAL,nx);
Ny := VAL(LONGREAL,ny);
Za := InvProbInt((1.0-Alpha),1.0E-07);
F := sqrt(Nx*Ny*(Nx + Ny + 1.0) / 12.0);
Ua := Za*F;
END;
IF (Ua = -1.0) THEN (* Kleine Stichprobe ohne Vorgabewert *)
iFehl := 2;
Fehlerflag:=
'Sehr kleine Probe ohne Vorgabewert - unsicher (StatLib.Wilcoxon) !';
Nx := VAL(LONGREAL,nx);
Ny := VAL(LONGREAL,ny);
Za := InvProbInt((1.0-Alpha),1.0E-07);
F := sqrt(Nx*Ny*(Nx + Ny + 1.0) / 12.0);
Ua := Za*F;
END;
DEALLOCATE(XY ,nx*ny*SIZELONGREAL);
DEALLOCATE(index,nx*ny*SIZECARDINAL);
RETURN Uxy <= Ua;
END Wilcoxon;
END StatLib.