<<<<<<<< FILE CONTAINS BOTH A FORTRAN 77 VERSION AND
<<<<<<<< A SHAR FILE CONTAINING A C VERSION
<<<<<<<< LOOK FOR <<<<<< AS A SEPARATOR
<<< C version also provided by Larry Weissman >>>
Subject: Re: TOMS algorithm 322
Date: Tue, 29 Mar 94 15:34:51 -0800
From: Larry Weissman 5-2011
To: ehg@research.att.com
In netlib.att.com:/netlib/toms/322.Z:
> Sorry, netlib does not have this algorithm online. If you locate it
> elsewhere, please email a copy to ehg@research.att.com.
Below is a Fortran '77 coding of this algorithm, followed by a test program.
Larry Weissman Center for Bioengineering, WD-12
larryw@nsr.bioeng.washington.edu Univ of Washington, Seattle WA 98195
Phone: (206)685-2011 FAX: (206)685-3300
---------------------------------------------------------------------------
FUNCTION PROB(LM,LN,X)
C
C Probabilities in several variance distributions
C
C References:
C Egon Dorrer, "Algorithm 322: F-Distribution [S14]",
C Communications of the ACM, Vol 11(2), 1968, p116-117.
C J.B.F. Field, "Certification of Algorithm 322",
C Communications of the ACM, Vol 12(1), 1969, p39.
C Hubert Tolman, "Remark on Algorithm 322 [S14]",
C Communications of the ACM, Vol 14(2), 1979, p117.
C
C Probabilities are returned as the integral of P(q)dq for q
C in the range zero to X (F-ratio and chi-square) or from -X to +X
C (Student's t and Normal).
C
C For F-ratios:
C M = numerator degrees of freedom
C N = denominator degrees of freedom
C X = F-ratio
C For Student's t (two tailed):
C M = 1
C N = Degrees of freedom
C X = The square of t
C For normal deviates (two tailed):
C M = 1
C N = 5000
C X = The square of the deviate
C For chi-square:
C M = Degrees of freedom
C N = 5000
C X = Chi-square/M
C
C The returned probability will be in the range 0 to 1, unless an error
C occurred, in which case -1.0 is returned.
C Recoded by LW, 28-Mar-94.
C
INTEGER M,N,LM,LN,A,B,I,J
REAL X,W,Y,Z,D,P,ZK
C
M = MIN(LM,300)
N = MIN(LN,5000)
IF(MIN(M,N) .LT. 1) THEN
PROB = -1.0
RETURN
ENDIF
C
A = 2*(M/2)-M+2
B = 2*(N/2)-N+2
W = X*FLOAT(M)/FLOAT(N)
Z = 1.0/(1.0+W)
IF(A .EQ. 1) THEN
IF(B .EQ. 1)THEN
P=SQRT(W)
Y=0.3183098862
D=Y*Z/P
P=2.0*Y*ATAN(P)
ELSE
P=SQRT(W*Z)
D=P*Z/(2.0*W)
ENDIF
ELSE
IF(B .EQ. 1)THEN
P=SQRT(Z)
D=Z*P/2.0
P=1.0-P
ELSE
D=Z*Z
P=W*Z
ENDIF
ENDIF
Y = 2.0*W/Z
IF(A .EQ. 1) THEN
DO 10 J=B+2,N,2
D = (1.0 + FLOAT(A)/FLOAT(J-2)) * D * Z
P = P + D * Y/FLOAT(J-1)
10 CONTINUE
ELSE
ZK = Z**((N-1)/2)
D = D * ZK * N/B
P = P * ZK + W * Z * (ZK-1.0)/(Z-1.0)
ENDIF
Y = W * Z
Z = 2.0/Z
B = N-2
DO 20 I=A+2,M,2
J = I + B
D = (Y*D*FLOAT(J))/FLOAT(I-2)
20 P = P - Z * D/FLOAT(J)
PROB = MAX(0.0,MIN(1.0,P))
RETURN
END
-------------
PROGRAM tstprob
c
c Test function prob() that generates critical values of standard
c probability functions
c
c L. Weissman 29-Mar-94 larryw@bioeng.washington.edu
c
PARAMETER (NZ=13)
REAL z(NZ),pz(NZ)
PARAMETER (NT=25)
REAL t(NT),pt(NT)
INTEGER nut(NT)
PARAMETER (NX=25)
REAL x(NX),px(NX)
INTEGER nux(NX)
PARAMETER (NF=100)
REAL f(NF),pf(NF)
INTEGER nf1(NF),nf2(NF)
c
DATA z/4.417, 3.891, 3.291, 2.576, 2.241,
1 1.960, 1.645, 1.150, 0.674, 0.319,
2 0.126, 0.063, 0.0125/
DATA pz/1.e-5, 1.e-4, 1.e-3, 1.e-2, 0.025,
1 0.050, 0.100, 0.250, 0.500, 0.750,
2 0.900, 0.950, 0.990/
c
DATA t/1.000, 0.817, 0.700, 0.687, 0.677,
1 6.314, 2.920, 1.813, 1.725, 1.658,
2 12.706, 4.303, 2.228, 2.086, 1.980,
3 63.657, 9.925, 3.169, 2.845, 2.617,
4 127.32,14.089, 3.581, 3.153, 2.860/
DATA nut/ 1, 2, 10, 20, 120,
1 1, 2, 10, 20, 120,
2 1, 2, 10, 20, 120,
3 1, 2, 10, 20, 120,
4 1, 2, 10, 20, 120/
DATA pt/0.500, 0.500, 0.500, 0.500, 0.500,
1 0.100, 0.100, 0.100, 0.100, 0.100,
2 0.050, 0.050, 0.050, 0.050, 0.050,
3 0.010, 0.010, 0.010, 0.010, 0.010,
4 0.005, 0.005, 0.005, 0.005, 0.005/
c
DATA x/ 0.45, 1.39, 9.34, 19.34, 99.33,
1 3.84, 5.99, 18.31, 31.41,124.34,
2 5.02, 7.38, 20.48, 34.17,129.56,
3 6.63, 9.21, 23.21, 37.57,135.81,
4 7.88, 10.60, 25.19, 40.00,140.17/
DATA nux/ 1, 2, 10, 20, 100,
1 1, 2, 10, 20, 100,
2 1, 2, 10, 20, 100,
3 1, 2, 10, 20, 100,
4 1, 2, 10, 20, 100/
DATA px/0.500, 0.500, 0.500, 0.500, 0.500,
1 0.050, 0.050, 0.050, 0.050, 0.050,
2 0.025, 0.025, 0.025, 0.025, 0.025,
3 0.010, 0.010, 0.010, 0.010, 0.010,
4 0.005, 0.005, 0.005, 0.005, 0.005/
c
DATA f/ 16211.,198.50,12.826,9.9439,8.1790,
1 20000.,199.00,9.4270,6.9865,5.5393,
2 24224.,199.40,5.8467,3.8470,2.7052,
3 24836.,199.45,5.2740,3.3178,2.1881,
4 25359.,199.49,4.7501,2.8058,1.6055,
5 4052.2,98.503,10.044,8.0960,6.8510,
6 4999.5,99.000,7.5594,5.8489,4.7865,
7 6055.8,99.399,4.8492,3.3682,2.4721,
8 6208.7,99.490,4.4054,2.9377,2.0346,
9 6339.4,99.491,3.9965,2.5168,1.5330,
A 161.45,18.513,4.9646,4.3513,3.9201,
B 199.50,19.000,4.1028,3.4928,3.0718,
C 241.88,19.396,2.9782,2.3479,1.9105,
D 248.01,19.446,2.7740,2.1242,1.6587,
E 253.25,19.487,2.5801,1.8963,1.3519,
F 1.0000,.66667,.48973,.47192,.45774,
G 1.5000,1.0000,.74349,.71773,.69717,
H 2.0419,1.3450,1.0000,.96626,.93943,
I 2.1190,1.3933,1.0349,1.0000,.97228,
J 2.1848,1.4344,1.0645,1.0285,1.0000/
DATA nf1/ 5* 1, 5 *2, 5* 10, 5* 20, 5*120,
1 5* 1, 5 *2, 5* 10, 5* 20, 5*120,
2 5* 1, 5 *2, 5* 10, 5* 20, 5*120,
3 5* 1, 5 *2, 5* 10, 5* 20, 5*120/
DATA nf2/1,2,10,20,120,1,2,10,20,120,1,2,10,20,120,1,2,10,20,120,
1 1,2,10,20,120,1,2,10,20,120,1,2,10,20,120,1,2,10,20,120,
2 1,2,10,20,120,1,2,10,20,120,1,2,10,20,120,1,2,10,20,120,
3 1,2,10,20,120,1,2,10,20,120,1,2,10,20,120,1,2,10,20,120,
4 1,2,10,20,120,1,2,10,20,120,1,2,10,20,120,1,2,10,20,120/
DATA pf/25*0.005, 25*0.01, 25*0.05, 25*0.50/
c
10 FORMAT(/15x,a//4x,
1'n1 n2 x p-exp p-calc e-diff e-rel'/
22x,
3'---- ---- ------ ---------- ---------- ---------- -----')
c
WRITE(*,10)'Testing two-tailed normal distribution'
DO 20 i=1,NZ
p=1.0-prob(1,5000,z(i)**2)
CALL report(1,5000,z(i),pz(i),p)
20 CONTINUE
c
WRITE(*,10)'Testing two-tailed Student''s t'
DO 30 i=1,NT
p=1.0-prob(1,nut(i),t(i)**2)
CALL report(1,nut(i),t(i),pt(i),p)
30 CONTINUE
c
WRITE(*,10)'Testing chi-square distribution'
DO 40 i=1,NX
p=1.0-prob(nux(i),5000,x(i)/nux(i))
CALL report(nux(i),5000,x(i),px(i),p)
40 CONTINUE
c
WRITE(*,10)'Testing F distribution'
DO 50 i=1,NF
p=1.0-prob(nf1(i),nf2(i),f(i))
CALL report(nf1(i),nf2(i),f(i),pf(i),p)
50 CONTINUE
STOP
END
SUBROUTINE report(n1,n2,x,pexp,pcalc)
c
c Report the difference between the expected probability pexp and
c the calculated probability pcalc. Also report the relative error.
c
pdiff = ABS(pexp-pcalc)
prel = 0.0
IF(ABS(pexp) .GT. 1.e-30) prel = pdiff/pexp
WRITE(*,10)n1,n2,x,pexp,pcalc,pdiff,prel
10 FORMAt(1x,i5,i5,f10.4,3g12.4,f8.3)
RETURN
END
<<<<<<< 322.sh
#! /bin/sh
# This is a shell archive, meaning:
# 1. Remove everything above the #! /bin/sh line.
# 2. Save the resulting text in a file.
# 3. Execute the file with /bin/sh (not csh) to create the files:
# Makefile
# fisher.c
# fisher.h
# ftest.c
# This archive created: Thu Oct 26 21:47:05 1995
export PATH; PATH=/bin:$PATH
if test -f 'Makefile'
then
echo shar: will not over-write existing file "'Makefile'"
else
cat << \SHAR_EOF > 'Makefile'
#Makefile for fisher.c and driver program.
#Copyright (c) 1995 by Matthew Belmonte .
#Permission for use and distribution is hereby granted, subject to the
#restrictions that this copyright notice be included in its entirety, and that
#any and all changes are clearly noted.
#
CC = /bin/cc
CFLAGS = -O2
ftest: ftest.o fisher.o
$(CC) $(CFLAGS) -o ftest ftest.o fisher.o -lm
ftest.o: ftest.c fisher.h
$(CC) $(CFLAGS) -c ftest.c
fisher.o: fisher.c
$(CC) $(CFLAGS) -c fisher.c
clean:
rm -f ftest.o fisher.o ftest
SHAR_EOF
fi # end of overwriting check
if test -f 'fisher.c'
then
echo shar: will not over-write existing file "'fisher.c'"
else
cat << \SHAR_EOF > 'fisher.c'
/*fisher.c - compute the two-tailed probability of correct rejection of the null
hypothesis with an F-ratio of x, for m degrees of freedom in the numerator and
n degrees of freedom in the denominator. In the special case of only two
populations, this is equivalent to Student's t-test with m=1 and x=t**2.
Coded by Matthew Belmonte , 28 September 1995. This
implementation Copyright (c) 1995 by Matthew Belmonte. Permission for use and
distribution is hereby granted, subject to the restrictions that this
copyright notice and reference list be included in its entirety, and that any
and all changes made to the program be clearly noted in the program text.
This software is provided 'as is', with no warranty, express or implied,
including but not limited to warranties of merchantability or fitness for a
particular purpose. The user of this software assumes liability for any and
all damages, whether direct or consequential, arising from its use. The
author of this implementation will not be liable for any such damages.
References:
Egon Dorrer, "Algorithm 322: F-Distribution [S14]", Communications of the
Association for Computing Machinery 11:2:116-117 (1968).
J.B.F. Field, "Certification of Algorithm 322 [S14] F-Distribution",
Communications of the Association for Computing Machinery 12:1:39 (1969).
Hubert Tolman, "Remark on Algorithm 322 [S14] F-Distribution", Communications
of the Association for Computing Machinery 14:2:117 (1971).
*/
#include
double fisher(m, n, x)
int m, n;
double x;
{
int a, b, i, j;
double w, y, z, zk, d, p;
a = 2*(m/2)-m+2;
b = 2*(n/2)-n+2;
w = (x*m)/n;
z = 1.0/(1.0+w);
if(a == 1)
{
if(b == 1)
{
p = sqrt(w);
y = 0.3183098862;
d = y*z/p;
p = 2.0*y*atan(p);
}
else
{
p = sqrt(w*z);
d = 0.5*p*z/w;
}
}
else if(b == 1)
{
p = sqrt(z);
d = 0.5*z*p;
p = 1.0-p;
}
else
{
d = z*z;
p = w*z;
}
y = 2.0*w/z;
if(a == 1)
for(j = b+2; j <= n; j += 2)
{
d *= (1.0+1.0/(j-2))*z;
p += d*y/(j-1);
}
else
{
zk = pow(z, (double)((n-1)/2));
d *= (zk*n)/b;
p = p*zk+w*z*(zk-1.0)/(z-1.0);
}
y = w*z;
z = 2.0/z;
b = n-2;
for(i = a+2; i <= m; i += 2)
{
j = i+b;
d *= (y*j)/(i-2);
p -= z*d/j;
}
return(p<0.0? 0.0: p>1.0? 1.0: p);
}
double student(df, t)
int df;
double t;
{
return(fisher(1, df, t*t));
}
SHAR_EOF
fi # end of overwriting check
if test -f 'fisher.h'
then
echo shar: will not over-write existing file "'fisher.h'"
else
cat << \SHAR_EOF > 'fisher.h'
/*fisher.h - header file for fisher.c.
Copyright (c) 1995 by Matthew Belmonte. Permission for use and distribution
is hereby granted, subject to the restrictions that this copyright notice be
included in its entirety, and that any and all changes made to the program be
clearly noted in the program text.
This software is provided 'as is', with no warranty, express or implied,
including but not limited to warranties of merchantability or fitness for a
particular purpose. The user of this software assumes liability for any and
all damages, whether direct or consequential, arising from its use. The
author of this software will not be liable for any such damages.
*/
extern double fisher(int, int, double);
extern double student(int, double);
SHAR_EOF
fi # end of overwriting check
if test -f 'ftest.c'
then
echo shar: will not over-write existing file "'ftest.c'"
else
cat << \SHAR_EOF > 'ftest.c'
/*ftest.c - driver program for fisher.c.
Copyright (c) 1995 by Matthew Belmonte. Permission for use and distribution
is hereby granted, subject to the restrictions that this copyright notice be
included in its entirety, and that any and all changes made to the program be
clearly noted in the program text.
This software is provided 'as is', with no warranty, express or implied,
including but not limited to warranties of merchantability or fitness for a
particular purpose. The user of this software assumes liability for any and
all damages, whether direct or consequential, arising from its use. The
author of this software will not be liable for any such damages.
*/
#include
#include
#include "fisher.h"
void main()
{
double f, t;
int df, df_num;
char test_type;
printf("Enter test type, degrees of freedom, and value, or 'q' to quit -\nfor example 'f 1 10 10.04' or 't 10 3.169'\n");
test_type = getchar();
while(test_type != 'q')
{
if((test_type == 't') || (test_type == 'f'))
{
if(test_type != 't')
scanf("%d", &df_num);
else
df_num = 1;
scanf("%d", &df);
if(test_type == 't')
{
scanf("%lf", &t);
f = t*t;
}
else
{
scanf("%lf", &f);
if(df_num == 1)
t = sqrt(f);
}
printf("\tF(%d,%d)=%f, ", df_num, df, f);
if(df_num == 1)
printf("t(%d)=%f, ", df, t);
printf("p=%f\n", 1.0-fisher(df_num, df, f));
}
test_type = getchar();
}
exit(0);
}
SHAR_EOF
fi # end of overwriting check
# End of shell archive
exit 0