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-distribution function - best version

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

  • Inverse F-distribution function - best version

    ' Inverse F-distribution function - best version.
    '
    ' 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 derived from a
    ' relationship with the incomplete beta function. A gamma function and an
    ' effective bisection method root finder is also used. The speed is not
    ' very impressive, but you can use very large degrees of freedom without any
    ' problems.
    '
    ' Thanks to Tony Burcham for a fine gamma function in this forum.
    '
    ' Erik Christensen ---- e.chr@email.dk ---- March 29, 2003
    '
    ' P.S. May 27, 2003: The incomplete beta function in the original code gave
    ' unpredictable errors - unfortunately. It has therefore been replaced by
    ' the present incomplete beta function, which appears to be reliable.
    Code:
    #COMPILE EXE
    #REGISTER NONE
    #DIM ALL
    '
    GLOBAL DN#,DD#,PP#
    '
    FUNCTION gammln(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
    '
    SUB IncomplBeta(BYVAL X AS DOUBLE,BYVAL P AS DOUBLE,BYVAL Q AS DOUBLE,BYREF betain AS DOUBLE,BYREF ifault AS LONG)
    '
    '     Derived from FORTRAN code based on:
    '
    '     algorithm AS 63  appl. statist. (1973), vol.22, no.3
    '
    '     computes incomplete beta function ratio for arguments
    '     x between zero and one, p and q positive.
    '
          LOCAL indx&,beta##
          LOCAL zero##,one##,acu##
          LOCAL psq##,cx##,xx##,pp##,qq##,term##,ai##,ns&,rx##,temp##
    '
    '     define accuracy and initialise
    '
          zero = 0.0E0## : one = 1.0E0## : acu = 0.1E-13##
          betain=x
    '
    '     test for admissibility of arguments
    '
          ifault=1&
          IF(p <= zero OR q <= zero) THEN EXIT SUB
          ifault=2&
          IF(x < zero OR x > one) THEN EXIT SUB
          ifault=0&
          IF(x = zero OR x = one) THEN EXIT SUB
    '
    '     calculate log of complete beta
    '
          beta = gammln(p)+gammln(q)-gammln(p+q)
    '
    '     change tail if necessary and determine s
    '
          psq=p+q
          cx=one-x
          IF(p >= psq*x) GOTO 1
          xx=cx
          cx=x
          pp=q
          qq=p
          indx = 1&
          GOTO 2
        1 xx=x
          pp=p
          qq=q
          indx = 0&
        2 term=one
          ai=one
          betain=one
          ns=qq+cx*psq
    '
    '     user Soper's reduction formulae.
    '
          rx=xx/cx
        3 temp=qq-ai
          IF(ns = 0&) THEN rx=xx
        4 term=term*temp*rx/(pp+ai)
          betain=betain+term
          temp=ABS(term)
          IF(temp <= acu AND temp <= acu*betain) THEN GOTO 5
          ai=ai+one
          ns=ns-1&
          IF(ns >= 0&) THEN GOTO 3
          temp=psq
          psq=psq+one
          GOTO 4
    '
    '     calculate result
    '
        5 betain=betain*EXP(pp*LOG(xx)+(qq-one)*LOG(cx)-beta)/pp
          IF indx = 1& THEN betain=one-betain
    END SUB
    '
    FUNCTION FCDFMP(BYVAL F AS DOUBLE) AS DOUBLE
    ' Calculate difference between original P and
    ' calculated P based on the current F-value.
          LOCAL a#,b#,x#,IER&,P#
          a=DD*.5# : b=DN*.5# : x = DD/(DD+DN*F)
          CALL IncomplBeta(x,a,b,P,IER)
          FUNCTION = P-PP
    END FUNCTION
    '
    SUB BRAK(BYREF A AS DOUBLE,BYREF B AS DOUBLE)
    '  Based on ZBRAC in Numerical Recipes by Press, Flannery,
    '  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,.00001) ' .00001 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 May 28, 2003).]

  • #2
    Just to inform you that this function has been improved - see above.

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

    Comment

    Working...
    X