自然正交函数分析(EOF)程序.docx
- 格式:docx
- 大小:63.26 KB
- 文档页数:7
5.4 主成分聚类与主成分回归5.4.1 变量聚类与样品分类主成分分析可用于聚类:变量聚类与样品聚类。
变量聚类:由主成分系数的差异,可将变量聚类。
例如例5.5中第2主成分中murder,rape, assult系数为负的, burglary,larceny, auto系数是正的。
按系数正负可把7个变量分为两类: murder, rape, assult属于暴力程度严重的一类;burglary,larceny,auto属于暴力程度较轻的一类。
按照这种方法,根据主成分系数的正负可以将变量聚类。
样品聚类:如果2个主成分能很好的概括随机向量的信息,计算每个样品的这两个主成分得分,把他们的散点图画出来,就能从图上将样品分类。
例5.5(续2)按照第一、第二主成分得分,画出散点图data crime; /*建立数据集crime*/input state $ 1-15 murder rape robbery assult burglary larceny auto;/*建立变量state murder rape robbery assult burglary larceny auto。
state $ 1-15表示前15列存州名。
murder rape robbery assult burglary larceny auto 表7种罪的犯罪率*/cards; /*以下为数据体*/Albama 14.2 25.2 96.8 278.3 1135.5 1881.9 280.7Alaska 10.8 51.6 96.8 284.0 1331.7 3369.8 753.3Arirona 9.5 34.2 138.2 312.3 2346.1 4467.4 439.5Arkansas 8.8 34.2 138.2 312.3 2346.1 4467.4 439.5Califonia 11.5 49.4 287.0 358.0 2139.4 3499.8 663.5Colorado 6.3 42.0 170.7 292.9 1935.2 3903.2 477.1Conecticat 4.2 16.8 129.5 131.8 1346.0 2620.7 593.2Delaware 6.0 24.9 157.0 194.2 1682.6 3678.4 467.0Florida 10.2 39.6 187.9 449.1 1859.9 3840.5 351.4Geogia 11.7 31.1 140.5 256.5 1351.1 2170.2 297.9Hawaii 7.2 25.5 128.0 64.1 1911.5 3920.4 489.4Idaho 5.5 19.4 39.6 172.5 1050.8 2599.6 237.6Illinois 9.9 21.8 211.3 209.0 1085.0 2828.5 528.6Indiana 7.4 26.5 123.2 153.5 1086.2 2498.7 377.4Iowa 2.3 10.6 41.2 89.8 812.5 2685.1 219.9Kansas 6.6 22.0 100.7 180.5 1270.4 2739.3 244.3Kentaky 10.1 19.1 81.1 123.3 872.2 1662.1 245.4Loisana 15.5 30.9 142.9 335.5 1165.5 2469.9 337.7Maine 2.4 13.5 38.7 170.0 1253.1 2350.7 246.9Maryland 8.0 34.8 292.1 358.9 1400.0 3177.7 428.5Masschusetts 3.1 20.8 169.1 231.6 1532.2 2311.3 1140.1Michigan 9.3 38.9 261.9 274.6 1522.7 3159.0 545.5Minnesota 2.7 19.5 85.9 85.8 1134.7 2559.3 343.1Mississippi 14.3 19.6 65.7 189.1 915.6 1239.9 144.4Missouri 9.6 28.3 189.0 233.5 1318.3 2424.2 378.4Montana 5.4 16.7 39.2 156.8 804.9 2773.2 309.3Nebraska 3.9 18.1 64.7 112.7 760.0 2316.1 249.1Nevada 15.8 49.1 323.1 355.0 2453.1 4212.6 559.2Mew Hampashare 3.2 10.7 23.2 76.0 1041.7 2343.9 293.4New Jersey 5.6 21.0 180.4 185.1 1435.8 2774.5 511.5New Maxico 8.8 39.1 109.6 343.4 1418.7 3008.6 259.5New York 10.7 29.4 472.6 319.1 1728.0 2782.0 745.8North Carolina 10.6 17.0 61.3 318.3 1154.1 2037.8 192.1North Dakoda 100.9 9.0 13.3 43.8 446.1 1843.0 144.7Ohio 7.8 27.3 190.5 181.1 1216.0 2696.8 400.4Oklahoma 8.6 29.2 73.8 205.0 1288.2 2228.1 326.8Oregan 4.9 39.9 124.1 286.9 1636.4 3506.1 388.9Pennsyvania 5.6 19.0 130.3 128.0 877.5 1624.1 333.2Rhode Island 3.6 10.5 86.5 201.0 1849.5 2844.1 791.4South Carolina 11.9 33.0 105.9 485.3 1613.6 2342.4 245.1South Dakoda 2.0 13.5 17.9 155.7 570.5 1704.4 147.5Tennessee 10.1 29.7 145.8 203.9 1259.7 1776.5 314.0Texas 13.3 33.8 152.4 208.2 1603.1 2988.7 397.6Utah 3.5 20.3 68.8 147.3 1171.6 3004.6 334.5Vermont 1.4 15.9 30.8 101.2 1348.2 2201.0 265.2Virginia 9.0 23.3 92.1 165.7 986.2 2521.2 226.7Wasinton 4.3 39.6 106.2 224.8 1605.6 3386.9 360.3West Viginia 6.0 13.2 42.2 90.9 597.4 1341.7 163.3Wiskonsin 2.8 12.9 52.2 63.7 846.9 2614.2 220.7Wyoming 5.4 21.9 39.7 173.9 811.6 2772.2 282.0;proc princomp out=crimprin n=2;var murder rape robbery assult burglary larceny auto;run;PROC PLOT data=crimprin;PLOT PRIN2*PRIN1=STATE/VPOS=31;TITLE2 ‘PLOT OF THE FIRST TWO PRINCIPAL COMPONENTS’;RUN;例5.7 (气温分析)本例的输入资料文件(TEMPERA T)是美国六十四个城市一月与七月的平均日温。
自然正交函数分析程序EOF方法是基于假设,即数据可以被表示为一系列正交函数的线性组合。
这些正交函数称为EOF模态,并代表了数据中的主要模式。
每个EOF 模态都具有相应的权重,称为EOF系数,用于描述该模态在总方差中的贡献程度。
EOF方法的步骤如下:1.数据预处理:首先,要对原始数据进行预处理。
这可以包括去除重复数据、去除异常值、进行数据平滑处理等。
2.协方差矩阵计算:接下来,需要计算数据的协方差矩阵。
协方差矩阵描述了数据中不同维度之间的相关性。
3.特征值分解:通过对协方差矩阵进行特征值分解,可以得到特征值和特征向量。
特征值表示了每个特征向量对总方差的贡献程度。
4.选择模态:根据特征值的大小,可以选择保留最重要的EOF模态,从而降低数据维度。
5.计算EOF系数:对于每个选定的EOF模态,可以计算其相应的EOF 系数,用于描述该模态在数据中的贡献程度。
6.重构数据:最后,通过将所有选定的EOF模态与相应的EOF系数进行线性组合,可以重构原始数据。
这样可以去除一些噪音和次要特征,从而提取出原始数据中的关键模式。
EOF方法有许多应用,特别是在气候学、地球物理学和图像处理等领域。
在气候学中,EOF可以帮助我们理解地球上不同地区的温度、降水和风向等变化模式。
在地球物理学中,EOF可以帮助我们分析地震数据、地磁数据和重力场数据等。
在图像处理中,EOF可以帮助我们提取图像中的关键特征,用于图像分类和识别。
总之,自然正交函数分析(EOF)是一种强大的数学工具,用于处理具有时间或空间相关性的数据。
通过对数据进行正交分解,EOF可以提取出关键的时间或空间模式,并帮助我们理解和分析数据中的重要特征。
eof经验正交函数正交设计是一种在实验设计中应用广泛的方法,它可以有效地降低实验的复杂性并提高实验结果的可信度。
在正交设计中,经验正交函数(EOF)是一种重要的工具,它可以帮助研究人员选择合适的实验因素和水平,以获得准确且可靠的实验结果。
经验正交函数是一组具有正交性质的基础函数,它们可以表示实验因素的不同水平。
通过使用经验正交函数,研究人员可以在有限的实验次数下,对多个因素进行全面的测试,从而节省时间和资源。
在实际应用中,经验正交函数可以用于设计和优化各种工程和科学实验。
例如,在材料科学领域,研究人员可以使用经验正交函数来优化材料的物理和化学性质。
在制造业中,经验正交函数可以用于优化生产过程中的参数设置,以提高产品质量和生产效率。
经验正交函数的选择和使用需要考虑因素的数量和水平,以及实验的目标和约束条件。
通常,研究人员可以使用经验正交函数表来选择合适的函数和水平组合。
经验正交函数表是经过统计和数学分析得出的,可以帮助研究人员快速准确地选择合适的函数。
除了选择适当的经验正交函数,研究人员还需要确定实验的因素和水平。
因素是指影响实验结果的变量,而水平则是指每个因素的具体取值。
通过选择合适的因素和水平,研究人员可以在实验设计中获得准确和可靠的结果。
在进行实验时,研究人员需要根据经验正交函数和选择的因素水平制定实验计划。
实验计划应考虑到实验次数、实验顺序和实验条件等因素,以确保实验结果的可靠性和可重复性。
经验正交函数的应用可以帮助研究人员在有限的实验条件下,获得准确和可靠的实验结果。
它不仅可以提高实验效率,还可以降低实验成本。
因此,经验正交函数在工程和科学领域中得到了广泛的应用和重视。
经验正交函数是一种重要的实验设计工具,可以帮助研究人员选择合适的实验因素和水平,以获得准确且可靠的实验结果。
它在各个领域的应用都取得了显著的成果,并为实验设计提供了有效的方法和策略。
通过合理地运用经验正交函数,研究人员可以在实验中更好地发现和理解因果关系,为工程和科学领域的发展做出贡献。
'*************************************' 全局变量,便于主函数调用。
' 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 i As 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)。
5・3自然正交函数分析(EOF)程序近年來,自然正交函数(乂称经验正交函数)展开在气象上应用比较广泛。
这种正交函数展开不彖三角函数展开、球函数展开那样有固定的展开形式。
它无固定的函数形式,不是事先人为地给定典型场函数,图形是由场木身来决定的,它具有收敛快又能更好地反映岀场的基木结构的特征。
它可以在有限的区域屮进行,既可以取空间不同站点进行分解,也可以对同一站点的不同吋间、不同高度的多种要素进行综和分析。
因此它在气彖中具有广泛的应用,可用于气象要素场分析、大气垂直结构分析、动力模型垂直分层等。
5. 3.1功能计算要素场的自然正交函数分解。
5. 3. 2方法说明口然止交函数分解是针对气彖要素场进行的,它的基本思想是把包含P个空间点(或P个变量)的n个时次的观测场随时间进行分解,即将某一区域的气象要素场序列Fq (i=l, 2,・・・,p;j=l,2,…,n,即p个空间点的n个时次的观测资料)分解成相互正交的时间函数与相互正交的空间函数的乘积Z和,常把空间函数VW看作典型场,时间函数看作典型场的权重系数,则不同时间的要素场是若干个典型场按不同权重线性叠加的结果,各个场之间的差别就在于各典型场的系数不同。
则气象耍素场可以表示为PEj =》%tkj = Vig+Vj2t2j+・・・+Viptpj (5. 3. 1)k=l英中Fq表示第i个场中的第j个测点的观测值。
可将(5.3.1)是写为矩阵的形式F =VT( 5 . 3 . 2 )式中F为pxn阶的均值为0的资料阵,V为pxp阶的空间函数阵,卩为pxn阶的时间函数阵。
由于V和0是根据场的资料阵F进行分解而得到的,分解的函数没有固定的函数形式,因而称为“经验”的,另外,我们还要求这种分解具有“正交”性,即要求满足下式PVk V, =X v ik v n =0 (kHl)i=1(5. 3. 3 )n兀齐-£,kjtij =0 (k H 1)冃事实上,我们对(5. 3. 2)式右乘厂可得FF =VTTV r( 5 . 3 . 4 )因FF'是pxp阶对称阵,其元素为距平变量的交义积。
根据实对称矩阵的分解定理有FF =VAV f(5. 3. 5 )其小A是FF'矩阵的特征值组成的对角阵,V是对应的特征向量为列向量组成的矩阵。
比较(5.3. 4)和(5. 3. 5 )式可知TT r = A(5 . 3. 6 )乂根据特征向量的性质有VV =W=I式中为I单位矩阵。
显然(5 . 3 . 6 )和(5 . 3. 7 )式满足(5 . 3 . 3 )式的要求。
由此町知空间函数矩阵可从FF'矩阵的特征向量求得,而时间函数则可利用(5. 3. 2)式左乘V'得到,即T =VF( 5. 3. 8)5. 3. 3子程序语句CALL EOF (X, P, N, XF)5. 3. 4哑元说明X——输入变量,二维实型数组,大小为PxN,存放原始观测值。
P——输入整型变量,空间格点数。
N——输入整型变量,序列的时间长度。
XF——输出变量,二维实型数组,人小为PxN,存放恢复值。
5. 3. 5子程序SUBROUTINE EOF(X,P,N,LW,XF)INTEGER::PINTEGER: :NINTEGER::LWREAL(8),DIMENSION(P,N)::X,XFREAL(8),DIMENSION(P,P):: A, V, V1REAL ⑻,DIMENSIONS,N)::TREAL ⑻,DIMENSION(P)::B,GM,GAREAL(8),DIMENSIONS,LW)::VFREAL(8),DIMENSION(LW,N)::TF! 求X乘以X的转置,即A=XX 'DO 1=1,PDO J=1,PA(I,J)=0DO K=1,NA(I,J)=A(I,J)+X(I,K)*X(J,K)END DOEND DOEND DO! 用Jacobi法求A的特征值和特征向屋!返回时B存放矩阵的全部特征值,V存放特征向量为列组成的矩阵CALL JCB(A,P, 1.0E・6,V,B,L)DO 1= 1 ,PGA(I)=()DO J=1,IGA(I)=GA(I)+B(J)END DOEND DODO 1=1,PGM(I)=GA(I)/GA(P)END DODO 1=1,PDO J=1,PV1(I,J)=V(J,I)END DOEND DOT=MATMUL(V1,X)WRITE(12;(H特征值”))WRITE( 12;(<P>I10)')(I,I= 1 ,P)WRITE( 12;(3X,<P>D 10.4)*)BWRITE(12;(H解释的方差(%)u nWRITE( 12;(<P>I7),)(I,I= 1 ,P)WRIT E( 12;(3X,<P>F7.2)')GM* 100WRITE(12,(“特征向量为列组成的矩阵,即空间函数V?)WRITE(12;(<P>F7.4),)((V(I,J)J= 1 ,P),I= 1 ,P)WRITE( 12,(“ 时间函数T")1)WRITE( 12;(<P>F10.4),)((T(I,J),J=1 ,P),I= 1 ,N)DO 1=1,PDO J=1,LWVF(I,J)=V(IJ)END DOEND DODO 1=1,LWDO J=1,NTF(I,J)=T(I,J)END DOEND DOXF=MATMUL(VF,TF)ENDSUBROUTINE JCB(A,N,EPS,V,B,L)! A:调用时存放实对称矩阵! B:返回时存放矩阵的全部特征值! V:存放特征向量,其中笫i列为与第i个特征值相对应的特征向屋! EPS:存放精度要求REAL(8),D1MENS1ON(N,N)::A,VREAL(8),DIMENSION(N)::BREAL ⑻::FM,CN,SN,OMEGA,X,YINTEGER:: P,QL=1V=0DO 1= 1 ,NV(I,I)=1END DO 10 FM=0DO 1=I,NDO J=1,NIF(I/=J.AND.ABS(A(I,J))>FM)THENFM=ABS(A(IJ))P=1Q=JEND IFEND DOEND DOIF(FMvEPS)THENL=1GOTO 15END IFIF(L>100)THENL=0WRITE(12,'(“L=”,I3,2X, “没有达到精度要求” ))L GOTO 15END IFL=L+1X=A(P,Q)Y=(A(Q,Q)・A(P,P))/2OMEGA二X/SQRT(X*X+Y*Y)IF(Y<0)OMEGA=-OMEGASN= 1 +SQRT(1 -OMEGA*OMEGA)SN=OMEGA/SQRT(2*SN)CN=SQRT(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)=0A(Q,P)=0DO J=1,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 IFEND DODO 1=1,NIF(I/=P.AND.1/=Q)THENFM=A(I,P)A(I,P)二FM*CN+A(I,Q)*SNA(I,Q)=-FM*SN+A(I,Q)*CNEND IFEND DODO 1=1,NFM=V(I,P)V(I,P)=FM*CN+V(I,Q)*SNV(I,Q)=-FM*SN+V(I,Q)*CNEND DODO 1=1,NB(I)=A(I,1)END DOGOTO 10!将特征值按人小排列15 M=N20 IF(M>0)THENJ=M-1M=0DO 1= 1 ,JIF(B(I).LT.B(I+1))THENB1=B(I)B(I)=B(I+1)B(I+1)=BIM=IDO K=1,NV1=V(K,I)V(K,I)=V(K,I+1)V(K,I+1)=V1END DOENDIFENDDOGOTO 20ENDIFEND5. 3. 6 例対1991年5月1日至8月31日123天某纬圈上20个格点的850hPa高度场进行自然正交分解。
PROGRAM MAININTEGER ‘PARAMETER::P=20 !资料的空间点数INTEGER,PARAMETER::N= 123 !资料的吋间长度REAL©),DIMENSIONS,N)::X,XF,ERRORREAL(8),DIMENSION(P)::XV !X 的平均值INTEGER: :LP=P;LW=3OPEN( 12,FILE-F58Z850.DAT')DO IT=1,123READ(12;(<LP>F8.1),)(X(I,IT)J= 1 ,P) ENDDOCLOSE(12) XV=O DO 1=1,P DO J=1,NXV(I)=XV(I)+X(I,J) END DO XV(I)=XV(I)/N END DO DO 1=I,P DO J=1,NX(I,J)=X(I,J)-XV(I) END DO END DOOPEN( 12,FILE-EOF.DAT) CALL EOF(X,P,N,LW,XF) ERROR=X-XF DO 1=1,P DO J=1,NXF(1,J)=XF(I,J)+XV(1) END DO END DOWRITE(12;(H 恢复场?)WRITE( 12;(<LP>F8.1 )*)((XF(I, J), J= 1 ,P),I 二 1 ,N) WRITE(12;(H 误差场WRITE( 12;(<LP>F8.1 )')((ERROR(I,J),J 二 1 ,P),I= 1 ,N) CLOSE(12) END输出结果为: 特征值:.39260+03 . 3108D+03 .2749D+03 . 2655D+03 . 2395D103 .19630)0317 18 19 20.1922D+03 ・ 1148D+03 . 1110D+03 . 9482D+02 解释的方差(%)1 23 4 5 6 7 89 10 78.018& 02 93. 19 95. 7797. 25 9& 2398.7199. 02 99. 1899.3111 1213 14 15 16 17 18 19 20 99. 44 99. 53 99. 62 99. 70 99. 78 99.8499. 90 99. 94 99. 97 100. 00可见第一主分量解释的方差达78.01%,前三个主分呈解释的方差达93.19%。