第三章多元统计分析.docx

资源描述

第三章多元统计分析.docx

《第三章多元统计分析.docx》由会员分享，可在线阅读，更多相关《第三章多元统计分析.docx（37页珍藏版）》请在冰豆网上搜索。

第三章多元统计分析.docx

第三章多元统计分析

3．２方差分析程序

3.2.1功能

对某种气象要素的时间序列进行方差分析。

3.2.2方法说明

3.2.3子程序语句

SUBROUTINEFANGCHA（X,N,XBAR,S,S2）

3.2.4哑元说明

X——输入参量，一维实型数组，大小为N，存放气象要素观测值。

N——输入参量，存放序列的长度。

XBAR——输出参数，存放x的平均值。

S——输出参数，存放x的标准差。

S2——输出参数，存放x的方差。

3.2.5子程序

SUBROUTINEFANGCHA（X,N,XBAR,S,S2）

IMPLICITNONE

INTEGER:

N,I

REAL（4）,DIMENSION（N）:

REAL（4）:

XBAR,S,S2

XBAR=0

DOI=1,N

XBAR=XBAR+X（I）

ENDDO

XBAR=XBAR/N

S2=0

DOI=1,N

S2=S2+（X（I）-XBAR）**2

ENDDO

S2=S2/N

S=SQRT（S2）

END

3.2.6例

以某海区22年逐月的海表面温度为例，计算如下。

PROGRAMMAIN

IMPLICITNONE

INTEGER,PARAMETER:

N=12*22

INTEGER:

REAL（4）,DIMENSION（N）:

REAL（4）:

XBAR,S,S2

OPEN（10,FILE='AA2.DAT'）

DOI=1,N

READ（10,'（F8.2）'）X（I）

ENDDO

CLOSE（10）

CALLFANGCHA（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.5718S=4.2888S2=18.3937

４．３回归分析程序

3.3.1一元线性回归分析

当考虑预报量只与某一个因子有关时，我们用一元线性回归，它是回归分析中最简单的一种。

3.3.1.1功能

　　对于给定的n个数据点（xi，yi）（i=1,2,3,…,n），用直线y=ax+b作回归分析。

3.3.1.2方法说明

3.3.1.3子程序语句

SUBROUTINEYXHG（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——实型变量，输出参数。

最小偏差。

U——实型变量，输出参数。

偏差的平均值。

3.3.1.5子程序（子程序名为：

YXHG）

SUBROUTINEYXHG（N,X,Y,A,B,Q,S,P,UMAX,UMIN,U）

REAL（4）,DIMENSION（N）:

X,Y

XV=0！

x的均值

YV=0！

y的均值

DOI=1,N

XV=XV+X（I）

YV=YV+Y（I）

ENDDO

XV=XV/N

YV=YV/N

DXX=0.0

DXY=0.0

DOI=1,N

Q=X（I）-XV

DXX=DXX+Q*Q

DXY=DXY+Q*（Y（I）-YV）

ENDDO

A=DXY/DXX

B=YV-A*XV

Q=0

U=0

P=0

UMAX=-1.0E10

UMIN=1.0E20

DOI=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

ENDDO

S=SQRT（Q/N）

END

3.3.1.6例

下面以某海区22年的1月平均气温作x，对应的海表温度作y，编写计算的主程序，并给出计算结果。

PROGRAMmain1

INTEGER,PARAMETER:

N=22

REAL（4）,DIMENSION（N）:

X,Y

DATAX/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/

DATAY/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/

CALLYXHG（N,X,Y,A,B,Q,S,P,UMAX,UMIN,U）

OPEN（9,FILE='RESULT.DAT'）

WRITE（9,10）A,B

10FORMAT（1X,'A=',E13.6,2X,'B=',E13.6）

WRITE（9,15）Q,S,P

15FORMAT（1X,'Q=',E13.6,2X,'S=',E13.6,2x,'P=',E13.6）

WRITE（9,20）UMAX,UMIN,U

20FORMAT（1X,'UMAX=',E13.6,2X,'UMIN=',E13.6,2x,'U=',E13.6）

CLOSE（9）

END

计算结果：

A=.680907E+00B=.511632E+01

Q=.372275E+02S=.130083E+01P=.343258E+02

UMAX=.343176E+01UMIN=.318956E-01U=.939135E+00

4．3．2多元线性回归分析

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子程序语句

SUBROUTINEDXHG（X,Y,M,N,A,Q,S,R,V,U,B）

3.3.2.4哑元说明

X——实型二维数组，大小为M×N，输入参数。

其中每一列存放m个自变量的一组观测值，即每一列为

，i=1,2,……,N

Y——实型一维数组，长度为Ｎ，输入参数。

存放ｙ的Ｎ个观测值。

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）

SUBROUTINEDXHG（M,N,X,Y,A,Q,S,R,V,U,B）

REAL（KIND=8）,DIMENSION（M,N）:

REAL（KIND=8）,DIMENSION（N）:

REAL（KIND=8）,DIMENSION（M+1）:

REAL（KIND=8）,DIMENSION（M+1,M+1）:

REAL（KIND=8）,DIMENSION（M）:

REAL（KIND=8）Q,S,R,U,YY,DYY,P,PP

MM=M+1

B（1,1）=N

DOJ=2,MM

B（1,J）=0

DOI=1,N

B（1,J）=B（1,J）+X（J-1,I）

ENDDO

B（J,1）=B（1,J）

ENDDO

DOI=2,MM

DOJ=I,MM

B（I,J）=0

DOK=1,N

B（I,J）=B（I,J）+X（I-1,K）*X（J-1,K）

ENDDO

B（J,I）=B（I,J）

ENDDO

（1）=0

DOI=1,N

（1）=A

（1）+Y（I）

ENDDO

DOI=2,MM

A（I）=0

DOJ=1,N

A（I）=A（I）+X（I-1,J）*Y（J）

ENDDO

CALLCHOLESKY（B,MM,1,A,L）

YY=0

DOI=1,N

YY=YY+Y（I）

ENDDO

YY=YY/N

Q=0

DYY=0

U=0

DOI=1,N

P=A

（1）

DOJ=1,M

P=P+A（J+1）*X（J,I）

ENDDO

Q=Q+（Y（I）-P）*（Y（I）-P）

DYY=DYY+（Y（I）-YY）*（Y（I）-YY）

U=U+（YY-P）*（YY-P）

ENDDO

S=SQRT（Q/N）

R=SQRT（1-Q/DYY）

DOJ=1,M

P=0

DOI=1,N

PP=A

（1）

DOK=1,M

IF（K/=J）PP=PP+A（K+1）*X（K,I）

ENDDO

P=P+（Y（I）-PP）*（Y（I）-PP）

ENDDO

V（J）=SQRT（1-Q/P）

ENDDO

END

SUBROUTINECHOLESKY（C,N,M,D,L）

REAL（KIND=8）,DIMENSION（N,N）:

REAL（KIND=8）,DIMENSION（N,M）:

L=1

IF（ABS（C（1,1））<1.0E-10）THEN

L=0

WRITE（*,'（"FAIL"）'）

RETURN

ENDIF

C（1,1）=SQRT（C（1,1））

DOJ=2,N

C（1,J）=C（1,J）/C（1,1）

ENDDO

DOI=2,N

DOJ=2,I

C（I,I）=C（I,I）-C（J-1,I）*C（J-1,I）

ENDDO

IF（ABS（C（I,I））<1.0E-10）THEN

L=0

WRITE（*,'（"FAIL"）'）

RETURN

ENDIF

C（I,I）=SQRT（C（I,I））

IF（I/=N）THEN

DOJ=I+1,N

DOK=2,I

C（I,J）=C（I,J）-C（K-1,I）*C（K-1,J）

ENDDO

C（I,J）=C（I,J）/C（I,I）

ENDDO

ENDIF

ENDDO

DOJ=1,M

D（1,J）=D（1,J）/C（1,1）

DOI=2,N

DOK=2,I

D（I,J）=D（I,J）-C（K-1,I）*D（K-1,J）

ENDDO

D（I,J）=D（I,J）/C（I,I）

ENDDO

DOJ=1,M

D（N,J）=D（N,J）/C（N,N）

DOK=N,2,-1

DOI=K,N

D（K-1,J）=D（K-1,J）-C（K-1,I）*D（I,J）

ENDDO

D（K-1,J）=D（K-1,J）/C（K-1,K-1）

ENDDO

END

3.3.2.6、例

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

根据50个观测资料（资料见下表）来作回归分析。

序号yx1x2x3x4

147695.91.902.8

26111910.42.327.0

36213210.22.357.3

450877.71.721.6

5571079.92.1010.9

634906.21.523.1

7481158.01.807.1

843849.01.8314.2

9371007.01.685.0

1054997.72.259.9

11611238.72.178.6

1232856.21.521.7

1343997.02.003.8

145512711.52.0911.1

1548959.22.138.8

16401047.71.785.3

1742907.91.594.1

1838856.21.773.1

19559710.02.0010.9

20391206.41.5411.4

2136907.51.604.6

2242806.01.754.4

2351927.71.866.8

2440869.21.984.9

255512410.22.3811.1

266512412.02.4512.4

27481159.92.007.1

2843957.51.933.6

2939976.41.614.0

3032753.91.661.4

315914512.82.208.8

32411337.91.543.1

33639812.32.6112.7

34521067.92.179.0

3538805.11.681.5

36491128.51.768.4

3733845.11.584.3

3856977.71.925.0

39431059.61.958.0

4044806.01.603.3

4148825.21.944.3

4237676.61.856.4

435510510.92.1515.2

4428633.01.161.4

4532545.01.34.8

4629805.21.544.3

47421087.71.925.7

48371007.71.534.4

49471469.61.8210.2

5042817.72.108.8

PROGRAMmain2

INTEGER,PARAMETER:

N=50

INTEGER,PARAMETER:

M=4

REAL（KIND=8）,DIMENSION（M,N）:

REAL（KIND=8）,DIMENSION（N）:

REAL（KIND=8）,DIMENSION（M+1）:

REAL（KIND=8）,DIMENSION（M+1,M+1）:

REAL（KIND=8）,DIMENSION（M）:

REAL（KIND=8）Q,S,R,U

OPEN（10,FILE=’DXHG.DAT’）

DOI=1,N

READ（10,*）I,Y（I）,X（1,I）,X（2,I）,X（3,I）,X（4,I）

ENDDO

CLOSE（10）

MM=M+1

CALLDXHG（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）

10FORMAT（1X,'A（',I2,'）=',E13.6）

WRITE（9,15）Q,S,R

15FORMAT（1X,'Q=',E13.6,2X,'S=',E13.6,2x,'R=',E13.6）

WRITE（9,20）（I,V（I）,I=1,M）

20FORMAT（1X,'V（',I2,'）=',E13.6）

CLOSE（9）

END

计算结果为：

（1）=-.699607E+01

（2）=.755634E-01

A（3）=.871720E+00

A（4）=.204831E+02

A（5）=-.396034E-01

Q=.887608E+03S=.421333E+01R=.896120E+00

（1）=.874766E+00

（2）=.860884E+00

V（3）=.994168E+00

V（4）=.694885E-01

U=.361871E+04

3．４判别分析程序

3.4.1功能

给出判别分析的程序。

3.4.2方法说明

3.4.2.1二级判别的Fisher准则概念

3.4.2.2多级判别的Fisher准则概念

3.4.2.3判别函数的显著性检验

3.4.3子程序语句

SUBROUTINEPBFX（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）,

SUBROUTINEPBFX（X,N,M,G,KG）

M个因子,序列总长度为N,分为G类,每类的个数为NG（G）,N=NG

（1）+NG

（2）+...+NG（G）

每类的均值为XV（M,G）,总的均值为XVV（M）,

IMPLICITNONE

INTEGER:

G,I,J,K,N,M,L,MG

REAL（8）,DIMENSION（N,M）:

REAL（8）,DIMENSION（M）:

XVV

INTEGER,DIMENSION（G）:

INTEGER,DIMENSION（N）:

REAL（8）,DIMENSION（M,G）:

REAL（8）,DIMENSION（M,M）:

T,S,B

REAL（8）,DIMENSION（M）:

ST,DT,XX1

REAL（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:

REAL（8）,DIMENSION（:

）,ALLOCATABLE:

INTEGER,DIMENSION（:

）,ALLOCATABLE:

KGJ

INTEGER:

L0,NP,NP1,N0,NN

REAL（8）:

AA1,XX2,XX3,DMIN,ALF,H

NG=0

DOI=1,N

NG（KG（I））=NG（KG（I））+1

ENDDO

MG=MAXVAL（NG）!

求数组的最大值

ALLOCATE（XG（M,G,MG））

XG=0

NG=0

DOI=1,N

NG（KG（I））=NG（KG（I））+1

DOJ=1,M

XG（J,KG（I）,NG（KG（I）））=X（I,J）

ENDDO

求每个因子总的均值

XVV=0

DOK=1,M

DOI=1,N

XVV（K）=XVV（K）+X（I,K）

ENDDO

XVV（K）=XVV（K）/N

ENDDO

WRITE（12,'（4X,"总的均值=",F8.2）'）XVV

求每个因子各组的均值

XV=0

DOK=1,M

DOI=1,G

DOJ=1,NG（I）

XV（K,I）=XV（K,I）+XG（K,I,J）

ENDDO

XV（K,I）=XV（K,I）/NG（I）

ENDDO

WRITE（12,*）

DOI=1,G

WRITE（12,'（"第",I2,"组的均值=",F8.2）'）I,（XV（J,I）,J=1,M）

ENDDO

求总的离差平方和阵T和组内间离差平方和阵S

T=0

DOK=1,M

DOL=1,M

DOI=1,G

DOJ=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,I））*（XG（L,I,J）-XV（L,I））

ENDDO

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=0

DOK=1,M

DOL=1,M

DOI=1,G

B（K,L）=B（K,L）+NG（I）*（XV（K,I）-XVV（K））*（XV（L,I）-XVV（L））

ENDDO

WRITE（12,*）

WRI

展开阅读全文