IncompleteBeta Function |
Unit
QESBPCSMath
Declaration
Function IncompleteBeta(X: Extended; P, Q: Extended): Extended;
Description
The Incomplete Beta function is the probability that a random variable from a Beta distribution having parameters P and Q will be less than or equal to X. Accuracy: Gives about 17 digits. Adapted From Collected Algorithms from CACM, Algorithm 179 - Incomplete Beta Function Ratios, Oliver G Ludwig.
Category
Arithmetic Routines for FloatsImplementation
function IncompleteBeta (X: Extended; P, Q: Extended): Extended; { Adapted From Collected Algorithms from CACM Algorithm 179 - Incomplete Beta Function Ratios Oliver G Ludwig } const Epsilon: Extended = 0.5E-18; MaxIterations = 1000; var FinSum, InfSum, Temp, Temp1, Term, Term1, QRecur, Index: Extended; I: Integer; Alter: Boolean; begin if not FloatIsPositive (P) or not FloatIsPositive (Q) or FloatIsNegative (X) or GreaterFloat (X, 1) then begin raise EMathError.Create (rsNotDefinedForValue) end; if (X = 0) or (X = 1) then Result := X else begin // Interchange arguments if necessary to get better convergence if X <= 0.5 then Alter := False else begin Alter := True; SwapXY (P, Q); X := 1 - X; end; // Recurs on the (effective) Q until the Power Series doesn't alternate FinSum := 0; Term := 1; Temp := 1 - X; QRecur := Q; Index := Q; repeat Index := Index - 1; if Index <= 0 then Break; QRecur := Index; Term := Term * (QRecur + 1) / (Temp * (P + QRecur)); FinSum := FinSum + Term; until False; // Sums a Power Series for non-integral effective Q and yields unity for integer Q InfSum := 1; Term := 1; for I := 1 to MaxIterations do begin if Term <= Epsilon then Break; Index := I; Term := Term * X * (Index - QRecur) * (P + Index - 1) / (Index * (P + Index)); InfSum := InfSum + Term; end; // Evaluates Gammas Temp := Gamma (QRecur); Temp1 := Temp; Term := Gamma (QRecur + P); Term1 := Term; Index := QRecur; repeat Temp1 := Temp1 * Index; Term1 := Term1 * (Index + P); Index := Index + 1; until Index >= Q - 0.5; Temp := XtoY (X, P) * (InfSum * Term / (P * Temp) + FinSum * Term1 * XtoY (1 - X, Q) / (Q * Temp1)) / Gamma (P); if Alter then Result := 1 - Temp else Result := Temp end; End; |
|