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

Two F-distribution functions

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

  • Two F-distribution functions

    ' Two F-distribution functions.
    '
    ' The first is based on series expansions to obtain the probability.

    ' The second is based on an incomplete beta function using a relation
    ' to the F-distribution. A separate gamma function is also used by this
    ' function.

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

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




    [This message has been edited by Erik Christensen (edited May 27, 2003).]

  • #2
    F-distribution version 1.
    Code:
    #COMPILE EXE
    #REGISTER NONE
    #DIM ALL
    
    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 necessary 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)) ' PI
            ' calculate ratio of factorials (= gamma of argument + 1) in 26.6.8
            ' Utilizes the principles in Functions 6.1.12, 6.1.8 and 6.1.1
            ' This calculation is accurate.
            X = 1 + (NY2 - 1) / 2
            LOCAL Ratio# : Ratio = 1 / SQR(4*ATN(1)) ' PI
            DO WHILE X>1.5
                DECR X : Ratio = Ratio * X / (X-.5)
            LOOP
            P5 = P5 * Ratio
    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 PBMAIN
        LOCAL Result$,P#,t$,F#,dft&,dfn&,i1&,i2&
        Result$ = INPUTBOX$("Enter F-value,dft,dfn","Get p-value")
        i1&=INSTR(Result$,",")
        i2&=INSTR(i1&+1,Result$,",")
        dft=VAL(MID$(Result$,i1&+1,i2&-i1&-1))
        dfn=VAL(MID$(Result$,i2&+1))
        P = ProbFromF(VAL(LEFT$(Result$,i1&-1)),dft,dfn)
        MSGBOX "P="+FORMAT$(P," ##.#####" )
    END FUNCTION
    ------------------


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

    Comment


    • #3
      F-distribution version 2. (code changed May 27, 2003) The
      previous version produced occational errors. Sorry. This version
      is much better. It allows a very large range of degrees of
      freedom. <IMG SRC="http://www.powerbasic.com/support/forums/smile.gif" ALT="smile">
      Code:
      #COMPILE EXE
      #REGISTER NONE
      #DIM ALL
      '
      FUNCTION gammln(num AS DOUBLE) AS DOUBLE 'logGamma 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 PfromF(BYVAL F AS DOUBLE,BYVAL dfn AS DOUBLE,BYVAL dfd AS DOUBLE) AS DOUBLE
          LOCAL ifault&,betain#
          CALL IncomplBeta(dfd/(dfd+dfn*F),dfd*.5#,dfn*.5#,betain,ifault)
          IF ifault = 0 THEN FUNCTION = betain ELSE FUNCTION = 1
      END FUNCTION
      '
      FUNCTION PBMAIN
          LOCAL Result$,P#,F#,dfn#,dfd#,i1&,i2&
          Result$ = INPUTBOX$("Enter F-value,dfn,dfd","Get p-value")
          i1&=INSTR(Result$,",")
          i2&=INSTR(i1&+1,Result$,",")
          F = VAL(LEFT$(Result$,i1&-1))
          dfn=VAL(MID$(Result$,i1&+1,i2&-i1&-1))
          dfd=VAL(MID$(Result$,i2&+1))
          P = PfromF(F,dfn,dfd)
          MSGBOX "P="+FORMAT$(P," ##.#####" ),,"Result:"
      END FUNCTION

      [This message has been edited by Erik Christensen (edited May 28, 2003).]

      Comment


      • #4
        The above code has been renewed. The present version is much
        better. It also surpasses that presented in NUMERICAL RECIPES,
        which is not precise.

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


        [This message has been edited by Erik Christensen (edited May 28, 2003).]

        Comment

        Working...
        X