www.gusucode.com > signal 工具箱matlab源码程序 > signal/private/cpsingle.m

    function [icp, residue] = cpsingle(y, statistic, Lmin)
%CPSINGLE single best changepoint between two equal distributions
%   This file is for internal use only and may be removed in a future
%   release.
%
%   References:
%   * Tony F. Chan, Gene H. Golub, Randall J. LeVeque, Algorithms for
%     computing the sample variance: analysis and recommendations"  The
%     American Statistician.  Vol 37, No. 3 (Aug., 1983) pp. 242-247.
%   * Philippe Pierre Pebay. "Formulas for robust one-pass parallel
%     computation of covariances and arbitrary-order statistical moments."
%     SAND2008-6212.  2008. Published through SciTech Connect.

%   Copyright 2015 The MathWorks, Inc.


% farm out to correct algorithm
if strcmp(statistic,'mean')
  [fwd, rev] = meanResidual(y);
elseif strcmp(statistic,'rms')
  [fwd, rev] = rmsResidual(y);
elseif strcmp(statistic,'std')
  [fwd, rev] = stdResidual(y);
elseif strcmp(statistic,'linear')
  [fwd, rev] = linearResidual(y);
end

[residue, icp] = min(fwd(Lmin:end-Lmin)+rev(Lmin+1:end-Lmin+1));
icp = icp + Lmin;

% ------------------------------------------------------------------------
function [fwd,rev] = meanResidual(y)
fwd = zeros(size(y));
rev = zeros(size(y));
n = numel(y);

ymean = 0;
Syy = 0;
for ix=1:n
  ydelta = y(ix) - ymean;
  npoints = ix;
  ymean = ymean + ydelta./npoints;
  Syy = Syy + ydelta * (y(ix) - ymean);
  fwd(ix) = Syy;
end

ymean = 0;
Syy = 0;
for ix=n:-1:1
  ydelta = y(ix) - ymean;
  npoints = n-ix+1;
  ymean = ymean + ydelta./npoints;
  Syy = Syy + ydelta * (y(ix) - ymean);
  rev(ix) = Syy;
end

% ------------------------------------------------------------------------
function [fwd,rev] = rmsResidual(y)
fwd = zeros(size(y));
rev = zeros(size(y));
n = numel(y);

logrealmin = log(realmin);

Syy = 0;
for ix=1:n
  npoints = ix;
  Syy = Syy + y(ix).^2;
  fwd(ix) = npoints*max(logrealmin,log(Syy./npoints));
end

Syy = 0;
for ix=n:-1:1
  npoints = n+1-ix;
  Syy = Syy + y(ix).^2;
  rev(ix) = npoints*max(logrealmin,log(Syy./npoints));
end

% ------------------------------------------------------------------------
function [fwd,rev] = stdResidual(y)
fwd = zeros(size(y));
rev = zeros(size(y));
n = numel(y);

logrealmin = log(realmin);

ymean = 0;
Syy = 0;
for ix=1:n
  npoints = ix;
  ydelta = y(ix) - ymean;
  ymean = ymean + ydelta./npoints;
  Syy = Syy + ydelta * (y(ix) - ymean);
  fwd(ix) = npoints*max(logrealmin,log(Syy./npoints));
end

ymean = 0;
Syy = 0;
for ix=n:-1:1
  ydelta = y(ix) - ymean;
  npoints = n+1-ix;
  ymean = ymean + ydelta./npoints;
  Syy = Syy + ydelta * (y(ix) - ymean);
  rev(ix) = npoints*max(logrealmin,log(Syy./npoints));
end

% ------------------------------------------------------------------------
function [fwd,rev] = linearResidual(y)
fwd = zeros(size(y));
rev = zeros(size(y));

xmean = 0;
ymean = 0;
Sxx = 0;
Syy = 0;
Sxy = 0;
SxxSSE = 0;

for ix=1:numel(y)
  npoints = ix;
  ydelta = y(ix) - ymean;
  xdelta = ix - xmean;
  ymean = ymean + ydelta/npoints;
  xmean = xmean + xdelta/npoints;
  dSyy = ydelta .* (y(ix) - ymean);
  dSxx = xdelta .* (ix - xmean);
  dSxy = xdelta .* ydelta .* (npoints - 1) ./ npoints;
  Syy = Syy + dSyy;
  dSxxSSE = dSxx.*Syy + dSyy.*Sxx - dSxy.*(2*Sxy+dSxy);
  Sxx = Sxx + dSxx;
  Sxy = Sxy + dSxy;
  SxxSSE = SxxSSE + dSxxSSE;
  fwd(ix) = SxxSSE ./ Sxx;
end

xmean = 0;
ymean = 0;
Sxx = 0;
Syy = 0;
Sxy = 0;
n = numel(y);
for ix=n:-1:1
  npoints = n+1-ix;
  ydelta = y(ix) - ymean;
  xdelta = ix - xmean;
  ymean = ymean + ydelta/npoints;
  xmean = xmean + xdelta/npoints;
  Syy = Syy + ydelta * (y(ix) - ymean);
  Sxx = Sxx + xdelta * (ix - xmean);
  Sxy = Sxy + xdelta .* ydelta .* (npoints - 1) ./ npoints;
  rev(ix) = Syy - Sxy.^2 ./ Sxx;
end