F distribution on the 41C

[Dieter]

(12-20-2015 11:15 PM)Namir Wrote:  I am implementing the equations in Abramowitz & Stegun's handbook of math.

OK, what about A&S 26.6.2 then?
Here the Fisher CDF is defined by means of the regularized Beta function. ;-)

(12-20-2015 11:15 PM)Namir Wrote:  Show me some VB code using the Beta function?

The following code uses the continued fraction method according to A&S 26.5.8. It is based on the modified Lentz algorithm as suggested in "Numerical Recipes in C", chapter 6.4, with some changes and simplifications. For fast convergence and improved accuracy x should be low, so either x = df2/(df2 + F*df1) is used, or 1–x = F*df1/(df2 + x*df1) with swapped parameters a and b.

I have split the functionality into four independent routines so that these can be used for other applications as well. For instance for Student's t or the Binomial distribution.

Note that here the upper tail integral Q(F) is calculated. For the lower tail P(F) simply exchange ibeta and 1–ibeta in the two if-then-else branches of FCDF.

Code:

Const Pi = 3.141592653589793

Function FCDFQ(ByVal F, ByVal df1, ByVal df2)

x1 = df2 / (df2 + F * df1)
x2 = F * df1 / (df2 + F * df1)

' calculate regularized Beta using min(x, 1-x)

If x1 < x2 Then
   FCDFQ = ibeta(x1, df2 / 2, df1 / 2)
Else
   FCDFQ = 1 - ibeta(x2, df1 / 2, df2 / 2)
End If

End Function


Function ibeta(ByVal x, ByVal a, ByVal b)

tiny = 1E-300   ' a value close to smallest floating point number
c = 1
d = 1 - (a + b) * x / (a + 1)
If (Abs(d) < tiny) Then d = tiny
d = 1 / d
h = d
m = 1
Do
   amm = a + m + m
   aa = m * (b - m) * x / ((amm - 1) * amm)
   d = 1 + aa * d
   If (Abs(d) < tiny) Then d = tiny
   c = 1 + aa / c
   If (Abs(c) < tiny) Then c = tiny
   d = 1 / d
   h = (d * c) * h
   aa = -(a + m) * (a + b + m) * x / (amm * (amm + 1))
   d = 1 + aa * d
   If (Abs(d) < tiny) Then d = tiny
   c = 1 + aa / c
   If (Abs(c) < tiny) Then c = tiny
   d = 1 / d
   adj = d * c
   h = h * adj
   m = m + 1
Loop Until (adj = 1) Or (m > 200)

If m > 200 Then
   ibeta = -1    ' flag error due to no convergence, should never be returned
Else
   ibeta = h * x ^ a * (1 - x) ^ b / a / betafn(a, b)
End If

End Function


Function hGamma(ByVal k)  ' evaluates Gamma for k=0.5, 1, 1.5, 2, 2.5, 3, ...

k = Int(2 * k + 0.5) / 2 ' round input to integer or half-integer
If k = Int(k) Then g = 1 Else g = Sqr(Pi)

Do While k > 1.25
  k = k - 1
  g = g * k
Loop

hGamma = g

End Function


Function betafn(ByVal v, ByVal w)

betafn = hGamma(v) * hGamma(w) / hGamma(v + w)

End Function

This is what I currently use in Excel for calculating the upper F integral.

Dieter