曼-肯德尔法matlab分享

/tpに泳ぐ鱼

          曼-肯德尔法的计算公式，在不同的参考资料中，略有不同。有些参考资料还描述的有些模糊和错误，请同学们在使用时注意。

                   在这里，我们用一个实例来说明曼-肯德尔法的具体做法。

         应用实例

                   我们现在用曼-肯德尔(Mann-Kendall)法检测1900-1990年上海年平均气温序列的突变。     

         计算步骤

                   在计算时，对于具有n个样本量的时间序列x,我们首先需要构造一个秩序列sk，秩序列sk是第i 时刻的数值大于j 时刻数值的个数的累计数。若想获得魏凤英老师《现代气候统计诊断与预测技术》一书中的曼-肯德尔统计图线图（也就是书中的图5.2），那么我们不能使用这本书中的计算公式。        

                   在这里，sk=sk-1+从i=1到k-1对ri求和，k=3,4,...,n。s2=r1。其中：当xk≥xi时，ri=1,当xk<xi时，ri=0。i=1,2,...,k-1。

                   然后，在时间序列随机独立的假定下，定义统计量UFk=(sk-E(Sk))/根号下var(Sk)，k=2,3,...,n

                            其中E(Sk)，var(Sk)是累计数Sk的均值和方差。在x1,x2,…,xn相互独立，而且有相同连续分布时，E(Sk)=k(k-1)/4,var(sk)=k(k-1)(2k+5)/72,k=2,3,...,n。

                            当k=1时，sk、E(Sk)，var(Sk)均无意义，我们将UF1赋值为0。

                   所构造的统计量UFk为标准正态分布，它是按时间序列x的顺序x1, x2,…,xn计算出的统计量序列。

                   给定显著性水平α，根据正态分布的对称性，从正态分布函数表上查出与α/2水平相应的数值，即可确定出临界值uα=u(0.025)=1.96 ，如果|UFk |>uα，则表明序列存在着明显的趋势变化。uα也可由MATLAB的norminv(1-α/2,0,1)计算得到。

                   除了正序统计量UFk外，我们还需要计算逆序统计量UBk。

                   首先将时间序列x 逆序排列，也就是xn, xn-1,…, x1，得到一个序列y1,y2,...,yn;

然后计算s2k=s2k-1+从i=1到k-1对ri求和，k=3,4,...,n。k=2时，s22=r1。其中：当yk≥yi时，ri=1,当yk<yi时，ri=0。i=1,2,...,k-1。

                   然后计算累计数S2k的均值和方差E(S2k)，var(S2k)。在y1, y2,…,yn相互独立，而且有相同连续分布时，E(S2k)=k(k-1)/4,var(S2k)=k(k-1)(2k+5)/72,k=2,3,...,n。均值和方差的计算公式跟前面是一样的。

                   然后，在时间序列随机独立的假定下，定义统计量UBk=-UFk=-(sk-E(Sk))/根号下var(Sk)，k=2,3,...,n。

                   这里得到的UBk表现的是逆序列在逆序时间上的趋势统计量。在与UFk做曼-肯德尔统计曲线图寻找突变点时，UFk和UBk这两条曲线应该具有相同的时间轴，因此，再将统计量UBk按时间序列逆序排列，得到时间正序的UBk2。

                   最后，将UFk和UBk2这两条曲线，以及临界值u0.05=±1.96这两条直线均绘在同一张图上，得到曼-肯德尔统计曲线图。

伽蓝鸟

这东西很早前就有人发过了。。还需要两个贡献啊。。

/tpに泳ぐ鱼

伽蓝鸟 发表于 2019-10-27 11:52
这东西很早前就有人发过了。。还需要两个贡献啊。。

是因为都没有跟书上的图一样的。这个是查阅大量的资料，才搞定的。

/tpに泳ぐ鱼

伽蓝鸟 发表于 2019-10-27 11:52
这东西很早前就有人发过了。。还需要两个贡献啊。。

是因为没有跟书上一样的图，书上的算法并不能得到书上的图。我也是查阅了大量的资料，画了很多功夫才弄清楚的。不知道为什么提交了多次，回复不了。

zhouxiaosky

多谢分享~~~

何佳田力

这个是我自己写的，可以拿魏凤英老师的数据跑一下，我记得自己测试的时候是完全一样的
function [UFk,UBk,je]=mk(D,T,ALPHA,k)
data=D;time=T;U=norminv(1-ALPHA/2,0,1);%显著水平
N=size(data,1);n=2;nb=1;nB=2;ri=zeros(N,1);rbi=zeros(N,1);s=0;datac=flipud(data);
SFk=zeros(N,1);SBk=zeros(N,1);UFk=zeros(N,1);UBk=zeros(N,1);%data be limited in column
%in inverse order(RbI,UBK,SBK) & natural order(RI,UFK,SFK)
for n=2:N
    ri(n)=sum(double(data(1:n-1)<=data(n)));
    rbi(n)=sum(double(datac(1:n-1)<=datac(n)));
    SFk(n)=SFk(n-1)+ri(n);
    SBk(n)=SBk(n-1)+rbi(n);
    UFk(n)=(SFk(n)-n*(n-1)/4)/sqrt(n*(n-1)*(2*n+5)/72);
    UBk(N-n+1)=(-(SBk(n)-n*(n-1)/4)/sqrt(n*(n-1)*(2*n+5)/72));
end