B | |
| bdtr [Ocephes] | Binomial distribution
SYNOPSIS:
int k, n;
double p, y, bdtr();
y = bdtr( k, n, p );
DESCRIPTION:
Returns the sum of the terms 0 through k of the Binomial
probability density:
k
-- ( n ) j n-j
> ( ) p (1-p)
-- ( j )
j=0
The terms are not summed directly; instead the incomplete
beta integral is employed, according to the formula
y = bdtr( k, n, p ) = incbet( n-k, k+1, 1-p ).
The arguments must be positive, with p ranging from 0 to 1.
ACCURACY:
Tested at random points (a,b,p), with p between 0 and 1.
a,b Relative error:
arithmetic domain # trials peak rms
For p between 0.001 and 1:
IEEE 0,100 100000 4.3e-15 2.6e-16
See also incbet.c.
ERROR MESSAGES:
message condition value returned
bdtr domain k < 0 0.0
n < k
x < 0, x > 1
|
| bdtrc [Ocephes] | Complemented binomial distribution
SYNOPSIS:
int k, n;
double p, y, bdtrc();
y = bdtrc( k, n, p );
DESCRIPTION:
Returns the sum of the terms k+1 through n of the Binomial
probability density:
n
-- ( n ) j n-j
> ( ) p (1-p)
-- ( j )
j=k+1
The terms are not summed directly; instead the incomplete
beta integral is employed, according to the formula
y = bdtrc( k, n, p ) = incbet( k+1, n-k, p ).
The arguments must be positive, with p ranging from 0 to 1.
ACCURACY:
Tested at random points (a,b,p).
a,b Relative error:
arithmetic domain # trials peak rms
For p between 0.001 and 1:
IEEE 0,100 100000 6.7e-15 8.2e-16
For p between 0 and .001:
IEEE 0,100 100000 1.5e-13 2.7e-15
ERROR MESSAGES:
message condition value returned
bdtrc domain x<0, x>1, n<k 0.0
|
| bdtri [Ocephes] | Inverse binomial distribution
SYNOPSIS:
int k, n;
double p, y, bdtri();
p = bdtr( k, n, y );
DESCRIPTION:
Finds the event probability p such that the sum of the
terms 0 through k of the Binomial probability density
is equal to the given cumulative probability y.
This is accomplished using the inverse beta integral
function and the relation
1 - p = incbi( n-k, k+1, y ).
ACCURACY:
Tested at random points (a,b,p).
a,b Relative error:
arithmetic domain # trials peak rms
For p between 0.001 and 1:
IEEE 0,100 100000 2.3e-14 6.4e-16
IEEE 0,10000 100000 6.6e-12 1.2e-13
For p between 10^-6 and 0.001:
IEEE 0,100 100000 2.0e-12 1.3e-14
IEEE 0,10000 100000 1.5e-12 3.2e-14
See also incbi.c.
ERROR MESSAGES:
message condition value returned
bdtri domain k < 0, n <= k 0.0
x < 0, x > 1
|
| btdtr [Ocephes] | Beta distribution
SYNOPSIS:
double a, b, x, y, btdtr();
y = btdtr( a, b, x );
DESCRIPTION:
Returns the area from zero to x under the beta density
function:
x
- -
| (a+b) | | a-1 b-1
P(x) = ---------- | t (1-t) dt
- - | |
| (a) | (b) -
0
This function is identical to the incomplete beta
integral function incbet(a, b, x).
The complemented function is
1 - P(1-x) = incbet( b, a, x );
ACCURACY:
See incbet.c.
|
C | |
| chbevl [Ocephes] | Evaluate Chebyshev series
SYNOPSIS:
int N;
double x, y, coef[N], chebevl();
y = chbevl( x, coef, N );
DESCRIPTION:
Evaluates the series
N-1
- '
y = > coef[i] T (x/2)
- i
i=0
of Chebyshev polynomials Ti at argument x/2.
Coefficients are stored in reverse order, i.e. the zero
order term is last in the array. Note N is the number of
coefficients, not the order.
If coefficients are for the interval a to b, x must
have been transformed to x -> 2(2x - b - a)/(b-a) before
entering the routine. This maps x from (a, b) to (-1, 1),
over which the Chebyshev polynomials are defined.
If the coefficients are for the inverted interval, in
which (a, b) is mapped to (1/b, 1/a), the transformation
required is x -> 2(2ab/x - b - a)/(b-a). If b is infinity,
this becomes x -> 4a/x - 1.
SPEED:
Taking advantage of the recurrence properties of the
Chebyshev polynomials, the routine requires one more
addition per loop than evaluating a nested polynomial of
the same degree.
|
| chdtr [Ocephes] | Chi-square distribution
SYNOPSIS:
double df, x, y, chdtr();
y = chdtr( df, x );
DESCRIPTION:
Returns the area under the left hand tail (from 0 to x)
of the Chi square probability density function with
v degrees of freedom.
inf.
-
1 | | v/2-1 -t/2
P( x | v ) = ----------- | t e dt
v/2 - | |
2 | (v/2) -
x
where x is the Chi-square variable.
The incomplete gamma integral is used, according to the
formula
y = chdtr( v, x ) = igam( v/2.0, x/2.0 ).
The arguments must both be positive.
ACCURACY:
See igam().
ERROR MESSAGES:
message condition value returned
chdtr domain x < 0 or v < 1 0.0
|
| chdtrc [Ocephes] | Complemented Chi-square distribution
SYNOPSIS:
double v, x, y, chdtrc();
y = chdtrc( v, x );
DESCRIPTION:
Returns the area under the right hand tail (from x to
infinity) of the Chi square probability density function
with v degrees of freedom:
inf.
-
1 | | v/2-1 -t/2
P( x | v ) = ----------- | t e dt
v/2 - | |
2 | (v/2) -
x
where x is the Chi-square variable.
The incomplete gamma integral is used, according to the
formula
y = chdtr( v, x ) = igamc( v/2.0, x/2.0 ).
The arguments must both be positive.
ACCURACY:
See igamc().
ERROR MESSAGES:
message condition value returned
chdtrc domain x < 0 or v < 1 0.0
|
| chdtri [Ocephes] | Inverse of complemented Chi-square distribution
SYNOPSIS:
double df, x, y, chdtri();
x = chdtri( df, y );
DESCRIPTION:
Finds the Chi-square argument x such that the integral
from x to infinity of the Chi-square density is equal
to the given cumulative probability y.
This is accomplished using the inverse gamma integral
function and the relation
x/2 = igami( df/2, y );
ACCURACY:
See igami.c.
ERROR MESSAGES:
message condition value returned
chdtri domain y < 0 or y > 1 0.0
v < 1
|
E | |
| erf [Ocephes] | Error function
SYNOPSIS:
double x, y, erf();
y = erf( x );
DESCRIPTION:
The integral is
x
-
2 | | 2
erf(x) = -------- | exp( - t ) dt.
sqrt(pi) | |
-
0
The magnitude of x is limited to 9.231948545 for DEC
arithmetic; 1 or -1 is returned outside this range.
For 0 <= |x| < 1, erf(x) = x * P4(x**2)/Q5(x**2); otherwise
erf(x) = 1 - erfc(x).
ACCURACY:
Relative error:
arithmetic domain # trials peak rms
DEC 0,1 14000 4.7e-17 1.5e-17
IEEE 0,1 30000 3.7e-16 1.0e-16
|
| erfc [Ocephes] | Complementary error function
SYNOPSIS:
double x, y, erfc();
y = erfc( x );
DESCRIPTION:
1 - erf(x) =
inf.
-
2 | | 2
erfc(x) = -------- | exp( - t ) dt
sqrt(pi) | |
-
x
For small x, erfc(x) = 1 - erf(x); otherwise rational
approximations are computed.
A special function expx2.c is used to suppress error amplification
in computing exp(-x^2).
ACCURACY:
Relative error:
arithmetic domain # trials peak rms
IEEE 0,26.6417 30000 1.3e-15 2.2e-16
ERROR MESSAGES:
message condition value returned
erfc underflow x > 9.231948545 (DEC) 0.0
.
|
| expx2 [Ocephes] | Exponential of squared argument
SYNOPSIS:
double x, y, expx2();
int sign;
y = expx2( x, sign );
DESCRIPTION:
Computes y = exp(x*x) while suppressing error amplification
that would ordinarily arise from the inexactness of the
exponential argument x*x.
If sign < 0, the result is inverted; i.e., y = exp(-x*x) .
ACCURACY:
Relative error:
arithmetic domain # trials peak rms
IEEE -26.6, 26.6 10^7 3.9e-16 8.9e-17
|
F | |
| fdtr [Ocephes] | F distribution
SYNOPSIS:
int df1, df2;
double x, y, fdtr();
y = fdtr( df1, df2, x );
DESCRIPTION:
Returns the area from zero to x under the F density
function (also known as Snedcor's density or the
variance ratio density). This is the density
of x = (u1/df1)/(u2/df2), where u1 and u2 are random
variables having Chi square distributions with df1
and df2 degrees of freedom, respectively.
The incomplete beta integral is used, according to the
formula
P(x) = incbet( df1/2, df2/2, (df1*x/(df2 + df1*x) ).
The arguments a and b are greater than zero, and x is
nonnegative.
ACCURACY:
Tested at random points (a,b,x).
x a,b Relative error:
arithmetic domain domain # trials peak rms
IEEE 0,1 0,100 100000 9.8e-15 1.7e-15
IEEE 1,5 0,100 100000 6.5e-15 3.5e-16
IEEE 0,1 1,10000 100000 2.2e-11 3.3e-12
IEEE 1,5 1,10000 100000 1.1e-11 1.7e-13
See also incbet.c.
ERROR MESSAGES:
message condition value returned
fdtr domain a<0, b<0, x<0 0.0
|
| fdtrc [Ocephes] | Complemented F distribution
SYNOPSIS:
int df1, df2;
double x, y, fdtrc();
y = fdtrc( df1, df2, x );
DESCRIPTION:
Returns the area from x to infinity under the F density
function (also known as Snedcor's density or the
variance ratio density).
inf.
-
1 | | a-1 b-1
1-P(x) = ------ | t (1-t) dt
B(a,b) | |
-
x
The incomplete beta integral is used, according to the
formula
P(x) = incbet( df2/2, df1/2, (df2/(df2 + df1*x) ).
ACCURACY:
Tested at random points (a,b,x) in the indicated intervals.
x a,b Relative error:
arithmetic domain domain # trials peak rms
IEEE 0,1 1,100 100000 3.7e-14 5.9e-16
IEEE 1,5 1,100 100000 8.0e-15 1.6e-15
IEEE 0,1 1,10000 100000 1.8e-11 3.5e-13
IEEE 1,5 1,10000 100000 2.0e-11 3.0e-12
See also incbet.c.
ERROR MESSAGES:
message condition value returned
fdtrc domain a<0, b<0, x<0 0.0
|
| fdtri [Ocephes] | Inverse of complemented F distribution
SYNOPSIS:
int df1, df2;
double x, p, fdtri();
x = fdtri( df1, df2, p );
DESCRIPTION:
Finds the F density argument x such that the integral
from x to infinity of the F density is equal to the
given probability p.
This is accomplished using the inverse beta integral
function and the relations
z = incbi( df2/2, df1/2, p )
x = df2 (1-z) / (df1 z).
Note: the following relations hold for the inverse of
the uncomplemented F distribution:
z = incbi( df1/2, df2/2, p )
x = df2 z / (df1 (1-z)).
ACCURACY:
Tested at random points (a,b,p).
a,b Relative error:
arithmetic domain # trials peak rms
For p between .001 and 1:
IEEE 1,100 100000 8.3e-15 4.7e-16
IEEE 1,10000 100000 2.1e-11 1.4e-13
For p between 10^-6 and 10^-3:
IEEE 1,100 50000 1.3e-12 8.4e-15
IEEE 1,10000 50000 3.0e-12 4.8e-14
See also fdtrc.c.
ERROR MESSAGES:
message condition value returned
fdtri domain p <= 0 or p > 1 0.0
v < 1
|
G | |
| gamma [Ocephes] | gamma x computes the gamma function for x, this is the extension of the
factorial function to real numbers.
|
| gdtr [Ocephes] | Gamma distribution function
SYNOPSIS:
double a, b, x, y, gdtr();
y = gdtr( a, b, x );
DESCRIPTION:
Returns the integral from zero to x of the gamma probability
density function:
x
b -
a | | b-1 -at
y = ----- | t e dt
- | |
| (b) -
0
The incomplete gamma integral is used, according to the
relation
y = igam( b, ax ).
ACCURACY:
See igam().
ERROR MESSAGES:
message condition value returned
gdtr domain x < 0 0.0
|
| gdtrc [Ocephes] | Complemented gamma distribution function
SYNOPSIS:
double a, b, x, y, gdtrc();
y = gdtrc( a, b, x );
DESCRIPTION:
Returns the integral from x to infinity of the gamma
probability density function:
inf.
b -
a | | b-1 -at
y = ----- | t e dt
- | |
| (b) -
x
The incomplete gamma integral is used, according to the
relation
y = igamc( b, ax ).
ACCURACY:
See igamc().
ERROR MESSAGES:
message condition value returned
gdtrc domain x < 0 0.0
|
I | |
| igam [Ocephes] | igam a x computes the regularized (divided by gamma a)
lower incomplete (integrate from 0 to x) gamma function.
|
| igamc [Ocephes] | igamc a x computes the regularized (divided by gamma a)
uppwer incomplete (integrate x to infinity) gamma function.
|
| igami [Ocephes] | Inverse of complemented imcomplete gamma integral
SYNOPSIS:
double a, x, p, igami();
x = igami( a, p );
DESCRIPTION:
Given p, the function finds x such that
igamc( a, x ) = p.
Starting with the approximate value
3
x = a t
where
t = 1 - d - ndtri(p) sqrt(d)
and
d = 1/9a,
the routine performs up to 10 Newton iterations to find the
root of igamc(a,x) - p = 0.
ACCURACY:
Tested at random a, p in the intervals indicated.
a p Relative error:
arithmetic domain domain # trials peak rms
IEEE 0.5,100 0,0.5 100000 1.0e-14 1.7e-15
IEEE 0.01,0.5 0,0.5 100000 9.0e-14 3.4e-15
IEEE 0.5,10000 0,0.5 20000 2.3e-13 3.8e-14
|
| incbet [Ocephes] | Incomplete beta integral
SYNOPSIS:
double a, b, x, y, incbet();
y = incbet( a, b, x );
DESCRIPTION:
Returns incomplete beta integral of the arguments, evaluated
from zero to x. The function is defined as
x
- -
| (a+b) | | a-1 b-1
----------- | t (1-t) dt.
- - | |
| (a) | (b) -
0
The domain of definition is 0 <= x <= 1. In this
implementation a and b are restricted to positive values.
The integral from x to 1 may be obtained by the symmetry
relation
1 - incbet( a, b, x ) = incbet( b, a, 1-x ).
The integral is evaluated by a continued fraction expansion
or, when b*x is small, by a power series.
ACCURACY:
Tested at uniformly distributed random points (a,b,x) with a and b
in "domain" and x between 0 and 1.
Relative error
arithmetic domain # trials peak rms
IEEE 0,5 10000 6.9e-15 4.5e-16
IEEE 0,85 250000 2.2e-13 1.7e-14
IEEE 0,1000 30000 5.3e-12 6.3e-13
IEEE 0,10000 250000 9.3e-11 7.1e-12
IEEE 0,100000 10000 8.7e-10 4.8e-11
Outputs smaller than the IEEE gradual underflow threshold
were excluded from these statistics.
ERROR MESSAGES:
message condition value returned
incbet domain x<0, x>1 0.0
incbet underflow 0.0
|
| incbi [Ocephes] | Inverse of incomplete beta integral
SYNOPSIS:
double a, b, x, y, incbi();
x = incbi( a, b, y );
DESCRIPTION:
Given y, the function finds x such that
incbet( a, b, x ) = y .
The routine performs interval halving or Newton iterations to find the
root of incbet(a,b,x) - y = 0.
ACCURACY:
Relative error:
x a,b
arithmetic domain domain # trials peak rms
IEEE 0,1 .5,10000 50000 5.8e-12 1.3e-13
IEEE 0,1 .25,100 100000 1.8e-13 3.9e-15
IEEE 0,1 0,5 50000 1.1e-12 5.5e-15
VAX 0,1 .5,100 25000 3.5e-14 1.1e-15
With a and b constrained to half-integer or integer values:
IEEE 0,1 .5,10000 50000 5.8e-12 1.1e-13
IEEE 0,1 .5,100 100000 1.7e-14 7.9e-16
With a = .5, b constrained to half-integer or integer values:
IEEE 0,1 .5,10000 10000 8.3e-11 1.0e-11
|
L | |
| lgam [Ocephes] | lgam x computes the logarithm of the gamma function.
|
N | |
| nbdtr [Ocephes] | Negative binomial distribution
SYNOPSIS:
int k, n;
double p, y, nbdtr();
y = nbdtr( k, n, p );
DESCRIPTION:
Returns the sum of the terms 0 through k of the negative
binomial distribution:
k
-- ( n+j-1 ) n j
> ( ) p (1-p)
-- ( j )
j=0
In a sequence of Bernoulli trials, this is the probability
that k or fewer failures precede the nth success.
The terms are not computed individually; instead the incomplete
beta integral is employed, according to the formula
y = nbdtr( k, n, p ) = incbet( n, k+1, p ).
The arguments must be positive, with p ranging from 0 to 1.
ACCURACY:
Tested at random points (a,b,p), with p between 0 and 1.
a,b Relative error:
arithmetic domain # trials peak rms
IEEE 0,100 100000 1.7e-13 8.8e-15
See also incbet.c.
|
| nbdtrc [Ocephes] | Complemented negative binomial distribution
SYNOPSIS:
int k, n;
double p, y, nbdtrc();
y = nbdtrc( k, n, p );
DESCRIPTION:
Returns the sum of the terms k+1 to infinity of the negative
binomial distribution:
inf
-- ( n+j-1 ) n j
> ( ) p (1-p)
-- ( j )
j=k+1
The terms are not computed individually; instead the incomplete
beta integral is employed, according to the formula
y = nbdtrc( k, n, p ) = incbet( k+1, n, 1-p ).
The arguments must be positive, with p ranging from 0 to 1.
ACCURACY:
Tested at random points (a,b,p), with p between 0 and 1.
a,b Relative error:
arithmetic domain # trials peak rms
IEEE 0,100 100000 1.7e-13 8.8e-15
See also incbet.c.
|
| nbdtri [Ocephes] | Functional inverse of negative binomial distribution
SYNOPSIS:
int k, n;
double p, y, nbdtri();
p = nbdtri( k, n, y );
DESCRIPTION:
Finds the argument p such that nbdtr(k,n,p) is equal to y.
ACCURACY:
Tested at random points (a,b,y), with y between 0 and 1.
a,b Relative error:
arithmetic domain # trials peak rms
IEEE 0,100 100000 1.5e-14 8.5e-16
See also incbi.c.
|
| ndtr [Ocephes] | Normal distribution function
SYNOPSIS:
double x, y, ndtr();
y = ndtr( x );
DESCRIPTION:
Returns the area under the Gaussian probability density
function, integrated from minus infinity to x:
x
-
1 | | 2
ndtr(x) = --------- | exp( - t /2 ) dt
sqrt(2pi) | |
-
-inf.
= ( 1 + erf(z) ) / 2
= erfc(z) / 2
where z = x/sqrt(2). Computation is via the functions
erf and erfc with care to avoid error amplification in computing exp(-x^2).
ACCURACY:
Relative error:
arithmetic domain # trials peak rms
IEEE -13,0 30000 1.3e-15 2.2e-16
ERROR MESSAGES:
message condition value returned
erfc underflow x > 37.519379347 0.0
|
| ndtri [Ocephes] | ndtri p computes the value x such that the integral of the normal
probability density function from neg_infinity to x is p.
|
P | |
| pdtr [Ocephes] | Poisson distribution SYNOPSIS: int k; double m, y, pdtr(); y = pdtr( k, m ); DESCRIPTION: Returns the sum of the first k terms of the Poisson distribution: k j -- -m m > e -- -- j! j=0 The terms are not summed directly; instead the incomplete gamma integral is employed, according to the relation y = pdtr( k, m ) = igamc( k+1, m ). The arguments must both be positive. ACCURACY: See igamc(). |
| pdtrc [Ocephes] | Complemented poisson distribution SYNOPSIS: int k; double m, y, pdtrc(); y = pdtrc( k, m ); DESCRIPTION: Returns the sum of the terms k+1 to infinity of the Poisson distribution: inf. j -- -m m > e -- -- j! j=k+1 The terms are not summed directly; instead the incomplete gamma integral is employed, according to the formula y = pdtrc( k, m ) = igam( k+1, m ). The arguments must both be positive. ACCURACY: See igam.c. |
| pdtri [Ocephes] | Inverse Poisson distribution
SYNOPSIS:
int k;
double m, y, pdtr();
m = pdtri( k, y );
DESCRIPTION:
Finds the Poisson variable x such that the integral
from 0 to x of the Poisson density is equal to the
given probability y.
This is accomplished using the inverse gamma integral
function and the relation
m = igami( k+1, y ).
ACCURACY:
See igami.c.
ERROR MESSAGES:
message condition value returned
pdtri domain y < 0 or y >= 1 0.0
k < 0
|
S | |
| stdtr [Ocephes] | Student's t distribution
SYNOPSIS:
double t, stdtr();
short k;
y = stdtr( k, t );
DESCRIPTION:
Computes the integral from minus infinity to t of the Student
t distribution with integer k > 0 degrees of freedom:
t
-
| |
- | 2 -(k+1)/2
| ( (k+1)/2 ) | ( x )
---------------------- | ( 1 + --- ) dx
- | ( k )
sqrt( k pi ) | ( k/2 ) |
| |
-
-inf.
Relation to incomplete beta integral:
1 - stdtr(k,t) = 0.5 * incbet( k/2, 1/2, z )
where
z = k/(k + t**2).
For t < -2, this is the method of computation. For higher t,
a direct method is derived from integration by parts.
Since the function is symmetric about t=0, the area under the
right tail of the density is found by calling the function
with -t instead of t.
ACCURACY:
Tested at random 1 <= k <= 25. The "domain" refers to t.
Relative error:
arithmetic domain # trials peak rms
IEEE -100,-2 50000 5.9e-15 1.4e-15
IEEE -2,100 500000 2.7e-15 4.9e-17
|
| stdtri [Ocephes] | stdtri k p computes the value x such that the integral
(from neg_infinity to x) of the Student's T probability density
function, with k degrees of freedom is p.
|