Attribute VB_Name = "modMethod"
Option Explicit
'相關分析
'x(1 To n)：變量，n為觀測次數，已知
'y(1 To n)：變量，n為觀測次數，已知
'R：相關系數，計算結果
't：t檢驗值，計算結果
Public Sub Relation(x() As Double, y() As Double, R As Double, t As Double)
    Dim Xa As Double, Ya As Double, Sxx As Double, Sxy As Double, Syy As Double
    Dim n As Integer, I As Double
    n = UBound(x, 1)
    For I = 1 To n
        Xa = Xa + x(I): Ya = Ya + y(I)
    Next I
    Xa = Xa / n: Ya = Ya / n                        '平均值
    For I = 1 To n
        Sxx = Sxx + (x(I) - Xa) ^ 2
        Sxy = Sxy + (x(I) - Xa) * (y(I) - Ya)
        Syy = Syy + (y(I) - Ya) ^ 2
    Next I
    R = Sxy / Sqr(Sxx * Syy)
'求t值
    t = (R / Sqr(1 - R ^ 2)) * Sqr(n - 2)
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

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

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

'求t分布的分位數
'n：自由度，已知
'Q：上側概率（<=0.5），已知
'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