Attribute VB_Name = "modMethod"
'方差的假設(shè)檢驗(yàn)
Option Explicit
'計(jì)算卡方分布的分位數(shù)
'n：自由度
'Q：上側(cè)概率
'xx：分位數(shù)
Public Sub PCX2(n As Integer, Q As Double, xx As Double)
    Dim I As Integer, x As Double, p As Double, W As Double
    Dim x0 As Double, pp As Double, d As Double
    If n = 1 Then
        PNorm Q / 2, x: xx = x * x
        Exit Sub
    End If
    If n = 2 Then
        xx = -2 * Log(Q)
        Exit Sub
    End If
    p = 1 - Q: PNorm Q, x: W = 2 / (9 * n)
    x0 = n * (1 - W + x * Sqr(W)) ^ 3
    For I = 1 To 30
        CX2 n, x0, pp, d
        If d = 0 Then
            xx = x0
            Exit Sub
        End If
        xx = x0 - (pp - p) / d
        If Abs(xx - x0) < 10 - 6 * Abs(xx) Then Exit Sub Else x0 = xx
    Next I
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

'計(jì)算卡方分布函數(shù)和概率密度
'n：自由度
'x2：卡方值
'F：下側(cè)概率
'd：概率密度
Public Sub CX2(n As Integer, X2 As Double, F As Double, d As Double)
    Dim PIS As Double, x As Double, CHS As Double, u As Double
    Dim IAI As Integer, pp As Double, n2 As Integer, I As Integer
    Const PI As Double = 3.14159265359
    If X2 = 0 Then
        F = 0: d = 0: Exit Sub
    End If
    PIS = Sqr(PI)
    x = X2 / 2
    CHS = Sqr(X2)
    If (n \ 2) * 2 = n Then                 'n為偶數(shù)
        u = x * Exp(-x)
        F = 1 - Exp(-x)
        IAI = 2
    Else                                    'n為奇數(shù)
        u = Sqr(x) * Exp(-x) / PIS
        Norm CHS, pp                        '調(diào)用正態(tài)分布函數(shù)計(jì)算過程
        F = 2 * (pp - 0.5)
        IAI = 1
    End If
    If IAI = n Then GoTo LL1 Else n2 = n - 2
    For I = IAI To n2 Step 2
        F = F - 2 * u / I
        u = X2 * u / I
    Next I
LL1:
    d = u / X2
End Sub

'計(jì)算正態(tài)分布函數(shù)
'x：正態(tài)偏離點(diǎn)
'F：下側(cè)概率
Public Sub Norm(x, F)
    Dim y As Double, ER As Double, Q As Double
    Dim A
    Const A1 As Double = 0.0705230784
    Const a2 As Double = 0.0422820123
    Const a3 As Double = 0.0092705272
    Const a4 As Double = 0.0001520143
    Const a5 As Double = 0.0002765672
    Const a6 As Double = 0.0000430638
    y = 0.707106781187 * Abs(x)
    A = a4 + y * (a5 + y * a6)
    A = a3 + y * A
    A = a2 + y * A
    A = A1 + y * A
    ER = 1 - (1 + y * A) ^ (-16)
    Q = 0.5 * ER
    If x < 0 Then F = 0.5 - Q Else F = 0.5 + Q
End Sub

'計(jì)算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

'求Gamma函數(shù)的對(duì)數(shù)LogGamma(x)
'x：自變量
'G：Gamma函數(shù)的對(duì)數(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

'計(jì)算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