MATLAB 进行MK检验和滑动T检验问题

liuwenyan150825

用MK检验与滑动T检验得到的结果不一致，怎么办？
MK检验程序

1. <p><p><p>MK检验程序
   
2. clc,clear;
   
3. [a,textdata]=xlsread('问题.xlsx');
   
4.    s=mean(reshape(a(:,2),12,[]));
   
5. s;
   
6. timeSeries=[1960:2010]';
   
7. dataSeries=s';  %
   
8. dataCount=length(dataSeries);
   
9. %% 正序列计算
   
10. SK=zeros(size(dataSeries));  % 定义SK
    
11. UFK=zeros(size(dataSeries));  % 定义UFK
    
12. s=0;
    
13. for i=2:dataCount
    
14.     for j=1:i
    
15.         if dataSeries(i)>dataSeries(j)
    
16.             s=s+1;
    
17.         else
    
18.             s=s+0;
    
19.         end
    
20.     end
    
21.     SK(i)=s;
    
22.     E=i*(i-1)/4;  % SK均值
    
23.     Var=i*(i-1)*(2*i+5)/72; % SK方差
    
24.     UFK(i)=(SK(i)-E)/sqrt(Var);
    
25. end
    
26. %% 逆序列计算
    
27. oppSK=zeros(size(dataSeries));
    
28. UBK=zeros(size(dataSeries));
    
29. s=0;
    
30. opphumiditySeries=flipud(dataSeries);
    
31. for i=2:dataCount
    
32.     for j=1:i
    
33.         if opphumiditySeries(i)>opphumiditySeries(j)
    
34.             s=s+1;
    
35.         else
    
36.             s=s+0;
    
37.         end
    
38.     end
    
39.     oppSK(i)=s;
    
40.     E=i*(i-1)/4;  % SK均值
    
41.     Var=i*(i-1)*(2*i+5)/72; % SK方差
    
42.     UBK(i)=0-(oppSK(i)-E)/sqrt(Var);
    
43. end
    
44. oppUBK=flipud(UBK);
    
45. %  xlswrite(outputdata,UFK,'Sheet2','A1');% 输出结果
    
46. %  xlswrite(outputdata,oppUBK,'Sheet2','B1');
    
47. %% 制图
    
48. figure(3);
    
49. plot(timeSeries,UFK,'r-','linewidth',1.5);
    
50. hold on
    
51. plot(timeSeries,oppUBK,'b-.','linewidth',1.5);
    
52. plot(timeSeries,1.96*ones(dataCount,1),':','linewidth',1);
    
53. axis([min(timeSeries),max(timeSeries),-5,5]);
    
54. legend('UF统计量','UB统计量','0.05显著水平');
    
55. xlabel('t (year)','FontName','TimesNewRoman','FontSize',12);
    
56. ylabel('统计量','FontName','TimesNewRoman','Fontsize',12);
    
57. hold on
    
58. plot(timeSeries,0*ones(dataCount,1),'-.','linewidth',1);
    
59. plot(timeSeries,1.96*ones(dataCount,1),':','linewidth',1);
    
60. plot(timeSeries,-1.96*ones(dataCount,1),':','linewidth',1);
    
61. 
    
62. p=[];u=[];
    
63. for n=1:length(UFK)
    
64.        if UFK(n)-oppUBK(n)>0
    
65.            p=[p,n];
    
66.        else
    
67.            u=[u,n];
    
68.        end
    
69. end
    
70. p
    
71. u
    
72. nian=timeSeries;
    
73. delete=nian(p)
    
74. 
    
75. 
    
76. 
    
77. 
    
78. 
    
79. 滑动T检验程序
    
80. 
    
81. clc,clear;
    
82. [a,textdata]=xlsread('问题.xlsx');
    
83.    s=mean(reshape(a(:,2),12,[]));
    
84. s;
    
85. timeSeries=[1960:2010]';
    
86. dataSeries=s';  %
    
87. dataCount=length(dataSeries);
    
88. 
    
89. %% 设置步长与检验值
    
90. step=10;  % 步长
    
91. v=step+step-2; % 计算自由度
    
92. ttest=1.734; % 查表得t检验值，修改
    
93. len1=step;
    
94. len2=step;
    
95. x=timeSeries(step:dataCount-step);
    
96. for i=step:dataCount-step
    
97.     n1=dataSeries(i-step+1:i);
    
98.     n2=dataSeries(i+1:i+step);
    
99.     mean1=mean(n1);
    
100.     mean2=mean(n2);
     
101.     c=(len1+len2)/(len1*len2);
     
102.     var1=1/len1*sum((n1-mean1).^2);
     
103.     var2=1/len2*sum((n2-mean2).^2);
     
104.     delta1=len1*var1+len2*var2;
     
105.     delta=delta1/(len1+len2-2);
     
106.     t(i-step+1)=(mean1-mean2)/sqrt(delta*c);
     
107. end
     
108. %% 制图
     
109. figure(1);
     
110. plot(x,t,'r-','linewidth',1.5);
     
111. xlabel('t (年)','FontName','TimesNewRoman','FontSize',12);
     
112. ylabel('统计量','FontName','TimesNewRoman','Fontsize',12);
     
113. axis([min(x),max(x),-4,4]);
     
114. hold on
     
115. plot(x,0*ones(i-step+1,1),'-.','linewidth',1);
     
116. plot(x,ttest*ones(i-step+1,1),':','linewidth',1);% 更改数字
     
117. plot(x,-ttest*ones(i-step+1,1),':','linewidth',1);% 更改数字
     
118. legend('t统计量','0.05显著水平');
     
119. 所用数据为51年的观测值，上传不上去啊</p><p>1028416496这是我的qq，知道的请加我一下，我可以把数据再发给您
     
120. </p>

复制代码

snowie

你可以尝试改变滑动t检验的步长试一下

liuwenyan150825

不好意思，刚看到，步长我基本上已经都尝试了，但是还是不行，要不就是MK检验与滑动T检验的不重合，要不就是有两个交点