DataMuseum.dk

top - download

    Length: 2121 (0x849)
    Types: TextFile
    Names: »f.c«

└─⟦a0efdde77⟧ Bits:30001252 EUUGD11 Tape, 1987 Spring Conference Helsinki
    └─ ⟦this⟧ »EUUGD11/stat-5.3/eu/stat/src/f.c« 

TextFile

/*  Copyright 1981 Gary Perlman */

/*LINTLIBRARY*/
#include "stat.h"
FUN(f,F-distribution functions,5.2,12/26/85)

/*
	Compute Probability of F ratio
	Adapted from Collected Algorithms of the CACM
	Algorithm 322
	Egon Dorrer
*/
double
pof (F, df1, df2) double F;
	{
	int	i, j;
	int	a, b;
	double	w, y, z, d, p;

	if (fzero (F) || 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 = I_PI; /* 1 / 3.14159 */
			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;
#ifdef	REMARK /* speedup modification suggested by Tolman (wrong answer!) */
	if (a == 1)
		for (j = b + 2; j <= df2; j += 2)
			{
			d *= (1.0 + a / (j - 2.0)) * z;
			p += d * y / (j - 1.0);
			}
	else
		{
		double	zk = 1.0;
		for (j = (df2 - 1) / 2; j; j--)
			zk *= z;
		d *= zk * df2/b;
		p *= zk + w * z * (zk - 1.0)/(z-1.0);
		}
#else /* original version */
	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);
		}
#endif	REMARK
	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;
		}
	/* correction for approximation errors suggested in certification */
	if (p < 0.0)
		p = 0.0;
	else if (p > 1.0)
		p = 1.0;
	return (1.0-p);
	}

double
critf (p, df1, df2)
double	p;
int 	df1;
int 	df2;
	{
	double	fval;
	double	fabs ();
	double	pof ();
	double	maxf = 99999.0;     /* maximum possible F ratio */
	double	minf = .000001;     /* minimum possible F ratio */

	if (p <= 0.0 || p >= 1.0)
		return (0.0);
	
	fval = 1.0 / p;             /* the smaller the p, the larger the F */

	while (fabs (maxf - minf) > .000001)
		{
		if (pof (fval, df1, df2) < p) /* F too large */
			maxf = fval;
		else /* F too small */
			minf = fval;
		fval = (maxf + minf) * 0.5;
		}

	return (fval);
	}