题目:时间序列的贝叶斯突变检测算法MATLAB源码
突变分为如下主要的几种:均值突变(最常见)、方差突变、线性回归突变(也称趋势突变)、概率突变、空间型突变、谱突变、模型参数突变,等等。贝叶斯突变检测属于概率突变检测方法,其特点是能给出突变点的概率分布图。以下源码实现了贝叶斯突变检测,提一点经验性的东西,就是在做累乘计算的时候,为防止CPU溢出,应该使用对数把乘法转化为加法,最后再使用指数计算把结果返回。本源码由GreenSim团队原创,转载请注明,有意购买源码或代写相关程序,请与GreenSim团队联系(主页http://blog.sina.com.cn/greensim)。
function P=BMDCP(Y,R,M)
%% Bayesian Method for Detecting Change Points
%% 输入参数列表
%
Y 时间序列原始数据
%
R 控制参数Var0/Var,Var为全序列方差,Var0为均值方差,R应大于等于4
%
M 端点间隔,取值3~5
%% 输入参数列表
%
P 突变点的概率分布
%% 第一步:数据归一化处理
%归一化处理的目的在于避免浮点计算时溢出,属程序设计时考虑的问题,它不影响结果
Y=Y./max(Y);
%% 第二步:计算后验分布的数字特征
N=length(Y);%序列的长度
Mu0=mean(Y);%全序列的均值
Var=var(Y);%全序列方差
MuA=zeros(N,1);
MuB=zeros(N,1);
VarB=zeros(N,1);
for i=M:(N-M)
YA=Y(1:i);%左子序列
YB=Y((i+1):N);%右子序列
VarA(i)=Var/(1/R+i);
VarB(i)=Var/(1/R+N-i);
end
MuA(1:(M-1))=MuA(M);
MuB(1:(M-1))=MuB(M);
MuA((N-M+1):N)=MuA(N-M);
MuB((N-M+1):N)=MuB(N-M);
%% 第三步:计算后验概率分布
P=zeros(N,1);
for k=M:(N-M)
A=0;
for
i=1:k
% 注意:这里对似然函数进行了取对数处理,也是为了防止CPU溢出
a=log((1/sqrt(2*pi*Var))*exp(-(Y(i)-MuA(i))^2/(2*Var)));
end
B=0;
for
j=(k+1):N
b=log((1/sqrt(2*pi*Var))*exp(-(Y(j)-MuB(j))^2/(2*Var)));
end
P(k)=A+B;
end
%% 第四步:计算突变点的概率分布图
P=exp(P);
P=P/sum(P);
maxP=max(P);
bar(1:N,P);
axis([1,N,0,1.2*maxP]);
hold on
p0=(1/N)*ones(N,1);
plot(1:N,p0,'--r');
legend('后验概率','先验概率');
xlabel('k (month)','FontName','TimesNewRoman','FontSize',12);
ylabel('P(k)','FontName','TimesNewRoman','Fontsize',12);