Attribute VB_Name = "modMethod"
'多項式逐步回歸
Option Explicit
'xMy(1 To n, 1 To m)：觀測數據，已知，n是觀測次數，m是多項式最高冪數
'F1：指定的F臨界值，用于引入，已知
'F2：指定的F臨界值，用于剔出，已知
'    要求F1>=F2。如果F1=F2=0，則引入除線性相關外的全部變量
'F：F檢驗值，計算結果
'L：選出的重要冪次的個數，計算結果
'b(0 To m)：各個冪次的回歸系數，計算結果
'Ti(1 To m)：各冪次的t檢驗值，計算結果
Public Sub StrdM(xMy() As Double, F1 As Double, F2 As Double, F As Double, _
            L As Integer, b() As Single, Ti() As Single)
    Dim I As Integer, J As Integer, K As Integer
    Dim n As Integer, m As Integer, y As Integer
    Dim Imax As Integer, Imin As Integer
    Dim Ry12m As Double, Sy As Double, Syy As Double, V As Double
    Dim F12 As Double, K12 As Integer
    Dim Mx(1 To 101) As Double, Vx(1 To 101) As Double, Vyx(1 To 101) As Double
    Dim R(1 To 101, 1 To 101) As Double, Ri(1 To 101) As Double
    Dim d As Double, Sp As Integer, Q As Double, Vmax As Double, Vmin As Double
    Dim Bi As Double
    n = UBound(xMy, 1)                      'N是觀測數據點數
    y = UBound(xMy, 2)                      'y是最高冪次+因變量數
    m = y - 1                               'm是最高冪次
'求平均值，保存于Mx()
    For I = 1 To y
        d = 0
        For K = 1 To n
            d = xMy(K, I) + d
        Next K
        Mx(I) = d / n
    Next I
'計算離差矩陣，放在R的下三角部分
    For K = 1 To n
        For I = 1 To y
            d = xMy(K, I) - Mx(I): Vx(I) = d
            For J = 1 To I
                R(I, J) = d * Vx(J) + R(I, J)
            Next J
        Next I
    Next K
    For I = 1 To y
        Syy = R(I, I)
        If Syy = 0 Then
            MsgBox "某變量為常數，無法計算相關系數！"
            Exit Sub
        Else
        Vx(I) = Sqr(Syy)
        End If
    Next I
'計算相關矩陣，放在R的上三角部分
    For I = 2 To y
        d = Vx(I)
        For J = 1 To I - 1
            R(J, I) = R(I, J) / (d * Vx(J))
        Next J
    Next I
    d = Sqr(1 / (n - 1))
    For I = 1 To y
        Vx(I) = d * Vx(I)
        Vyx(I) = R(I, y)
    Next I
    For I = 1 To y
        R(I, I) = 1: Vyx(I) = Vx(y) / Vx(I)
        For J = I + 1 To y
            R(J, I) = R(I, J)
        Next J
    Next I
'法方程已建立，下面進入逐步計算
'計算各變量的貢獻V，從已入選的變量中找出最小的V，從未選量中找出最大的V
L2:
    L = 0: Sp = 0: Q = 1
LStep:
    Sp = Sp + 1: Vmax = 0: Vmin = 10
    For I = 1 To m
        Ti(I) = 0: d = R(I, I)
        If d > 0.00000001 Then
            V = (R(y, I) / d) * R(I, y)
            If V < 0 Then
                Ti(I) = d
                If -V < Vmin Then
                    Vmin = -V: Imin = I
                End If
            Else
                If V > Vmax Then
                    Vmax = V: Imax = I
                End If
            End If
        End If
    Next I
    If L <> 0 Then
        d = 0
        For I = 1 To m
            If Ti(I) = 0 Then
                b(I) = 0: Ri(I) = 0
            Else
                Bi = R(I, y): b(I) = Vyx(I) * Bi
                d = d + b(I) * Mx(I)
                Ri(I) = Bi / Sqr(Ti(I) * Q + Bi ^ 2)
                Ti(I) = Bi / Sqr(Ti(I) * Q / (n - L - 1))
            End If
        Next I
        b(0) = Mx(y) - d
    End If
    F12 = (n - L - 1) * Vmin / Q
    If F12 < F2 Then
        L = L - 1: K = Imin: K12 = -K
    Else
        F12 = (n - L - 2) * Vmax / (Q - Vmax)
        If F12 <= (F1 + 0.00000001) Then
            GoTo L3
        Else
            L = L + 1: K = Imax: K12 = K
        End If
    End If
'下面對R矩陣的第K列作消去計算
    d = 1 / R(K, K): R(K, K) = 1
    For J = 1 To y
        R(K, J) = R(K, J) * d
    Next J
    For I = 1 To y
        If I = K Then
        Else
            d = R(I, K): R(I, K) = 0
            For J = 1 To y
                R(I, J) = R(I, J) - d * R(K, J)
            Next J
        End If
    Next I
    Q = R(y, y): Ry12m = Sqr(1 - Q)
    F = (n - L - 1) * (1 - Q) / (L * Q)
    Sy = Sqr(Syy * Q / (n - L - 1))
    GoTo LStep
L3:
    If L = 0 Then
        MsgBox "在當前的引入F、剔出F下，不能選出重要的變量！"
        Exit Sub
    End If
    d = 0
    For I = 1 To m
        If Ti(I) = 0 Then
            b(I) = 0: Ri(I) = 0
        Else
            Bi = R(I, y): b(I) = Vyx(I) * Bi
            d = d + b(I) * Mx(I)
            Ri(I) = Bi / Sqr(Abs(Ti(I) * Q + Bi ^ 2))
            Ti(I) = Bi / Sqr(Abs(Ti(I) * Q / (n - L - 1)))
            Ti(I) = Abs(Ti(I))
        End If
    Next I
    b(0) = Mx(y) - d
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

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