Attribute VB_Name = "modMethod"
't分布
Option Explicit
'計算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

'計算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