| 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; |
|
|