'******************************************************** '* 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