Implementation of the gamma and beta functions, and their integrals.
Stephen L. Moshier (original C code). Conversion to D by Don Clugston
- real MAXGAMMA [manifest] ¶
-
The maximum value of x for which gamma(x) < real.infinity.
- real sgnGamma(real x) ¶
-
Returns -1 if Γ(x) < 0, +1 if Γ(x) > 0,
NAN if sign is indeterminate.
- real gamma(real 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∞ t
z-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) ¶
-
The beta function is defined as
beta(x, y) = (Γ(x) Γ(y))/Γ(x + y)
- real betaIncomplete(real aa, real bb, real xx) ¶
-
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 t
a-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 t
a-1 dt )/ Γ(a)
gammaIncompleteCompl(a,x) = 1 - gammaIncomplete(a,x)
= (
∫x∞ e
-t t
a-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 t
3, 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) ¶
-
The digamma function is the logarithmic derivative of the gamma function.
digamma(x) = d/dx logGamma(x)