【多元线性回归方程】算法求解~救命呀呀呀呀哎呀~~

哈哈哈

我在写多元线性回归方程~可是他怎么也算不对…哪位大神帮帮忙来帮小女子看看呀
我用的是高斯消元法求的逆矩阵

[Fortran] 纯文本查看 复制代码

查看代码

打印?

01	program multiple_linear_regression

02	implicit none

03	integer :: n,num, i,j

04	real, allocatable :: x(:,:),y(:),xtrans(:,:),xx(:,:),aa(:,:), b(:)

05	write(*,*) "请输入自变量的个数："

06	read(*,*) n

07	write(*,*) "请输入次数："

08	read(*,*) num

09	allocate( x(num,n+1) )

10	allocate( xtrans(n+1,num) )

11	allocate( xx(num,num) )

12	allocate( aa(num,num) )!aa是用来存放求逆后的矩阵

13	allocate( y(num) )

14	allocate( b(n+1) )

15

x(:,1)=1

16	write(*,*) "请输入自变量X矩阵："

17	do i=1, num

18	read(*,*) (x(i,j),j=2,n+1)

19

end do

20	write(*,*) "请输入因变量y的值："

21	read(*,*) y(:)

22

!转置矩阵

23	xtrans=transpose(x)

24

!矩阵相乘

25	xx=matmul(x,xtrans)

26	!（X'X）矩阵求逆

27	call nijuzhen(xx,num,aa)

28	!求B=（X'X）的逆矩阵 * X'* Y

29	b=matmul( matmul(aa,xtrans),y )

30	write(*,*) b

31

32

end

33

34

35	subroutine nijuzhen(a,n,aa)!高斯消去法求逆矩阵

36	implicit none

37	integer :: i,k,j,n

38	real :: ae(n,n+n),a(n,n),aa(n,n)!a矩阵右侧添加一个单位矩阵

39	ae(:,:)=0!!!!

40

do i=1,n

41	ae(i,1:n)=a(i,:)

42	ae(i,n+i)=1

43

end do

44	!ae化为上三角矩阵

45

do k=2,n

46

do i=k,n

47	ae(i,:)=ae(i,:)-ae(k-1,:)*ae(i,k-1)/ae(k-1,k-1)

48

end do

49

end do

50	!将第i行第i列元素化成1

51

do i=1,n

52	a(i,:)=a(i,:)/a(i,i)

53

end do

54	!矩阵变换将左侧化为单位矩阵

55

do j=2,n

56	do i=1,j-1

57	if( a(i,j)/=0 ) then

58	a(i,:)=a(i,:)-a(i,j)*a(j,:)

59

end if

60

end do

61

end do

62	!矩阵右侧即得到逆矩阵

63

do i=1,n

64	aa(i,:)=ae(i,(n+1):(n+n))

65

end do

66

end

kerb

我觉得你可能是想表达：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