【多元线性回归方程】算法求解~救命呀呀呀呀哎呀~~
我在写多元线性回归方程~可是他怎么也算不对…哪位大神帮帮忙来帮小女子看看呀{:4_117:}我用的是高斯消元法求的逆矩阵
program multiple_linear_regression
implicit none
integer :: n,num, i,j
real, allocatable :: x(:,:),y(:),xtrans(:,:),xx(:,:),aa(:,:), b(:)
write(*,*) "请输入自变量的个数:"
read(*,*) n
write(*,*) "请输入次数:"
read(*,*) num
allocate( x(num,n+1) )
allocate( xtrans(n+1,num) )
allocate( xx(num,num) )
allocate( aa(num,num) )!aa是用来存放求逆后的矩阵
allocate( y(num) )
allocate( b(n+1) )
x(:,1)=1
write(*,*) "请输入自变量X矩阵:"
do i=1, num
read(*,*) (x(i,j),j=2,n+1)
end do
write(*,*) "请输入因变量y的值:"
read(*,*) y(:)
!转置矩阵
xtrans=transpose(x)
!矩阵相乘
xx=matmul(x,xtrans)
!(X'X)矩阵求逆
call nijuzhen(xx,num,aa)
!求B=(X'X)的逆矩阵 * X'* Y
b=matmul( matmul(aa,xtrans),y )
write(*,*) b
end
subroutine nijuzhen(a,n,aa)!高斯消去法求逆矩阵
implicit none
integer :: i,k,j,n
real :: ae(n,n+n),a(n,n),aa(n,n)!a矩阵右侧添加一个单位矩阵
ae(:,:)=0!!!!
do i=1,n
ae(i,1:n)=a(i,:)
ae(i,n+i)=1
end do
!ae化为上三角矩阵
do k=2,n
do i=k,n
ae(i,:)=ae(i,:)-ae(k-1,:)*ae(i,k-1)/ae(k-1,k-1)
end do
end do
!将第i行第i列元素化成1
do i=1,n
a(i,:)=a(i,:)/a(i,i)
end do
!矩阵变换将左侧化为单位矩阵
do j=2,n
do i=1,j-1
if( a(i,j)/=0 ) then
a(i,:)=a(i,:)-a(i,j)*a(j,:)
end if
end do
end do
!矩阵右侧即得到逆矩阵
do i=1,n
aa(i,:)=ae(i,(n+1):(n+n))
end do
end
我觉得你可能是想表达:AX=B,然后用最小二乘:A^TAX=A^TB,然后求出X=(A^TA)^{-1}(A^TB),但是要记住最小二乘法,方程多余自变量才有解,方程少于自变量,解是不确定的 kerb 发表于 2016-9-21 21:38
我觉得你可能是想表达:AX=B,然后用最小二乘:A^TAX=A^TB,然后求出X=(A^TA)^{-1}(A^TB),但是要记住最小 ...
我在测试的时候确实是自变量n小于方程数num{:4_117:}
页:
[1]