二元多项式逐步回归m2.bas

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

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 Double, Ti() As Double)
    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是多项式的项数+因变量数（=1）
    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

'计算X和Y的乘幂
Public Sub Power(X, Y, I0, E, F)
    Dim a As Integer, C As Integer
    a = 0
SS1:
    If ((a + 1) * (a + 2)) / 2 >= I0 Then GoTo SS2:
    a = a + 1: GoTo SS1
SS2:
    C = (a + 1) * (a + 2) / 2 - I0
    If X <> 0 Or C <> 0 Then E = X ^ C Else E = 1
    If Y <> 0 Or a - C <> 0 Then F = Y ^ (a - C) Else F = 1
End Sub

'求预测值
'X0：预测点的X坐标值
'Y0：预测点的Y坐标值
'NT：总项数
'Z：预测值，结果
Public Sub PreValue(X0, Y0, NT, w)
    Dim I As Integer, E As Double, F As Double
    w = 0
    For I = 1 To NT                   '多项式的项数
        Call Power(X0, Y0, I, E, F)
        w = H(I) * E * F + w
    Next I
End Sub

'近N点按距离加权平均法作曲面插值
'X：数据点X坐标数组
'Y：数据点Y坐标数组
'Z：数据点Z坐标(函数值)数组
'A：插值点X坐标
'B：插值点Y坐标
'F：插值点函数值
'N：取离插值点最近的数据点个数
Public Sub NDS(X() As Double, Y() As Double, Z() As Double, _
            a As Double, b As Double, F As Double, N As Integer)
    Dim I As Integer, J As Integer, IC As Integer, M As Integer
    Dim S1 As Double, S2 As Double, D As Double
    Dim DIS(1000) As Double
    On Error GoTo errL
    M = UBound(X, 1)                        '数据点数
'求各数据点与插值点的距离
    For I = 1 To M
        DIS(I) = (X(I) - a) ^ 2 + (Y(I) - b) ^ 2
    Next I
    S1 = 0#: S2 = 0#
    For I = 1 To N
        IC = 1
        For J = 1 To M
            If DIS(J) < DIS(IC) Then IC = J
        Next J
        If DIS(IC) < 0.0001 Then
'当距离很近时数据点函数值即为插值点函数值
            F = Z(IC)
            Exit Sub
        End If
        D = Sqr(DIS(IC))
        S1 = S1 + Z(IC) / D
        S2 = S2 + 1 / D
        DIS(IC) = 10000000
    Next I
    F = S1 / S2
    Exit Sub
errL:
    MsgBox "不同的数据点有相同的X、Y坐标，造成除数为0"
End Sub

'按所求得的趋势面多项式计算网格点的趋势值
'X：数组，观测数据的X坐标
'Y：数组，观测数据的Y坐标
'H：数组，趋势面多项式的系数
'B：数组，保存网格点的趋势值
'G：数组，保存网格点的残差值
Public Sub GRID(X() As Double, Y() As Double, H() As Double, _
            b() As Double, G() As Double)
    Dim n1 As Integer, NT As Integer, M As Integer, N As Integer
    Dim I As Integer, J As Integer, I0 As Integer
    Dim XX As Double, YY As Double
    Dim miX As Double, maX As Double, miY As Double, maY As Double
    Dim DX As Double, DY As Double, E As Double, F As Double
    n1 = UBound(X, 1)                       'N1为观测点个数
    NT = UBound(H, 1)                       '多项式的项数
    M = UBound(b, 1)                        '网格的行数
    N = UBound(b, 2)                        '网格的列数
    miX = X(1): miY = Y(1): maX = X(1): maY = Y(1)
    For I = 1 To n1
        If X(I) < miX Then miX = X(I)       '求观测值X坐标最小值
        If Y(I) < miY Then miY = Y(I)       '求观测值Y坐标最小值
        If X(I) > maX Then maX = X(I)       '求观测值X坐标最大值
        If Y(I) > maY Then maY = Y(I)       '求观测值Y坐标最大值
    Next I
    DX = (maX - miX) / (N - 1)              '网格在X方向上的增量
    DY = (maY - miY) / (M - 1)              '网格在Y方向上的增量
    For I = 1 To M
        YY = miY + DY * (I - 1)
        For J = 1 To N
            XX = miX + DX * (J - 1)
            For I0 = 1 To NT
                Call Power(XX, YY, I0, E, F)
'b为网格值
                b(I, J) = H(I0) * E * F + b(I, J)
            Next I0
        Next J
    Next I
    For I = 1 To M