Attribute VB_Name = "modFCheck"
'模糊識別
'F檢驗方法模塊
Option Explicit
'求F檢驗值
'X：試驗數(shù)據(jù)
'IJ：分類結(jié)果
'F：F檢驗值
Public Sub F_Check(X() As Double, IJ() As Integer, F As Double)
    Dim R As Integer, Nj As Integer
    Dim QA As Double, QB As Double, A As Double, B As Double
    Dim UA As Double, UB As Double, S As Double
    Dim SA2 As Double, SB2 As Double
    Dim N As Integer, M As Integer
    Dim Xka(1 To 100) As Double, Xkja(1 To 100, 1 To 100) As Double
    Dim Xkia(1 To 100) As Double
    Dim I As Integer, J As Integer, K As Integer
'N：樣本數(shù)。M：指標(biāo)數(shù)。
    N = UBound(X, 1): M = UBound(X, 2)
'求總體樣本的中心向量Xka
    For K = 1 To M
        Xka(K) = 0
        For I = 1 To N
            Xka(K) = Xka(K) + X(I, K)
        Next I
        Xka(K) = Xka(K) / N
    Next K
'求在某一個“入”下的分類數(shù)R
    R = 0                               '類計數(shù)器
    For I = 1 To 100
        If IJ(I, 1) <> 0 Then R = R + 1
    Next I
'求第J類Nj個元素第K個指標(biāo)的平均值Xkja
    For J = 1 To R
        Nj = 0
        For I = 1 To 100                '求第J類的元素個數(shù)Nj
            If IJ(J, I) <> 0 Then Nj = Nj + 1
        Next I
        For K = 1 To M
            For I = 1 To Nj
                Xkja(J, K) = Xkja(J, K) + X(IJ(J, I), K)
            Next I
            Xkja(J, K) = Xkja(J, K) / Nj
        Next K
    Next J
    QA = 0: QB = 0
    For J = 1 To R                      '對每種分類的循環(huán)
        Nj = 0                          '第j類的樣本數(shù)
        For I = 1 To 100
            If IJ(J, I) <> 0 Then Nj = Nj + 1
        Next I
'求類間的平方和
        A = 0
        For K = 1 To M
            A = A + (Xkja(J, K) - Xka(K)) ^ 2
        Next K
        A = Nj * Sqr(A)
        QA = QA + A                     '“QA/自由度”表征類間距離
'求類內(nèi)平方和
        B = 0
        For I = 1 To Nj
            S = 0
            For K = 1 To M
                S = S + (X(IJ(J, I), K) - Xkja(J, K)) ^ 2
            Next K
            B = B + Sqr(S)
        Next I
        QB = QB + B                     '“QB/自由度”表征類內(nèi)距離
    Next J
    UA = R - 1: UB = N - R              '類間自由度：類內(nèi)自由度
'SA2表征類間距離；SB2表征類內(nèi)距離
    If UA = 0 Or UB = 0 Then
        F = 9999
    Else
        SA2 = QA / UA: SB2 = QB / UB: F = SA2 / SB2
    End If
End Sub

'以下各公有過程為計算F值使用
'求Gamma函數(shù)的對數(shù)LogGamma(x)
'x：自變量
'G：Gamma函數(shù)的對數(shù)
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ài)分布的分位數(shù)
'Q：上側(cè)概率
'x：分位數(shù)
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分布的分布函數(shù)
'n1：自由度，已知
'n2：自由度，已知
'F：F值，已知
'p：下側(cè)概率，所求
'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分布的分位數(shù)
'n1：自由度，已知
'n2：自由度，已知
'Q：上側(cè)概率，已知
'F：分位數(shù)，所求
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