趋势面分析m2.bas

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

Attribute VB_Name = "modMethod"
'趋势面分析
Option Explicit

'计算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

'计算趋势面多项式的系数，并求观测点的趋势值和残差值
'X：数组，观测数据的X坐标
'Y：数组，观测数据的Y坐标
'Z：数组，观测数据的Z坐标
'N0：趋势面的次数
'H：数组，保存趋势面多项式的系数
'C：保存正规方程系数
'T：数组，保存观测点的趋势值
'D：数组，保存观测点的残差值
Public Sub TREND(X() As Double, Y() As Double, Z() As Double, _
        N0 As Integer, H() As Double, C() As Double, _
        T() As Double, D() As Double)
    Dim NT As Integer, I As Integer, J As Integer, K As Integer
    Dim I0 As Integer, J0 As Integer, I9 As Integer, J9 As Integer
    Dim XX As Double, YY As Double, ZZ As Double
    Dim E As Double, F As Double, M9 As Double, w As Double
    Dim b1 As Double, n1 As Integer
'如果多项式的次数超过19则需要加大H、S的下标上界
    Dim S(1 To 200, 1 To 200) As Double
    Dim miX As Double, miY As Double, maX As Double, maY As Double
    n1 = UBound(X, 1)                       'N1为观测点个数
    NT = UBound(H, 1)                       '多项式的项数
    For I = 1 To n1
        XX = X(I): YY = Y(I): ZZ = Z(I)
        For J = 1 To NT
            Call Power(XX, YY, J, E, F)
            H(J) = E * F
        Next J
        For I0 = 1 To NT
            For J0 = 1 To NT
                S(I0, J0) = H(I0) * H(J0) + S(I0, J0)   '系数矩阵
            Next J0
            S(I0, NT + 1) = H(I0) * ZZ + S(I0, NT + 1)  '右端向量
        Next I0
    Next I
'将系数矩阵保存到C（）
    For I = 1 To NT
        For J = 1 To NT
            C(I, J) = S(I, J)
        Next J
    Next I
'***************************************************************
'解线性代数方程组
    For K = 1 To NT
        M9 = 0
        For I = K To NT
            For J = 1 To NT
                If Abs(M9) < Abs(S(I, J)) Then
                    M9 = S(I, J): I9 = I: J9 = J
                End If
            Next J
        Next I
        For J = 1 To NT + 1
            b1 = S(I9, J) / M9
            S(I9, J) = S(K, J)
            S(K, J) = b1
        Next J
        For I = 1 To NT
            If I <> K Then
                b1 = -S(I, J9)
                For J = 1 To NT + 1
                    S(I, J) = S(I, J) + b1 * S(K, J)
                Next J
            End If
        Next I
    Next K
    For J = 1 To NT
        For I = 1 To NT
'H为计算出的多项式系数
            If S(I, J) >= 0.5 Then H(J) = S(I, NT + 1)
        Next I
    Next J
'**********************************************************
'回代进行平滑处理
    For I = 1 To n1                         'N1为观测点个数
        w = 0
        XX = X(I): YY = Y(I)
        For I0 = 1 To NT                    '多项式的项数
            Call Power(XX, YY, I0, E, F)
            w = H(I0) * E * F + w
        Next I0
        T(I) = w                            '趋势值
        D(I) = Z(I) - w                     '残差值
    Next I
End Sub

'全主元高斯-约当消去法求逆矩阵
'A(1 To m, 1 To m)：开始存放欲求逆的矩阵，最终存求逆的结果矩阵，m是自变量个数
Public Sub Invert(a() As Double)
    Dim N As Integer, ep As Double
    Dim I As Integer, J As Integer, K As Integer
    Dim I0 As Integer, J0 As Integer
    Dim w As Double, Z As Double
    Dim b(1 To 100) As Double, C(1 To 100) As Double
    Dim p(1 To 100) As Double, Q(1 To 100) As Double
    N = UBound(a, 1)
    ep = 0.0000000001
    For K = 1 To N
        w = 0#
        For I = K To N
            For J = K To N
                If Abs(a(I, J)) > Abs(w) Then
                    w = a(I, J): I0 = I: J0 = J
                End If
            Next J
        Next I
        If Abs(w) < ep Then
            MsgBox "全主元素的绝对值小于0.0000000001，矩阵是奇异的！"
            Exit Sub
        End If
        If I0 <> K Then
            For J = 1 To N
                Z = a(I0, J): a(I0, J) = a(K, J): a(K, J) = Z
            Next J
        End If
        If J0 <> K Then
            For I = 1 To N
                Z = a(I, J0): a(I, J0) = a(I, K): a(I, K) = Z
            Next I
        End If
        p(K) = I0: Q(K) = J0
        For J = 1 To N
            If J = K Then
                b(J) = 1 / w: C(J) = 1
            Else
                b(J) = -a(K, J) / w: C(J) = a(J, K)
            End If
            a(K, J) = 0#: a(J, K) = 0#
        Next J
        For I = 1 To N
            For J = 1 To N
                a(I, J) = a(I, J) + C(I) * b(J)
            Next J
        Next I
    Next K
    For K = N To 1 Step -1
        I0 = p(K): J0 = Q(K)
        If I0 <> K Then
            For I = 1 To N
                Z = a(I, I0): a(I, I0) = a(I, K): a(I, K) = Z
            Next I
        End If
        If J0 <> K Then
            For J = 1 To N
                Z = a(J0, J): a(J0, J) = a(K, J): a(K, J) = Z
            Next J
        End If
    Next K
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
        YY = miY + DY * (I - 1)