visualise and select subset

F2_cpsingle.m

+function [icp, residue] = F2_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-2017 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)
+m = size(y,1);
+n = size(y,2);
+fwd = zeros(1,n,'like',y);
+rev = zeros(1,n,'like',y);
+
+ymean = zeros(m,1,'like',y);
+Syy = zeros(m,1,'like',y);
+for ix=1:n
+  ydelta = y(:,ix) - ymean;
+  npoints = ix;
+  ymean = ymean + ydelta./npoints;
+  Syy = Syy + ydelta .* (y(:,ix) - ymean);
+  fwd(ix) = sum(Syy);
+end
+
+ymean = zeros(m,1,'like',y);
+Syy = zeros(m,1,'like',y);
+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) = sum(Syy);
+end
+
+% ------------------------------------------------------------------------
+function [fwd,rev] = rmsResidual(y)
+m = size(y,1);
+n = size(y,2);
+fwd = zeros(1,n,'like',y);
+rev = zeros(1,n,'like',y);
+
+logrealmin = repmat(log(realmin),m,1);
+
+Syy = zeros(m,1,'like',y);
+for ix=1:n
+  npoints = ix;
+  Syy = Syy + y(:,ix).^2;
+  fwd(ix) = npoints*sum(max(logrealmin-log(npoints),log(Syy./npoints)));
+end
+
+Syy = zeros(m,1,'like',y);
+for ix=n:-1:1
+  npoints = n+1-ix;
+  Syy = Syy + y(:,ix).^2;
+  rev(ix) = npoints*sum(max(logrealmin-log(npoints),log(Syy./npoints)));
+end
+
+% ------------------------------------------------------------------------
+function [fwd,rev] = stdResidual(y)
+m = size(y,1);
+n = size(y,2);
+fwd = zeros(1,n);
+rev = zeros(1,n);
+
+logrealmin = repmat(log(realmin),m,1);
+
+ymean = zeros(m,1,'like',y);
+Syy = zeros(m,1,'like',y);
+for ix=1:n
+  npoints = ix;
+  ydelta = y(:,ix) - ymean;
+  ymean = ymean + ydelta./npoints;
+  Syy = Syy + ydelta .* (y(:,ix) - ymean);
+  fwd(ix) = npoints*sum(max(logrealmin-log(npoints),log(Syy./npoints)));
+end
+
+ymean = zeros(m,1,'like',y);
+Syy = zeros(m,1,'like',y);
+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*sum(max(logrealmin-log(npoints),log(Syy./npoints)));
+end
+
+% ------------------------------------------------------------------------
+function [fwd,rev] = linearResidual(y)
+m = size(y,1);
+n = size(y,2);
+fwd = zeros(1,n,'like',y);
+rev = zeros(1,n,'like',y);
+
+xmean = zeros(m,1,'like',y);
+ymean = zeros(m,1,'like',y);
+Sxx = zeros(m,1,'like',y);
+Syy = zeros(m,1,'like',y);
+Sxy = zeros(m,1,'like',y);
+SxxSSE = zeros(m,1,'like',y);
+
+for ix=1:n
+  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) = sum(SxxSSE ./ Sxx);
+end
+
+xmean = zeros(m,1,'like',y);
+ymean = zeros(m,1,'like',y);
+Sxx = zeros(m,1,'like',y);
+Syy = zeros(m,1,'like',y);
+Sxy = zeros(m,1,'like',y);
+SxxSSE = zeros(m,1,'like',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;
+  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;
+  rev(ix) = sum(SxxSSE ./ Sxx);
+end

S10_penalty_postprocess.asv

@@ -132,6 +132,8 @@ for i = 1:length(fields_)
     
 end
     
+clear ans days endb field_ i randindx start selday
+
 % plot both time series
 for i = 1:length(fields_)
     
@@ -156,13 +158,31 @@ for i = 1:length(fields_)
 end
 
 % summarize and calculate the metrics of the CP analysis
-%    for "mean" ==> residue = n*sum(var(x,1,2))
+%    for "mean" ==> residue = n*sum(var(x,1,2));
 %    for "std" ==> residue = sum(n*log(var(x,1,2)));
 %    for "rms" ==> residue = sum(n*log(sum(x.^2,2)/n));
 % calculate residu when only one CP
 %    
 
-cpsums = table(fields_,'VariableNames','cows');
+cpsums = table((1:length(fields_))','VariableNames',{'No'});
+for i = 1:height(cpsums)
+    % cowid
+    cpsums.cow(i,:) = seldata.(fields_{i}).object_name(1);
+    
+    % residus for each of the stats - Z data
+    z = seldata.(fields_{i}).avg_z_sm(~isnan(seldata.(fields_{i}).avg_z_sm))';
+    p = seldata.(fields_{i}).numtime(~isnan(seldata.(fields_{i}).avg_z_sm))';
+    n = length(z);
+    cpsums.n(i) = n;
+    cpsums.ZmeanR(i) = n*sum(var(z,1,2));
+    cpsums.ZstdR(i) = sum(n*log(var(z,1,2)));
+    cpsums.ZrmsR(i) = sum(n*log(sum(z.^2,2)/n));
+    
+    % residus for each of the stats with one CP - Z data
+    [icp,residu] = F2_cpsingle(z,'mean',300)
+    
+end
+
 
 
 

S10_penalty_postprocess.m

@@ -132,12 +132,15 @@ for i = 1:length(fields_)
     
 end
     
+clear ans days endb field_ i randindx start selday
+
 % plot both time series
+close all
 for i = 1:length(fields_)
     
     % figure;
     figure('Units','centimeters','OuterPosition',[1 1 28 20]);
-    subplot(2,1,1); hold on;box on; title('Z position'); 
+    subplot(2,1,1); hold on;box on; title([fields_{i} ', Z position']); 
        xlabel('day [h]'); ylabel('Z-position [m]')
        plot(seldata.(fields_{i}).numtime-min(seldata.(fields_{i}).numtime),...
             seldata.(fields_{i}).avg_z_sm,'Color',[178/255 34/255 34/255],...
@@ -156,16 +159,56 @@ for i = 1:length(fields_)
 end
 
 % summarize and calculate the metrics of the CP analysis
-%    for "mean" ==> residue = n*sum(var(x,1,2))
+%    for "mean" ==> residue = n*sum(var(x,1,2));
 %    for "std" ==> residue = sum(n*log(var(x,1,2)));
 %    for "rms" ==> residue = sum(n*log(sum(x.^2,2)/n));
 % calculate residu when only one CP
 %    
 
-cpsums = table(fields_,'VariableNames','cows');
-
+cpsums = table((1:length(fields_))','VariableNames',{'No'});
+for i = 1:height(cpsums)
+    % cowid
+    cpsums.cow(i,:) = seldata.(fields_{i}).object_name(1);
+    
+    % residus for each of the stats - Z data
+    z = seldata.(fields_{i}).avg_z_sm(~isnan(seldata.(fields_{i}).avg_z_sm))';
+    p = seldata.(fields_{i}).numtime(~isnan(seldata.(fields_{i}).avg_z_sm))-...
+            min(seldata.(fields_{i}).numtime);
+    n = length(z);
+    cpsums.n(i) = n;
+    cpsums.ZmeanR(i) = n*sum(var(z,1,2));
+    cpsums.ZstdR(i) = sum(n*log(var(z,1,2)));
+    cpsums.ZrmsR(i) = sum(n*log(sum(z.^2,2)/n));
+    
+    % residus for each of the stats with one CP - Z data
+    [~,cpsums.ZmeanR1(i)] = F2_cpsingle(z,'mean',300);
+    [~,cpsums.ZstdR1(i)] = F2_cpsingle(z,'std',300);
+    [~,cpsums.ZrmsR1(i)] = F2_cpsingle(z,'rms',300);
+    
+    %    
+    cdd = seldata.(fields_{i}).centerdist(~isnan(seldata.(fields_{i}).centerdist))';
+    p = seldata.(fields_{i}).numtime(~isnan(seldata.(fields_{i}).centerdist))-...
+            min(seldata.(fields_{i}).numtime);
+    n = length(cdd);
+    cpsums.n(i) = n;
+    cpsums.CDmeanR(i) = n*sum(var(cdd,1,2));
+    cpsums.CDstdR(i) = sum(n*log(var(cdd,1,2)));
+    cpsums.CDrmsR(i) = sum(n*log(sum(cdd.^2,2)/n));
+    
+    % residus for each of the stats with one CP - Z data
+    [~,cpsums.CDmeanR1(i)] = F2_cpsingle(cdd,'mean',300);
+    [~,cpsums.CDstdR1(i)] = F2_cpsingle(cdd,'std',300);
+    [~,cpsums.CDrmsR1(i)] = F2_cpsingle(cdd,'rms',300);
 
+end
 
+% difference between zero and one cp
+cpsums.ZmeanDif = cpsums.ZmeanR-cpsums.ZmeanR1;
+cpsums.ZstdDif = cpsums.ZstdR-cpsums.ZstdR1;
+cpsums.ZrmsDif = cpsums.ZrmsR-cpsums.ZrmsR1;
+cpsums.CDmeanDif = cpsums.CDmeanR-cpsums.CDmeanR1;
+cpsums.CDstdDif = cpsums.CDstdR-cpsums.CDstdR1;
+cpsums.CDrmsDif = cpsums.CDrmsR-cpsums.CDrmsR1;