MATLAB,M-K和滑动T检验分析气候突变，帮我看看分析的正确吗？

落侠孤鹜

数据：

1957	7.4 （℃）
1958	8.1
1959	7.9
1960	8.0
1961	7.9
1962	7.6
1963	7.7
1964	7.5
1965	7.7
1966	8.0
1967	6.9
1968	7.6
1969	7.8
1970	7.1
1971	8.0
1972	7.6
1973	8.2
1974	7.6
1975	7.5
1976	6.9
1977	7.5
1978	7.9
1979	7.7
1980	7.6
1981	7.5
1982	7.7
1983	7.1
1984	6.8
1985	7.5
1986	7.7
1987	8.3
1988	7.6
1989	7.5
1990	7.9
1991	7.9
1992	7.5
1993	7.6
1994	8.1
1995	8.1
1996	7.7
1997	8.6
1998	9.1
1999	8.6
2000	8.5
2001	8.6
2002	8.9
2003	8.5
2004	8.5
2005	8.3
2006	9.4
2007	8.9
2008	8.3
2009	8.7
2010	8.8
2011	8.4

1.MK检验
程序（显著性=0.05）
% Mann-Kendall突变检测
% 数据序列y
% 结果序列UFk，UBk2
%--------------------------------------------
%读取excel中的数据，赋给矩阵y
%获取y的样本数
%A为时间和径流数据列
A=xlsread('F:\DATA.xls')
x=A(:,1);%时间序列
y=A(:,2);%径流数据列
N=length(y);
n=length(y);
% 正序列计算---------------------------------
% 定义累计量序列Sk，长度=y，初始值=0
Sk=zeros(size(y));
% 定义统计量UFk，长度=y，初始值=0
UFk=zeros(size(y));
% 定义Sk序列元素s
s = 0;
% i从2开始，因为根据统计量UFk公式，i=1时，Sk(1)、E(1)、Var(1)均为0
% 此时UFk无意义，因此公式中，令UFk(1)=0
for i=2:n
    for j=1:i
          if y(i)>y(j)
            s=s+1;
          else
            s=s+0;
          end;
    end;
    Sk(i)=s;
    E=i*(i-1)/4; % Sk(i)的均值
  Var=i*(i-1)*(2*i+5)/72; % Sk(i)的方差
  UFk(i)=(Sk(i)-E)/sqrt(Var);
end;
% ------------------------------正序列计算end
% 逆序列计算---------------------------------
% 构造逆序列y2，长度=y，初始值=0
y2=zeros(size(y));
% 定义逆序累计量序列Sk2，长度=y，初始值=0
Sk2=zeros(size(y));
% 定义逆序统计量UBk，长度=y，初始值=0
UBk=zeros(size(y));
% s归0
s=0;
% 按时间序列逆转样本y
% 也可以使用y2=flipud(y);或者y2=flipdim(y,1);
for i=1:n
     y2(i)=y(n-i+1);
end;
% i从2开始，因为根据统计量UBk公式，i=1时，Sk2(1)、E(1)、Var(1)均为0
% 此时UBk无意义，因此公式中，令UBk(1)=0
for i=2:n
    for j=1:i
          if y2(i)>y2(j)
            s=s+1;
          else
            s=s+0;
          end;
    end;
    Sk2(i)=s;
    E=i*(i-1)/4; % Sk2(i)的均值
  Var=i*(i-1)*(2*i+5)/72; % Sk2(i)的方差
% 由于对逆序序列的累计量Sk2的构建中，依然用的是累加法，即后者大于前者时s加1，
% 则s的大小表征了一种上升的趋势的大小，而序列逆序以后，应当表现出与原序列相反
% 的趋势表现，因此，用累加法统计Sk2序列，统计量公式(S(i)-E(i))/sqrt(Var(i))
% 也不应改变，但统计量UBk应取相反数以表征正确的逆序序列的趋势
  UBk(i)=0-(Sk2(i)-E)/sqrt(Var);
end;
% ------------------------------逆序列计算end
% 此时上一步的到UBk表现的是逆序列在逆序时间上的趋势统计量
% 与UFk做图寻找突变点时，2条曲线应具有同样的时间轴，因此
% 再按时间序列逆转结果统计量UBk，得到时间正序的UBk2，做图用
UBk2=zeros(size(y));
% 也可以使用UBk2=flipud(UBk);或者UBk2=flipdim(UBk,1);
for i=1:n
    UBk2(i)=UBk(n-i+1);
end;
% 做突变检测图时，使用UFk和UBk2
% 写入目标xlsx文件：C:\Users\chenfei\Desktop\test2.xlsx
% 目标表单：Sheet1
% 目标区域：UFk从A1开始，UBk2从B1开始
xlswrite('F:\test.xls',UFk,'Sheet1','A1');
xlswrite('F:\test.xls',UBk2,'Sheet1','B1');
figure(3)%画图
  figure1 = figure('Color',[1 1 1]);
plot(x,UFk,'r-','linewidth',1.5);
hold on
plot(x,UBk2,'b-.','linewidth',1.5);
plot(x,1.96*ones(N,1),':','linewidth',1);
axis([min(x),max(x),-6,6]);
legend('UF统计量','UB统计量','0.05显著水平');
xlabel('t (year)','FontName','TimesNewRoman','FontSize',12);
ylabel('统计量','FontName','TimesNewRoman','Fontsize',12);
%grid on
hold on
plot(x,0*ones(N,1),'-.','linewidth',1);
plot(x,1.96*ones(N,1),':','linewidth',1);
plot(x,-1.96*ones(N,1),':','linewidth',1);
图1

2.滑动T检验
（1）程序（步长为5,0.05）
clc,clear;
%% 读取数据  
excelFile='F:\DATA.XLS';% 文件路径  
myData=xlsread(excelFile);% 读取数据
timeSeries=myData(:,1);  % 时间序列数据
dataSeries=myData(:,2);  %  
dataCount=length(dataSeries);
%% 设置步长与检验值
step=5;  % 步长  
v=step+step-2; % 计算自由度  
ttest=2.306; % 查表得t检验值，修改
len1=step;
len2=step;  
x=timeSeries(step:dataCount-step);
for i=step:dataCount-step      
    n1=dataSeries(i-step+1:i);   
    n2=dataSeries(i+1:i+step);     
    mean1=mean(n1);     
    mean2=mean(n2);      
    c=(len1+len2)/(len1*len2);      
    var1=1/len1*sum((n1-mean1).^2);     
    var2=1/len2*sum((n2-mean2).^2);     
    delta1=len1*var1+len2*var2;     
    delta=delta1/(len1+len2-2);      
    t(i-step+1)=(mean1-mean2)/sqrt(delta*c);
end
%% 制图
figure(1);
figure1 = figure('Color',[1 1 1]);
plot(x,t,'r-','linewidth',1.5);  
xlabel('t (year)','FontName','TimesNewRoman','FontSize',12);
ylabel('统计量','FontName','TimesNewRoman','Fontsize',12);
axis([min(x),max(x),-10,10]);
hold on  
plot(x,0*ones(i-step+1,1),'-.','linewidth',1);  
plot(x,ttest*ones(i-step+1,1),':','linewidth',1);% 更改数字
plot(x,-ttest*ones(i-step+1,1),':','linewidth',1);% 更改数字
legend('t统计量','0.05显著水平');
图2

（2）程序（步长为10,0.05）
ttest=2.101
图3


3.分析
（1）MK有一个交叉点，初步判断1995年前后发生了突变；
（2）T检验步长为5的时候在1994-1995年后发生突变；
（3）T检验步长为10的时候有两处超过显著水平，其中1987年之后的突变明显。
4.结论
1995年发生过一次明显的突变。

jiangsaiping

15682865660 发表于 2019-4-10 11:25
等号的后面加一个负号

为什么要加负号呢？我见定义是前端序列均值减去后端序列均值，所得差值再除以一堆数。还是不明白

落侠孤鹜

顶起来，大神们帮我分析分析

落侠孤鹜

数据分析已完成，经分析突变发生在1997年。程序中t(i-step+1)=(mean1-mean2)/sqrt(delta*c); 这里应是-号。

Zoe_zuo

分析突变发生在1997年

落侠孤鹜

Zoe_zuo 发表于 2018-1-9 17:30
分析突变发生在1997年

是啊，帖子发了几天没人回，自己弄了好久终于研究懂了。

jiangsaiping

落侠孤鹜 发表于 2018-1-9 17:21
数据分析已完成，经分析突变发生在1997年。程序中t(i-step+1)=(mean1-mean2)/sqrt(delta*c); 这里应是-号。

楼主，您好！请问程序哪里改成-号？按照您的数据，我现在画出的还是您最初画出来的结果。不知道改哪里？

15682865660

jiangsaiping 发表于 2019-4-9 17:53
楼主，您好！请问程序哪里改成-号？按照您的数据，我现在画出的还是您最初画出来的结果。不知道改哪里？

等号的后面加一个负号

编程渣渣

楼主程序可以分享一下嘛

jiangsaiping

15682865660 发表于 2019-4-10 11:25
等号的后面加一个负号

为什么要加负号呀？还是不懂。书上定义是前端子序列均值减去后端子序列均值，所得结果再除以一堆数。