% M-K 趋势检定

clear, close all,clc % clear:清变数 close all:清图面 clc:清画面

% define and assign the full file path using "file open" dialog

[filename filepath]=uigetfile('data1.xls');

full_filepath=[filepath filename];

[X,TXTX,RAWX]=xlsread(full_filepath,1); % 数据全部读入,数据缺失不影响结果

x=X(:,1); % x 时间

y=X(:,2); % y 数据

% 计算 S

n=size(y,1); % 数据个数

S=0;

for i=1:n-1

S = S + sum(sign(y(i+1:n) - y(i))); % S 计算式

end

% 计算 VarS

VarS=n*(n-1)*(2*n+5)/18;

%计算 Z

if S>0

Z=(S-1)/sqrt(VarS);

else

Z=(S+1)/sqrt(VarS);

end

% 计算 Zabs

alpha1=0.05; % 信度 95% 的显著水平

alpha2=0.01; % 信度 99% 的显著水平

PZ1=norminv(1-alpha1/2,0,1);

PZ2=norminv(1-alpha2/2,0,1);

H=0; % 虚无假设

Zabs=abs(Z);

if Zabs >= PZ1

H=1;

else

H=0;

end

P_value=2*(1-normcdf(abs(Z),0,1)); % 若 P_value 比 alpha1 小,则否定虚无假设

% 计算倾斜度

ndash=n*(n -1)/2; % 对称矩阵上半部

slope1= zeros( ndash, 1 ); % 起始归零

m=0;

for k = 1:n-1,

for j = k+1:n,

m=m+1;

slope1(m) = ( y(j) - y(k) ) / (x(j) - x(k) ) ;% 分母非 (j-k)

end;

end;

slope= median( slope1 ); % 中位数

% 历线绘图

yd=max(y)-min(y);

figure

plot(x,y,'b-o','linewidth',1.5);

axis([min(x),max(x),min(y)-0.2*yd,max(y)+0.2*yd]); % 全距外扩 20%

xlabel('时间','FontName','TimesNewRoman','FontSize',12);

ylabel('数据','FontName','TimesNewRoman','Fontsize',12);

title('数据历线图') %添加标题

grid on

output='数据历线图';

saveas(gcf, output, 'jpg')

% M-K 突变检定

Sk=zeros(size(y)); % 起始归零

UFk=zeros(size(y)); % 起始归零

s1=0;

for i=2:n

for j=1:i

if y(i)>y(j)

s1=s1+1;

else

s1=s1+0;

end;

end;

Sk(i)=s1;

E=i*(i-1)/4; % 均值

Var=i*(i-1)*(2*i+5)/72; % 方差

UFk(i)=(Sk(i)-E)/sqrt(Var);

end;

% 起始归零

y2=zeros(size(y));

Sk2=zeros(size(y));

UBk=zeros(size(y));

s2=0;

for i=1:n

y2(i)=y(n-i+1); % 逆序

end;

for i=2:n

for j=1:i

if y2(i)>y2(j)

s2=s2+1;

else

s2=s2+0;

end;

end;

Sk2(i)=s2;

E=i*(i