多元回归fortran程序，预估值y哪里有错

jojingwj

各位大神，求指教！
我用fortran做二元回归（回归程序就是气象家园上的），回归系数是没问题的，但是在算预估值的时候，发现与实际值差别很大，个人觉得是常数项y0的问题，但是又不知道哪里有错，请大神们看一下，谢谢！
module LinearAlgebra
implicit none
contains
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!                     !!!!!!!!!
!!!!!!!! Gauss_Jordan法 !!!!!!!!!
!!!!!!!! !!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine Gauss_Jordan(A,S,ANS)
implicit none
real :: A(:,:)
real :: S(:)
real :: ANS(:)
real, allocatable :: B(:,:)
integer :: i, N
N = size(A,1)
allocate(B(N,N))
! 保存原先的矩阵A,及数组S
B=A
ANS=S
! 把B化成对角矩阵
call Upper(B,ANS,N) ! 把B化成上三角矩阵
call Lower(B,ANS,N) ! 把B化成下三角矩
! 求解
forall(i=1:N)
ANS(i)=ANS(i)/B(i,i)
end forall
return
end subroutine Gauss_Jordan
! 输出等式子程序
subroutine output(M,S)
implicit none
real :: M(:,:), S(:)
integer :: N,i,j
N = size(M,1)
! write中加上advance="no",可以中止断行发生,使下一次的
! write接续在同一行当中.
do i=1,N
!write(*,"(1x,f5.2,a1)", advance="NO") M(i,1),'A'
do j=2,N
if ( M(i,j) < 0 ) then
!write(*,"('-',f5.2,a1)",advance="NO") -M(i,j),char(64+j)
else
!write(*,"('+',f5.2,a1)",advance="NO") M(i,j),char(64+j)
end if
end do
!write(*,"('=',f8.4)") S(i)
end do
return
end subroutine output
! 求上三角矩阵的子程序
subroutine Upper(M,S,N)
implicit none
integer :: N
real :: M(N,N)
real :: S(N)
integer :: I,J
real :: E
do I=1,N-1
do J=I+1,N
E=M(J,I)/M(I,I)
M(J,I:N)=M(J,I:N)-M(I,I:N)*E
S(J)=S(J)-S(I)*E
end do
end do
return
end subroutine Upper
! 求下三角矩阵的子程序
subroutine Lower(M,S,N)
implicit none
integer :: N
real :: M(N,N)
real :: S(N)
integer :: I,J
real :: E
do I=N,2,-1
do J=I-1,1,-1
E=M(J,I)/M(I,I)
M(J,1:N)=M(J,1:N)-M(I,1:N)*E
S(J)=S(J)-S(I)*E
end do
end do
return
end subroutine Lower
end module
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!! !!!!!!!!!
!!!!!!!!! 求解齐次线性方程组 !!!!!!!!!
!!!  created by LIU Qian , 23 Sep 2013  ,QingDao
!!!  courtesy lanxizhishui
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
program multi_reg
use LinearAlgebra
implicit none
integer,parameter:: N=2 , it=51 ! (N:线性方程组矩阵大小/多元回归‘元’的个数 , it:时间序列长度)      !!!
real::b(it,n),y(144,73,it),p(n),A(n,n),S(n),ans(n) ,hcy(144,73)
real::q,r,u,y0
integer :: i,j,k,ii,jj,kk
!读入要进行多元线性回归的数据，以二元为例，建立x1、x2到y的线性回归
open(11,file='x1.grd',form='binary')
open(12,file='x2.grd',form='binary')
open(13,file='y.grd',form='binary')
open(15,file='a.grd',form='binary')
open(16,file='b.grd',form='binary')
open(17,file='reg_y.grd',form='binary')
open(18,file='reg_y_new.txt')  !写这两个txt只是为了看结果
open(19,file='reg_y_old.txt')

do i=1,it
do jj=1,73
do ii=1,144
read(13) y(ii,jj,i)
enddo;enddo;enddo
do jj=1,73
do ii=1,144
do i=1,it
write(19,*)  y(ii,jj,i)
enddo;enddo;enddo
do i=1,it
read(11)b(i,1)
read(12)b(i,2)
enddo


do jj=1,73
do ii=1,144
p=0;q=0
do i=1,it

p(1)=p(1)+b(i,1)

p(2)=p(2)+b(i,2)

q=q+y(ii,jj,i)
end do
!close(11)
!close(12)
!close(13)
!计算回归系数阵
do i=1,n
do j=1,n
r=0;u=0
do k=1,it
r=r+b(k,i)*b(k,j)
u=u+b(k,i)*y(ii,jj,k)
end do
a(i,j)=r-p(i)*p(j)/it
s(i)=u-p(i)*q/it
end do
end do
print*,s(:)
!解线性方程组
!write(*,*) 'Equation:'
call output(A,S)
call Gauss_Jordan(A,S,ANS)
y0=(q-s(1)*p(1)-s(2)*p(2))/it
!write(*,*) 'Ans:'
do i=1,N
!write(*,"(1x,a1,'=',F8.4)") char(64+i),ANS(i)  !A、B依次为第一、二个元的回归系数
end do
!write(15) ANS(1)
!write(16) ANS(2)
!write(17) y0
!write(25,*) y0
!print*,ANS(1),ANS(2),ii,jj
!!print*,y0 !此处输出的是常数项
write(15) ANS(1)
write(16) ANS(2)
do kk=1,it
write(17) ANS(1)*b(kk,1)+ANS(2)*b(kk,2)+y0   !!!!!!!!!!!!!!!!!!!!!这里的y0有大几百，加上之后差距很大，不知道是哪里的错
write(18,*) ANS(1)*b(kk,1)+ANS(2)*b(kk,2)+y0   
enddo
enddo;enddo
end program

无心恋唱

你的问题现在解决了吗？

zx1234

这个是最想要的，积分积分啊