Attribute VB_Name = "modMethod"
Option Explicit
'计算F分布的分位数
'n1：自由度，已知
'n2：自由度，已知
'Q：上侧概率，已知
'F：分位数，所求
Public Sub PF_DIST(n1 As Integer, n2 As Integer, _
Q As Double, F As Double)
Dim DF12 As Double, DF22 As Double, A As Double, B As Double
Dim A1 As Double, B1 As Double, p As Double, YQ As Double
Dim E As Double, FO As Double, pp As Double, d As Double
Dim GA1 As Double, GA2 As Double, GA3 As Double
Dim K As Integer
DF12 = n1 / 2: DF22 = n2 / 2
A = 2 / (9 * n1): A1 = 1 - A
B = 2 / (9 * n2): B1 = 1 - B
p = 1 - Q: PNorm Q, YQ
E = B1 * B1 - B * YQ * YQ
If E > 0.8 Then
FO = ((A1 * B1 + YQ * Sqr(A1 * A1 * B + A * E)) / E) ^ 3
Else
lnGamma DF12 + DF22, GA1
lnGamma DF12, GA2
lnGamma DF22, GA3
FO = (2 / n2) * (GA1 - GA2 - GA3 + 0.69315 + (DF22 - 1) * Log(n2) _
- DF22 * Log(n1) - Log(Q))
FO = Exp(FO)
End If
For K = 1 To 30
F_DIST n1, n2, FO, pp, d
If d = 0 Then
F = FO: Exit Sub
End If
F = FO - (pp - p) / d
If Abs(FO - F) < 0.000001 * Abs(F) Then Exit Sub Else FO = F
Next K
End Sub
'求Gamma函数的对数LogGamma(x)
'x：自变量
'G：Gamma函数的对数
Public Sub lnGamma(x As Double, G As Double)
Dim y As Double, z As Double, A As Double
Dim B As Double, B1 As Double, n As Integer
Dim I As Integer
If x < 8 Then
y = x + 8: n = -1
Else
y = x: n = 1
End If
z = 1 / (y * y)
A = (y - 0.5) * Log(y) - y + 0.9189385
B1 = (0.0007663452 * z - 0.0005940956) * z
B1 = (B1 + 0.0007936431) * z
B1 = (B1 - 0.002777778) * z
B = (B1 + 0.0833333) / y
G = A + B
If n >= 0 Then Exit Sub
y = y - 1: A = y
For I = 1 To 7
A = A * (y - I)
Next I
G = G - Log(A)
End Sub
'求正态分布的分位数
'Q：上侧概率
'x：分位数
Public Sub PNorm(Q, x)
Dim p As Double, y As Double, z As Double
Dim B0 As Double, B1 As Double, B2 As Double
Dim B3 As Double, B4 As Double, B5 As Double
Dim B6 As Double, B7 As Double, B8 As Double
Dim B9 As Double, B10 As Double, B As Double
B0 = 1.570796288: B1 = 0.03706987906
B2 = -0.0008364353589: B3 = -0.0002250947176
B4 = 0.000006841218299: B5 = 0.000005824238515
B6 = -0.00000104527497: B7 = 8.360937017E-08
B8 = -3.231081277E-09: B9 = 3.657763036E-11
B10 = 6.936233982E-13
If Q = 0.5 Then
x = 0: GoTo PN01
End If
If Q > 0.5 Then p = 1 - Q Else p = Q
y = -Log(4 * p * (1 - p))
B = y * (B9 + y * B10)
B = y * (B8 + B): B = y * (B7 + B)
B = y * (B6 + B): B = y * (B5 + B)
B = y * (B4 + B): B = y * (B3 + B)
B = y * (B2 + B): B = y * (B1 + B)
z = y * (B0 + B): x = Sqr(z)
If Q > 0.5 Then x = -x
PN01:
End Sub
'计算F分布的分布函数
'n1：自由度，已知
'n2：自由度，已知
'F：F值，已知
'p：下侧概率，所求
'd：概率密度，所求
Public Sub F_DIST(n1 As Integer, n2 As Integer, F As Double, _
p As Double, d As Double)
Dim x As Double, u As Double, Lu As Double
Dim IAI As Integer, IBI As Integer, nn1 As Integer, nn2 As Integer
Dim I As Integer
Const PI As Double = 3.14159265359
If F = 0 Then
p = 0: d = 0: Exit Sub
End If
x = n1 * F / (n2 + n1 * F)
If (n1 \ 2) * 2 = n1 Then
If (n2 \ 2) * 2 = n2 Then
u = x * (1 - x): p = x: IAI = 2: IBI = 2
Else
u = x * Sqr(1 - x) / 2: p = 1 - Sqr(1 - x): IAI = 2: IBI = 1
End If
Else
If (n2 \ 2) * 2 = n2 Then
p = Sqr(x): u = p * (1 - x) / 2: IAI = 1: IBI = 2
Else
u = Sqr(x * (1 - x)) / PI
p = 1 - 2 * Atn(Sqr((1 - x) / x)) / PI: IAI = 1: IBI = 1
End If
End If
nn1 = n1 - 2: nn2 = n2 - 2
If u = 0 Then
d = u / F
Exit Sub
Else
Lu = Log(u)
End If
If IAI = n1 Then GoTo LL1
For I = IAI To nn1 Step 2
p = p - 2 * u / I
Lu = Lu + Log((1 + IBI / I) * x)
u = Exp(Lu)
Next I
LL1:
If IBI = n2 Then
d = u / F: Exit Sub
End If
For I = IBI To nn2 Step 2
p = p + 2 * u / I
Lu = Lu + Log((1 + n1 / I) * (1 - x))
u = Exp(Lu)
Next I
d = u / F
End Sub
