BG算法

xqtn

matlab小白一个，论文第一次接触，BG算法里定义变量x,y,P0,L0具体要怎么写，x为时间，y为数据列，P0=0.95，L0=6

xqtn

程序如下
function [ FLAGS ] = BGSA( x,y,P0,L0 )
% [ FLAGS ] = BGSA( x,y,P0,L0 )
%  x: 序列的x坐标（仅用于绘图，如果不使用则设置为[]）
%  y: 降雨量序列
%  P0: 显著性水平门限值，低于此值的不再分割
%  L0: 最小分割尺度，子段长度小于此值的不再分割
%  FLAGS: 返回的逻辑向量与y大小相同，值为1表示相应的位置定位分割点
% 通过Bernaola-Galvan分割算法（BGSA）找出并绘制突变点及其最大密度间隔。
%   如何解释这个数字？
%   蓝色细曲线是原始序列; 而粗一个表示每个段中原始序列的间隔平均值t; 红色垂直线表示分割点所在的位置;阴影指示最大突变密度间隔。\
n = length(y);
if isempty(x)
    x = 1:n;
end
FLAGS = zeros(1,n);
FLAGS(n) = 1;

% 通过标记1在FLAGS中找到分割点的位置。
num_sp = sum(FLAGS);
while 1
    I = find(FLAGS==1);
    numsp_bf = num_sp;
    for k = 1:num_sp
        if k==1
            subx = x(1:I(k));
            suby = y(1:I(k));
        else
            subx = x(I(k-1)+1:I(k));
            suby = y(I(k-1)+1:I(k));
        end
        if length(subx)>L0
            [ Xpos,PTmax ] = calc_TP( subx,suby );
            if PTmax>=P0
                FLAGS(x==min(Xpos)) = 1;
                num_sp = num_sp+1;
            end
        end
    end
    if num_sp==numsp_bf
        break
    end
end

% 计算每个分割的平均序列。
for k = 1:num_sp
    if k==1
        my(1:I(k)) = mean(y(1:I(k)));
    else
        my(I(k-1)+1:I(k)) = mean(y(I(k-1)+1:I(k)));
    end
end

% 调整一些不合理的突变点。
FLAGS(n) = 0;
for j = 2:n
    if FLAGS(j)==1 && FLAGS(j-1)==1
        FLAGS(j) = 0;
    end
end

%return
figure

% 计算并找出最大突变密度的间隔。
NDT = L0;
for j = 1:n-NDT
    n = sum(FLAGS(j:j+NDT));
    eta(j) = n/NDT;
end
Ieta = find((eta==max(eta))+(eta~=0)==2);

% 影响间隔
ya = [(max(y)+1)*1.1 (max(y)+1)*1.1];
bv = (min(y)-1)*1.1;
xamin = x(Ieta);
xamax = x(Ieta+NDT);
if ~isempty(Ieta)
xa(1,:) = [xamin(1) xamax(1)];
kk = 1;
for l = 2:length(Ieta)
    if xamax(l-1)>=xamin(l)
        xa(kk,:) = [min(xa(kk,:)) xamax(l)];
    else
        kk = kk+1;
        xa(kk,:) = [xamin(l) xamax(l)];
    end
end
for l = 1:kk
    area(xa(l,:),ya,bv,'LineStyle','none','FaceColor',[.6 .6 .6],'FaceAlpha',.6,'ShowBaseLine','off');
    hold on
end
end

% 在分割点绘制曲线和垂直线。
plot(x,y,'color',[0 .447 .741]);
hold on
plot(x,my,'color',[0 .447 .741],'linewidth',2);
hold on
I = find(FLAGS==1);
for k = 1:sum(FLAGS)
    plot([x(I(k)) x(I(k))],[bv max(ya)],'r-');
    hold on
end
axis([min(x) max(x) (min(y)-1)*1.05 (max(y)+1)*1.05]);
title('Bernaola-Galvan Segmentation Algorithm')
end

function [ Xpos,PTmax ] = calc_TP( X,Y )
% 子函数计算统计量T，并找出与X系列对应的位置（Xpos）的Tmax，然后计算显着性检验的统计PTmax。
N = length(Y);
T = zeros(1,N);
for i=2:N-1
    nl = length(Y(1:i));
    nr = length(Y(i:N));
    ml = mean(Y(1:i));
    mr = mean(Y(i:N));
    vl = var(Y(1:i));
    vr = var(Y(i:N));
    T(i) = abs((ml-mr)/sqrt(1/nl+1/nr)/sqrt(((nl-1)*vl+(nr-1)*vr)/...
        (nl+nr-2)));
end
Tmax = max(T);
Xpos = X(T==Tmax);
gamma = 4.19*log(N)-11.54;
delta = .4;
v = N-2;
PTmax = (1-betainc(v/(v+Tmax^2),delta*v,delta))^gamma;
end

名字是在不好整

我去试试BG，自己琢磨好久，老是有问题

名字是在不好整

请问下 楼主做出来了吗 能运行吗？

kingjames

请问 楼主你捉出来了吗？

世界很大啊

楼主您好，想麻烦问一下，这个L0该怎么写呢，不能运行