/*
**  prob.c	BSD 4.3 Unix	910922	SIO/ODF	fmd
**
**  Functions for calculating various probability distributions.
*/
#include <math.h>
#include <odflib.h>

#define	INV_PI		0.31830988618379067153	/* 1.0/PI */
#define EPS		1.0e-23
#define	fzero(x)	(fabs(x) < EPS)

/*
**  Compute the probability of a normal distribution z value
**  Adapted from a polynomial approximation in:
**  Ibbetson D
**  Algorithm 209
**  Collected Algorithms of the CACM 1963 p. 616
*/
double				/* <- probability of z */
ProbOfZ(z)
	double	z;		/* -> normal distribution z variate */
{
	double	y, x, w;
	if (z == 0.0)
	    x = 0.0;
	else {
	    y = 0.5 * fabs(z);
	    if (y >= 3.0)
		x = 1.0;
	    else if (y < 1.0) {
		w = y*y;
		x = (((((((( 0.000124818987  * w
			    -0.001075204047) * w
			    +0.005198775019) * w
			    -0.019198292004) * w
			    +0.059054035642) * w
			    -0.151968751364) * w
			    +0.319152932694) * w
			    -0.531923007300) * w
			    +0.797884560593) * y * 2.0;
	    } else {
		y -= 2.0;
		x = (((((((((((((-0.000045255659 * y
				+0.000152529290) * y
				-0.000019538132) * y
				-0.000676904986) * y
				+0.001390604284) * y
				-0.000794620820) * y
				-0.002034254874) * y
				+0.006549791214) * y
				-0.010557625006) * y
				+0.011630447319) * y
				-0.009279453341) * y
				+0.005353579108) * y
				-0.002141268741) * y
				+0.000535310849) * y
				+0.999936657524;
	    }
	}
	return (z > 0.0 ? ((x + 1.0) / 2.0) : ((1.0 - x) / 2.0));
}

double				/* <- standard distribution z variate */
ZOfProb(p)
	double	p;		/* -> probability of z */
{
	double	minz = -6.0;
	double	maxz =  6.0;
	double	zval = 0.0;
	double	ProbOfZ(), pval;

	if (p <= 0.0 || p >= 1.0)
	    return (0.0);
	
	while ((maxz - minz) > 0.000001) {
	    pval = ProbOfZ(zval);
	    if (pval > p)
		maxz = zval;
	    else
		minz = zval;
	    zval = (maxz + minz) * 0.5;
	}
	return (zval);
}

/*
**	Compute Probability of F ratio
**	Adapted from Collected Algorithms of the CACM
**	Algorithm 322
**	Egon Dorrer
*/
double					/* <- probabilty of F */
ProbOfF(f, df1, df2)
	double	f;			/* -> F ratio */
	int	df1;			/* -> degrees of freedom for pop 1 */
	int	df2;			/* -> degrees of freedom for pop 2 */
{
	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 = INV_PI;
		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;
	}
	/* 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				/* <- F ratio */
FOfProb(p, df1, df2)
	double	p;		/* -> probability of F */
	int	df1;		/* -> degrees of freedom for pop 1 */
	int	df2;		/* -> degrees of freedom for pop 2 */
{
	double	fval;
	double	fabs();
	double	ProbOfF();
	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 (ProbOfF(fval, df1, df2) < p) /* F too large */
		maxf = fval;
	    else /* F too small */
		minf = fval;
	    fval = (maxf + minf) * 0.5;
	}

	return (fval);
}

/*
**  Compute probability of Pearson r.
**  Adapted from "Numerical Recipies in C",
**  1988, Cambridge University Press, pg. 506
*/
double				/* <- probability of R */
ProbOfR(r, n)
	double		r;	/* -> Pearson correlation coefficient */
	int		n;	/* -> sample size */
{
	int		df;
	double		t, betai();

	df		= n-2;
	t		= r*sqrt(df/((1.0-r+EPS)*(1.0+r+EPS)));
	return (1.0-betai(0.5*df, 0.5, df/(df+t*t)));
}

double
gammln(xx)
	double		xx;
{
	double		x, tmp, ser;
	static double	c[6] = {
			76.18009173,-86.50532033,24.01409822,
			-1.231739516,0.120858003e-2,-0.536382e-5};
	int		j;

	x		= xx-1.0;
	tmp		= x+5.5;
	tmp	       -= (x+0.5)*log(tmp);
	ser		= 1.0;
	for (j=0; j<=5; j++) {
	    x	+= 1.0;
	    ser	+= c[j]/x;
	}
	return (-tmp+log(2.50662827465*ser));
}

#define	ITMAX	100
#define	BEPS	3.0e-7

double
betacf(a, b, x)
	double		a;
	double		b;
	double		x;
{
	double		qap, qam, qab, em, tem, d;
	double		bz, bm=1.0, bp, bpp;
	double		az=1.0, am=1.0, ap, app, aold;
	int		m;

	qab		= a+b;
	qap		= a+1.0;
	qam		= a-1.0;
	bz		= 1.0-qab*x/qap;
	for (m=1; m <= ITMAX; m++) {
	    em		= (double)m;
	    tem		= em+em;
	    d		= em*(b-em)*x/((qam+tem)*(a+tem));
	    ap		= az+d*am;
	    bp		= bz+d*bm;
	    d		= -(a+em)*(qab+em)*x/((qap+tem)*(a+tem));
	    app		= ap+d*az;
	    bpp		= bp+d*bz;
	    aold	= az;
	    am		= ap/bpp;
	    bm		= bp/bpp;
	    az		= app/bpp;
	    bz		= 1.0;
	    if (fabs(az-aold) < (BEPS*fabs(az)))
		return (az);
	}
	return (quiet_nan());
}

double
betai(a, b, x)
	double		a;
	double		b;
	double		x;
{
	double		bt;

	if (x < 0.0 || x > 1.0)
	    return (quiet_nan());
	if (x == 0.0 || x == 1.0)
	    bt		= 0.0;
	else
	    bt		= exp(gammln(a+b)-gammln(a)-gammln(b)+
			  a*log(x)+b*log(1.0-x));
	if (x < (a+1.0)/(a+b+2.0))
	    return (bt*betacf(a,b,x)/a);
	else
	    return (1.0-bt*betacf(b,a,1.0-x)/b);
}