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子程序语句 
SUBROUTINE DXHG(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） 
SUBROUTINE DXHG(M,N,X,Y,A,Q,S,R,V,U,B) 
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,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 
   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) 
 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,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)THEN 
   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-1,I)*D(I,J) 
END DO 
D(K-1,J)=D(K-1,J)/C(K-1,K-1) 
 END DO 
END DO 
END
[color=00FF00]上面是段多元回归的算法，大家帮忙读下或给点读的思路，本人还是新手，难以入手，正急着查资料解决，路过的请发表点意见，先谢谢了[/color]