MK突变检验（Matlab）

function [ UF,UB ] = MannKendall( x,y,p )

% x表示时间磨拍如1982-2015

% y表示对应时间的结果

% p表咐游毕示显著水平衡芹

N = length(y)

UF = SMK(y)

yy = reshape(y,1,length(y))

yy = fliplr(yy)

UB = -fliplr(SMK(yy))

function U = SMK( Y )

N = length(Y)

s = zeros(1,N)

U(1) = 0

for k=2:N

r = 0

s(k) = 0

for j=1:k-1

if Y(k)>Y(j)

r = r+1

end

s(k) = s(k-1)+r

end

E = k (k-1)/4

VAR = k (k-1) (2 k+5)/72

U(k) = (s(k)-E)/sqrt(VAR)

end

分别运用分层聚类、K均值聚类以及高斯念散混仔衫合模型来进行分析，然后比较三者的结果仔戚氏

生成随机二维分布图形，三个中心

% 使用高斯分布（正态分布）

% 随机生成3个中心以及标准差

s = rng(5,'v5normal')

mu = round((rand(3,2)-0.5)*19)+1

sigma = round(rand(3,2)*40)/10+1

X = [mvnrnd(mu(1,:),sigma(1,:),200)...

mvnrnd(mu(2,:),sigma(2,:),300)...

mvnrnd(mu(3,:),sigma(3,:),400)]

UF,UB]=MannKendall(Y)

% % Y即为待检测的时间序列

figure,hold on

plot(UF,'-b')

plot(UB,'-r')

function [UF,UB]=MannKendall(timeseries)

N=length(timeseries)

UF=SMK(timeseries,N)

for i=1:N

YY(i)=timeseries(N+1-i)

end

u_res=SMK(YY,N)

for i=1:N

UB(i)=-u_res(N+1-i)

end

function u_res=SMK(Y,N)

m_res=zeros(N,1)md_res=zeros(N,1)u_res=zeros(N,1)

m_res(1)=0

for i=2:N

m_res(i)=0

md_res(i)=0

for j=1:i-1

if Y(i)<Y(j)

m_res(i)=m_res(i)+0

else

m_res(i)=m_res(i)+1

end

md_res(i)=md_res(i-1)+m_res(i)

end

end

u_res(1)=0

for i=2:N

E=i*(i-1)/4

VAR=i*(i-1)*(2*i+5)/72

u_res(i)=(md_res(i)-E)/sqrt(VAR)

end

---