From: Jim Weigang <jimw@chilton.com>
Newsgroups: comp.lang.apl
Date: 19 Nov 95 18:01:22 GMT
Subject: Re: F-distribution

Stephen Brooks wrote on 16 Nov 1995:

I am looking for an APL function that will return critical values of the F-distribution based on degrees of freedom in numerator and denominator and degree of confidence. Can anyone help?

Roger Hui replied on 17 Nov 1995:

What I used to do (and it worked quite well), is to find in Abramowitz and Stegun's "Handbook of Mathematical Functions" a suitable series for computing the probability given the F value, then invert using Newton iteration.

Beware that some of the F distribution approximations in the Handbook of Mathematical Functions aren't very accurate. For example, formula 26.6.15 is not very good when the degrees of freedom are small, and it can cause large errors when computing the inverse function. (E.g., P(F|3,5)=.999 being solved as F=44.15 instead of F=33.20.) A more accurate approach is to express the F distribution in terms of the incomplete beta function (formula 26.6.2) and use the continued fraction expansion in formula 26.5.8 for the inc. beta. The INTFD function below does this.

INVFD is the function Stephen wants: it computes the inverse of INTFD, taking a probability (confidence) P and finding the X value such that P = DF INTFD X.

Jim

(I provided some code, which quickly became obsolete.)


From: Jim Weigang <jimw@chilton.com>
Newsgroups: comp.lang.apl
Date: 20 Nov 95 16:52:11 GMT
Subject: Re: F-distribution

Roger Hui wrote on 19 Nov 1995:

First, Formula 2.6.15 may be approximate, but it should provide a better initial estimate for your Newton iteration function than the initial 1 that it currently uses.

Good point. I've put an approximation formula in function INVFDX below and used it to pick the starting point. For generating random numbers with an F distribution, this approximation alone may be adequate. (But beware of the limitations described on line [2].)

Second, There is some advantage to making the FDIST function work on an array of F values instead of scalars.

Yes, of course. The modified code appended below allows this. (And it has some additional changes to avoid overflow with large DF.)

The next step in this stone soup is for someone to point out how useful it would be if the program allowed the two degrees of freedom to be arrays as well, with all three arguments combined in outer product fashion. Left as an exercise for the reader!

Jim

Function calling tree:


INVFD        - Inverse F distribution by Newton's method
|
|-INVFDX     - Approximate inverse to the F distribution
|
|-INTFD      - Integral of F distribution from 0 to {omega}, with
| |            {alpha} = df1,df2
| '-INCBETA  - Incomplete beta function Ix(a,b), with {alpha} = a,b
|   |          and {omega} = x (array)
|   |-TRIM   - Replaces small values with TINY
|   |
|   '-LGAMMA - Log (ln) of the gamma function
|
'-FDIST      - F distribution (non-cumulative) with {alpha} = df1,df2
  |            and {omega} = x
  '-LGAMMA   - (see above)


     {del} X{<-}DF INVFD P;D;E;I;L;R;S;Y
[1]    @Inverse F distribution by Newton's method
[2]    @ Arguments:
[3]    @    P   - array of probability (confidence) values
[4]    @    DF  - numerator and denominator degrees of freedom
[5]    @ Result:
[6]    @    X   - value such that DF INTFD X is P
[7]    @
[8]    R{<-}{rho}P & P{<-},P      @ work with vector
[9]    L{<-}X{<-}(2{max}DF)INVFDX P{min}.999 @ approximate answer
[10]   E{<-}1E{neg}8              @ stopping precision
[11]  L1:Y{<-}(DF INTFD X)-P      @   g(x) = F(x)-P
[12]   D{<-}DF FDIST X            @   g'(x)
[13]   X{<-}X-S{<-}Y{divide}D     @   next approximation
[14]   I{<-}(X{<=}0)/{iota}{rho}X @   If any X{<=}0,
[15]   X[I]{<-}L[I]{divide}10     @     just move closer to 0
[16]   L{<-}X                     @   remember previous X
[17]   {->}({or}/E<|S)/L1         @ repeat if any steps were significant
[18]   X{<-}R{rho}X
     {del}

     {del} F{<-}DF INVFDX P;H;L;T;W;YP
[1]    @Approximate inverse to the F distribution
[2]    @ Requirements:  ^/DF>1  and  P<1
[3]    @ From formulas 26.6.16, 26.5.22, and 26.2.23 of the
[4]    @   Handbook of Mathematical Functions
[5]    T{<-}({ln}(1-P)*{neg}2)*.5
[6]    YP{<-}2.515517+T{times}0.802853+T{times}0.010328
[7]    YP{<-}T-YP{divide}1+T{times}1.432788+T{times}0.189269+T{times}0.001308
[8]    H{<-}2{divide}+/{divide}DF-1
[9]    L{<-}({neg}3+YP*2){divide}6
[10]   W{<-}(YP{times}(0{max}H+L)*.5){divide}H
[11]   W{<-}W-(-/{divide}DF-1){times}L+(5{divide}6)-2{divide}3{times}H
[12]   F{<-}*2{times}W
     {del}

     {del} Z{<-}D INTFD X
[1]    @Integral of F distribution from 0 to {omega}, with {alpha} = df1,df2
[2]    Z{<-}1-(D[2 1]{divide}2)INCBETA D[2]{divide}D[2]+D[1]{times}X
     {del}

     {del} I{<-}AB INCBETA X;A;B;C;D;F;H;J;K;M;N;P;R;S;Z;TINY
[1]    @Incomplete beta function Ix(a,b), with {alpha} = a,b and {omega} = {+
   +}x (array)
[2]    @ From the Handbook of Mathematical Functions, formula 26.6.8
[3]    @
[4]    TINY{<-}1E{neg}30           @ tiny value used by TRIM
[5]    P{<-}1E{neg}15              @ precision desired
[6]    H{<-}{rho}X & X{<-},X       @ work with vector X
[7]    R{<-}X>(AB[1]+1){divide}2++/AB @ 1s mark "large" Xs
[8]    K{<-}R/{iota}{rho}R
[9]    X[K]{<-}1-X[K]              @ use symmetry relation for faster{+
                                     +} convergence
[10]   A{<-}AB[1+R]                @ AB[1] for small X, AB[2] for large
[11]   B{<-}AB[1+~R]               @ the reverse
[12]   Z{<-}X{/=}0                 @ special-case zeros
[13]   X{<-}Z/X & A{<-}Z/A & B{<-}Z/B
[14]
[15]   @ Use modified Lentz algorithm to evaluate the continued fraction:
[16]   @  (As described in "Numerical Recipes in C", 2nd ed, pg. 171)
[17]   C{<-}{divide}TINY
[18]   D{<-}F{<-}1
[19]   M{<-}J{<-}0
[20]  L1:J{<-}J+1                  @ Loop for each term, j=1 2 3...
[21]   {->}(2|J)/L2
[22]   M{<-}M+1                    @   bump the half-J counter
[23]   N{<-}M{times}B-M            @   even c.f. term
[24]   {->}L3
[25]  L2:N{<-}-(A+M){times}A+B+M   @   odd term
[26]  L3:N{<-}(N{times}X){divide}(A+J-1){times}A+J @   numerator for c.f. term
[27]   D{<-}{divide}TRIM 1+N{times}D
[28]   C{<-}TRIM 1+N{divide}C
[29]   F{<-}F{times}S{<-}C{times}D @   next approximation
[30]   {->}({or}/P{<=}|S-1)/L1     @ repeat if any changes are significant
[31]
[32]   C{<-}*(B{times}{ln}1-X)+(A{times}{ln}X)+(LGAMMA A+B)-(LGAMMA A)+{+
   +}LGAMMA B
[33]   I{<-}F{times}C{divide}A
[34]   I{<-}Z\I                    @ insert zeros we deleted
[35]   I[K]{<-}1-I[K]              @ undo symmetry transformation
[36]   I{<-}H{rho}I                @ restore original shape
     {del}

     {del} Z{<-}TRIM A;B
[1]    @Replaces small values with TINY
[2]    B{<-}TINY<|A
[3]    Z{<-}(A{times}B)+TINY{times}~B
     {del}

     {del} Z{<-}LGAMMA B;I;S;T
[1]    @Log (ln) of the gamma function
[2]    S{<-}{rho}B & Z{<-}B{<-},B  @ work with vector
[3]    I{<-}(B{<=}50)/{iota}{rho}B@
[4]    Z[I]{<-}{ln}!B[I]-1         @ use APL for small B
[5]    @ Use series (which is very accurate for B{>=}19) for large B:
[6]    I{<-}(B>50)/{iota}{rho}B
[7]    B{<-}B[I]
[8]    T{<-}-/(B{jot}.*-1 3 5 7){divide}(({rho}B),4){rho}12 360 1260 1680
[9]    Z[I]{<-}T+(0.5{times}{ln}{pitimes}2)+(({ln}B){times}B-0.5)-B
[10]   Z{<-}S{rho}Z
     {del}

     {del} Z{<-}DF FDIST X;A;B;C;D
[1]    @F distribution (non-cumulative) with {alpha} = df1,df2 and {omega} = x
[2]    A{<-}(LGAMMA .5{times}+/DF)-+/LGAMMA DF{divide}2
[3]    C{<-}{divide}/DF
[4]    A{<-}A+({ln}C){times}DF[1]{divide}2
[5]    B{<-}{neg}1+DF[1]{divide}2
[6]    D{<-}{neg}.5{times}+/DF
[7]    @ {take} All the above could be hoisted out of a root-finding loop
[8]    Z{<-}*A+(B{times}{ln}X)+D{times}{ln}1+C{times}X
     {del}

Home Page