IncompleteBeta Function
Returns the Incomplete Beta Ix(P, Q), where 0 <= X <= 1 and P and Q are positive.

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 Floats

Implementation

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;


HTML generated by Time2HELP
http://www.time2help.com