tango.math.GammaFunction

License:

BSD style: see license.txt

Authors:

Stephen L. Moshier (original C code). Conversion to D by Don Clugston

real MAXGAMMA [const] ¶#

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) = ∫₀^∞ 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) ¶#

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)) * ∫₀^x 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 = ( ∫₀^x 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³, 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)