M-K方法存疑

Jasonluo

魏老师书后边的M-K检验方法代码如下：
——————————————————————————
C        THIS IS A PROGRAM FOR DETECTING ABRUPT CLIMATIC CHANGE
C        BY USING MANN-KENDALL TECHNIQUE
        PROGRAM MK
        DIMENSION Y(1000),YY(1000),U(1000),UF(1000),UB(1000),
     &        M(1000),MD(1000)
        WRITE(*,10)
  10        FORMAT(2X,'N=?,NYEAR=?')
        READ(*,*) N,NYEAR
C        ***************************************************
C        * N:     SAMPLE SIZE                              *
C        * NYEAR: FIRST YEAR OF THE TIME SERIES            *
C        * Y(N):  ORIGINAL TIME SERIES                     *
C        * UF(N): ORIGINAL SERIES OF U(LN)                 *
C        * UB(N): COUNTER SERIES OF U(LN)                  *
C        * A,B:   CRITICAL VALUE 1.96 AND -1.96            *
C        ***************************************************          
        OPEN(2,FILE='T.txt')
        READ(2,*)(Y(I),I=1,N)
        CALL SMK(Y,M,MD,UF,N)
        DO 20 I=1,N
  20        YY(I)=Y(N+1-I)
        CALL SMK(YY,M,MD,U,N)
        DO 30 I=1,N
  30        UB(I)=-U(N+1-I)
        OPEN(3,FILE='T_0522a.txt',STATUS='NEW')
        A=1.96
        B=-1.96
        DO 40 I=1,N
        WRITE(3,50)NYEAR+I-1,UF(I),UB(I),A,B
  50        FORMAT(1X,I4,4F8.2)
  40        CONTINUE
        CLOSE(3)
        STOP
        END
C***********************************************************
        SUBROUTINE SMK(Y,M,MD,U,N)
        DIMENSION Y(N),M(N),MD(N),U(N)
        M(1)=0
        DO 10 I=2,N
        M(I)=0
        MD(I)=0
        DO 20 J=1,I-1
        IF(Y(I).LT.Y(J))GOTO 20
        M(I)=M(I)+1
  20        CONTINUE
        MD(I)=MD(I-1)+M(I)
  10        CONTINUE
        U(1)=0.0
        DO 30 I=2,N
        E=I*(I-1)/4.00
        VAR=I*(I-1)*(2*I+5)/72.00
        U(I)=(MD(I)-E)/SQRT(VAR)
  30        CONTINUE
        RETURN
        END
————————————————————————————
其中，红色字体醒目标注的这一行代码，在有的程序中为J=1,I。

计算的检验结果差别还是比较大的，见下图：


究竟哪一种对了呢？从数据原始序列判断，还是J=1，I比较靠谱一些。
这差别在哪里呢？

river

没仔细研究过，可以看一下书上的原始公式吧。还有就是这两个结果的差别就在于沿着Y轴有一个为1的位移，导致了交叉点水平的位移

言深深

按源程序走，貌似是I吧

zihan6688

用多个方法验证一下，看哪个比较好。

Jasonluo

言深深 发表于 2013-5-23 08:49

                               
登录/注册后可看大图

按源程序走，貌似是I吧

魏老师的源程序是i-1，比较了一下原序列，发现i-1应该靠谱一些。
魏老师的程序应该是对的吧。但之前有网友给魏老师的M-K方法跳出了好几处错误，导致我对这方法发生了疑惑。
小木虫上的网友说咱中国人用M-K方法错误百出，唉，作为初学者，有点无所适从啊

言深深

Jasonluo 发表于 2013-5-23 09:28

                               
登录/注册后可看大图

魏老师的源程序是i-1，比较了一下原序列，发现i-1应该靠谱一些。
魏老师的程序应该是对的吧。但之前有网 ...

mk方法本身就是有缺陷的，李建平先生在90年代的时候就提出来了，但是很多人一直沿用，所以也没有办法···突变检测方法很多的，不下于30种···用这个方法，不在于会跑程序，而应该看看这个方法的本身有何要义，关键在什么地方，是否通过改进适合于自己的数据。因此出现多种mk检测的程序也是正常的···
瞬变说一下，不妨自己阅读以下mk的计算方法，很多中文书籍都有介绍的，或者可以阅读以下mk的原始文献【Goossens C and Berger A 1986 Annals of Geophysics 4 385】

言深深

Jasonluo 发表于 2013-5-23 09:28

                               
登录/注册后可看大图

魏老师的源程序是i-1，比较了一下原序列，发现i-1应该靠谱一些。
魏老师的程序应该是对的吧。但之前有网 ...

mk方法本身就是有缺陷的，李建平先生在90年代的时候就提出来了，但是很多人一直沿用，所以也没有办法···突变检测方法很多的，不下于30种···用这个方法，不在于会跑程序，而应该看看这个方法的本身有何要义，关键在什么地方，是否通过改进适合于自己的数据。因此出现多种mk检测的程序也是正常的···
瞬变说一下，不妨自己阅读以下mk的计算方法，很多中文书籍都有介绍的，或者可以阅读以下mk的原始文献【Goossens C and Berger A 1986 Annals of Geophysics 4 385】

Jasonluo

言深深 发表于 2013-5-23 12:26

                               
登录/注册后可看大图

mk方法本身就是有缺陷的，李建平先生在90年代的时候就提出来了，但是很多人一直沿用，所以也没有办法·· ...

谢谢深深，你说的这篇文献，能给出题目吗？我在google学术里没搜到。如果有原文的话，方便共享一下吗？

言深深

Jasonluo 发表于 2013-5-23 14:43

                               
登录/注册后可看大图

谢谢深深，你说的这篇文献，能给出题目吗？我在google学术里没搜到。如果有原文的话，方便共享一下吗？

抱歉，我也没有这篇文章，年代比较久远，该杂志没有收录到那边远的。不过帮你找到一个详细的文献结果：Goossens, C. and A. Berger, 1986. Annual and seasonal climatic variations of northern hemisphere and europe over the last century. Ann. Geo-Phys., 4: 385-400.

或许可以在史料馆里面能找到

星图

很好，学习了