DataMuseum.dk

top - download

    Length: 1486 (0x5ce)
    Types: TextFile
    Names: »z.c«

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

TextFile

/*  Copyright 1986 Gary Perlman */

/*LINTLIBRARY*/
#include "stat.h"
FUN(z,normal z-distribution functions,5.1,12/26/85)

/*
	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
poz (z)
double	z;
	{
	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
critz (p)
double	p;
	{
	double	minz = -6.0;
	double	maxz = 6.0;
	double	zval = 0.0;
	double	poz (), pval;

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