correlation_dimension.cls

修改过去混沌算法中的参数选择问题,使算法变得更准确.

VERSION 1.0 CLASS
BEGIN
  MultiUse = -1  'True
  Persistable = 0  'NotPersistable
  DataBindingBehavior = 0  'vbNone
  DataSourceBehavior  = 0  'vbNone
  MTSTransactionMode  = 0  'NotAnMTSObject
END
Attribute VB_Name = "correlation_dimension"
Attribute VB_GlobalNameSpace = False
Attribute VB_Creatable = True
Attribute VB_PredeclaredId = False
Attribute VB_Exposed = False

Option Explicit
Const tao As Long = 8
Const min_d As Long = 5
Const max_d As Long = 40
Const KK As Integer = 100
Const JJ As Integer = 80
Dim A As Double
Dim B As Double

Public Sub computing_correlation_dimension(x_series() As Double)
    Dim i, j, k As Long
    Dim N As Long
    N = UBound(x_series())
    Dim d As Long
    
    Dim xx() As Double
    Dim rr() As Double
    Dim h() As Double
    Dim x_temp_i() As Double
    Dim x_temp_j() As Double
    
    Dim h_temp(KK) As Double
    Dim r_temp(KK) As Double
    Dim h_log(JJ) As Double
    Dim r_log(JJ) As Double
    
    Dim SIG(KK) As Double
    
    Dim Dimen(max_d) As Double
    Dim max_r As Double
    Dim min_r As Double
    
    
    For d = min_d To max_d
        ReDim xx(N - (d - 1) * tao, d) As Double
        ReDim rr(N - (d - 1) * tao, N - (d - 1) * tao) As Double
        ReDim h(N - (d - 1) * tao, N - (d - 1) * tao) As Double
        ReDim x_temp_i(d) As Double
        ReDim x_temp_j(d) As Double
        
        '重构d维相空间，共有N - (d - 1) * tao个点
        For i = 1 To N - (d - 1) * tao
            For k = 1 To d
                 xx(i, k) = x_series(i + (k - 1) * tao)
            Next k
        Next i
        '计算Euclid范数矩阵
        min_r = 10000
        max_r = 0
        For i = 1 To N - (d - 1) * tao
            For j = 1 To N - (d - 1) * tao
                For k = 1 To d
                   x_temp_i(k) = xx(i, k)
                   x_temp_j(k) = xx(j, k)
                Next k
                rr(i, j) = fnorm(x_temp_i(), x_temp_j())
                If rr(i, j) > max_r Then
                   max_r = rr(i, j)
                End If
                If rr(i, j) < min_r Then
                   min_r = rr(i, j)
                End If
            Next j
        Next i
        
        '计算heaviside阶跃函数
        For k = 1 To KK
            h_temp(k) = 0
            For i = 1 To N - (d - 1) * tao
                For j = 1 To N - (d - 1) * tao
                    h(i, j) = heaviside(min_r + k * (max_r - min_r) / KK, rr(i, j))
                    h_temp(k) = h_temp(k) + h(i, j)
                Next j
            Next i
            '计算关联积分函数
            h_temp(k) = (h_temp(k) - (N - (d - 1) * tao)) / ((N - (d - 1) * tao) * ((N - (d - 1) * tao) - 1))
            r_temp(k) = k * (max_r - min_r) / KK
            SIG(k) = 0.5
        Next k
        For j = 1 To JJ
            h_log(j) = h_temp(j + 10)
            r_log(j) = r_temp(j + 10)
        Next j
        '直线拟合
        Call FIT(r_log(), h_log(), JJ, SIG(), 0)
        '计算混沌关联维数

        Dimen(d) = B
        Debug.Print B
        'If Abs(Dimen(d) - Dimen(d - 1)) < 0.001 Then
        '    GoTo flai_end
        'End If
    Next d
flai_end:
End Sub

Private Function fnorm(ri() As Double, rj() As Double) As Double '计算Euclid范数
    Dim M As Long
    M = UBound(ri())
    Dim dist As Double
    dist = 0
    Dim k As Long
    For k = 1 To M
        dist = dist + (ri(k) - rj(k)) ^ 2
    Next k
    fnorm = Sqr(dist)
End Function

Private Function heaviside(number1 As Double, number2 As Double) As Long 'heaviside阶跃函数
    If number1 >= number2 Then
       heaviside = 1
    Else
       heaviside = 0
    End If
End Function

Sub FIT(X() As Double, Y() As Double, NDATA As Integer, SIG() As Double, MWT As Integer) '直线拟合函数
    
    Dim SIGA As Double
    Dim SIGB As Double
    Dim CHI2 As Double
    Dim Q As Double
    Dim SX As Double
    Dim SY As Double
    Dim ST2 As Double
    Dim SS As Double
    Dim SXOSS As Double
    Dim SIGDAT As Double
    Dim i As Integer
    Dim WT As Double
    Dim T As Double
    
    
    SX = 0#
    SY = 0#
    ST2 = 0#
    B = 0#
    If MWT <> 0 Then
        SS = 0#
        For i = 1 To NDATA
            WT = 1# / (SIG(i) ^ 2)
            SS = SS + WT
            SX = SX + X(i) * WT
            SY = SY + Y(i) * WT
        Next i
    Else
        For i = 1 To NDATA
            SX = SX + X(i)
            SY = SY + Y(i)
        Next i
        SS = NDATA
    End If
    SXOSS = SX / SS
    If MWT <> 0 Then
        For i = 1 To NDATA
            T = (X(i) - SXOSS) / SIG(i)
            ST2 = ST2 + T * T
            B = B + T * Y(i) / SIG(i)
        Next i
    Else
        For i = 1 To NDATA
            T = X(i) - SXOSS
            ST2 = ST2 + T * T
            B = B + T * Y(i)
    Next i
    End If
    B = B / ST2
    A = (SY - SX * B) / SS
    SIGA = Sqr((1# + SX * SX / (SS * ST2)) / SS)
    SIGB = Sqr(1# / ST2)
    CHI2 = 0#
    If MWT = 0 Then
        For i = 1 To NDATA
            CHI2 = CHI2 + (Y(i) - A - B * X(i)) ^ 2
        Next i
        Q = 1#
        SIGDAT = Sqr(CHI2 / (NDATA - 2))
        SIGA = SIGA * SIGDAT
        SIGB = SIGB * SIGDAT
    Else
        For i = 1 To NDATA
            CHI2 = CHI2 + ((Y(i) - A - B * X(i)) / SIG(i)) ^ 2
        Next i
        Q = 0
    End If
End Sub