第三章多元统计分析.docx

1、第三章多元统计分析第三章 多元统计分析3方差分析程序3.2.1功能对某种气象要素的时间序列进行方差分析。3.2.2方法说明3.2.3子程序语句SUBROUTINE FANGCHA(X,N,XBAR,S,S2)3.2.4哑元说明X输入参量，一维实型数组，大小为N，存放气象要素观测值。N输入参量，存放序列的长度。XBAR输出参数，存放x的平均值。S输出参数，存放x的标准差。S2输出参数，存放x的方差。3.2.5子程序SUBROUTINE FANGCHA(X,N,XBAR,S,S2) IMPLICIT NONE INTEGER:N,I REAL(4),DIMENSION(N):X REAL(4):X

2、BAR,S,S2 XBAR=0 DO I=1,N XBAR=XBAR+X(I) END DO XBAR=XBAR/N S2=0 DO I=1,N S2=S2+(X(I)-XBAR)*2 END DOS2=S2/N S=SQRT(S2) END3.2.6例以某海区22年逐月的海表面温度为例，计算如下。 PROGRAM MAIN IMPLICIT NONE INTEGER,PARAMETER:N=12*22 INTEGER:I REAL(4),DIMENSION(N):X REAL(4):XBAR,S,S2 OPEN(10,FILE=AA2.DAT) DO I=1,N READ(10,(F8.2)

3、X(I) END DO CLOSE(10) CALL FANGCHA(X,N,XBAR,S,S2) OPEN(12,FILE=FANGCHA.DAT) WRITE(12,(4X,XBAR,7X,S,9X,S2) WRITE(12,(3F10.4)XBAR,S,S2 CLOSE(12) END计算结果如下：XBAR= 22.5718 S= 4.2888 S2= 18.3937回归分析程序3.3.1一元线性回归分析 当考虑预报量只与某一个因子有关时，我们用一元线性回归，它是回归分析中最简单的一种。3.3.1.1功能对于给定的n个数据点(xi，yi)(i=1,2,3,n)，用直线y=ax+b作回归分

4、析。3.3.1.2方法说明3.3.1.3子程序语句SUBROUTINE YXHG(N,X,Y, A,B,Q,S,P,UMAX,UMIN,U)3.3.1.4哑元说明N整型变量，输入参数。观测数列长度X长度为N的一维实型数组，输入参数。存放自变量x的N个值。Y长度为N的一维实型数组，输入参数。存放与自变量x相对应的N个Y观测值。A实型变量，输出参数。回归系数，即线性回归方程的一次项系数。B实型变量，输出参数。回归系数，即线性回归方程的常数项。Q实型变量，输出参数。偏差平方和。S实型变量，输出参数。平均标准偏差。P实型变量，输出参数。回归平方和。UMAX实型变量，输出参数。最大偏差。UMIN实型变量

5、，输出参数。最小偏差。U实型变量，输出参数。偏差的平均值。3.3.1.5子程序（子程序名为：YXHG） SUBROUTINE YXHG(N,X,Y,A,B,Q,S,P,UMAX,UMIN,U) REAL(4),DIMENSION(N):X,Y XV=0 ！x的均值 YV=0 ！y的均值 DO I=1,N XV=XV+X(I) YV=YV+Y(I) END DO XV=XV/N YV=YV/N DXX=0.0 DXY=0.0 DO I=1,N Q=X(I)-XV DXX=DXX+Q*Q DXY=DXY+Q*(Y(I)-YV) END DO A=DXY/DXX B=YV-A*XV Q=0 U=0

6、P=0 UMAX=-1.0E10 UMIN=1.0E20 DO I=1,N S=A*X(I)+B Q=Q+(Y(I)-S)*2 P=P+(S-YV)*2 DX=ABS(Y(I)-S) IF(DX.GT.UMAX)THEN UMAX=DX ENDIF IF(DX.LT.UMIN)THEN UMIN=DX ENDIF U=U+DX/N END DO S=SQRT(Q/N) END3.3.1.6例下面以某海区22年的1月平均气温作x，对应的海表温度作y，编写计算的主程序，并给出计算结果。 PROGRAM main1 INTEGER,PARAMETER:N=22 REAL(4),DIMENSION(N

7、):X,Y DATA X/16.80,15.60,15.60,14.00,18.95,18.50,16.50,15.35,17.50,17.00,17.30,& 14.30,13.75,14.00,12.00,16.30,14.40,12.90,13.00,15.00,13.70,14.00/ DATA Y/16.20,15.80,16.20,18.00,17.70,18.00,16.00,19.00,18.40,15.35,17.00,& 14.00,14.65,13.30,14.35,15.00,14.35,12.50,14.00,14.50,13.25,14.10/ CALL YXHG(

8、N,X,Y,A,B,Q,S,P,UMAX,UMIN,U) OPEN(9,FILE=RESULT.DAT) WRITE(9,10)A,B10 FORMAT(1X,A=,E13.6,2X,B=,E13.6) WRITE(9,15)Q,S,P15 FORMAT(1X,Q=,E13.6,2X,S=,E13.6,2x,P=,E13.6) WRITE(9,20)UMAX,UMIN,U20 FORMAT(1X,UMAX=,E13.6,2X,UMIN=,E13.6,2x,U=,E13.6) CLOSE(9) END计算结果： A= .680907E+00 B= .511632E+01 Q= .372275E+

9、02 S= .130083E+01 P= .343258E+02 UMAX= .343176E+01 UMIN= .318956E-01 U= .939135E+00432 多元线性回归分析3.3.2.1功能 假定预报量y与m个因子(x1，x2，xm)的关系是线性的，对于它的n组观测值（x1i，x2i，xmi）(i=1，2，n)作线性回归分析。3.3.2.2方法说明3.3.2.3子程序语句SUBROUTINE DXHG(X,Y,M,N,A,Q,S,R,V,U,B)3.3.2.4哑元说明X实型二维数组，大小为MN，输入参数。其中每一列存放m个自变量的一组观测值，即每一列为，i=1,2,NY实型一

10、维数组，长度为，输入参数。存放的个观测值。M整型变量，输入参数。自变量个数。N整型变量，输入参数。观测数据的组数。A实型一维数组，长度为M+1，输出参数。存放明M+1个回归系数a1，a2，am+1。Q实型变量，输出参数。偏差平方和。S实型变量，输出参数。平均标准偏差。R实型变量，输出参数。复相关系数。V实型一维数组，长度为M，输出参数。M个自变量的偏相关系数。U实型变量，输出参数。回归平方和。B实型二维数组，大小为(M+1)(M+1)。工作数组，存放CCT。3.3.2.5子程序（子程序名为：DXHG） SUBROUTINE DXHG(M,N,X,Y,A,Q,S,R,V,U,B) REAL(KI

11、ND=8),DIMENSION(M,N):X REAL(KIND=8),DIMENSION(N):Y REAL(KIND=8),DIMENSION(M+1):A REAL(KIND=8),DIMENSION(M+1,M+1):B REAL(KIND=8),DIMENSION(M):V REAL(KIND=8) Q,S,R,U,YY,DYY,P,PP MM=M+1 B(1,1)=N DO J=2,MM B(1,J)=0 DO I=1,N B(1,J)=B(1,J)+X(J-1,I) END DO B(J,1)=B(1,J) END DO DO I=2,MM DO J=I,MM B(I,J)=0

12、DO K=1,N B(I,J)=B(I,J)+X(I-1,K)*X(J-1,K) END DO B(J,I)=B(I,J) END DO END DO A(1)=0 DO I=1,N A(1)=A(1)+Y(I) END DO DO I=2,MM A(I)=0 DO J=1,N A(I)=A(I)+X(I-1,J)*Y(J) END DO END DO CALL CHOLESKY(B,MM,1,A,L) YY=0 DO I=1,N YY=YY+Y(I) END DO YY=YY/N Q=0 DYY=0 U=0 DO I=1,N P=A(1) DO J=1,M P=P+A(J+1)*X(J,I)

13、 END DO Q=Q+(Y(I)-P)*(Y(I)-P) DYY=DYY+(Y(I)-YY)*(Y(I)-YY) U=U+(YY-P)*(YY-P) END DO S=SQRT(Q/N) R=SQRT(1-Q/DYY) DO J=1,M P=0 DO I=1,N PP=A(1) DO K=1,M IF(K/=J)PP=PP+A(K+1)*X(K,I) END DO P=P+(Y(I)-PP)*(Y(I)-PP) END DO V(J)=SQRT(1-Q/P) END DO END SUBROUTINE CHOLESKY(C,N,M,D,L) REAL(KIND=8),DIMENSION(N,

14、N):C REAL(KIND=8),DIMENSION(N,M):D L=1 IF(ABS(C(1,1)1.0E-10)THEN L=0 WRITE(*,( FAIL) RETURN END IF C(1,1)=SQRT(C(1,1) DO J=2,N C(1,J)=C(1,J)/C(1,1) END DO DO I=2,N DO J=2,I C(I,I)=C(I,I)-C(J-1,I)*C(J-1,I) END DO IF(ABS(C(I,I)1.0E-10)THEN L=0 WRITE(*,( FAIL) RETURN END IF C(I,I)=SQRT(C(I,I) IF(I/=N)T

15、HEN DO J=I+1,N DO K=2,I C(I,J)=C(I,J)-C(K-1,I)*C(K-1,J) END DO C(I,J)=C(I,J)/C(I,I) END DO END IF END DO DO J=1,M D(1,J)=D(1,J)/C(1,1) DO I=2,N DO K=2,I D(I,J)=D(I,J)-C(K-1,I)*D(K-1,J) END DO D(I,J)=D(I,J)/C(I,I) END DO END DO DO J=1,M D(N,J)=D(N,J)/C(N,N) DO K=N,2,-1 DO I=K,N D(K-1,J)=D(K-1,J)-C(K-

16、1,I)*D(I,J) END DO D(K-1,J)=D(K-1,J)/C(K-1,K-1) END DO END DO END3.3.2.6、例以南京冬季西路冷锋移速Y（公里/小时）为预报对象，以西路冷锋地区（过银川站）时700hPa图上东胜与锡林浩特两站位势高度差为x1，它反映了高空引导气流的强弱；x2为850hPa图上地面锋后风速的垂直分量；x3为冷锋过银川站时的锋后气压梯度，取自锋线至锋后冷高压中心（hPa/纬距），它反映了锋后高压的强度对冷锋移速的影响；x4为700hPa图上东胜与锡林浩特的位势高度差（即x1）乘以700hPa图上太原与二连浩特站的温度差和等温线与等高线间的夹角的正

17、弦。根据50个观测资料（资料见下表）来作回归分析。 序号 y x1 x2 x3 x4 1 47 69 5.9 1.90 2.8 2 61 119 10.4 2.32 7.0 3 62 132 10.2 2.35 7.3 4 50 87 7.7 1.72 1.6 5 57 107 9.9 2.10 10.9 6 34 90 6.2 1.52 3.1 7 48 115 8.0 1.80 7.1 8 43 84 9.0 1.83 14.2 9 37 100 7.0 1.68 5.0 10 54 99 7.7 2.25 9.9 11 61 123 8.7 2.17 8.6 12 32 85 6.2 1

18、.52 1.7 13 43 99 7.0 2.00 3.8 14 55 127 11.5 2.09 11.1 15 48 95 9.2 2.13 8.8 16 40 104 7.7 1.78 5.3 17 42 90 7.9 1.59 4.1 18 38 85 6.2 1.77 3.1 19 55 97 10.0 2.00 10.9 20 39 120 6.4 1.54 11.4 21 36 90 7.5 1.60 4.6 22 42 80 6.0 1.75 4.4 23 51 92 7.7 1.86 6.8 24 40 86 9.2 1.98 4.9 25 55 124 10.2 2.38

19、11.1 26 65 124 12.0 2.45 12.4 27 48 115 9.9 2.00 7.1 28 43 95 7.5 1.93 3.6 29 39 97 6.4 1.61 4.0 30 32 75 3.9 1.66 1.4 31 59 145 12.8 2.20 8.8 32 41 133 7.9 1.54 3.1 33 63 98 12.3 2.61 12.7 34 52 106 7.9 2.17 9.0 35 38 80 5.1 1.68 1.5 36 49 112 8.5 1.76 8.4 37 33 84 5.1 1.58 4.3 38 56 97 7.7 1.92 5.

20、0 39 43 105 9.6 1.95 8.0 40 44 80 6.0 1.60 3.3 41 48 82 5.2 1.94 4.3 42 37 67 6.6 1.85 6.4 43 55 105 10.9 2.15 15.2 44 28 63 3.0 1.16 1.4 45 32 54 5.0 1.34 .8 46 29 80 5.2 1.54 4.3 47 42 108 7.7 1.92 5.7 48 37 100 7.7 1.53 4.4 49 47 146 9.6 1.82 10.2 50 42 81 7.7 2.10 8.8 PROGRAM main2 INTEGER,PARAM

21、ETER:N=50 INTEGER,PARAMETER:M=4 REAL(KIND=8),DIMENSION(M,N):X REAL(KIND=8),DIMENSION(N):Y REAL(KIND=8),DIMENSION(M+1):A REAL(KIND=8),DIMENSION(M+1,M+1):B REAL(KIND=8),DIMENSION(M):V REAL(KIND=8) Q,S,R,U OPEN(10,FILE=DXHG.DAT) DO I=1,N READ(10,*)I,Y(I),X(1,I),X(2,I),X(3,I),X(4,I) END DO CLOSE(10) MM=

22、M+1 CALL DXHG(M,N,X,Y,A,Q,S,R,V,U,B) OPEN(9,FILE=RESULT2.DAT) WRITE(9,10)(I,A(I),I=1,MM)10 FORMAT(1X,A(,I2,)=,E13.6) WRITE(9,15)Q,S,R15 FORMAT(1X,Q=,E13.6,2X,S=,E13.6,2x,R=,E13.6) WRITE(9,20)(I,V(I),I=1,M)20 FORMAT(1X,V(,I2,)=,E13.6) CLOSE(9) END计算结果为： A( 1)= -.699607E+01 A( 2)= .755634E-01 A( 3)= .

23、871720E+00 A( 4)= .204831E+02 A( 5)= -.396034E-01 Q= .887608E+03 S= .421333E+01 R= .896120E+00 V( 1)= .874766E+00 V( 2)= .860884E+00 V( 3)= .994168E+00 V( 4)= .694885E-01 U= .361871E+043判别分析程序3.4.1功能给出判别分析的程序。3.4.2方法说明3.4.2.1二级判别的Fisher准则概念3.4.2.2多级判别的Fisher准则概念3.4.2.3判别函数的显著性检验3.4.3子程序语句 SUBROUTINE

24、 PBFX(X,N,M,G,KG)3.4.4哑元说明X输入参数，二维实型数组，大小为N*M，存放M个要素的N次观测值。N输入参数，整型变量，存放样本的长度。M输入参数，整型变量，存放要素的个数。NG输入参数，一维整型数组，大小为G，存放分类的数目。KG输入参数，一维整型数组，大小为N，存放样本所属类数。3.4.5子程序（子程序名为：PBFX）! 判别分析程序 ! M个因子,序列总长度为N,分为G类,每类的个数为NG(G),N=NG(1)+NG(2)+.+NG(G) ! 每类的均值为XV(M,G),总的均值为XVV(M),SUBROUTINE PBFX(X,N,M,G,KG) ! M个因子,序列

25、总长度为N,分为G类,每类的个数为NG(G),N=NG(1)+NG(2)+.+NG(G) ! 每类的均值为XV(M,G),总的均值为XVV(M), IMPLICIT NONE INTEGER:G,I,J,K,N,M,L,MG REAL(8),DIMENSION(N,M):X REAL(8),DIMENSION(M):XVV INTEGER,DIMENSION(G):NG INTEGER,DIMENSION(N):KG REAL(8),DIMENSION(M,G):XV REAL(8),DIMENSION(M,M):T,S,B REAL(8),DIMENSION(M):ST,DT,XX1 REA

26、L(8),DIMENSION(M,M):ST12,ST_12,SS_12,SS,V,S12,S_12,VD,D,VS REAL(8),DIMENSION(:),ALLOCATABLE:AA,X2,X22 REAL(8),DIMENSION(:,:,:),ALLOCATABLE:XG REAL(8),DIMENSION(:,:,:),ALLOCATABLE:DS INTEGER,DIMENSION(:,:),ALLOCATABLE:KGJ INTEGER:L0,NP,NP1,N0,NN REAL(8):AA1,XX2,XX3,DMIN,ALF,H NG=0 DO I=1,N NG(KG(I)=N

27、G(KG(I)+1 END DO MG=MAXVAL(NG) !求数组的最大值 ALLOCATE(XG(M,G,MG) XG=0 NG=0 DO I=1,N NG(KG(I)=NG(KG(I)+1 DO J=1,M XG(J,KG(I),NG(KG(I)=X(I,J) END DO END DO! 求每个因子总的均值 XVV=0 DO K=1,M DO I=1,N XVV(K)=XVV(K)+X(I,K) END DO XVV(K)=XVV(K)/N END DO WRITE(12,(4X,总的均值=,F8.2)XVV! 求每个因子各组的均值 XV=0 DO K=1,M DO I=1,G DO

28、 J=1,NG(I) XV(K,I)=XV(K,I)+XG(K,I,J) END DO XV(K,I)=XV(K,I)/NG(I) END DO END DO WRITE(12,*) DO I=1,G WRITE(12,(第,I2,组的均值=,F8.2)I,(XV(J,I),J=1,M) END DO ! 求总的离差平方和阵T和组内间离差平方和阵S T=0 DO K=1,M DO L=1,M DO I=1,G DO J=1,NG(I) T(K,L)=T(K,L)+(XG(K,I,J)-XVV(K)*(XG(L,I,J)-XVV(L) S(K,L)=S(K,L)+(XG(K,I,J)-XV(K,

29、I)*(XG(L,I,J)-XV(L,I) END DO END DO END DO END DO WRITE(12,*) WRITE(12,(总的离差阵T) WRITE(12,(F9.2)T WRITE(12,*) WRITE(12,(组内的离差阵S) WRITE(12,(F9.2)S! 求组间离差平方和阵B，可以采用两种方法：1. B=T-S；2. 直接计算法。两种方法计算结果完全一致。! 第一种方法:B=T-S B=T-S! 第二种方法:直接计算法! B=0! DO K=1,M! DO L=1,M! DO I=1,G! B(K,L)=B(K,L)+NG(I)*(XV(K,I)-XVV(K)*(XV(L,I)-XVV(L)! END DO! END DO! END DO WRITE(12,*) WRI