Announcement

Collapse

Forum Guidelines

This forum is for finished source code that is working properly. If you have questions about this or any other source code, please post it in one of the Discussion Forums, not here.
See more
See less

Inverse F-function

Collapse
X
  • Filter
  • Time
  • Show
Clear All
new posts

  • Inverse F-function

    ' 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
    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).]

  • #2
    Hi, Eric,

    I have not read my e-mail in two days, but since the last time I responded to your posting on the f-distribution, I have had a chance to try out the current one -- inverse f-distribution.

    In a few words, it is remarkably accurate across the range of degrees of freedom, from 1,1 to 24,120, which I have tested.

    Great work, great contribution.

    George.

    ------------------
    GCO

    Comment

    Working...
    X