module Ocephes:sig..end
Documentation is organized alphabetically by
function/original-source.
val bdtr : k:int -> n:int -> float -> float 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
val bdtrc : k:int -> n:int -> float -> float 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
val bdtri : k:int -> n:int -> float -> float 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
val btdtr : a:float -> b:float -> float -> float 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.
val chbevl : ?n:int -> float -> float array -> float 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.
val chdtr : df:float -> float -> float 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
val chdtrc : v:float -> float -> float 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
val chdtri : df:float -> float -> float 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
val expx2 : float -> int -> float 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
val fdtr : df1:int -> df2:int -> float -> float 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
val fdtrc : df1:int -> df2:int -> float -> float 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
val fdtri : df1:int -> df2:int -> float -> float 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
val gamma : float -> floatgamma x computes the gamma function for x, this is the extension of the
factorial function to real numbers.
Gamma function
SYNOPSIS:
double x, y, gamma();
extern int sgngam;
y = gamma( x );
DESCRIPTION:
Returns gamma function of the argument. The result is
correctly signed, and the sign (+1 or -1) is also
returned in a global (extern) variable named sgngam.
This variable is also filled in by the logarithmic gamma
function lgam().
Arguments |x| <= 34 are reduced by recurrence and the function
approximated by a rational function of degree 6/7 in the
interval (2,3). Large arguments are handled by Stirling's
formula. Large negative arguments are made positive using
a reflection formula.
ACCURACY:
Relative error:
arithmetic domain # trials peak rms
DEC -34, 34 10000 1.3e-16 2.5e-17
IEEE -170,-33 20000 2.3e-15 3.3e-16
IEEE -33, 33 20000 9.4e-16 2.2e-16
IEEE 33, 171.6 20000 2.3e-15 3.2e-16
Error for arguments outside the test range will be larger
owing to error amplification by the exponential function.
val lgam : float -> floatlgam x computes the logarithm of the gamma function.
Natural logarithm of gamma function
SYNOPSIS:
double x, y, lgam();
extern int sgngam;
y = lgam( x );
DESCRIPTION:
Returns the base e (2.718...) logarithm of the absolute
value of the gamma function of the argument.
The sign (+1 or -1) of the gamma function is returned in a
global (extern) variable named sgngam.
For arguments greater than 13, the logarithm of the gamma
function is approximated by the logarithmic version of
Stirling's formula using a polynomial approximation of
degree 4. Arguments between -33 and +33 are reduced by
recurrence to the interval [2,3] of a rational approximation.
The cosecant reflection formula is employed for arguments
less than -33.
Arguments greater than MAXLGM return MAXNUM and an error
message. MAXLGM = 2.035093e36 for DEC
arithmetic or 2.556348e305 for IEEE arithmetic.
ACCURACY:
arithmetic domain # trials peak rms
DEC 0, 3 7000 5.2e-17 1.3e-17
DEC 2.718, 2.035e36 5000 3.9e-17 9.9e-18
IEEE 0, 3 28000 5.4e-16 1.1e-16
IEEE 2.718, 2.556e305 40000 3.5e-16 8.3e-17
The error criterion was relative when the function magnitude
was greater than one but absolute when it was less than one.
The following test used the relative error criterion, though
at certain points the relative error could be much higher than
indicated.
IEEE -200, -4 10000 4.8e-16 1.3e-16
val gdtr : a:float -> b:float -> float -> float 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
val gdtrc : a:float -> b:float -> float -> float 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
val igam : a:float -> float -> floatigam a x computes the regularized (divided by gamma a)
lower incomplete (integrate from 0 to x) gamma function.
Incomplete gamma integral
SYNOPSIS:
double a, x, y, igam();
y = igam( a, x );
DESCRIPTION:
The function is defined by
x
-
1 | | -t a-1
igam(a,x) = ----- | e t dt.
- | |
| (a) -
0
In this implementation both arguments must be positive.
The integral is evaluated by either a power series or
continued fraction expansion, depending on the relative
values of a and x.
ACCURACY:
Relative error:
arithmetic domain # trials peak rms
IEEE 0,30 200000 3.6e-14 2.9e-15
IEEE 0,100 300000 9.9e-14 1.5e-14
val igamc : a:float -> float -> floatigamc a x computes the regularized (divided by gamma a)
uppwer incomplete (integrate x to infinity) gamma function.
Complemented incomplete gamma integral
SYNOPSIS:
double a, x, y, igamc();
y = igamc( a, x );
DESCRIPTION:
The function is defined by
igamc(a,x) = 1 - igam(a,x)
inf.
-
1 | | -t a-1
= ----- | e t dt.
- | |
| (a) -
x
In this implementation both arguments must be positive.
The integral is evaluated by either a power series or
continued fraction expansion, depending on the relative
values of a and x.
ACCURACY:
Tested at random a, x.
a x Relative error:
arithmetic domain domain # trials peak rms
IEEE 0.5,100 0,100 200000 1.9e-14 1.7e-15
IEEE 0.01,0.5 0,100 200000 1.4e-13 1.6e-15
val igami : a:float -> float -> float 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
val incbet : a:float -> b:float -> float -> float 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
val incbi : a:float -> b:float -> float -> float 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
val nbdtr : k:int -> n:int -> float -> float 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.
val nbdtrc : k:int -> n:int -> float -> float 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.
val nbdtri : k:int -> n:int -> float -> float 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.
val ndtr : float -> float 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
val erf : float -> float 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
val erfc : float -> float 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
.val ndtri : float -> floatndtri p computes the value x such that the integral of the normal
probability density function from neg_infinity to x is p.
Inverse of Normal distribution function
SYNOPSIS:
double x, y, ndtri();
x = ndtri( y );
DESCRIPTION:
Returns the argument, x, for which the area under the
Gaussian probability density function (integrated from
minus infinity to x) is equal to y.
For small arguments 0 < y < exp(-2), the program computes
z = sqrt( -2.0 * log(y) ); then the approximation is
x = z - log(z)/z - (1/z) P(1/z) / Q(1/z).
There are two rational functions P/Q, one for 0 < y < exp(-32)
and the other for y up to exp(-2). For larger arguments,
w = y - 0.5, and x/sqrt(2pi) = w + w**3 R(w**2)/S(w**2)).
ACCURACY:
Relative error:
arithmetic domain # trials peak rms
DEC 0.125, 1 5500 9.5e-17 2.1e-17
DEC 6e-39, 0.135 3500 5.7e-17 1.3e-17
IEEE 0.125, 1 20000 7.2e-16 1.3e-16
IEEE 3e-308, 0.135 50000 4.6e-16 9.8e-17
ERROR MESSAGES:
message condition value returned
ndtri domain x <= 0 -MAXNUM
ndtri domain x >= 1 MAXNUM
val pdtr : k:int -> float -> floatPoisson 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().
val pdtrc : k:int -> float -> floatComplemented 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.
val pdtri : k:int -> float -> float 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
val stdtr : k:int -> float -> float 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
val stdtri : k:int -> float -> floatstdtri 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.
Functional inverse of Student's t distribution
SYNOPSIS:
double p, t, stdtri();
int k;
t = stdtri( k, p );
DESCRIPTION:
Given probability p, finds the argument t such that stdtr(k,t)
is equal to p.
ACCURACY:
Tested at random 1 <= k <= 100. The "domain" refers to p:
Relative error:
arithmetic domain # trials peak rms
IEEE .001,.999 25000 5.7e-15 8.0e-16
IEEE 10^-6,.001 25000 2.0e-12 2.9e-14