正交函数
-
2023年3月17日发(作者：阴阳互根)'*************************************
'全局变量，便于主函数调用。
'VB6.0的函数返回的参数偏少，
'使用全局变量在一定程度可以解决这个问题。
'****************************************
PublicA()AsSingle'协方差/相关系数矩阵A
PublicV()AsSingle'特征向量为列组成的矩阵，即空间函数V（EOF）
PublicT()AsSingle'时间系数矩阵T（PC）
PublicB()AsSingle'特征值λ（E）,按从大到小排列
PublicGM()AsSingle'解释的方差(%)（特征向量对X场的累积贡献率）P
PublicGA()AsSingle
PublicGB()AsSingle'个体i特征向量对X场的贡献率ρ
PublicXF()AsSingle'模拟结果
'********************************************************
'函数名：CovarMat
'函数用途:计算协方差（相关系数）矩阵
'参数说明：矩阵下标为1：N,从1开始；
'X,存放原始观测值，二维实型数组，X(P,P)。
'返回：计算协方差（相关系数）矩阵。
'*******************************************************
FunctionCovarMat(X()AsSingle)AsSingle()
DimXX()AsSingle
DimPAsInteger,NAsInteger
DimpxAsSingle
P=UBound(X,1)
N=UBound(X,2)
px=IIf(N>0,1/N,1)
ReDimPreserveXX(1ToP,1ToP)
DimiAsInteger,jAsInteger,kAsInteger
'求X乘以X的转置,即A=XXˊ
Fori=1ToP
Forj=1ToP
XX(i,j)=0
Fork=1ToN
XX(i,j)=XX(i,j)+X(i,k)*X(j,k)
Nextk
XX(i,j)=XX(i,j)*px
Nextj
Nexti
CovarMat=XX
EndFunction
'********************************************************
'函数名：EOF
'函数用途:EOF,经验正交分解法(EOF)
'参数说明：矩阵下标为1：N,从1开始；
'X,存放原始观测值，二维实型数组，X(P,N)。
'P,整型变量，空间格点数。
'N,整型变量，序列的时间长度。
'XF,返回时存放恢复值，二维实型数组，XF(P,N)。
'*******************************************************
SubEOF(X()AsSingle,ByRefS()AsString)
DimV1()AsSingle
DimVF()AsSingle,TF()AsSingle
DimPAsInteger,NAsInteger
P=UBound(X,1)
N=UBound(X,2)
ReDimXF(1ToP,1ToN)
ReDimT(1ToP,1ToN)
ReDimA(1ToP,1ToP)
ReDimV(1ToP,1ToP)
ReDimV1(1ToP,1ToP)
ReDimB(1ToP)
ReDimGM(1ToP)
ReDimGA(1ToP)
ReDimGB(1ToP)
ReDimVF(1ToP,1ToP)
ReDimTF(1ToP,1ToN)
'求X的协方差（相关系数）矩阵
A=CovarMat(X)
DimiAsInteger,jAsInteger,kAsInteger,LAsInteger
'用Jacobi法求A的特征值和特征向量
'返回时B存放矩阵的全部特征值，V存放特征向量为列组成的矩阵
CallJacobi(A,P,0.000001,V,B,L)
Fori=1ToP
GA(i)=0
Forj=1Toi
GA(i)=GA(i)+B(j)
Nextj
Nexti
Fori=1ToP
GM(i)=GA(i)/GA(P)
GB(i)=B(i)/GA(P)
Nexti
Fori=1ToP
Forj=1ToP
V1(i,j)=V(j,i)
Nextj
Nexti
T=MATMUL(V1,X)'时间系数
'输出结果字符串
DimlsAsInteger
ls=UBound(S)
ReDimPreserveS(ls+2)
S(ls)=MATtoString(B,1,,"特征值λ（E）")
S(ls+1)=MATtoString(GB,1,100,"个体i特征向量对X场的贡献率ρ")
S(ls+2)=MATtoString(GM,1,100,"解释的方差(%)（特征向量对X场的累积贡献率）P")
Fori=1ToP
Forj=1Tolw
VF(i,j)=V(i,j)
Nextj
Nexti
Fori=1Tolw
Forj=1ToN
TF(i,j)=T(i,j)
Nextj
Nexti
XF=MATMUL(VF,TF)'模拟值
EndSub
'*******************
'函数名：MATtoString
'函数作用：矩阵转成字符串输出
'参数说明：
'MAT，用以存储矩阵数值，最多为2维数组
'retS，返回时存储字符串，以换行符（vbcrlf）结尾
'nDim，矩阵维数，最大为2
'strMatNote,矩阵备注信息，默认为空字符串
'********************
FunctionMATtoString(MAT()AsSingle,_
nDimAsInteger,_
OptionalZoomCoefAsSingle=1,_
OptionalMatNoteStringAsString="")AsString
DimretStringAsString
DimNAsInteger,iAsInteger
'如果数组说明为非空串，则添加换行符
retString=IIf(Len(MatNoteString)>0,MatNoteString&vbCrLf,MatNoteString)
SelectCasenDim
Case1
N=UBound(MAT)
Fori=1ToN
retString=retString&Format(MAT(i)*ZoomCoef,"#0.##")&vbTab
Nexti
retString=retString&vbCrLf
Case2
DimmAsInteger,jAsInteger
N=UBound(MAT,1)
m=UBound(MAT,2)
Fori=1ToN
Forj=1Tom
retString=retString&Format(MAT(i,j)*ZoomCoef,"#0.##")&vbTab
Next
retString=retString&vbCrLf
Nexti
EndSelect
MATtoString=retString
EndFunction
'*************************
'函数名：MATMUL
'函数作用：矩阵相乘
'参数说明：
'Mat1，用以存储矩阵1数值
'Mat2，用以存储矩阵1数据
'返回：矩阵相乘结果
'**************************
FunctionMATMUL(MAT1()AsSingle,MAT2()AsSingle)AsSingle()
DimMatXX()AsSingle
DimNAsInteger,mAsInteger,LAsInteger,l2AsInteger
DimiAsInteger,jAsInteger,kAsInteger
N=UBound(MAT1,1)
m=UBound(MAT2,2)
L=UBound(MAT1,2)
l2=UBound(MAT2,1)
IfLl2ThenEnd
ReDimMatXX(1ToN,1Tom)
Fori=1ToN
Forj=1Tom
MatXX(i,j)=0
Fork=1ToL
MatXX(i,j)=MatXX(i,j)+MAT1(i,k)*MAT2(k,j)
Nextk
Nextj
Nexti
MATMUL=MatXX
EndFunction
'********************************************************
'函数名：Jacobi
'函数用途:用Jacobi法求A的特征值和特征向量
'参数说明：矩阵下标为1：N,从1开始；
'A:调用时存放实对称矩阵,A(N,N)
'N:矩阵长度
'Bx:返回时存放矩阵的全部特征值,B(N)
'Vx:特征向量为列组成的矩阵,V(N,N),其中第i列为与第i个特征值相对应的特征向量
'EPS:存放精度要求
'*********************************************************
PrivateSubJacobi(A()AsSingle,_
NAsInteger,_
EPSAsSingle,_
ByRefVx()AsSingle,_
ByRefBx()AsSingle,_
ByRefLAsInteger)
'A:调用时存放实对称矩阵
'Bx:返回时存放矩阵的全部特征值
'Vx:存放特征向量,其中第i列为与第i个特征值相对应的特征向量
'EPS:存放精度要求
ReDimVx(1ToN,1ToN)
ReDimBx(1ToN)
DimFMAsSingle,CNAsSingle,SNAsSingle
DimOmegaAsSingle,XAsSingle,YAsSingle
DimPAsInteger,QAsInteger
DimiAsInteger,jAsInteger,kAsInteger
L=1
'初始化特征向量矩阵
Fori=1ToN
Forj=1ToN
Vx(i,j)=0
Nextj
Vx(i,i)=1
Nexti
'如果计算次数大于给定次数，则跳出循环
DoWhile(L<100)
FM=0
'查找矩阵中非对角元素的最大值，并记录其位置（P,Q)
Fori=1ToN
Forj=1ToN
If(ijAndAbs(A(i,j))>FM)Then
FM=Abs(A(i,j))
P=i
Q=j
EndIf
Nextj
Nexti
'如果最大值小于给定最小值，则跳出循环
If(FM L=-1 ExitDo EndIf L=L+1 X=-A(P,Q) Y=(A(Q,Q)-A(P,P))/2 Omega=X/(X*X+Y*Y) If(Y<0)ThenOmega=-Omega SN=1+(1-Omega*Omega) SN=Omega/(2*SN) CN=(1-SN*SN) FM=A(P,P) A(P,P)=FM*CN*CN+A(Q,Q)*SN*SN+A(P,Q)*Omega A(Q,Q)=FM*SN*SN+A(Q,Q)*CN*CN-A(P,Q)*Omega A(P,Q)=0'处理完当前最大非对角元素值，赋值0，下一循环不再计算 A(Q,P)=0'处理完当前最大非对角元素值，赋值0，下一循环不再计算 Forj=1ToN If(jPAndjQ)Then FM=A(P,j) A(P,j)=FM*CN+A(Q,j)*SN A(Q,j)=-FM*SN+A(Q,j)*CN EndIf Nextj Fori=1ToN If(iPAndiQ)Then FM=A(i,P) A(i,P)=FM*CN+A(i,Q)*SN A(i,Q)=-FM*SN+A(i,Q)*CN EndIf Nexti Fori=1ToN FM=Vx(i,P) Vx(i,P)=FM*CN+Vx(i,Q)*SN Vx(i,Q)=-FM*SN+Vx(i,Q)*CN Nexti Fori=1ToN Bx(i)=A(i,i) Nexti Loop '将特征值按大小排列 m=N DoWhile(m>0) j=m-1 m=0 Fori=1Toj If(Bx(i) B1=Bx(i) Bx(i)=Bx(i+1) Bx(i+1)=B1 m=i Fork=1ToN V1=Vx(k,i) Vx(k,i)=Vx(k,i+1) Vx(k,i+1)=V1 Nextk EndIf Nexti Loop EndSub