'********************************************************
'*      Calculate Incomplete Beta Function Ix(a,b)      *
'* ---------------------------------------------------- *
'* SAMPLE RUN:                                          *
'*                                                      *
'*  Incomplete Beta Function Ix(a,b)                    *
'*                                                      *
'*  A = .5                                              *
'*  B = 5                                               *
'*  x = .2                                              *
'*                                                      *
'*  Y = .8550723997280966                               *
'*                                                      *
'* ---------------------------------------------------- *
'* Reference:                                           *
'* "Numerical Recipes, By W.H. Press, B.P. Flannery,    *
'*  S.A. Teukolsky and T. Vetterling, Cambridge         *
'*  University Press, 1986" [BIBLI 08].                 *
'*                                                      *
'*                  Basic Release By J-P Moreau, Paris. *
'*                          (www.jpmoreau.fr)           *
'********************************************************
'PROGRAM TEST_BETAI

DEFDBL A-H, O-Z
DEFINT I-N

' A,B,X: Double

  A = .5#: B = 5#: X = .2#

  GOSUB 1000 'call BETAI(A,B,X)

  CLS
  PRINT
  PRINT " Incomplete Beta Function Ix(A,B)"
  PRINT
  PRINT " A = "; A
  PRINT " B = "; B
  PRINT " X = "; X
  PRINT
  PRINT " Y = "; BETAI
  PRINT

END 'of main program


1000 'Function BETAI(A,B,X)
' returns the incompleteBeta function Ix(a,b)
' BT: Double
  IF X < 0# OR X > 1# THEN
    PRINT " BETAI: Bad argument X (must be 0<=X<=1)."
    RETURN
  END IF

  IF X = 0# OR X = 1# THEN
    BT = 0#
  ELSE
    XX = A + B: GOSUB 2000: GAB = GAMMLN
    XX = A: GOSUB 2000: GA = GAMMLN
    XX = B: GOSUB 2000: GB = GAMMLN
    BT = EXP(GAB - GA - GB + A * LOG(X) + B * LOG(1# - X))
  END IF

  IF X < (A + 1#) / (A + B + 2#) THEN
    GOSUB 1500
    BETAI = BT * BETACF / A
    RETURN
  ELSE
    AS1 = A: BS = B: XS = X  'save current A, B, X
    TMP = A: A = B: B = TMP: X = 1# - X: GOSUB 1500
    A = AS1: B = BS: X = XS  'restore current A, B, X
    BETAI = 1# - BT * BETACF / B
    RETURN
  END IF
RETURN

1500 'Function BETACF(A,B,X)
' continued fraction for incomplete Beta function, used by BETAI
  ITMAX = 100
  EPS = .0000003#
' AM,BM,AZ,QAB,QAP,QAM,BZ,EM,TEM,D,AP,BP,APP,BPP,AOLD: Double

  AM = 1#
  BM = 1#
  AZ = 1#
  QAB = A + B
  QAP = A + 1#
  QAM = A - 1#
  BZ = 1# - QAB * X / QAP
  FOR M = 1 TO ITMAX
    EM = M
    TEM = EM + EM
    D = EM * (B - M) * X / ((QAM + TEM) * (A + TEM))
    AP = AZ + D * AM
    BP = BZ + D * BM
    D = -(A + EM) * (QAB + EM) * X / ((A + TEM) * (QAP + TEM))
    APP = AP + D * AZ
    BPP = BP + D * BZ
    AOLD = AZ
    AM = AP / BPP
    BM = BP / BPP
    AZ = APP / BPP
    BZ = 1#
    IF ABS(AZ - AOLD) < EPS * ABS(AZ) THEN GOTO 1
  NEXT M
  PRINT " BETACF: A or B too big, or ITMAX too small."
  RETURN
1 BETACF = AZ
RETURN
  
2000 'Function GAMMLN(XX)
' returns the value Ln(Gamma(XX) for XX > 0
' STP,HALF,ONE,FPF,XTMP,TMP,SER: Double
DIM COF(6)

  COF(1) = 76.18009173#: COF(2) = -86.50532033#: COF(3) = 24.01409822#
  COF(4) = -1.231739516#: COF(5) = .00120858003#: COF(6) = -.00000536382#
  STP = 2.50662827465#
  HALF = .5: ONE = 1#: FPF = 5.5
  XTMP = XX - ONE
  TMP = XTMP + FPF
  TMP = (XTMP + HALF) * LOG(TMP) - TMP
  SER = ONE
  FOR J = 1 TO 6
    XTMP = XTMP + ONE
    SER = SER + COF(J) / XTMP
  NEXT J
  GAMMLN = TMP + LOG(STP * SER)
RETURN

' end of file ibeta.bas