Attribute VB_Name = "modMethod"
Option Explicit
'求t分布的分位数
'n：自由度，已知
'Q：上侧概率（ 'T：分位数，所求
Public Sub PT_DIST(n As Integer, Q As Double, T As Double)
Dim PIS As Double, DFR2 As Double, C As Double
Dim Q2 As Double, p As Double, YQ As Double, E As Double
Dim GA1 As Double, GA2 As Double, GA3 As Double
Dim T0 As Double, pp As Double, d As Double
Dim K As Integer
Const PI As Double = 3.14159265359
PIS = Sqr(PI): DFR2 = n / 2
If n = 1 Then
T = Tan(PI * (0.5 - Q)): Exit Sub
End If
If n = 2 Then
If Q > 0.5 Then C = -1 Else C = 1
Q2 = (1 - 2 * Q) ^ 2
T = Sqr(2 * Q2 / (1 - Q2)) * C
Exit Sub
End If
p = 1 - Q: PNorm Q, YQ '正态分布分位数
E = (1 - 1 / (4 * n)) ^ 2 - YQ * YQ / (2 * n)
If E > 0.5 Then
T0 = YQ / Sqr(E)
Else
lnGamma DFR2, GA1: lnGamma DFR2 + 0.5, GA2
GA3 = Exp((GA1 - GA2) / n)
T0 = Sqr(n) / (PIS * Q * n) ^ (1 / n) / GA3
End If
For K = 1 To 30
T_Dist n, T0, pp, d
If d = 0 Then
T = T0: Exit Sub
End If
T = T0 - (pp - p) / d
If Abs(T0 - T) < 0.000001 * Abs(T) Then _
Exit Sub Else T0 = T
Next K
End Sub
'计算GAMMA函数
'x：自变量
'z：GAMMA函数值
Public Sub GAMMA(x As Double, z As Double)
Dim H As Double, y As Double, y1 As Double
H = 1: y = x
LL1:
If y = 2 Then
z = H
Exit Sub
ElseIf y < 2 Then
H = H / y: y = y + 1: GoTo LL1
ElseIf y >= 3 Then
y = y - 1: H = H * y: GoTo LL1
End If
y = y - 2
y1 = y * (0.005159 + y * 0.001606)
y1 = y * (0.004451 + y1)
y1 = y * (0.07211 + y1)
y1 = y * (0.082112 + y1)
y1 = y * (0.41174 + y1)
y1 = y * (0.422787 + y1)
H = H * (0.999999 + y1)
z = H
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
'计算t分布的分布函数
'n：自由度，已知
'T：t值，已知
'pp：下侧概率，所求
'dd：概率密度，所求
Public Sub T_Dist(n As Integer, T As Double, pp As Double, dd As Double)
Dim Sign As Integer, TT As Double, x As Double
Dim p As Double, u As Double, GA1 As Double, GA2 As Double
Dim IBI As Integer, N2 As Integer, I As Integer
Const PI As Double = 3.14159265359
If T = 0 Then
Call GAMMA(n / 2, GA1): Call GAMMA(n / 2 + 0.5, GA2): pp = 0.5
dd = GA2 / (Sqr(n * PI) * GA1): Exit Sub
End If
If T < 0 Then Sign = -1 Else Sign = 1
TT = T * T: x = TT / (n + TT)
If (n \ 2) * 2 = n Then 'n为偶数
p = Sqr(x): u = p * (1 - x) / 2
IBI = 2
Else 'n为奇数
u = Sqr(x * (1 - x)) / PI
p = 1 - 2 * Atn(Sqr((1 - x) / x)) / PI
IBI = 1
End If
If IBI = n Then GoTo LL1 Else N2 = n - 2
For I = IBI To N2 Step 2
p = p + 2 * u / I
u = u * (1 + I) / I * (1 - x)
Next I
LL1:
dd = u / Abs(T)
pp = 0.5 + Sign * p / 2
End Sub
