' Inverse F-function.

'

' This code calculates F from the p-value (the probability) and the degrees

' of freedom for numerator and denominator.

'

' The program incorporates code for: the usual F-distribution, the gamma

' function and an effective bisection method root finder. The speed is not

' bad, but it may be increased somewhat using the Newton-Raphson method.

' This I will leave to you.

'

' Thanks to Roger L. Berger for his FORTRAN code and to Tony Burcham for

' a fine gamma function in this forum.

'

' Erik Christensen ---- e.chr@email.dk ---- March 8, 2003

------------------

[This message has been edited by Erik Christensen (edited March 08, 2003).]

'

' This code calculates F from the p-value (the probability) and the degrees

' of freedom for numerator and denominator.

'

' The program incorporates code for: the usual F-distribution, the gamma

' function and an effective bisection method root finder. The speed is not

' bad, but it may be increased somewhat using the Newton-Raphson method.

' This I will leave to you.

'

' Thanks to Roger L. Berger for his FORTRAN code and to Tony Burcham for

' a fine gamma function in this forum.

'

' Erik Christensen ---- e.chr@email.dk ---- March 8, 2003

Code:

#COMPILE EXE #REGISTER NONE #DIM ALL GLOBAL DN#,DD#,PP# FUNCTION LogGam(num AS DOUBLE) AS DOUBLE 'LogE-Gamma function ' Thanks to Tony Burcham who provided this fine function. LOCAL fl AS EXT LOCAL j AS LONG LOCAL k AS EXT LOCAL z AS EXT LOCAL zpwr AS EXT 'power of z LOCAL zsqr AS EXT 'z squared DIM cf(5) AS EXT cf(0) = 1## / 12## cf(1) = 1## / -360## cf(2) = 1## / 1260## cf(3) = 1## / -1680## cf(4) = 1## / 1188## FOR k = 0## TO 10## fl = fl - LOG(num + k ) NEXT k z = num + 11## zpwr = z 'z^1 zsqr = z * z FOR j = 0 TO 4 fl = fl + ( cf(j) / zpwr) zpwr = zpwr * zsqr NEXT j fl = fl + ((z - 0.5##) * LOG(z)) + 0.9189385332046727## - z FUNCTION = fl ' Gam = EXP(fl) END FUNCTION ' FUNCTION ProbFromF(F AS DOUBLE, DFnum AS LONG, DFdenom AS LONG) AS DOUBLE ' CALCULATES THE P-VALUE FROM F AND THE DEGREES OF FREEDOM. ' The probability is calculated using series expansions 26.6.4, 26.6.5, ' and 26.6.8. The reflexive relation 26.6.9 is also used. ' In: Handbook of mathematical Functions. Editors: Abramowitz M & Stegun IA, ' Dover, New York, 1972. ' Original code by Erik Christensen. LOCAL NY1#,NY2#,SUM#,PROD#,PIPI#,PNI#,XX#,NY1mod2#,NY2mod2#,T1#,T2#,IG&,FA#,I#,X#,lx# LOCAL PHI#,SINPHI#,COSPHI#,COS2#,SIN2#,BB1#,AA#,D5#,CC#,P5#,TT#,TS#,P# LOCAL flag& NY1 = DFnum: NY2 = DFdenom ' F = F-value ' NY1 = degrees of freedom for numerator ' NY2 = degrees of freedom for denominator PIPI = 1 IF NY1 < 1 OR NY1 <> INT(NY1) THEN GOTO 3720 ' Error IF NY2 < 1 OR NY2 <> INT(NY2) THEN GOTO 3720 ' Error IF F < 1E-11 THEN 3490 SUM = 1: PROD = 1: PNI = 2: XX = NY2 / (NY2 + NY1 * F) NY1mod2 = NY1 MOD 2: NY2mod2 = NY2 MOD 2 IF NY1mod2 + NY2mod2 > 1.5 THEN 3500 ' both df's odd IF NY2mod2 > NY1mod2 + .2 THEN 3420 IF NY1mod2 - NY2mod2 < .1 AND NY2 > NY1 THEN 3420 T1 = NY1: T2 = XX: IG = INT(NY2 / 2) - 1: FA = (1 - XX) ^ (NY1 / 2): FA = -FA: GOTO 3430 3420 T1 = NY2: T2 = 1 - XX: IG = INT(NY1 / 2) - 1: FA = XX ^ (NY2 / 2) 3430 IF IG < .5 THEN 3470 FOR I = 1 TO IG PROD = PROD * T1 * T2 / PNI: SUM = SUM + PROD: T1 = T1 + 2: PNI = PNI + 2 NEXT I 3470 PIPI = SUM * FA IF PIPI < 0 THEN PIPI = 1 + PIPI ' An important safeguard. 3490 GOTO 3720 ' ' This swapping is done to avoid NY2 becoming 1 while NY1 is >= 3 3500 IF NY1 > NY2 THEN SWAP NY1,NY2 : F = 1 / F : flag = 1 ' PHI = ATN(SQR(NY1 * F / NY2)): SINPHI = SIN(PHI): COSPHI = COS(PHI): PROD = COSPHI: SUM = COSPHI COS2 = COSPHI * COSPHI: SIN2 = SINPHI * SINPHI: BB1 = 0 IF NY2 = 1 THEN 3580 IF NY2 = 3 THEN 3570 FOR I = 5 TO NY2 STEP 2 PROD = PROD * COS2 * (I - 3) / (I - 2): SUM = SUM + PROD NEXT I 3570 BB1 = SUM * SINPHI 3580 AA = (2 / (4 * ATN(1))) * (PHI + BB1): PROD = 1: SUM = 1 IF NY1 = 1 THEN D5 = 0: GOTO 3700 IF NY1 = 3 THEN GOTO 3640 FOR I = 5 TO NY1 STEP 2 PROD = PROD * SIN2 * (I + NY2 - 4) / (I - 2): SUM = SUM + PROD NEXT I 3640 CC = SUM * SINPHI * COSPHI ^ NY2 : P5 = 2 / SQR(4 * ATN(1)) ' calculate (ratio of) factorials (gamma + 1) X = (NY2 - 1) / 2 : lx = LogGam(X+1) X = (NY2 - 2) / 2 : lx = lx - LogGam(X+1) P5 = P5 * EXP(lx) 3690 D5 = CC * P5 3700 PIPI = 1 - AA + D5 IF PIPI < 0 THEN PIPI = 1 + PIPI ' An important safeguard. IF flag = 1 THEN PIPI = 1 - PIPI 3720 P = PIPI FUNCTION = P END FUNCTION ' FUNCTION FCDFMP(BYVAL F AS DOUBLE) AS DOUBLE ' Calculate difference between original P and ' calculated P based on the current F-value. LOCAL DNU&,DDE& DNU = DN# : DDE = DD# FUNCTION = ProbFromF(F,DNU,DDE)-PP END FUNCTION ' SUB BRAK(BYREF A AS DOUBLE,BYREF B AS DOUBLE) 'C Based on ZBRAC in Numerical Recipes by Press, Flannery, 'C Teukolsky and Vetterling, New York: Cambridge. LOCAL ZERO#,F1#,F2#,TP#,FAC#,NTRY&,I& FAC=1.6E0 NTRY=60 ZERO=0.0E0 F1 = FCDFMP(A) F2 = FCDFMP(B) FOR I=1 TO NTRY IF(F1*F2 <= ZERO) THEN EXIT SUB IF(ABS(F1) < ABS(F2)) THEN TP = A A = MAX(ZERO,A + FAC*(A-B)) F1 = FCDFMP(A) B = TP ELSE TP = B B = MAX(ZERO,B + FAC*(B-A)) F2 = FCDFMP(B) A = TP END IF NEXT END SUB ' FUNCTION BISECT(BYVAL F AS DOUBLE,BYVAL A AS DOUBLE,BYVAL B AS DOUBLE,BYVAL XACC AS DOUBLE) AS DOUBLE ' Bisection method root finder. Value returned will be within XACC of a root. LOCAL JMAX&,J& LOCAL ZERO#,HALF#,TWO#,DX#,FHI#,FLO#,FMID#,XHI#,XLO#,XMID# ZERO=0.0E0 HALF=0.5E0 TWO=2.0E0 FLO = FCDFMP(A) IF(FLO = ZERO) THEN BISECT = A : EXIT FUNCTION END IF FHI = FCDFMP(B) IF(FHI = ZERO) THEN BISECT = B : EXIT FUNCTION END IF IF(FLO*FHI > ZERO) THEN MSGBOX "A ="+STR$(A)+" F(A) ="+STR$(FLO)+$CRLF+"B ="+STR$(B)+" F(B) ="+STR$(FHI) _ ,,"BISECT ERROR: Endpoints have same sign." EXIT FUNCTION END IF IF(FLO < ZERO) THEN XLO = A : XHI = B ELSE XLO = B : XHI = A FMID = FLO : FLO = FHI : FHI = FMID END IF DX = XHI-XLO IF(ABS(DX) <= XACC) GOTO 20 ' JMAX is number of bisections needed to get XACC accuracy. JMAX = 1 + FIX(LOG(ABS(DX)/XACC)/LOG(TWO)) FOR J = 1 TO JMAX DX = DX*HALF XMID = XLO+DX FMID = FCDFMP(XMID) IF(FMID = ZERO) THEN BISECT = XMID : EXIT FUNCTION ELSEIF(FMID > ZERO) THEN XHI = XMID : FHI = FMID ELSE XLO = XMID : FLO = FMID END IF NEXT 20 BISECT = XLO IF(FHI < ABS(FLO)) THEN BISECT = XHI END FUNCTION ' FUNCTION FINV(BYVAL P AS DOUBLE,BYVAL DFN AS DOUBLE,BYVAL DFD AS DOUBLE) AS DOUBLE ' Inverse F function calculated from P, DFN and DFD ' DFN and DFD being the numerator and denominator degrees of freedom. ' Based on FORTRAN code by Roger L. Berger 06/23/94 LOCAL A#,B#,FF# FUNCTION = 0 IF (P > 0) THEN DN = DFN : DD = DFD : PP = P A = 1 : B = 20 ' Starting range - may be widened if necessary. CALL BRAK(A,B) FUNCTION = BISECT(FF,A,B,.00005) ' .00005 is the current accuracy - may be changed. END IF END FUNCTION ' FUNCTION PBMAIN LOCAL Result$,i1&,i2&,p#,dfnum#,dfDenom#,F# Result$ = INPUTBOX$("Enter p, df-num, df-denom","Get F from: p, df-numerator and df-denominator") i1&=INSTR(Result$,",") ' get position i2&=INSTR(i1&+1,Result$,",") ' of commas ' extract parameters p = VAL(LEFT$(Result$,i1&-1)) dfNum = VAL(MID$(Result$,i1&+1,i2&-i1&-1)) dfDenom = VAL(MID$(Result$,i2&+1)) DN#=dfNum : DD#=dfDenom : PP#=p ' Initialize global variables. ' Get F F = FINV(p,dfNum,dfDenom) MSGBOX "F="+FORMAT$(F,"#######.#######"),,"Inverse F from: p="+STR$(p)+" dfn="+STR$(dfNum)+" dfd="+STR$(dfDenom) END FUNCTION

[This message has been edited by Erik Christensen (edited March 08, 2003).]

## Comment