--- a
+++ b/StatLib.mod.m2pp
@@ -0,0 +1,513 @@
+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.