This is your lucky day!
Here is a Visual Basic implementation for the cumulative distribution function for the F statistics. It is based on the Handbook of Mathematical Function. I just coded it today for my upcoming Console Programmable RPN Calculation app for Windows. It handles odd and even degrees of freedom. I checked the results against the HP Prime and the function did well. The function uses an implementation of Gamma function, which I also include.
Code:
Function Fpdf(ByVal F As Double, ByVal Df1 As Integer, ByVal Df2 As Integer) As Double
' Use equations on page 946 (sections 26.6.4 and 26/6/5) of "Handbook of Mathematical Functions", 9th printing, 1970,
' by Milton Abramowitz And Irene Stegun
Dim ddf1, ddf2 As Double
Dim Prod, Sum As Double
Dim Theta, X0 As Double
Dim X, A, B As Double
Dim CosTheta, SinTheta As Double
Dim CosPwr, SinPwr As Double
X = Df2 / (Df2 + Df1 * F)
Sum = 0
If Df1 Mod 2 = 0 Then
If Df1 = 2 Then
Sum = 1
ElseIf Df1 = 4 Then
Sum = 1 + Df2 / 2 * (1 - X)
Else
Prod = Df2 * (Df2 + 2) / 2 / 4
Sum = 1 + Df2 / 2 * (1 - X) + Prod * (1 - X) ^ 2
For ddf1 = 8 To Df1 Step 2
Prod *= (Df2 + ddf1 - 4) / (ddf1 - 2)
Sum += Prod * (1 - X) ^ ((ddf1 - 2) \ 2)
Next
End If
F = X ^ (Df2 / 2) * Sum
ElseIf Df2 Mod 2 = 0 Then
If Df2 = 2 Then
Sum = 1
ElseIf Df2 = 4 Then
Sum = 1 + Df1 / 2 * X
Else
Prod = Df1 * (Df1 + 2) / 2 / 4
Sum = 1 + Df1 / 2 * X + Prod * X ^ 2
For ddf2 = 8 To Df2 Step 2
Prod *= (Df1 + ddf2 - 4) / (ddf2 - 2)
Sum += Prod * X ^ ((ddf2 - 2) \ 2)
Next
End If
F = 1 - (1 - X) ^ (Df1 / 2) * Sum
Else ' df1 and df2 are odd
Theta = Math.Atan(Math.Sqrt(Df1 / Df2 * F))
CosTheta = Math.Cos(Theta)
SinTheta = Math.Sin(Theta)
' Calculate A
If Df2 > 1 Then
Sum = CosTheta
Prod = 1
CosPwr = CosTheta
For ddf2 = 5 To Df2 Step 2
Prod *= (ddf2 - 3) / (ddf2 - 2)
CosPwr *= CosTheta * CosTheta
Sum += Prod * CosPwr
Next
A = 2 / Math.PI * (Theta + SinTheta * Sum)
Else
A = 2 * Theta / Math.PI
End If
' Calculate B
If Df2 > 1 And Df1 > 1 Then
Sum = 1
Prod = 1
SinPwr = 1
For ddf1 = 5 To Df1 Step 2
Prod *= (Df2 + ddf1 - 4) / (ddf1 - 2)
SinPwr *= SinTheta * SinTheta
Sum += Prod * SinPwr ' SinTheta ^ (ddf1 - 3)
Next
B = 2 / Math.Sqrt(Math.PI) * Gamma((Df2 - 1) / 2 + 1) / Gamma((Df2 - 2) / 2 + 1) * SinTheta * (CosTheta ^ Df2) * Sum
ElseIf Df1 = 1 Then
B = 0
Else
B = 0
End If
' calculate result
F = 1 - A + B
End If
Return 1 - F
End Function
Function Gamma(ByVal X As Double)
Const MAX = 26
Dim Sum, XPower, Prod As Double
Static C(MAX) As Double, DataFlag As Integer
Dim I As Integer
If DataFlag <> 123 Then
DataFlag = 123
C(1) = 1
C(2) = 0.577215664901533
C(3) = -0.655878071520254
C(4) = -0.0420026350340952
C(5) = 0.166538611382292
C(6) = -0.0421977345555443
C(7) = -0.009621971527877
C(8) = 0.007218943246663
C(9) = -0.0011651675918591
C(10) = -0.0002152416741149
C(11) = 0.0001280502823882
C(12) = -0.0000201348547807
C(13) = -0.00000125049348311
C(14) = 0.00000113302723202
C(15) = -0.00000020563384173
C(16) = 0.000000006116095
C(17) = 0.00000000500200753
C(18) = -0.0000000011812746
C(19) = 0.00000000010434271
C(20) = 0.00000000000778234
C(21) = -0.0000000000036968
C(22) = 0.00000000000051
C(23) = -0.0000000000000206
C(24) = -0.0000000000000054
C(25) = 0.0000000000000014
C(26) = 0.0000000000000001
End If
Prod = 1
Do While X > 1
X -= 1
Prod *= X
Loop
Sum = 0
XPower = X
For I = 1 To MAX
Sum += C(I) * XPower
XPower *= X
Next I
Return Prod / Sum
End Function
I know you want HP-41C code, but the above VB code gives you an idea what the HP-41C code must do!!! Not a small task!
Namir
Search
Member List
Calendar
Help
