### Program to calculate one- and two-way probabilities for t-tests

DEF FNt(x) = (1 + x^2/v)^(exponent)
DIM n9(16), t(136)
INPUT "enter n (v=n-2 will be computed) = ",n
v = n-2
x = (v+1)/2
GOSUB gamma
g1 = gamma
x = v/2
GOSUB gamma
g2 = gamma
g2 = SQR(v*pi)*g2
constant = g1/g2
exponent = -(v+1)/2
GOSUB t:
gamma:
pi = 3.14159
y = x + 2
g = SQR(2*pi/y) * EXP(y * LOG(y) + (1 - 1/(30*y^2))/(12*y) - y)
gamma = g/(x * (x+1))
RETURN
t:
lower = 0
INPUT "t = ",tstat
upper = ABS(tstat)
tol = .001
GOSUB integration
PRINT
PRINT "probability of"
PRINT
IF tstat < 0 THEN
PRINT " t <=";
PRINT tstat;
PRINT USING " = #.######"; .5 - sum*constant
ELSE
TABLE 10-4 (CONTINUED)
PRINT " t >=";
PRINT tstat;
PRINT USING " = #.######"; .5 - sum*constant
END IF
PRINT
PRINT " ";:
PRINT -tstat;
PRINT " < t < ";
PRINT tstat;
PRINT USING " = #.######"; (.5 - sum*constant)*2
a$ = INPUT$(1)
END
integration:
panels = 1
n9(1) = 1
delta = (upper - lower)/panels
c = (FNt(lower) + FNt(upper))/2
t(1) = delta * c
n = 1
n8 = 2
sum = c
again:
n = n+1
f3 = 4
n9(n) = n8
panels = panels*2
l = panels - 1
delta = (upper - lower)/panels
FOR ii = 1 TO (l+1)/2
i = ii*2 - 1
x = lower + i*delta
sum = sum + FNt(x)
NEXT ii
t(n8) = delta*sum
n7 = n9(n-1)
k = n-1
FOR m = 1 TO k
j = n8 + m
n6 = n9(n-1) + m - 1
t(j) = (f3*t(j-1) - t(n6))/(f3 - 1)
f3 = f3*4
NEXT m
IF n<=4 THEN more
IF t(n8 + 1) = 0 THEN more
IF ABS(t(n7 + 1) - t(n8 + 1)) <= ABS(t(n8+1)*tol) THEN out
IF ABS(t(n8-1) - t(j)) = ABS(t(j)*tol) THEN out
IF n > 15 THEN fail
TABLE 10-4 (CONTINUED)
more: n8 = j + 1
GOTO again
out: sum = t(j)
RETURN
fail: PRINT "no convergence"
RETURN

Parts of this program were adapted from Alan R. Miller,
*BASIC Programs for Scientists and Engineers*, Berkeley: Sybex, 1981.