求读程序高手，专家赐教

[jinwangwan]

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

[jinwangwan]

真的很急，懂的朋友给点意见，哪怕一点点的提示，顺便问下这是什么语言

[yiyaoyao58958]

FORTRAN语言。。。。。

[yiyaoyao58958]

vb形式如下（几个if语句的条件可能不对~）(不知道你会不会vb。。。。但是因为是vb社区，所以转成了vb方式。。。)


不知道楼主要这个方法是用来学习还是泡女朋友.....

FORTRAN语言是用来专门解决计算问题的语言,你是研究数学的?

DXHG(X,Y,M,N,A,Q,S,R,V,U,B) 中的一个参数就是一个变量


REAL(KIND=8),DIMENSION(M,N)::X 是你们特有的,不明白!

但是你可以从观察每个变量作了什么来判断他的原方程,或者说原算式的形态~

如果你已经知道原算式,只是看不懂程序的话,建议你从CHOLESKY(...)开始看,看它作了什么,能对应到算式的那个小部分~

如果不知道算式的话,我也帮不了你...没研究过纯算式问题....不好意思..

Dim i As Integer = 0
	Dim j As Integer = 0
	Sub DXHG(ByVal M, ByVal N, ByVal X, ByVal Y, ByVal A, ByVal Q, ByVal S, ByVal R, ByVal V, ByVal U, ByVal 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
		For j = 2 To MM
			B(1, J) = 0

			For i = 1 To N
				B(1, J) = B(1, J) + X(J - 1, i)
			Next

			B(J, 1) = B(1, J)
		Next

		For i = 2 To MM
			For J = i To MM
				B(i, J) = 0
				For K = 1 To N
					B(i, J) = B(i, J) + X(i - 1, K) * X(J - 1, K)
				Next
				B(J, i) = B(i, J)
			Next
		Next

		A(1) = 0

		For i = 1 To N
			A(1) = A(1) + Y(i)
		Next

		For i = 2 To MM
			A(i) = 0
			For J = 1 To N
				A(i) = A(i) + X(i - 1, J) * Y(J)
			Next
		Next

		Call CHOLESKY(B, MM, 1, A, L)

		YY = 0

		For i = 1 To N
			YY = YY + Y(i)
		Next

		YY = YY / N
		Q = 0
		DYY = 0
		U = 0

		For i = 1 To N
			P = A(1)
			For J = 1 To M
				P = P + A(J + 1) * X(J, i)
			Next
			Q = Q + (Y(i) - P) * (Y(i) - P)
			DYY = DYY + (Y(i) - YY) * (Y(i) - YY)
			U = U + (YY - P) * (YY - P)
		Next

		S = SQRT(Q / N)
		R = SQRT(1 - Q / DYY)

		For J = 1 To M
			P = 0
			For i = 1 To N
				PP = A(1)
				For K = 1 To M
					If K / j Then
						PP = PP + A(K + 1) * X(K, i)
					End If
				Next
				P = P + (Y(i) - PP) * (Y(i) - PP)
			Next
			V(j) = SQRT(1 - Q / P)
		Next
	End Sub

	Sub CHOLESKY(ByVal C, ByVal N, ByVal M, ByVal D, ByVal L)
		'REAL(KIND=8),DIMENSION(N,N)::C  
		'REAL(KIND=8),DIMENSION(N,M)::D  
		L = 1

		If ABS(C(1, 1)) < 0.0000000001 Then
			L = 0
			MsgBox(" FAIL")
			Return
		End If

		C(1, 1) = SQRT(C(1, 1))
		For J = 2 To N
			C(1, J) = C(1, J) / C(1, 1)
		Next
		For i = 2 To N
			For J = 2 To i
				C(i, i) = C(i, i) - C(J - 1, i) * C(J - 1, i)
			Next
			If ABS(C(i, i)) < 0.0000000001 Then
				L = 0
				MsgBox(" FAIL")
				Return
			End If

			C(i, i) = SQRT(C(i, i))

			If i / N Then
				For j = i + 1 To N
					For K = 2 To i
						C(i, j) = C(i, j) - C(K - 1, i) * C(K - 1, j)
					Next
					C(i, j) = C(i, j) / C(i, i)
				Next
			End If
		Next

		For J = 1 To M
			D(1, J) = D(1, J) / C(1, 1)
			For i = 2 To N
				For K = 2 To i
					D(i, J) = D(i, J) - C(K - 1, i) * D(K - 1, J)
				Next
				D(i, J) = D(i, J) / C(i, i)
			Next
		Next

		For J = 1 To M
			D(N, J) = D(N, J) / C(N, N)
			For K = N To 2 Step -1
				For i = K To N
					D(K - 1, J) = D(K - 1, J) - C(K - 1, i) * D(i, J)
				Next
				D(K - 1, J) = D(K - 1, J) / C(K - 1, K - 1)
			Next
		Next

	End Sub

[yiyaoyao58958]

啊~~写错了~~


FORTRAN语言的“/=”是vb的“<>”。。。。。


这次的是等价的vb程序

Dim i As Integer = 0
	Dim j As Integer = 0
	Sub DXHG(ByVal M, ByVal N, ByVal X, ByVal Y, ByVal A, ByVal Q, ByVal S, ByVal R, ByVal V, ByVal U, ByVal 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
		For j = 2 To MM
			B(1, J) = 0

			For i = 1 To N
				B(1, J) = B(1, J) + X(J - 1, i)
			Next

			B(J, 1) = B(1, J)
		Next

		For i = 2 To MM
			For J = i To MM
				B(i, J) = 0
				For K = 1 To N
					B(i, J) = B(i, J) + X(i - 1, K) * X(J - 1, K)
				Next
				B(J, i) = B(i, J)
			Next
		Next

		A(1) = 0

		For i = 1 To N
			A(1) = A(1) + Y(i)
		Next

		For i = 2 To MM
			A(i) = 0
			For J = 1 To N
				A(i) = A(i) + X(i - 1, J) * Y(J)
			Next
		Next

		Call CHOLESKY(B, MM, 1, A, L)

		YY = 0

		For i = 1 To N
			YY = YY + Y(i)
		Next

		YY = YY / N
		Q = 0
		DYY = 0
		U = 0

		For i = 1 To N
			P = A(1)
			For J = 1 To M
				P = P + A(J + 1) * X(J, i)
			Next
			Q = Q + (Y(i) - P) * (Y(i) - P)
			DYY = DYY + (Y(i) - YY) * (Y(i) - YY)
			U = U + (YY - P) * (YY - P)
		Next

		S = SQRT(Q / N)
		R = SQRT(1 - Q / DYY)

		For J = 1 To M
			P = 0
			For i = 1 To N
				PP = A(1)
				For K = 1 To M
					If K <> j Then
						PP = PP + A(K + 1) * X(K, i)
					End If
				Next
				P = P + (Y(i) - PP) * (Y(i) - PP)
			Next
			V(j) = SQRT(1 - Q / P)
		Next
	End Sub

	Sub CHOLESKY(ByVal C, ByVal N, ByVal M, ByVal D, ByVal L)
		REAL(KIND=8),DIMENSION(N,N)::C  
		REAL(KIND=8),DIMENSION(N,M)::D  
		L = 1

		If ABS(C(1, 1)) < 0.0000000001 Then
			L = 0
			MsgBox(" FAIL")
			Return
		End If

		C(1, 1) = SQRT(C(1, 1))
		For J = 2 To N
			C(1, J) = C(1, J) / C(1, 1)
		Next
		For i = 2 To N
			For J = 2 To i
				C(i, i) = C(i, i) - C(J - 1, i) * C(J - 1, i)
			Next
			If ABS(C(i, i)) < 0.0000000001 Then
				L = 0
				MsgBox(" FAIL")
				Return
			End If

			C(i, i) = SQRT(C(i, i))

			If i <> N Then
				For j = i + 1 To N
					For K = 2 To i
						C(i, j) = C(i, j) - C(K - 1, i) * C(K - 1, j)
					Next
					C(i, j) = C(i, j) / C(i, i)
				Next
			End If
		Next

		For J = 1 To M
			D(1, J) = D(1, J) / C(1, 1)
			For i = 2 To N
				For K = 2 To i
					D(i, J) = D(i, J) - C(K - 1, i) * D(K - 1, J)
				Next
				D(i, J) = D(i, J) / C(i, i)
			Next
		Next

		For J = 1 To M
			D(N, J) = D(N, J) / C(N, N)
			For K = N To 2 Step -1
				For i = K To N
					D(K - 1, J) = D(K - 1, J) - C(K - 1, i) * D(i, J)
				Next
				D(K - 1, J) = D(K - 1, J) / C(K - 1, K - 1)
			Next
		Next

	End Sub

[khan45]

楼上太吊了。