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 using Newton-Raphson root finding

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

  • Inverse F-distribution using Newton-Raphson root finding

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