Hallo,
ich habe gerade versucht, die Funktionen für den t-Test aus dem Buch „Numerical Recipes in Pascal“ auszuprobieren. Dummerweise bricht er aber immer mit einem Fehler ab, und zwar in der Procedure betacf (-> ‚a or b too big, or itmax too small‘, s.u.). Ich kann nicht erkennen, was schief läuft, brauche aber einen funktionierenden T-Test (eigentlich auch den „paired t-test“, aber soweit bin ich noch nicht).
Ich bin sehr dankbar für jede Hilfestellung - entweder ein Tipp, was hier nicht stimmt (s.u.) oder vielleicht sogar ein Freeware-Code würden sehr helfen.
Die Prozedur avevar funktioniert und liefert Mittelwert und Varianz der Werteliste.
Der t-Wert für data1 = [1,2,3,4] und data2 = [5,6,7,8] wird zu -4 berechnet. Der p-Wert wird dann nicht mehr berechnet.
aus betai wird betacf mit a=4, b=0.5 und x=0.33… aufgerufen. Dort wird in der Schleife beim ersten Durchlauf az=1.4329… berechnet und ändert sich dann nicht mehr, woraufhin die Schleife im 2. Durchlauf mit dem Fehler abgebrochen wird.
Hier mal der Code aus „Numerical Recipes“, etwas modifiziert, und zwar verwende ich den Typ PFltArray als ein Zeiger auf ein Array of Float mit den Indizes 0…(n-1).
FUNCTION gammln(xx: real): Float;
CONST
stp = 2.50662827465;
VAR
x,tmp,ser: double;
BEGIN
x := xx-1.0;
tmp := x+5.5;
tmp := (x+0.5)\*ln(tmp)-tmp;
ser := 1.0+76.18009173/(x+1.0)-86.50532033/(x+2.0)+24.01409822/(x+3.0)
-1.231739516/(x+4.0)+0.120858003e-2/(x+5.0)-0.536382e-5/(x+6.0);
gammln := tmp+ln(stp\*ser)
END;
//---------------------------------------------------------------------
PROCEDURE avevar(data: PFltArray;
n: integer;
VAR ave,svar: Float);
VAR
j: integer;
s: real;
BEGIN
ave := 0.0;
svar := 0.0;
FOR j := 0 TO n-1 DO ave := ave+data[j];
ave := ave/n;
FOR j := 0 TO n-1 DO BEGIN
s := data[j]-ave;
svar := svar+s\*s
END;
svar := svar/(n-1)
END;
//---------------------------------------------------------------------
FUNCTION betacf(a,b,x: Float): Float;
CONST
itmax = 100;
eps = 3.0e-7;
VAR
tem,qap,qam,qab,em,d: Float;
bz,bpp,bp,bm,az,app: Float;
am,aold,ap: Float;
m: integer;
BEGIN
am := 1.0;
bm := 1.0;
az := 1.0;
qab := a+b;
qap := a+1.0;
qam := a-1.0;
bz := 1.0-qab\*x/qap;
FOR m := 1 TO itmax DO BEGIN
em := m;
tem := em+em;
d := em\*(b-m)\*x/((qam+tem)\*(a+tem));
ap := az+d\*am;
bp := bz+d\*bm;
d := -(a+em)\*(qab+em)\*x/((a+tem)\*(qap+tem));
app := ap+d\*az;
bpp := bp+d\*bz;
aold := az;
am := ap/bpp;
bm := bp/bpp;
az := app/bpp;
bz := 1.0;
IF abs(az-aold) 1.0) THEN
raise EMathError.create('pause in routine BETAI');
IF (x = 0.0) OR (x = 1.0) THEN bt := 0.0
ELSE bt := exp(gammln(a+b)-gammln(a)-gammln(b)
+a\*ln(x)+b\*ln(1.0-x));
IF x
Gruß
Jochen