[轉貼] t分佈

laigs218

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

laigs218

抱歉漏了一個
'計算t分佈函數的下側概率和密度函數
Sub Main()
    Dim I As Integer, TT As Double, p As Double, d As Double
    Dim sngP As Single, sngD As Single
    For I = 10 To 20 Step 5
        For TT = 1 To 6
            T_Dist I, TT, p, d
            sngP = p: sngD = d
            Debug.Print "n="; I; "   t="; TT; _
                "   下側概率="; sngP; "   密度函數="; sngD
        Next TT
    Next I
End Sub