'*************************************' 全局变量，便于主函数调用。

' VB 6.0 的函数返回的参数偏少，' 使用全局变量在一定程度可以解决这个问题。

'****************************************Public A() As Single ' 协方差/相关系数矩阵APublic V() As Single '特征向量为列组成的矩阵，即空间函数V （EOF）Public T() As Single '时间系数矩阵T（PC）Public B() As Single '特征值λ（E）,按从大到小排列Public GM() As Single '解释的方差(%)（特征向量对X场的累积贡献率）P Public GA() As SinglePublic GB() As Single '个体i特征向量对X场的贡献率ρPublic XF() As Single '模拟结果'********************************************************' 函数名：CovarMat' 函数用途: 计算协方差（相关系数）矩阵' 参数说明：矩阵下标为1：N,从1开始；' X,存放原始观测值，二维实型数组，X(P,P)。

' 返回：计算协方差（相关系数）矩阵。

'*******************************************************Function CovarMat(X() As Single) As Single()Dim XX() As SingleDim P As Integer, N As IntegerDim px As SingleP = UBound(X, 1)N = UBound(X, 2)px = IIf(N > 0, 1 / N, 1)ReDim Preserve XX(1 To P, 1 To P)Dim iAs Integer, j As Integer, k As Integer' 求X乘以X的转置,即A=XXˊFor i = 1 To PFor j = 1 To PXX(i, j) = 0For k = 1 To NXX(i, j) = XX(i, j) + X(i, k) * X(j, k)Next kXX(i, j) = XX(i, j) * pxNext jNext iCovarMat = XXEnd Function'********************************************************' 函数名：EOF' 函数用途: EOF,经验正交分解法(EOF)' 参数说明：矩阵下标为1：N,从1开始；' X,存放原始观测值，二维实型数组，X(P,N)。

' P,整型变量，空间格点数。

' N,整型变量，序列的时间长度。

' XF,返回时存放恢复值，二维实型数组，XF(P,N)。

'*******************************************************Sub EOF(X() As Single, ByRef S() As String)Dim V1() As SingleDim VF() As Single, TF() As SingleDim P As Integer, N As IntegerP = UBound(X, 1)N = UBound(X, 2)ReDimXF(1 To P, 1 To N)ReDimT(1 To P, 1 To N)ReDimA(1 To P, 1 To P)ReDimV(1 To P, 1 To P)ReDimV1(1 To P, 1 To P)ReDimB(1 To P)ReDimGM(1 To P)ReDimGA(1 To P)ReDimGB(1 To P)ReDimVF(1 To P, 1 To P)ReDimTF(1 To P, 1 To N)' 求X的协方差（相关系数）矩阵A = CovarMat(X)Dim iAs Integer, j As Integer, k As Integer, L As Integer' 用Jacobi法求A的特征值和特征向量' 返回时B存放矩阵的全部特征值，V存放特征向量为列组成的矩阵Call Jacobi(A, P, 0.000001, V, B, L)For i = 1 To PGA(i) = 0For j = 1 ToiGA(i) = GA(i) + B(j)Next jNext iFor i = 1 To PGM(i) = GA(i) / GA(P)GB(i) = B(i) / GA(P)Next iFor i = 1 To PFor j = 1 To PV1(i, j) = V(j, i)Next jNext iT = MATMUL(V1, X) '时间系数'输出结果字符串Dim lsAs Integerls = UBound(S)ReDim Preserve S(ls + 2)S(ls) = MA TtoString(B, 1, , "特征值λ（E）")S(ls + 1) = MA TtoString(GB, 1, 100, "个体i特征向量对X场的贡献率ρ")S(ls + 2) = MA TtoString(GM, 1, 100, "解释的方差(%)（特征向量对X场的累积贡献率）P")For i = 1 To PFor j = 1 TolwVF(i, j) = V(i, j)Next jNext iFor i = 1 TolwFor j = 1 To NTF(i, j) = T(i, j)Next jNext iXF = MATMUL(VF, TF) '模拟值End Sub'*******************'函数名：MATtoString'函数作用：矩阵转成字符串输出'参数说明：' MAT，用以存储矩阵数值，最多为2维数组' retS，返回时存储字符串，以换行符（vbcrlf）结尾' nDim，矩阵维数，最大为2' strMatNote,矩阵备注信息，默认为空字符串'********************Function MATtoString(MAT() As Single, _nDim As Integer, _Optional ZoomCoefAs Single = 1, _Optional MatNoteStringAs String = "") As StringDim retString As StringDim N As Integer, i As Integer'如果数组说明为非空串，则添加换行符retString = IIf(Len(MatNoteString) > 0, MatNoteString&vbCrLf, MatNoteString)Select Case nDimCase 1N = UBound(MAT)For i = 1 To NretString = retString& Format(MAT(i) * ZoomCoef, "#0.##") &vbTabNext iretString = retString&vbCrLfCase 2Dim m As Integer, j As IntegerN = UBound(MA T, 1)m = UBound(MA T, 2)For i = 1 To NFor j = 1 To mretString = retString& Format(MAT(i, j) * ZoomCoef, "#0.##") &vbTabNextretString = retString&vbCrLfNext iEnd SelectMATtoString = retStringEnd Function'*************************'函数名：MATMUL'函数作用：矩阵相乘'参数说明：' Mat1，用以存储矩阵1数值' Mat2，用以存储矩阵1数据' 返回：矩阵相乘结果'**************************Function MATMUL(MAT1() As Single, MA T2() As Single) As Single()Dim MatXX() As SingleDim N As Integer, m As Integer, L As Integer, l2 As IntegerDim iAs Integer, j As Integer, k As IntegerN = UBound(MAT1, 1)m = UBound(MAT2, 2)L = UBound(MAT1, 2)l2 = UBound(MAT2, 1)If L <> l2 Then EndReDimMatXX(1 To N, 1 To m)For i = 1 To NFor j = 1 To mMatXX(i, j) = 0For k = 1 To LMatXX(i, j) = MatXX(i, j) + MAT1(i, k) * MAT2(k, j)Next kNext jNext iMATMUL = MatXXEnd Function'********************************************************' 函数名：Jacobi' 函数用途: 用Jacobi法求A的特征值和特征向量' 参数说明：矩阵下标为1：N,从1开始；' A:调用时存放实对称矩阵,A(N,N)' N:矩阵长度' Bx:返回时存放矩阵的全部特征值,B(N)' Vx:特征向量为列组成的矩阵,V(N,N),其中第i列为与第i个特征值相对应的特征向量' EPS:存放精度要求'*********************************************************Private Sub Jacobi(A() As Single, _N As Integer, _EPS As Single, _ByRefVx() As Single, _ByRefBx() As Single, _ByRef L As Integer)' A:调用时存放实对称矩阵' Bx:返回时存放矩阵的全部特征值' Vx:存放特征向量,其中第i列为与第i个特征值相对应的特征向量' EPS:存放精度要求ReDimVx(1 To N, 1 To N)ReDimBx(1 To N)Dim FM As Single, CN As Single, SN As SingleDim Omega As Single, X As Single, Y AsSingleDim P As Integer, Q As IntegerDim iAs Integer, j As Integer, k As IntegerL = 1'初始化特征向量矩阵For i = 1 To NFor j = 1 To NVx(i, j) = 0Next jVx(i, i) = 1Next i'如果计算次数大于给定次数，则跳出循环Do While (L < 100)FM = 0'查找矩阵中非对角元素的最大值，并记录其位置（P,Q)For i = 1 To NFor j = 1 To NIf (i<> j And Abs(A(i, j)) > FM) ThenFM = Abs(A(i, j))P = iQ = jEnd IfNext jNext i'如果最大值小于给定最小值，则跳出循环If (FM < EPS) ThenL = -1Exit DoEnd IfL = L + 1X = -A(P, Q)Y = (A(Q, Q) - A(P, P)) / 2Omega = X / VBA.Sqr(X * X + Y * Y)If (Y < 0) Then Omega = -OmegaSN = 1 + VBA.Sqr(1 - Omega * Omega)SN = Omega / VBA.Sqr(2 * SN)CN = VBA.Sqr(1 - SN * SN)FM = A(P, P)A(P, P) = FM * CN * CN + A(Q, Q) * SN * SN + A(P, Q) * OmegaA(Q, Q) = FM * SN * SN + A(Q, Q) * CN * CN - A(P, Q) * OmegaA(P, Q) = 0 '处理完当前最大非对角元素值，赋值0，下一循环不再计算A(Q, P) = 0 '处理完当前最大非对角元素值，赋值0，下一循环不再计算For j = 1 To NIf (j <> P And j <> Q) ThenFM = A(P, j)A(P, j) = FM * CN + A(Q, j) * SNA(Q, j) = -FM * SN + A(Q, j) * CNEnd IfNext jFor i = 1 To NIf (i<> P Andi<> Q) ThenFM = A(i, P)A(i, P) = FM * CN + A(i, Q) * SNA(i, Q) = -FM * SN + A(i, Q) * CNEnd IfNext iFor i = 1 To NFM = Vx(i, P)Vx(i, P) = FM * CN + Vx(i, Q) * SNVx(i, Q) = -FM * SN + Vx(i, Q) * CNNext iFor i = 1 To NBx(i) = A(i, i)Next iLoop' 将特征值按大小排列m = NDo While (m > 0)j = m - 1m = 0For i = 1 To jIf (Bx(i) <Bx(i + 1)) ThenB1 = Bx(i)Bx(i) = Bx(i + 1)Bx(i + 1) = B1m = iFor k = 1 To NV1 = Vx(k, i) Vx(k, i) = Vx(k, i + 1)Vx(k, i + 1) = V1Next kEnd IfNext iLoopEnd Sub。