【已解决】相关系数求出来怎么是错的，各位帮帮忙啊

穿云搅浪

大家帮帮忙啊，检查了好长时间，不知道哪里出问题了，各位好心人帮帮我吧！！


program main
parameter(nx=30,ny=22)
external corr
real slpyear(nx,ny,23),indexncep(23),corrslp(nx,ny)
!读入要素场
open(100,file='E:/project/ncep/slpyear.grd',form='binary')
      read(100) (((slpyear(i,j,n),i=1,nx),j=1,ny),n=1,23)
close(100)
!读入季风指数
open(101,file='E:/project/ncep/indexncep.txt')   
         read(101,*) (indexncep(n),n=1,23)
close(101)
!得到相关系数
do j=1,ny
    do i=1,nx      
              corrslp(i,j)=corr(indexncep(:),slpyear(i,j,:),23)         
    end do
end do

!输出指数
open(200,file='E:/project/ncep/corrslp.grd',form='binary')   
        write(200) ((corrslp(i,j),i=1,nx),j=1,ny)
close(200)


end program

!计算相关系数函数
function corr(x,y,m)
    implicit none
    integer i,m
        real x(m),y(m),corr,sum1,sum2,sum3,u,d        
        sum1=0
        sum2=0
        sum3=0
        do i=1,m
                sum1=sum1+x(i)**2
                sum2=sum2+y(i)**2
                sum3=sum3+x(i)*y(i)
        end do
        u=sum3-m*(sum(x)/m)*(sum(y)/m)
        d=sqrt(abs((sum1-(sum(x)**2)/m)*(sum2-(sum(y)**2)/m)))
        corr=u/d
end function corr

等待中幸福

这问题确实很怪呀，程序好像没看出有什么错误呀？

river

没太看明白相关系数用的哪个公式，还是试试一步一步来计算的吧
subroutine corra(x,y,m,cc)
        real x(m),y(m),cc
        real up,down1,down2,down,xx,yy
        integer m
        xx=0
        yy=0
        up=0
        down1=0
        down2=0       
        do k=1,m
        xx=xx+x(k)/m
        yy=yy+y(k)/m
        enddo
        do k=1,m
        up=up+(x(k)-xx)*(y(k)-yy)
        down1=down1+(x(k)-xx)**2
        down2=down2+(y(k)-yy)**2
        enddo
      
        down=sqrt(down1*down2)
        cc=up/down

穿云搅浪

river 发表于 2014-4-18 09:10
没太看明白相关系数用的哪个公式，还是试试一步一步来计算的吧
subroutine corra(x,y,m,cc)
        real x(m),y ...

非常感谢啊，试了一下跟excel算的一样！应该是我以前公式搞错了

river

穿云搅浪 发表于 2014-4-18 13:00
非常感谢啊，试了一下跟excel算的一样！应该是我以前公式搞错了

一样了就好，本来还想学一下你这个简单的公式呢O(∩_∩)O哈哈~

lqouc

楼主解决了没有，解决了记得修改下标题。

来福

博主写的太精炼了，顶顶！