求教大家缺测值计算程序里怎么的处理？

ARS09

在求两向量相关系数的程序基础上，编的求向量x(n)和三维数组y(l,m,n)在第三维方向上向量的二者的相关系数，但是不知道缺测值部分怎么写，下面的程序出去缺测处理那部分其它计算结果是对的，x（n），y是输了几个数测试。求教大家处理缺测该怎么写？    program mean
    integer,parameter::l=3,m=4,n=5
    real,parameter::undef=-9.99e33
    integer p,q,i,num
    real ax,vx,sx
    real y(l,m,n),ay(l,m),vy(l,m),sy(l,m),sxy(l,m),r(l,m)
    real::x(n)=(/1,2,3,4,5/),xw(n)
    open(12,file='st.txt')
    read(12,*)(((y(p,q,i),q=1,4),p=1,3),i=1,5)
    close(12)

    !处理缺测
    num=0
    do 10 p=1,l
        do q=1,m
            do i=1,n
                if(x(i).eq.undef.or.y(p,q,i).eq.undef)goto 10
                num=num+1
                xw(num)=x(i)
                y(p,q,num)=y(p,q,i)
            end do
        end do
10  continue
    if(num.le.5)then
        do p=1,l
            do q=1,m
                r(p,q)=undef
            end do
        end do
    else


!计算meanvar  y(l,m,n)
    do p=1,l
        do q=1,m
            ay(p,q)=0.0
            vy(p,q)=0.0
            sy(p,q)=0.0
        end do
    end do
    !求和
    do q=1,m
        do p=1,l
            do i=1,n
                ay(p,q)=ay(p,q)+y(p,q,i)      
            end do
        end do
    end do
    do p=1,l
        do q=1,m
            ay(p,q)=ay(p,q)/float(n)         
        end do
    end do

    !求平均
    do p=1,l
        do q=1,m
            do i=1,n
                vy(p,q)=vy(p,q)+(y(p,q,i)-ay(p,q))**2
            end do
        end do
    end do

    do p=1,l
        do q=1,m
            vy(p,q)=vy(p,q)/float(n)
            sy(p,q)=sqrt(vy(p,q))
        end do
    end do
! 计算meanvar  x(n)
    call meanvar(n,x,ax,sx,vx)

!求x(n)和y(l,m,n)在N方向上相关系数
    do p=1,l
        do q=1,m
            sxy(p,q)=0.0
        end do
    end do
    do p=1,l
        do q=1,m
            do i=1,n
                sxy(p,q)=sxy(p,q)+(x(i)-ax)*(y(p,q,i)-ay(p,q))
            end do
        end do
    end do
    do p=1,l
        do q=1,m
            sxy(p,q)=sxy(p,q)/float(n)
        end do
    end do
    do p=1,l
        do q=1,m
            r(p,q)=sxy(p,q)/(sx*sy(p,q))
        end do
    end do
    write(*,*)((r(p,q),p=1,l),q=1,m)

    end if

    end




    subroutine meanvar(n,x,ax,sx,vx)
    dimension x(n)
    ax=0.0
    vx=0.0
    sx=0.0
    do 10 i=1,n
        ax=ax+x(i)
10  continue
    ax=ax/float(n)
    do 20 i=1,n
        vx=vx+(x(i)-ax)**2
20  continue
    vx=vx/float(n)
    sx=sqrt(vx)
    return
    end

lqouc

把每一个点上两个值都不为缺省的提取为新序列，然后对这个新序列求相关。
以上为个人看法，没实现过，也不知道别的带缺测的程序是怎么算的，你可以多看看别人的算法。

ARS09

lqouc 发表于 2014-11-2 17:43
把每一个点上两个值都不为缺省的提取为新序列，然后对这个新序列求相关。
以上为个人看法，没实现过，也不 ...

谢谢版主！我是想让带有缺测值得序列计算的出的相关系数也为缺测值，但是写的程序运行有问题，不知怎么写。

lqouc

ARS09 发表于 2014-11-2 18:16
谢谢版主！我是想让带有缺测值得序列计算的出的相关系数也为缺测值，但是写的程序运行有问题，不知怎么写 ...

序列中只要有缺省值就给结果赋缺省么？这个很简单的啊，用一个数组函数any就可以判断了。
至于自己写的程序有什么问题说的具体一点，报错信息贴上来，光看代码我也看不出问题。
顺便说一下，上面那个程序实在是太麻烦了，多用数组函数你的程序能缩减十几行。

ARS09

lqouc 发表于 2014-11-2 18:34
序列中只要有缺省值就给结果赋缺省么？这个很简单的啊，用一个数组函数any就可以判断了。
至于自己写的 ...

    num=0
    do 10 p=1,l
        do q=1,m
            do i=1,n
                if(x(i).eq.undef.or.y(p,q,i).eq.undef)goto 10
                num=num+1
                xw(num)=x(i)
                y(p,q,num)=y(p,q,i)
            end do
        end do
10  continue
    if(num.le.5)then
        do p=1,l
            do q=1,m
                r(p,q)=undef
            end do
        end do
    else就是这部分选择判断，判断完了把r赋值成缺省值，后面计算没问题，这部分运行后就中断了，虽然没有语法错误。

ARS09

这部分写的有问题，就是不知道怎么判断完，再赋值为缺测后，若是没缺测再执行后面的计算过程。

lqouc

不明白你想怎么去检测缺测，把你对缺测判断的具体要求说一下，我看看逻辑有没有问题，至于程序怎么不对，还是那句话，把中断时候的报错给出来，如果是秒退，自己先确定一下断点位置。

ARS09

lqouc 发表于 2014-11-2 19:34
不明白你想怎么去检测缺测，把你对缺测判断的具体要求说一下，我看看逻辑有没有问题，至于程序怎么不对，还 ...

恩，我用的是一个点的月时间序列，和海温每年5月月序列，因为陆地都是用缺测值表示的，所以就想把陆地那部分参与计算的缺测值不参与计算或者把结果直接设为某个值，方便后面作图。
运行后出现的结果：C:\Users\Administrator\Desktop\QQ截图20141102193856.png

ARS09

ARS09 发表于 2014-11-2 19:45
恩，我用的是一个点的月时间序列，和海温每年5月月序列，因为陆地都是用缺测值表示的，所以就想把陆地那 ...

没有逻辑判断那部分，程序运行结果是正确的。

lqouc

大致看了下你的程序，可以理解你的l、m是海温的空间，然后n是时间，你那么写的话前面不管怎么改都不可能输出正确结果，因为后面的相关性计算是对全部空间场，这样前面被你赋的缺省值又被覆盖了。
建议你先对所有点计算相关，然后再判断是否为陆地赋缺省值。也就是把你的程序倒过来用。这里需要注意的是关于缺省值的判断，我觉得你goto跳出的位置不对，跳出第一层即可，这个你自己考虑下吧。
ps：你的图没有显示，要先把图传到附件中再贴过来，提问前一定要阅读‘提问的智慧’。