求高手、大虾赐教
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
上面是段多元回归的算法，大家帮忙读下或给点读的思路，本人还是新手，难以入手，正急着查资料解决，路过的请发表点意见，先谢谢了
