Implementation of the gamma and beta functions, and their integrals.
License:
BSD style: see license.txt
Authors:
Stephen L. Moshier (original C code). Conversion to D by Don Clugston
- const real
MAXGAMMA
;
- The maximum value of x for which gamma(x) < real.infinity.
- real
sgnGamma
(real x);
- The sign of Γ(x).
Returns -1 if Γ(x) < 0, +1 if Γ(x) > 0,
NAN if sign is indeterminate.
- real
gamma
(real x);
- The Gamma function, Γ(x)
Γ(x) is a generalisation of the factorial function
to real and complex numbers.
Like x!, Γ(x+1) = x*Γ(x).
Mathematically, if z.re > 0 then
Γ(z) = ∫0∞ tz-1e-t dt
Special Values
| x | Γ(x) |
|---|
| NAN | NAN |
| ±0.0 | ±∞ |
| integer > 0 | (x-1)! |
| integer < 0 | NAN |
| +∞ | +∞ |
| -∞ | NAN |
- real
logGamma
(real x);
- Natural logarithm of gamma function.
Returns the base e (2.718...) logarithm of the absolute
value of the gamma function of the argument.
For reals,
logGamma
is equivalent to log(fabs(gamma(x))).
Special Values
| x |
logGamma
(x) |
|---|
| NAN | NAN |
| integer <= 0 | +∞ |
| ±∞ | +∞ |
- real
beta
(real x, real y);
- Beta function
The
beta
function is defined as
beta
(x, y) = (Γ(x) Γ(y))/Γ(x + y)
- real
betaIncomplete
(real aa, real bb, real xx);
- Incomplete beta integral
Returns incomplete beta integral of the arguments, evaluated
from zero to x. The regularized incomplete beta function is defined as
betaIncomplete
(a, b, x) = Γ(a+b)/(Γ(a) Γ(b)) *
∫0x ta-1(1-t)b-1 dt
and is the same as the the cumulative distribution function.
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
betaIncompleteCompl(a, b, x ) =
betaIncomplete
( b, a, 1-x )
The integral is evaluated by a continued fraction expansion
or, when b*x is small, by a power series.
- real
betaIncompleteInv
(real aa, real bb, real yy0);
- Inverse of incomplete beta integral
Given y, the function finds x such that
betaIncomplete(a, b, x) == y
Newton iterations or interval halving is used.
- real
gammaIncomplete
(real a, real x);
real
gammaIncompleteCompl
(real a, real x);
- Incomplete gamma integral and its complement
These functions are defined by
gammaIncomplete
= ( ∫0x e-t ta-1 dt )/ Γ(a)
gammaIncompleteCompl(a,x) = 1 -
gammaIncomplete
(a,x)
= (∫x∞ e-t ta-1 dt )/ Γ(a)
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.
- real
gammaIncompleteComplInv
(real a, real p);
- Inverse of complemented incomplete gamma integral
Given a and y, the function finds x such that
gammaIncompleteCompl( a, x ) = p.
Starting with the approximate value x = a t3, where
t = 1 - d - normalDistributionInv(p) sqrt(d),
and d = 1/9a,
the routine performs up to 10 Newton iterations to find the
root of incompleteGammaCompl(a,x) - p = 0.
- real
digamma
(real x);
- Digamma function
The
digamma
function is the logarithmic derivative of the gamma function.
digamma
(x) = d/dx logGamma(x)
|
|