greensim
2010-03-06, 17:16
本源码实现了下面参考文献中的算法，并对该文献中的实例进行了仿真，属于GreenSim团队原创作品，转载请注明。
封国林,龚志强,董文杰等.基于启发式分割算法的气候突变检测研究[J].物理学报,2005,54(11):5494-5499
function [FLAG,AllT,AllTmax,AllPTmax]=BGA(X,P0,L0)
%% 非平稳时间序列突变检测的启发式分割算法——BG算法
% GreenSim团队原创作品，转载请注明
% 欢迎访问GreenSim——算法仿真团队→http://blog.sina.com.cn/greensim
%% 输入参数列表
% X 待检测的数据，列向量存储
% P0 显著性水平门限值，低于此值的不再分割
% L0 最小分割尺度，子段长度小于此值的不再分割
%% 输出参数列表
% FLAG 分割点标记，列向量存储，长度与X相同
% AllT 与分割点对应的全部t检验序列，其首位数字为起点坐标
% AllTmax 与分割点对应的全部t检验序列的最大值
% AllPTmax 与分割点对应的全部t检验序列对应的统计显著性
%% 第一步：变量初始化
N=length(X);
FLAG=zeros(N,1);
FLAG(1)=0.1;
FLAG(N)=0.1;
AllT=cell(0,0);
AllTmax=cell(0,0);
AllPTmax=cell(0,0);
%% 第二步：产生第一个突变点，并对序列进行分割
[T,Tmax,p,PTmax]=Tseries(X);
T=[1;T];
if PTmax<P0
flag=FLAG;
flag(1)=0;
flag(N)=0;
pos3=flag(pos2);
return
end
%记录输出数据
FLAG(p)=1;
AllT=[AllT;T];
AllTmax=[AllTmax;Tmax];
AllPTmax=[AllPTmax;PTmax];
%以下为两个控制计数器
counter=2;%下一个突变点的序号
TC=0;%临时计数器
%%
while 1%设置死循环
%% 第三步：对每个段进行突变检测，能分割则分割，直到不能分割为止
pos=find(FLAG>0);
M=length(pos)-1;%当前子段数目
for m=1:M
s=pos(m);
t=pos(m+1);
L=length(SubX);
if L>=L0
[T,Tmax,p,PTmax]=Tseries(SubX);
T=[s;T];
if PTmax>=P0
TC=TC+1;
FLAG(s+p-1)=counter;
AllT=[AllT;T];
AllPTmax=[AllPTmax;PTmax];
counter=counter+1;
end
end
end
%% 第四步：返回输出数据
if TC==0
flag=FLAG;
flag(1)=0;
flag(N)=0;
pos3=flag(pos2);
FLAG=[pos2,pos3];
return
end
%%
TC=0;
%%
end
function [T,Tmax,p,PTmax]=Tseries(x)
%% 计算t检验统计序列的子函数
% GreenSim团队原创作品，转载请注明
% Email:[email protected]
% GreenSim团队主页：http://blog.sina.com.cn/greensim
% 欢迎访问GreenSim——算法仿真团队→http://blog.sina.com.cn/greensim
%% 参数列表
% x 时间序列，N×1列向量
% T t检验序列，N×1列向量
% Tmax t检验序列的最大值
% p t检验序列最大值对应的下标
% PTmax Tmax对应的统计显著性
%% 参数初始化
N=length(x);
T=zeros(N,1);
%% 以下是主循环，用于创建t检验序列
for i=3:(N-2)%最左边以及最右边的两个点没有对应的t检验值（或者说，其值初始化为0）
x1=x(1:i);%序列左边部分
N1=length(x1);%左边序列的长度
x2=x(i:N);%序列右边部分
N2=length(x2);%右边序列的长度
mean_x2=mean(x2);%右边部分的均值
std_x2=std(x2);%右边部分的标准差
%下面是计算合并偏差的公式，中英文文献里的这个公式略有不同，此处以英文文献为准
SD=sqrt(1/N1+1/N2)*sqrt(((N1-1)*std_x1^2+(N2-1)*std_x2^2)/(N1+N2-2));
T(i)=abs((mean_x1-mean_x2)/SD);
end
%% 计算其它三个输出参数
Tmax=max(T);%t检验序列的最大值
pos=find(T==Tmax);
Eta=4.19*log(N)-11.54;%计算PTmax用的参数
Delta=0.40;%计算PTmax用的参数
v=N-2;%计算PTmax用的参数
c=v/(v+Tmax^2);%不完全B函数的下标
PTmax=(1-betainc(c,Delta*Eta,Delta));%调用不完全beta函数
封国林,龚志强,董文杰等.基于启发式分割算法的气候突变检测研究[J].物理学报,2005,54(11):5494-5499
function [FLAG,AllT,AllTmax,AllPTmax]=BGA(X,P0,L0)
%% 非平稳时间序列突变检测的启发式分割算法——BG算法
% GreenSim团队原创作品，转载请注明
% 欢迎访问GreenSim——算法仿真团队→http://blog.sina.com.cn/greensim
%% 输入参数列表
% X 待检测的数据，列向量存储
% P0 显著性水平门限值，低于此值的不再分割
% L0 最小分割尺度，子段长度小于此值的不再分割
%% 输出参数列表
% FLAG 分割点标记，列向量存储，长度与X相同
% AllT 与分割点对应的全部t检验序列，其首位数字为起点坐标
% AllTmax 与分割点对应的全部t检验序列的最大值
% AllPTmax 与分割点对应的全部t检验序列对应的统计显著性
%% 第一步：变量初始化
N=length(X);
FLAG=zeros(N,1);
FLAG(1)=0.1;
FLAG(N)=0.1;
AllT=cell(0,0);
AllTmax=cell(0,0);
AllPTmax=cell(0,0);
%% 第二步：产生第一个突变点，并对序列进行分割
[T,Tmax,p,PTmax]=Tseries(X);
T=[1;T];
if PTmax<P0
flag=FLAG;
flag(1)=0;
flag(N)=0;
pos3=flag(pos2);
return
end
%记录输出数据
FLAG(p)=1;
AllT=[AllT;T];
AllTmax=[AllTmax;Tmax];
AllPTmax=[AllPTmax;PTmax];
%以下为两个控制计数器
counter=2;%下一个突变点的序号
TC=0;%临时计数器
%%
while 1%设置死循环
%% 第三步：对每个段进行突变检测，能分割则分割，直到不能分割为止
pos=find(FLAG>0);
M=length(pos)-1;%当前子段数目
for m=1:M
s=pos(m);
t=pos(m+1);
L=length(SubX);
if L>=L0
[T,Tmax,p,PTmax]=Tseries(SubX);
T=[s;T];
if PTmax>=P0
TC=TC+1;
FLAG(s+p-1)=counter;
AllT=[AllT;T];
AllPTmax=[AllPTmax;PTmax];
counter=counter+1;
end
end
end
%% 第四步：返回输出数据
if TC==0
flag=FLAG;
flag(1)=0;
flag(N)=0;
pos3=flag(pos2);
FLAG=[pos2,pos3];
return
end
%%
TC=0;
%%
end
function [T,Tmax,p,PTmax]=Tseries(x)
%% 计算t检验统计序列的子函数
% GreenSim团队原创作品，转载请注明
% Email:[email protected]
% GreenSim团队主页：http://blog.sina.com.cn/greensim
% 欢迎访问GreenSim——算法仿真团队→http://blog.sina.com.cn/greensim
%% 参数列表
% x 时间序列，N×1列向量
% T t检验序列，N×1列向量
% Tmax t检验序列的最大值
% p t检验序列最大值对应的下标
% PTmax Tmax对应的统计显著性
%% 参数初始化
N=length(x);
T=zeros(N,1);
%% 以下是主循环，用于创建t检验序列
for i=3:(N-2)%最左边以及最右边的两个点没有对应的t检验值（或者说，其值初始化为0）
x1=x(1:i);%序列左边部分
N1=length(x1);%左边序列的长度
x2=x(i:N);%序列右边部分
N2=length(x2);%右边序列的长度
mean_x2=mean(x2);%右边部分的均值
std_x2=std(x2);%右边部分的标准差
%下面是计算合并偏差的公式，中英文文献里的这个公式略有不同，此处以英文文献为准
SD=sqrt(1/N1+1/N2)*sqrt(((N1-1)*std_x1^2+(N2-1)*std_x2^2)/(N1+N2-2));
T(i)=abs((mean_x1-mean_x2)/SD);
end
%% 计算其它三个输出参数
Tmax=max(T);%t检验序列的最大值
pos=find(T==Tmax);
Eta=4.19*log(N)-11.54;%计算PTmax用的参数
Delta=0.40;%计算PTmax用的参数
v=N-2;%计算PTmax用的参数
c=v/(v+Tmax^2);%不完全B函数的下标
PTmax=(1-betainc(c,Delta*Eta,Delta));%调用不完全beta函数