' Inverse F-distribution using Newton-Raphson root finding.

'

' 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

' using a gamma function and for the Newton-Raphson method,

' which is a highly effective iterative, derivative based, root

' finding method. Although numerical calculation of the

' derivative may be used, the method works much better for

' functions whose derivatives can be evaluated explicitly. This

' is the case for the cumulative F-distribution. The

' Newton-Raphson method uses the intersection of the tangent of

' the current function point with the X-axis as the next root

' approximation. Since for given degrees of freedom the

' cumulative F-distribution function is monotonous, the

' Newton-Raphson method will have no difficulty in finding the root.

' However, for very large degrees of freedom (> 1400) the derivative

' can become zero. If so, no convergence will take place.

'

' Thanks to Tony Burcham for a fine gamma function in this forum.

'

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

'

' P.S. Small improvements made March 29, 2003.

[This message has been edited by Erik Christensen (edited March 29, 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

' using a gamma function and for the Newton-Raphson method,

' which is a highly effective iterative, derivative based, root

' finding method. Although numerical calculation of the

' derivative may be used, the method works much better for

' functions whose derivatives can be evaluated explicitly. This

' is the case for the cumulative F-distribution. The

' Newton-Raphson method uses the intersection of the tangent of

' the current function point with the X-axis as the next root

' approximation. Since for given degrees of freedom the

' cumulative F-distribution function is monotonous, the

' Newton-Raphson method will have no difficulty in finding the root.

' However, for very large degrees of freedom (> 1400) the derivative

' can become zero. If so, no convergence will take place.

'

' Thanks to Tony Burcham for a fine gamma function in this forum.

'

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

'

' P.S. Small improvements made March 29, 2003.

Code:

#COMPILE EXE #REGISTER NONE #DIM ALL ' #INCLUDE "win32api.inc" ' FUNCTION LGam(BYVAL 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 EXT, NY1 AS EXT, NY2 AS EXT) AS EXT ' 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 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& ' F = F-value ' NY1 = degrees of freedom for numerator ' NY2 = degrees of freedom for denominator PIPI = 1 IF NY1 < 1 OR NY1 <> FIX(NY1) THEN GOTO 3720 ' Error IF NY2 < 1 OR NY2 <> FIX(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 = LGam(X+1) X = (NY2 - 2) / 2 : lx = lx - LGam(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 ' SUB funcd(BYVAL p AS EXT, BYVAL DN AS EXT, BYVAL DD AS EXT, _ BYVAL F AS EXT, BYREF fn AS EXT, BYREF df AS EXT, BYREF FirstFl AS LONG) STATIC fac## ' Calculate derivative of 26.6.1 IF ISTRUE FirstFl& THEN ' calculate first constant part, also using 6.2.2. fac = DN^(.5##*DN)*DD^(.5##*DD) / EXP(LGam(.5##*DN)+LGam(.5##*DD)-LGam(.5##*DN+.5##*DD)) FirstFL& = %FALSE END IF df = F^(.5##*(DN-2##)) * (DD+DN*F)^(-.5##*(DN+DD)) * fac * (-1##) ' = full derivative. ' Function value: Difference between currently calculated P and original P. fn = ProbFromF(F,CLNG(DN),CLNG(DD)) - P END SUB ' FUNCTION InverseF(BYVAL p AS EXT, BYVAL dfNum AS EXT, BYVAL dfDenom AS EXT, _ BYVAL Prec AS EXT, BYVAL FirstFl AS LONG) AS DOUBLE LOCAL jmax&,j&,df##,dx##,fn##,F##,fnPrev## IF p>0.5 THEN MSGBOX "P must be less than 0.5",%MB_ICONERROR,"Problem" : EXIT FUNCTION jmax = 20 ' Max number of iterations FUNCTION = 1 : F = 1 ' Starting value ' Newton Raphson iteration - very fast convergence. FOR j = 1 TO jmax fnPrev = fn CALL funcd(p, dfNum, dfDenom, F, fn, df, FirstFl) IF fn = fnPrev THEN MSGBOX "Does not converge.",%MB_ICONERROR,"Problem" : EXIT FUNCTION IF df = 0 THEN MSGBOX "Derivative is zero - cannot converge.",%MB_ICONERROR,"Problem" : EXIT FUNCTION dx = fn/df : F = F-dx : FUNCTION = F IF ABS(dx)<Prec THEN EXIT FUNCTION NEXT MSGBOX "Exceeded maximum iteration number",%MB_ICONERROR,"Problem" END FUNCTION ' FUNCTION PBMAIN LOCAL Result$,i1&,i2&,p##,dfnum##,dfDenom##,F##,Prec##,FirstFl& 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)) Prec = .00001 ' precision FirstFl = %TRUE F = InverseF(p,dfNum,dfDenom,Prec,FirstFl) MSGBOX "F="+FORMAT$(F,"#######.#####"),%MB_ICONINFORMATION,"Inverse F from: p="+STR$(p)+" dfn="+STR$(dfNum)+" dfd="+STR$(dfDenom) END FUNCTION

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