DataMuseum.dk

top - download

    Length: 1740 (0x6cc)
    Types: TextFile
    Names: »pof.c«

└─⟦87ddcff64⟧ Bits:30001253 CPHDIST85 Tape, 1985 Autumn Conference Copenhagen
    └─ ⟦this⟧ »cph85dist/stat/src/pof.c« 

TextFile

/*LINTLIBRARY*/
#include "unixstat.h"
FUN(pof,probability of F-ratio,5.0,1985)
/* Copyright (c) 1982 Gary Perlman (see Copyright file) */

double
pof (F, df1, df2) double F;
	{
	int	i, j;
	int	a, b;
	double	w, y, z, d, p;
	if (F <= 0.0 || df1 <= 0 || df2 <= 0) return (1.0);
	a = df1%2 ? 1 : 2;
	b = df2%2 ? 1 : 2;
	w = (F * df1) / df2;
	z = 1.0 / (1.0 + w);
	if (a == 1)
		if (b == 1)
			{
			p = sqrt (w);
			y = 1.0 / acos (-1.0);
			d = y * z / p;
			p = 2 * 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;
	for (j = b + 2; j <= df2; j += 2)
		{
		d *= (1.0 + a / (j - 2.0)) * z;
		p = (a == 1 ? p + d * y / (j - 1.0) : (p + w) * z);
		}
	y = w * z;
	z = 2.0 / z;
	b = df2 - 2;
	for (i = a + 2; i <= df1; i += 2)
		{
		j = i + b;
		d *= y * j / (i - 2.0);
		p -= z * d / j;
		}
	return (1.0-p);
	}

#ifdef STANDALONE
#include "unixstat.h"
PGM(POF,F-Ratio to Probability Conversion,5.0,3/5/85)

main (argc, argv) int argc; char **argv;
	{
	double	F;
	double	p;
	double	pof ();
	double	atof ();
	int 	df1;
	int 	df2;
	ARGV0;
	if (argc != 4)
		USAGE (F df1 df2)
	if (!number (argv[1]))
		ERRNUM (argv[1],probability of F-ratio)
	if (!INTEGER (argv[2]))
		ERRNUM (argv[2],numerator degrees of freedom)
	if (!INTEGER (argv[3]))
		ERRNUM (argv[3],denominator degrees of freedom)
	if (fzero (F = atof (argv[1])))
		ERRMSG0 (F-ratio must be positive)
	df1 = atoi (argv[2]);
	df2 = atoi (argv[3]);
	if (df1 <= 0 || df2 <= 0)
		ERRMSG0 (degrees of freedom must be positive)
	p = pof  (F, df1, df2);
	printf ("F(%d,%d) = %4.3f, p = %4.3f\n", df1, df2, F, p);
	exit (0);
	}
#endif