双因素交错m2.bas

<VB数理统计实用算法>书中的算法源程序

Attribute VB_Name = "modMethod"
Option Explicit
'双因素交错方差分析
'x：试验数据
'FA：行因素的F检验值
'FB：列因素的F检验值
'FAB：因素交错的F检验值
Public Sub DoubleAB(x() As Double, FA As Double, FB As Double, FAB As Double)
'p：垂直因素变化个数。q：水平因素变化个数。r：重复试验次数。n：试验的总次数
    Dim p As Integer, q As Integer, r As Integer, n As Integer
    Dim I As Integer, J As Integer, K As Integer
    Dim A As Double, B As Double, C As Double, D As Double
    Dim Qa As Double, Qb As Double, W As Double
    Dim Sa As Double, Sb As Double, Sab As Double, Se As Double
    p = UBound(x, 1): q = UBound(x, 2): r = UBound(x, 3)
    n = p * q * r                           'n：试验的总次数
    For I = 1 To p
        For J = 1 To q
            For K = 1 To r
                A = A + x(I, J, K)
                B = B + x(I, J, K)
                C = C + x(I, J, K)
                W = W + x(I, J, K) ^ 2
            Next K
            D = D + C * C
            C = 0
        Next J
        Qa = Qa + B * B
        B = 0
    Next I
    B = A ^ 2 / n
    Qa = Qa / (q * r)
    For J = 1 To q
        For I = 1 To p
            For K = 1 To r
                C = C + x(I, J, K)
            Next K
        Next I
        Qb = Qb + C * C
        C = 0
    Next J
    Qb = Qb / (p * r)
    D = D / r
    Sa = (Qa - B) / (p - 1)     '行因素引起的平均离差平方和
    Sb = (Qb - B) / (q - 1)     '列因素引起的平均离差平方和
'因素交互作用引起的平均离差平方和
    Sab = (D - Qa - Qb + B) / (p - 1) / (q - 1)
'误差引起的平均离差平方和
    Se = (W - D) / (p * q * (r - 1))
    FA = Sa / Se                '行因素的检验值
    FB = Sb / Se                '列因素的检验值
    FAB = Sab / Se              '交错因素的检验值
Debug.Print "行平方和", (Qa - B), "自由度", p - 1, "均方", Sa
Debug.Print "列平方和", (Qb - B), "自由度", q - 1, "均方", Sb
Debug.Print "交错平方和", (D - Qa - Qb + B), "自由度", (p - 1) * (q - 1), "均方", Sab
Debug.Print "误差平方和", (W - D), "自由度", (p * q * (r - 1)), "均方", Se
Debug.Print "行检验值", FA, "列检验值", FB, "交错检验值", FAB
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

'求正态分布的分位数
'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