MK检验突变分析 matlab

24444

有人分享吗

jiangph1986

同求啊，这个有没有好码啊！

小英

1. % Mann-Kendall突变检测
   
2. % 数据序列y
   
3. % 结果序列UFk，UBk2
   
4. %读取excel中的数据，赋给矩阵y
   
5. %获取y的样本数
   
6. %A为时间和降水数据列
   
7. x=降水(:,1);%时间序列
   
8. y=降水(:,2);%降水数据列
   
9. N=length(x);
   
10. n=length(y);
    
11. % 正序列计算---------------------------------
    
12. % 定义累计量序列Sk，长度=y，初始值=0
    
13. Sk=zeros(size(y));
    
14. % 定义统计量UFk，长度=y，初始值=0
    
15. UFk=zeros(size(y));
    
16. % 定义Sk序列元素s
    
17. s = 0;
    
18. % i从2开始，因为根据统计量UFk公式，i=1时，Sk(1)、E(1)、Var(1)均为0
    
19. % 此时UFk无意义，因此公式中，令UFk(1)=0
    
20. for i=2:n
    
21. for j=1:i
    
22. if y(i)>y(j)
    
23. s=s+1;
    
24. else
    
25. s=s+0;
    
26. end;
    
27. end;
    
28. Sk(i)=s;
    
29. E=i*(i-1)/4; % Sk(i)的均值
    
30. Var=i*(i-1)*(2*i+5)/72; % Sk(i)的方差
    
31. UFk(i)=(Sk(i)-E)/sqrt(Var);
    
32. end;
    
33. % ------------------------------正序列计算end
    
34. % 逆序列计算---------------------------------
    
35. % 构造逆序列y2，长度=y，初始值=0
    
36. y2=zeros(size(y));
    
37. % 定义逆序累计量序列Sk2，长度=y，初始值=0
    
38. Sk2=zeros(size(y));
    
39. % 定义逆序统计量UBk，长度=y，初始值=0
    
40. UBk=zeros(size(y));
    
41. % s归0
    
42. s=0;
    
43. % 按时间序列逆转样本y
    
44. % 也可以使用y2=flipud(y);或者y2=flipdim(y,1);
    
45. for i=1:n
    
46. y2(i)=y(n-i+1);
    
47. end;
    
48. % i从2开始，因为根据统计量UBk公式，i=1时，Sk2(1)、E(1)、Var(1)均为0
    
49. % 此时UBk无意义，因此公式中，令UBk(1)=0
    
50. for i=2:n
    
51. for j=1:i
    
52. if y2(i)>y2(j)
    
53. s=s+1;
    
54. else
    
55. s=s+0;
    
56. end;
    
57. end;
    
58. Sk2(i)=s;
    
59. E=i*(i-1)/4; % Sk2(i)的均值
    
60. Var=i*(i-1)*(2*i+5)/72; % Sk2(i)的方差
    
61. % 由于对逆序序列的累计量Sk2的构建中，依然用的是累加法，即后者大于前者时s加1，
    
62. % 则s的大小表征了一种上升的趋势的大小，而序列逆序以后，应当表现出与原序列相反
    
63. % 的趋势表现，因此，用累加法统计Sk2序列，统计量公式(S(i)-E(i))/sqrt(Var(i))
    
64. % 也不应改变，但统计量UBk应取相反数以表征正确的逆序序列的趋势
    
65. UBk(i)=0-(Sk2(i)-E)/sqrt(Var);
    
66. end;
    
67. % ------------------------------逆序列计算end
    
68. % 此时上一步的到UBk表现的是逆序列在逆序时间上的趋势统计量
    
69. % 与UFk做图寻找突变点时，2条曲线应具有同样的时间轴，因此
    
70. % 再按时间序列逆转结果统计量UBk，得到时间正序的UBk2，做图用
    
71. UBk2=zeros(size(y));
    
72. % 也可以使用UBk2=flipud(UBk);或者UBk2=flipdim(UBk,1);
    
73. for i=1:n
    
74. UBk2(i)=UBk(n-i+1);
    
75. end;
    
76. % 做突变检测图时，使用UFk和UBk2
    
77. figure(3)%画图
    
78. plot(x,UFk,'r-','linewidth',1);
    
79. hold on
    
80. plot(x,UBk2,'b-.','linewidth',1);
    
81. plot(x,1.96*ones(N,1),':','linewidth',0.5);
    
82. plot(x,2.56*ones(N,1),':','linewidth',0.5);
    
83. axis([min(x),max(x),-5,5]);
    
84. legend('UF统计量','UB统计量');
    
85. xlabel('年份','FontName','TimesNewRoman','FontSize',9);
    
86. ylabel('统计量','FontName','TimesNewRoman','Fontsize',9);
    
87. %grid on
    
88. hold on
    
89. plot(x,0*ones(N,1),'-.','linewidth',0.5);
    
90. plot(x,1.96*ones(N,1),':','linewidth',0.5);
    
91. plot(x,-1.96*ones(N,1),':','linewidth',0.5);
    
92. plot(x,2.56*ones(N,1),':','linewidth',0.5);
    
93. plot(x,-2.56*ones(N,1),':','linewidth',0.5);
    
94. time=1951:6:2014
    
95. Xlabel('年份');
    
96. Ylabel('统计量');

复制代码

我也是下载别人的，你试试看可不可以

24444

求帮助

我是小猪

LZ善用搜索功能

guoqu510

学习一下先

ttttt

金币不够使劲凑

huangzeqin

{:eb502:}{:eb502:}{:eb502:}学习了

curtainsky

学习了学习了学习了