问题很严重！MK检验代码的正确性

胆小车迷

1. 
   
2. <P>        SUBROUTINE SMK(Y,M,MD,U,N)
   
3.         DIMENSION Y(N),M(N),MD(N),U(N)
   
4.         M(1)=0
   
5.         DO 10 I=2,N
   
6.         M(I)=0
   
7.         MD(I)=0
   
8.         DO 20 J=1,I-1
   
9.         IF(Y(I).LT.Y(J))GOTO 20
   
10.         M(I)=M(I)+1
    
11.   20   CONTINUE</P>
    
12. <P>        MD(I)=MD(I-1)+M(I)
    
13. 10   CONTINUE
    
14.         U(1)=0.0
    
15.         DO 30 I=2,N
    
16.         E=I*(I-1)/4.00
    
17.         VAR=I*(I-1)*(2*I+5)/72.00
    
18.         U(I)=(MD(I)-E)/SQRT(VAR)
    
19.   30        CONTINUE
    
20.         RETURN
    
21.         END
    
22. </P>

复制代码

做过MK检验工作的人，对这段代码肯定熟悉，这是南气院程序集MK法的fortran代码中的一个子程序，功能就是计算UF,UB序列的。Y为我们要检验的序列。
重点来看上面红色标注的代码，从程序的上来看，意思应该是当第i时刻数值小于j时刻数值时，累计数不做任何处理，当第i时刻数值大于或者等于j时刻数值时，累计数加1。

然而，在魏老师现代气候统计诊断与预测技术（第2版）中的P63-64概述中指明秩序列sk是第i时刻数值大于j时刻数值个数的累计数。

程序中,当第i时刻数值等于j时刻数值时，累计数也加1了，这样和魏老师书上的说法不一样啊。难道程序有错啊？

想想觉得不应该啊，南气院的程序流传已久，多少老师、学生都用过，不可能那么多大虾没发现问题，让我发现了。于是反复看了算法和程序，还是没想通。
因本人也是最近才接触MK法，fortran也没完全掌握，有此疑问，请各位指点。

PS：Matlab板块现在就有两种MK算法的程序，一种是直接从fortran翻译过来的（记得是一个叫伽蓝的网友）；另一种就消除了上面说的i，j时刻相等，累计数加1的情况。真不知道谁是对滴。

        

胆小车迷

晕倒，我指的代码是这个
DO 20 J=1,I-1
        IF(Y(I).LT.Y(J))GOTO 20
        M(I)=M(I)+1
  20   CONTINUE

mofangbao

有人讨论过，楼主有耐心的话可以查查历史的帖子哈

胆小车迷

mofangbao 发表于 2016-6-22 08:33
有人讨论过，楼主有耐心的话可以查查历史的帖子哈

谢谢管理员提醒，可能之前没找到，我翻翻

胆小车迷

mofangbao 发表于 2016-6-22 08:33
有人讨论过，楼主有耐心的话可以查查历史的帖子哈

已经找到了“作者：yutian2999  关于《现代气候统计诊断与预测技术》书中Mann-Kendall突变检验法”，结论是“wei做图的时候采用xi>=xj，得出来的突变点是错误的，应该采用xi>xj来计算”。。。
这么说，坛子几乎所有MK检验的Fortran和部分Matlab源码都有问题。。。

fhzf4321

谢谢楼主分享!