第三章多元统计分析Word文档格式.docx
- 文档编号:17637924
- 上传时间:2022-12-07
- 格式:DOCX
- 页数:37
- 大小:28.69KB
第三章多元统计分析Word文档格式.docx
《第三章多元统计分析Word文档格式.docx》由会员分享,可在线阅读,更多相关《第三章多元统计分析Word文档格式.docx(37页珍藏版)》请在冰豆网上搜索。
CLOSE(10)
CALLFANGCHA(X,N,XBAR,S,S2)
OPEN(12,FILE='
FANGCHA.DAT'
WRITE(12,'
(4X,"
XBAR"
7X,"
S"
9X,"
S2"
)'
(3F10.4)'
)XBAR,S,S2
CLOSE(12)
计算结果如下:
XBAR=22.5718S=4.2888S2=18.3937
4.3回归分析程序
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)
X,Y
XV=0!
x的均值
YV=0!
y的均值
XV=XV+X(I)
YV=YV+Y(I)
XV=XV/N
YV=YV/N
DXX=0.0
DXY=0.0
Q=X(I)-XV
DXX=DXX+Q*Q
DXY=DXY+Q*(Y(I)-YV)
A=DXY/DXX
B=YV-A*XV
Q=0
U=0
P=0
UMAX=-1.0E10
UMIN=1.0E20
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
U=U+DX/N
S=SQRT(Q/N)
3.3.1.6例
下面以某海区22年的1月平均气温作x,对应的海表温度作y,编写计算的主程序,并给出计算结果。
PROGRAMmain1
N=22
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='
S='
E13.6,2x,'
P='
WRITE(9,20)UMAX,UMIN,U
20FORMAT(1X,'
UMAX='
UMIN='
U='
CLOSE(9)
计算结果:
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——实型一维数组,长度为N,输入参数。
存放y的N个观测值。
M——整型变量,输入参数。
自变量个数。
观测数据的组数。
A——实型一维数组,长度为M+1,输出参数。
存放明M+1个回归系数a1,a2,……am+1。
R——实型变量,输出参数。
复相关系数。
V——实型一维数组,长度为M,输出参数。
M个自变量的偏相关系数。
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):
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
DOJ=2,MM
B(1,J)=0
B(1,J)=B(1,J)+X(J-1,I)
B(J,1)=B(1,J)
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)
A
(1)=0
A
(1)=A
(1)+Y(I)
A(I)=0
DOJ=1,N
A(I)=A(I)+X(I-1,J)*Y(J)
CALLCHOLESKY(B,MM,1,A,L)
YY=0
YY=YY+Y(I)
YY=YY/N
DYY=0
P=A
(1)
DOJ=1,M
P=P+A(J+1)*X(J,I)
Q=Q+(Y(I)-P)*(Y(I)-P)
DYY=DYY+(Y(I)-YY)*(Y(I)-YY)
U=U+(YY-P)*(YY-P)
R=SQRT(1-Q/DYY)
PP=A
(1)
DOK=1,M
IF(K/=J)PP=PP+A(K+1)*X(K,I)
P=P+(Y(I)-PP)*(Y(I)-PP)
V(J)=SQRT(1-Q/P)
SUBROUTINECHOLESKY(C,N,M,D,L)
REAL(KIND=8),DIMENSION(N,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
ENDIF
C(1,1)=SQRT(C(1,1))
DOJ=2,N
C(1,J)=C(1,J)/C(1,1)
DOI=2,N
DOJ=2,I
C(I,I)=C(I,I)-C(J-1,I)*C(J-1,I)
IF(ABS(C(I,I))<
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)
C(I,J)=C(I,J)/C(I,I)
D(1,J)=D(1,J)/C(1,1)
DOK=2,I
D(I,J)=D(I,J)-C(K-1,I)*D(K-1,J)
D(I,J)=D(I,J)/C(I,I)
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)
D(K-1,J)=D(K-1,J)/C(K-1,K-1)
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
N=50
M=4
REAL(KIND=8)Q,S,R,U
OPEN(10,FILE=’DXHG.DAT’)
READ(10,*)I,Y(I),X(1,I),X(2,I),X(3,I),X(4,I)
CALLDXHG(M,N,X,Y,A,Q,S,R,V,U,B)
RESULT2.DAT'
WRITE(9,10)(I,A(I),I=1,MM)
A('
I2,'
)='
WRITE(9,15)Q,S,R
R='
WRITE(9,20)(I,V(I),I=1,M)
V('
计算结果为:
A
(1)=-.699607E+01
A
(2)=.755634E-01
A(3)=.871720E+00
A(4)=.204831E+02
A(5)=-.396034E-01
Q=.887608E+03S=.421333E+01R=.896120E+00
V
(1)=.874766E+00
V
(2)=.860884E+00
V(3)=.994168E+00
V(4)=.694885E-01
U=.361871E+04
3.4判别分析程序
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)
G,I,J,K,N,M,L,MG
REAL(8),DIMENSION(N,M):
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
ST,DT,XX1
ST12,ST_12,SS_12,SS,V,S12,S_12,VD,D,VS
REAL(8),DIMENSION(:
),ALLOCATABLE:
AA,X2,X22
:
XG
DS
INTEGER,DIMENSION(:
KGJ
L0,NP,NP1,N0,NN
REAL(8):
AA1,XX2,XX3,DMIN,ALF,H
NG=0
NG(KG(I))=NG(KG(I))+1
MG=MAXVAL(NG)!
求数组的最大值
ALLOCATE(XG(M,G,MG))
XG=0
XG(J,KG(I),NG(KG(I)))=X(I,J)
求每个因子总的均值
XVV=0
XVV(K)=XVV(K)+X(I,K)
XVV(K)=XVV(K)/N
总的均值="
<
M>
F8.2)'
)XVV
求每个因子各组的均值
XV=0
DOI=1,G
DOJ=1,NG(I)
XV(K,I)=XV(K,I)+XG(K,I,J)
XV(K,I)=XV(K,I)/NG(I)
WRITE(12,*)
第"
I2,"
组的均值="
)I,(XV(J,I),J=1,M)
求总的离差平方和阵T和组内间离差平方和阵S
T=0
DOL=1,M
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))
WRITE(12,'
总的离差阵T"
(<
F9.2)'
)T
组内的离差阵S"
)S
求组间离差平方和阵B,可以采用两种方法:
1.B=T-S;
2.直接计算法。
两种方法计算结果完全一致。
第一种方法:
B=T-S
B=T-S
第二种方法:
直接计算法
B=0
B(K,L)=B(K,L)+NG(I)*(XV(K,I)-XVV(K))*(XV(L,I)-XVV(L))
WRI
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 第三 多元 统计分析