update: test文件已上传，方便大家理解代码，希望对需要的同志有所帮助！

(1)  累计距平均

data=xlsread('test.xlsx');

data(:,1)=[];

mean_data=mean(data);

diff_streamflow=data(:,1)-mean_data(1);

diff_baseflow=data(:,2)-mean_data(2);

diff_surfq=data(:,3)-mean_data(3);

Sum_streamflow_anomaly=zeros(55,1);

Sum_baseflow_anomaly=zeros(55,1);

Sum_surfq_anomaly=zeros(55,1);

for i=1:55;

Sum_streamflow_anomaly(i)=sum(diff_streamflow(1:i));

Sum_baseflow_anomaly(i)=sum(diff_baseflow(1:i));

Sum_surfq_anomaly(i)=sum(diff_surfq(1:i));

end

years=(1960:2014)';

xlswrite('cumulative_anomaly.xlsx',[years,Sum_streamflow_anomaly,Sum_baseflow_anomaly,Sum_surfq_anomaly])

通过绘图找到分界点即得到break point

(2)   Pettitt检测

% This code is used to find the change point in a univariate continuous time series

% using Pettitt Test.

%

%

% The test here assumed is two-tailed test. The hypothesis are as follow:

%  H (Null Hypothesis): There is no change point in the series

%  H(Alternative Hypothesis): There is a change point in the series

%

% Input: univariate data series

% Output:

% The output of the answer in row wise respectively,

% loc: location of the change point in the series, index value in

% the data set

% K: Pettitt Test Statistic for two tail test

% pvalue: p-value of the test

%

%Reference: Pohlert, Thorsten. "Non-Parametric Trend Tests and Change-Point Detection." (2016).

%

function a=pettitt(data)

[m n]=size(data);

for t=2:1:m

for j=1:1:m

v(t-1,j)=sign(data(t-1,1)-data(j,1));

V(t-1)=sum(v(t-1,:));

end

end

U=cumsum(V);

loc=find(abs(U)==max(abs(U)));

K=max(abs(U));

pvalue=2*exp((-6*K^2)/(m^3+m^2));

a=[loc; K ;pvalue];

return

-------------------------------------

下面是函数的调用：

clc;

data=xlsread('test.xlsx');

y=data(1:end,2);%输入数据

a=pettitt(y);

参考文献：

(1)王随继等, 皇甫川流域降水和人类活动对径流量变化的贡献率分析——累积量斜率变化率比较方法的提出及应用. 地理学报, 2012. 67(3): 第388-397页

(2)  Pohlert, Thorsten. "Non-Parametric Trend Tests and Change-Point Detection." (2016).

转载本文请联系原作者获取授权，同时请注明本文来自张凌科学网博客。

链接地址：http://blog.sciencenet.cn/blog-922140-1116580.html

上一篇：ArcGIS Licience过期的解决方法

下一篇：DHSVM模型计算流程