m_k.m

时间序列趋势分析中

clear
clc

load 'Murrumbidgee_Balranald.dat';

%x=mean(hankou');%计算年平均序列
x=Murrumbidgee_Balranald;
e=1900;
d=2006;
[a,b]=size(x);

%*****************逆序生成序列y*************************

y=zeros(a,1);
for i=a:-1:1
    j=a-i+1;
    y(j)=x(i);
end

%*******************计算累计数m*************************

%m=zeros(a,1);
%d=zeros(a,1);
%m(1,1)=0;
%d(1,1)=0;
for i=2:a
    s1_1=0;s1_2=0;
    for j=1:i-1
        if x(i,1)>x(j,1)
            s1_1=s1_1+1;
        end
        if y(i,1)>y(j,1)
            s1_2=s1_2+1;
        end
        m1(i,1)=s1_1;
        m2(i,1)=s1_2;
    end
end

%u_d为正值表示增加趋势，负值表示减少趋势. u_d值在大于等于　∣1.28∣、∣1.64∣、
%∣2.32∣时分别表示通过了90%、95%、99%置信度显著性检验.
%E=zeros(a,1);
%var_d=zeros(a,1);
%u_d=zeros(a,1);
for k=2:a
    s2_1=0;s2_2=0;
    for i=1:k
        s2_1=s2_1+m1(i,1);
        s2_2=s2_2+m2(i,1);
    end
    d1(k,1)=s2_1;d2(k,1)=s2_2;
    E(k,1)=k*(k+1)/4;
    var_d(k,1)=k*(k-1)*(2*k+5)/72;
end

c=sqrt(var_d);
ufk=(d1-E)./c;
ubk0=-(d2-E)./c;

ubk=zeros(a,1);
for i=a:-1:1
    j=a-i+1;
    ubk(j)=ubk0(i);
end

ufkk(1:a-1,1)=ufk(2:a,1);
ubkk(1:a-1,1)=ubk(1:a-1,1);
xlswrite('pnp_MK',[ufk ubk]);

xx=e:d;
for i=1:a
    r1(i)=-1.96;
    r2(i)=1.96;
end
plot(xx,ufk,'k',xx,ubk,':k',xx,r1,'-.k',xx,r2,'-.k')
xlim([e d])
y1=[min(ufk),min(ubk),-1.96];
y3=[max(ufk),max(ubk),1.96];
ylim([1.2*min(y1) 1.2*max(y3)])
%legend('UFK','UBK')
xlabel('时间/年');
ylabel('MK统计值');