matlab怎么做M-K突变，求大神好友

279461970 · 发表于 2016-5-16 11:51:53

登录后查看更多精彩内容~

您需要登录才可以下载或查看，没有帐号？立即注册

x

本人小白，现不知怎么用matlab导入excel数据做M-K突变分析。
说明：1.XX地区近60年年平均气温突变分析，数据已经处理完了（excel里第一列为年份，第二列为年平均温度）；
2.运算公式程序也有；
3.如何导入（erxcel中格式怎么处理，行排序还是列），如何处理运算。
求大神指导帮忙，刚接触matlab1个星期，望指教。

jiajian0512 · 发表于 2016-5-16 14:48:54

% Mann-Kendall突变检验UF UB
% 数据序列y
% 结果序列UFk，UBk2
%读取excel中的数据，赋给矩阵y,获取y的样本数
%A为时间和径流数据列
clc;
clear all;
t=xlsread('h:\极端气温变化\年最高气温1.xls');
year=t(:,1);%时间序列
y=t(:,2);%极端最高温度
%tmin=t(:,4);极端最低温度
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
plot(year,UFk,'r-','linewidth',1.5);
hold on
plot(year,UBk2,'b-.','linewidth',1.5);
plot(year,1.96*ones(N,1),':','linewidth',1);
axis([min(year),max(year),-4,8]);
legend('UF统计量','UB统计量','0.05显著水平');
xlabel('t (year)','FontName','TimesNewRoman','FontSize',12);
ylabel('统计量','FontName','TimesNewRoman','Fontsize',12);
%grid on
hold on
plot(year,0*ones(N,1),'-.','linewidth',1);
plot(year,1.96*ones(N,1),':','linewidth',1);
plot(year,-1.96*ones(N,1),':','linewidth',1);

复制代码

jiajian0512 · 发表于 2016-5-16 14:48:55

% Mann-Kendall突变检验UF UB
% 数据序列y
% 结果序列UFk，UBk2
%读取excel中的数据，赋给矩阵y,获取y的样本数
%A为时间和径流数据列
clc;
clear all;
t=xlsread('h:\极端气温变化\年最高气温1.xls');
year=t(:,1);%时间序列
y=t(:,2);%极端最高温度
%tmin=t(:,4);极端最低温度
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
plot(year,UFk,'r-','linewidth',1.5);
hold on
plot(year,UBk2,'b-.','linewidth',1.5);
plot(year,1.96*ones(N,1),':','linewidth',1);
axis([min(year),max(year),-4,8]);
legend('UF统计量','UB统计量','0.05显著水平');
xlabel('t (year)','FontName','TimesNewRoman','FontSize',12);
ylabel('统计量','FontName','TimesNewRoman','Fontsize',12);
%grid on
hold on
plot(year,0*ones(N,1),'-.','linewidth',1);
plot(year,1.96*ones(N,1),':','linewidth',1);
plot(year,-1.96*ones(N,1),':','linewidth',1);

复制代码

279461970 · 发表于 2016-5-16 12:07:08

C:\Users\Administrator\Desktop\QQ截图20160516120659.jpg

279461970 · 发表于 2016-5-16 12:10:55

excel做出的图

279461970 · 发表于 2016-5-16 12:13:55

学习matlab是想验证下，觉得软件做出来应该准确些，大神帮帮忙，谢谢了

jiajian0512 · 发表于 2016-5-16 14:49:21

把这个程序改改应该就可以出图了

279461970 · 发表于 2016-5-16 19:37:26

jiajian0512 发表于 2016-5-16 14:49
把这个程序改改应该就可以出图了

非常感谢。我一会试试。

wsjqxjy · 发表于 2016-5-17 08:36:17

{:eb502:}{:eb502:}

sha · 发表于 2016-5-17 08:44:14

感谢指导，正巧有同样问题!

{:eb301:}{:eb301:}

		自动登录	找回密码
密码			立即注册

[讨论] matlab怎么做M-K突变，求大神好友

登录后查看更多精彩内容~