爱气象,爱气象家园! 

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

新浪微博登陆

只需一步, 快速开始

搜索
查看: 3429|回复: 2

[求助] 多元线性回归疑问

[复制链接]

新浪微博达人勋

发表于 2017-2-17 10:28:19 | 显示全部楼层 |阅读模式

登录后查看更多精彩内容~

您需要 登录 才可以下载或查看,没有帐号?立即注册 新浪微博登陆

x
根据帖子里面的方法,只是改变了我的时间序列n,但是就显示错误提示array  bounds exceeded ,不明白错在哪里,输入文件格式与测试数据格式一样
只是时间序列长度不一样,请大神相助http://bbs.06climate.com/forum.php?mod=viewthread&tid=26772&extra=&page=1
program multi_reg
implicit none
integer,parameter:: m=3,mm=4,n=14 ! (m:线性方程组矩阵大小/多元回归‘元’的个数 , n:时间序列长度,mm为m+1)
real::X(n,m),y(n),A(mm),B(mm,mm),BA(mm,mm+1)
integer :: i,j,k

!读入要进行多元线性回归的数据,建立y=a(1)+a(2)*x1+a(3)*x2+.....+a(m+1)*xm的线性回归
open(11,file='1.txt')
open(12,file='y.txt')

read(11,*) ((X(i,j),j=1,m),i=1,n)
print*,X(i,j)
read(12,*) (y(i),i=1,n)
b(1,1)=n
!!!!!计算回归系数阵!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
do j=2,mm
  b(1,j)=0.0
  do i=1,n
    b(1,j)=b(1,j)+x(i,j-1)
  enddo
  b(j,1)=b(1,j)
enddo
do i=2,mm
  do j=i,mm
     b(i,j)=0.0
  do k=1,n
     b(i,j)=b(i,j)+x(k,i-1)*x(k,j-1)
  enddo
  b(j,i)=b(i,j)
  enddo
enddo
a(1)=0.0
do i=1,n
  a(1)=a(1)+y(i)
enddo
do i=2,mm
  a(i)=0.0
  do j=1,n
    a(i)=a(i)+x(j,i-1)*y(j)
  enddo
enddo
ba(:,1:mm)=b
ba(:,mm+1)=a
!WRITE(*,*) (a(i),i=1,mm)
!write(*,100) ((ba(i,j),j=1,mm+1),i=1,mm)
!100 format(1x,5f10.3)
close(11)
close(12)
!!!!-----------------------------!!!!!!!!!!
!!!!!!调用子程序解方程组!!!!!!!
call xiaoyuan(ba,a,mm)
!!!!!!!!---------------------------!!!!!!!!!
!!!!输出回归系数a1,a2,a3,a4
write(*,"(1x,'a(',i1,')=',f10.6,2x)") (i,a(i),i=1,mm)
end program

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!! xiaoyuan法 !!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine xiaoyuan(ba,key,mm)
implicit none
integer :: i,j,k,mm
real ba(mm,mm+1),temp(mm+1),tmp,d(mm),key(mm)
do i=1,mm-1
if(ba(i,i)==0.0) then
  do j=i+1,mm
   if(ba(j,i)/=0.0) then
    temp=ba(i,1:mm+1)
    ba(i,1:mm+1)=ba(j,1:mm+1)
    ba(j,1:mm+1)=temp
   endif
  enddo
endif

do j=i+1,mm
  tmp=ba(j,i)/ba(i,i)
  ba(j,i:mm+1)=ba(j,i:mm+1)-tmp*ba(i,i:mm+1)
enddo
enddo
d=0
key(mm)=ba(mm,mm+1)/ba(mm,mm)
!print *, key(mm)
do i=mm-1,1,-1
do j=mm,i+1,-1
  d(i)=d(i)+key(j)*ba(i,j)
enddo
key(i)=(ba(i,mm+1)-d(i))/ba(i,i)
enddo
end
  

密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2017-2-17 11:10:47 | 显示全部楼层
这个好,不用下载直接用
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2017-4-11 19:33:47 | 显示全部楼层
不知道楼主解决了没有。用你贴的程序试了我自己的数据,没有报错。。。所以,程序本身没问题,楼主不妨检查自己的数据读入使用是否正确
密码修改失败请联系微信:mofangbao
您需要登录后才可以回帖 登录 | 立即注册 新浪微博登陆

本版积分规则

Copyright ©2011-2014 bbs.06climate.com All Rights Reserved.  Powered by Discuz! (京ICP-10201084)

本站信息均由会员发表,不代表气象家园立场,禁止在本站发表与国家法律相抵触言论

快速回复 返回顶部 返回列表