Skip to content
Snippets Groups Projects
Commit fffd8cd1 authored by Adriaens, Ines's avatar Adriaens, Ines :v_tone2:
Browse files

madd testscript MinDifference in improvemen

parent 67d2bc32
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
manuscript/20210908_Adriaens_UWB_.log
+ 29
− 22
View file @ fffd8cd1
This is pdfTeX, Version 3.141592653-2.6-1.40.22 (MiKTeX 21.2) (preloaded format=pdflatex 2021.5.7) 8 SEP 2021 15:01
This is pdfTeX, Version 3.141592653-2.6-1.40.22 (MiKTeX 21.2) (preloaded format=pdflatex 2021.5.7) 9 SEP 2021 13:25
entering extended mode
**./20210908_Adriaens_UWB_.tex
(20210908_Adriaens_UWB_.tex
......@@ -236,35 +236,42 @@ LaTeX Font Info: Trying to load font information for U+msb on input line 30.
(C:\Users\adria036\AppData\Local\Programs\MiKTeX\tex/latex/amsfonts\umsb.fd
File: umsb.fd 2013/01/14 v3.01 AMS symbols B
)
[2] [3] [4] [5] [6] (20210908_Adriaens_UWB_.aux) )
[2]
Overfull \hbox (25.41856pt too wide) detected at line 83
\OML/cmm/m/it/10 var\OT1/cmr/m/n/10 ([\OML/cmm/m/it/10 x[] [] x[]\OT1/cmr/m/n/1
0 ]) = [] [] []|
[]
[3] [4] [5] [6] (20210908_Adriaens_UWB_.aux) )
Here is how much of TeX's memory you used:
3818 strings out of 479334
50800 string characters out of 2855660
361535 words of memory out of 3000000
360535 words of memory out of 3000000
21231 multiletter control sequences out of 15000+200000
416185 words of font info for 62 fonts, out of 3000000 for 9000
1141 hyphenation exceptions out of 8191
68i,6n,74p,2272b,284s stack positions out of 5000i,500n,10000p,200000b,50000s
<C:\Users\adria036\AppData\L
ocal\MiKTeX\fonts/pk/ljfour/jknappen/ec/dpi600\ecti1000.pk> <C:\Users\adria036\
AppData\Local\MiKTeX\fonts/pk/ljfour/jknappen/ec/dpi600\ecbx1000.pk> <C:\Users\
adria036\AppData\Local\MiKTeX\fonts/pk/ljfour/jknappen/ec/dpi600\tcrm1000.pk> <
C:\Users\adria036\AppData\Local\MiKTeX\fonts/pk/ljfour/jknappen/ec/dpi600\ecbx1
440.pk> <C:\Users\adria036\AppData\Local\MiKTeX\fonts/pk/ljfour/jknappen/ec/dpi
600\ecbx1200.pk> <C:\Users\adria036\AppData\Local\MiKTeX\fonts/pk/ljfour/jknapp
en/ec/dpi600\ecrm1200.pk> <C:\Users\adria036\AppData\Local\MiKTeX\fonts/pk/ljfo
ur/jknappen/ec/dpi600\ecrm1728.pk> <C:\Users\adria036\AppData\Local\MiKTeX\font
s/pk/ljfour/jknappen/ec/dpi600\ecrm1000.pk><C:/Users/adria036/AppData/Local/Pro
grams/MiKTeX/fonts/type1/public/amsfonts/cm/cmex10.pfb><C:/Users/adria036/AppDa
ta/Local/Programs/MiKTeX/fonts/type1/public/amsfonts/cm/cmmi10.pfb><C:/Users/ad
ria036/AppData/Local/Programs/MiKTeX/fonts/type1/public/amsfonts/cm/cmmi7.pfb><
C:/Users/adria036/AppData/Local/Programs/MiKTeX/fonts/type1/public/amsfonts/cm/
cmr10.pfb><C:/Users/adria036/AppData/Local/Programs/MiKTeX/fonts/type1/public/a
msfonts/cm/cmr7.pfb><C:/Users/adria036/AppData/Local/Programs/MiKTeX/fonts/type
1/public/amsfonts/cm/cmsy10.pfb>
Output written on 20210908_Adriaens_UWB_.pdf (6 pages, 121055 bytes).
<C:\Users\adria036\AppData\Local
\MiKTeX\fonts/pk/ljfour/jknappen/ec/dpi600\ecti1000.pk> <C:\Users\adria036\AppD
ata\Local\MiKTeX\fonts/pk/ljfour/jknappen/ec/dpi600\ecbx1000.pk> <C:\Users\adri
a036\AppData\Local\MiKTeX\fonts/pk/ljfour/jknappen/ec/dpi600\tcrm1000.pk> <C:\U
sers\adria036\AppData\Local\MiKTeX\fonts/pk/ljfour/jknappen/ec/dpi600\ecbx1440.
pk> <C:\Users\adria036\AppData\Local\MiKTeX\fonts/pk/ljfour/jknappen/ec/dpi600\
ecbx1200.pk> <C:\Users\adria036\AppData\Local\MiKTeX\fonts/pk/ljfour/jknappen/e
c/dpi600\ecrm1200.pk> <C:\Users\adria036\AppData\Local\MiKTeX\fonts/pk/ljfour/j
knappen/ec/dpi600\ecrm1728.pk> <C:\Users\adria036\AppData\Local\MiKTeX\fonts/pk
/ljfour/jknappen/ec/dpi600\ecrm1000.pk><C:/Users/adria036/AppData/Local/Program
s/MiKTeX/fonts/type1/public/amsfonts/cm/cmex10.pfb><C:/Users/adria036/AppData/L
ocal/Programs/MiKTeX/fonts/type1/public/amsfonts/cm/cmmi10.pfb><C:/Users/adria0
36/AppData/Local/Programs/MiKTeX/fonts/type1/public/amsfonts/cm/cmmi7.pfb><C:/U
sers/adria036/AppData/Local/Programs/MiKTeX/fonts/type1/public/amsfonts/cm/cmr1
0.pfb><C:/Users/adria036/AppData/Local/Programs/MiKTeX/fonts/type1/public/amsfo
nts/cm/cmr7.pfb><C:/Users/adria036/AppData/Local/Programs/MiKTeX/fonts/type1/pu
blic/amsfonts/cm/cmsy10.pfb><C:/Users/adria036/AppData/Local/Programs/MiKTeX/fo
nts/type1/public/amsfonts/cm/cmsy7.pfb>
Output written on 20210908_Adriaens_UWB_.pdf (6 pages, 134688 bytes).
PDF statistics:
260 PDF objects out of 1000 (max. 8388607)
273 PDF objects out of 1000 (max. 8388607)
0 named destinations out of 1000 (max. 500000)
1 words of extra memory for PDF output out of 10000 (max. 10000000)
manuscript/20210908_Adriaens_UWB_.pdf
+ 0
− 0
View file @ fffd8cd1
No preview for this file type
manuscript/20210908_Adriaens_UWB_.synctex.gz
+ 0
− 0
View file @ fffd8cd1
No preview for this file type
manuscript/20210908_Adriaens_UWB_.tex
+ 18
− 5
View file @ fffd8cd1
......@@ -75,31 +75,44 @@
\subsubsection{Data selection}
\subsection{Changepoint analysis}
Changepoints are time instants or samples in which the statistical properties (i.e. distribution) of a (time) series change abruptly. In this study, we detected and combined the individual changepoints per cow per day in two time series of (x,y,z)-coordinate positioning data. Intuitively, one could argue to mainly rely on the position in the vertical (z) direction (height), as a cow that lies down is expected to remain in a lower and more stable position compared to when she is up. However, the z-position has been shown to be the most unreliable and noisy (range, variability,..) of all three coordinates, with its inaccuracy variable in time and space, and depending on e.g. the position in the barn, the behaviour and speed of the animals, the collar attachment and their interactions, etc. Similarly, relying on detection of a stable position in the (x,y)-direction (which is unmistakably true during lying bouts) is imprecise as well, as cow activity varies over the day, and oftentimes animals stand still for a longer period of time apart from their lying bouts, for example when grooming other animals, feeding, drinking or ruminating. These periods of “standing” inactivity might additionally depend on accessibility of cubicles and lying places, hierarchy, climate of the barn etc. In this study, we chose to work on a combination of two position derived time series. The first is the z-coordinate (height) of the animals, as this is the most straightforward one. The second time series is the ‘centre distance’ (CD), i.e. the position relative to the centre of the barn. The main advantage of using CD and not the raw (x,y)-position is that it summarizes position and movement of the animals in a single signal, less dependent on the actual direction of movement, and with a lower variability. Should a cow move in a perfect circle around the centre of the barn, however, CD remains constant (as is the case when a cow stands still or lies down). We assumed that this would be extremely rare, and when it would happen for a short period of time, this would not impair the analysis because movement as such causes the signal to be more variable, which also changes the statistical properties of the time series.
Changepoints are time instants or samples in which the statistical properties (i.e. distribution) of a (time) series change abruptly. In this study, we detected and combined the individual changepoints per cow per day in two time series of (x,y,z)-coordinate positioning data. Intuitively, one could argue to mainly rely on the position in the vertical (z) direction (height), as a cow that lies down is expected to remain in a lower and more stable position compared to when she is up. However, the z-position has been shown to be the most unreliable and noisy (range, variability,..) of all three coordinates, with its inaccuracy variable in time and space, and depending on e.g. the position in the barn, the behaviour and speed of the animals, the collar attachment and their interactions, etc. Similarly, relying on detection of a stable position in the (x,y)-direction (which is unmistakably true during lying bouts) is imprecise as well, as cow activity varies over the day, and oftentimes animals stand still for a longer period of time apart from their lying bouts, for example when grooming other animals, feeding, drinking or ruminating. These periods of “standing” inactivity might additionally depend on accessibility of cubicles and lying places, hierarchy, climate of the barn etc. In this study, we chose to work on a combination of two position derived time series. The first is the z-coordinate (height) of the animals, as this is the most straightforward one. The second time series is the ‘center distance’ (CD), i.e. the position relative to the center of the barn. The main advantage of using CD and not the raw (x,y)-position is that it summarizes position and movement of the animals in a single signal, less dependent on the actual direction of movement, and with a lower variability. Should a cow move in a perfect circle around the center of the barn, however, CD remains constant (as is the case when a cow stands still or lies down). We assumed that this would be extremely rare, and when it would happen for a short period of time, this would not impair the analysis because movement as such causes the signal to be more variable, which also changes the statistical properties of the time series.
The changepoint analysis relies on a parametric method that separates the time series in segments and computes either the level (average value, z-coordinate), or level and variability (mean and variance) of each segment from point \textit{m} to point \textit{n} as follows:
\[\dfrac{1}{n-m+1}\sum_{r=m}^{n}x_{r}\]
\[\dfrac{1}{n-m+1}\sum_{r=m}^{n}{(x_{r}-mean([x_{m}...x_{n}])^2}\]
\[mean([x_{m} \cdots x_{n}]) = \dfrac{1}{n-m+1}\sum_{r=m}^{n}x_{r}\]
\[var([x_{m} \cdots x_{n}]) = \dfrac{1}{n-m+1}\sum_{r=m}^{n}{(x_{r}-mean([x_{m} \cdots x_{n}]))^2(y_{r}-mean([y_{m} \cdots y_{n}]))^2}\]
%TODO = generalized equations with log likelihood
Based on the minimization of the sums of squares of the residual error (i.e. the maximum log likelihood), it finds the number of changepoints \textit{k} for which:
\[J = \sum_{i=1}^{k-1}{(x_{i}-mean([x_{1} \cdots x_{k-1}]))^2} + \sum_{i=k}^{N}{(x_{i}-mean([x_{k} \cdots x_{N}]))^2} \]
is smallest.
Changepoint detection in the \textit{z}-data relied on the level (i.e. the mean) alone, while for the CD, both mean and variance were taken into consideration. Based on the
\color{r}{TODO: rules to interpret the direction of the change in terms of level -- if higher --- higher : no change}
Is smallest.
Additionally, to deal with outliers in the CD time series when the variability was small, a penalty was implemented for detection of changepoints where the absolute difference in level and variance of the CD segment was below 0.5m. A maximum number of changepoints was set to aGiven the number of
\subsection{Statistical analysis}
%----------------------------------------------- results -------------------------------------------%
\pagebreak
\section{Results \& discussion}
\begin{itemize}
\item expert knowledge for settings of the CP analysis
\end{itemize}
......
test_chptsmin.m 0 → 100644
+ 464
− 0
View file @ fffd8cd1
%% S8_cpnts
% visualisation of the lying bouts + changepoints for development purposes
clear variables
close all
clc
% results directory
init_.resdir = ['C:\Users\adria036\OneDrive - Wageningen '...
'University & Research\iAdriaens_doc\Projects\' ...
'eEllen\B4F_indTracking_cattle\B4F_results\'];
% load data
load([init_.resdir 'D3_sdata_IQ_UWB.mat'])
load([init_.resdir 'D2_data_meta.mat'])
%% set constants
% delete cow 280
sdata = rmfield(sdata,'cow_280');
fields_ = fieldnames(sdata);
for f = 1:length(fields_)
% set cow
cow = fields_{f};
% cow = 'cow_99';
% set technique
cpts.Zstat = 'mean'; % technique based on which cpts are found for z
cpts.CDstat = 'std'; % technique based on which cpts are found for CD
% gapsize for new segment for changepoint analysis
cpts.gapsize = 300; % 5 minutes
% segment length
cpts.segmentlength = 600; % 10 minutes
% set windowsize in days
cpts.window = 600/(60*60*24); % windowsize in which cpts are detected
% set settings
cpts.MinDistance = 600; % 10 min
% cpts.NumChanges = 1; % number of changes per hour of data
cpts.MinThreshold = 1000;
% figure settings
cpts.Units = 'centimeters'; % set units
cpts.figPos = [1 1 30 22]; % set position
%% select data and plot lying bouts
% select data of one day per cow
clear data
close all
% data contains "z" and "CD" and runs from first lying to last lying bout
data = sdata.(cow)(:,[2 3 7 13 17 18 15]);
% reference date for numeric time
refdate = datenum(2019,7,3);
% summary of the lying bouts
data.islying_prev(1) = 0;
data.islying_prev(2:end) = data.islying(1:end-1);
data.lychange = ((data.islying == 1 & data.islying_prev == 0) | ...
(data.islying == 0 & data.islying_prev == 1));
data_ld = data(data.lychange == 1,:);
data_ld.dateup(:,1) = NaT;
data_ld.dateup(1:2:end-1) = data_ld.date(2:2:end);
data_ld(isnat(data_ld.dateup),:) = [];
data_ld = movevars(data_ld,'dateup','After','date');
data_ld.numdatedown = datenum(data_ld.date) - refdate;
data_ld.numdateup = datenum(data_ld.dateup) - refdate;
% visualise the data and the lying bouts (data_meta.selection)
h = figure('Units',cpts.Units,...
'OuterPosition',cpts.figPos); % figure
subplot(2,1,1); hold on; box on;title('Z'); % prepare z plot
ylim([0 2.5])
subplot(2,1,2); hold on; box on;title('CD'); % prepare CD plot
ylim([0 12])
% plot lying bouts
for i = 1:2
subplot(2,1,i);
for j = 1:height(data_ld)
x = [datenum(data_ld.date)-refdate ...
datenum(data_ld.dateup)-refdate]; % numeric time
y = [0.1 0.1]; % y variable = thick line
plot(x,y,'LineWidth',4.5,'Color',[204/255 0 102/255])
x = [datenum(data_ld.date)-refdate ... % numeric time
datenum(data_ld.date)-refdate];
y = [0 15]; % vertical line begin
plot(x,y,'LineWidth',0.8,'Color',[204/255 0 102/255])
x = [datenum(data_ld.dateup)-refdate ...
datenum(data_ld.dateup)-refdate]; % vertical line end
plot(x,y,'LineWidth',0.8,'Color',[204/255 0 102/255])
end
end
% plot z data and implement changepoint analysis
% stap 1: determine data with gaps not too long and plot gap
% step 2: split data according to gaps
% step 3: determine max amount of CP based on data
% step 4: CP analysis
% step 5: visualise result
data2 = data(isnan(data.avg_z_sm)~= 1 & isnan(data.centerdist) ~= 1,:);
data2.numdate = datenum(data2.date)-refdate;
% plot z data
subplot(2,1,1); plot(data2.numdate,data2.avg_z_sm, ...
'LineWidth',0.9,'Color',[60/255 179/255 113/255]);
xlim([min(data2.numdate) max(data2.numdate)])
% plot cd data
subplot(2,1,2); plot(data2.numdate,data2.centerdist, ...
'LineWidth',0.9,'Color',[60/255 179/255 113/255]);
xlim([min(data2.numdate) max(data2.numdate)])
% calculate gapsize
data2.gapsize(end) = 0;
data2.gapsize(1:end-1) = diff(data2.numdate)*60*60*24;
% select data segments and characterize the segments
ind = find(data2.gapsize > cpts.gapsize);
segments = data2(ind,:); % segments
segments.segstart(1) = data2.numdate(1); % segment start
segments.segstart(2:end) = data2.numdate(ind(1:end-1)+1); % segment start
segments.dur = (segments.numdate - segments.segstart)*60*60*24; % dur in s
% find the segment where transitions happens
for i = 1:height(data_ld)
% find whether the lie down event is in a segment
ind = find(segments.segstart < data_ld.numdatedown(i) & ...
segments.numdate > data_ld.numdatedown(i));
if ~isempty(ind)
data_ld.segdown(i) = ind;
else
data_ld.segdown(i) = NaN;
end
% find whether the get up event is in a segment
ind = find(segments.segstart < data_ld.numdateup(i) & ...
segments.numdate > data_ld.numdateup(i));
if ~isempty(ind)
data_ld.segup(i) = ind;
else
data_ld.segup(i) = NaN;
end
end
% implement changepoint analysis on the segments that are long enough, and
% that have at least one getup/liedown event
% segment number
segments.number = (1:height(segments))';
% segments with an event get up or lie down
seglist = unique([data_ld.segdown;data_ld.segup]);
seglist(isnan(seglist)) = [];
% add to segments whether there is a get up or lie down event
segments.liedown = ismember(segments.number, seglist);
%%
%---------------------------------------------------------------------
% changepoint analysis - ATTEMPT 2 - take full day independent of segment
% length and a fixed number of max changepoints
% --> evaluate validity of detection only
% if there is a getup or liedown event
days = unique(floor(data2.numdate));
for i = 1:length(days)
% select daydata
data_day = data2(floor(data2.numdate) == days(i),:);
if height(data_day) > 12*60*60
% prepare results
day = sprintf('day_%d',days(i));
% prepare figure
k = figure('Units',cpts.Units,...
'OuterPosition',cpts.figPos); % figure
subplot(2,1,1); hold on; box on; % prepare z plot
title([cow ', day = ' num2str(days(i)) ', Z'], ...
'Interpreter','none');
ylim([0 2.5])
subplot(2,1,2); hold on; box on;title('CD'); % prepare CD plot
ylim([0 12])
% plot z and cd of the segment
subplot(2,1,1);
plot(data_day.numdate, data_day.avg_z_sm, ...
'LineWidth',0.9,'Color',[60/255 179/255 113/255])
subplot(2,1,2);
plot(data_day.numdate, data_day.centerdist, ...
'LineWidth',0.9,'Color',[60/255 179/255 113/255])
% select events
downevents = data_ld(floor(data_ld.numdatedown) == days(i),:); % all down
upevents = data_ld(floor(data_ld.numdatedown) == days(i),:); % all up
% add detection window
downevents.winstart = downevents.numdatedown - cpts.window/2;
downevents.winend = downevents.numdatedown + cpts.window/2;
upevents.winstart = upevents.numdateup - cpts.window/2;
upevents.winend = upevents.numdateup + cpts.window/2;
% plot events as ground thruth + window
for j = 1:height(downevents)
% z subplot
subplot(2,1,1);
plot([downevents.numdatedown(j) downevents.numdatedown(j)], ...
[0 2.5],'LineWidth',0.8,'Color',[204/255 0 102/255])
plot([downevents.winstart(j) downevents.winend(j)],...
[1 1],'LineWidth',4.5,'Color',[204/255 0 102/255])
% cd subplot
subplot(2,1,2);
plot([downevents.numdatedown(j) downevents.numdatedown(j)], ...
[0 15],'LineWidth',0.8,'Color',[204/255 0 102/255])
plot([downevents.winstart(j) downevents.winend(j)],...
[6 6],'LineWidth',4.5,'Color',[204/255 0 102/255])
end
for j = 1:height(upevents)
% z subplot
subplot(2,1,1); xlim([days(i) days(i)+1])
plot([upevents.numdateup(j) upevents.numdateup(j)], ...
[0 2.5],'LineWidth',0.8,'Color',[204/255 0 102/255])
plot([upevents.winstart(j) upevents.winend(j)],...
[1 1],'LineWidth',4.5,'Color',[204/255 0 102/255])
% cd subplot
subplot(2,1,2);
xlim([days(i) days(i)+1])
plot([upevents.numdateup(j) upevents.numdateup(j)], ...
[0 15],'LineWidth',0.8,'Color',[204/255 0 102/255])
plot([upevents.winstart(j) upevents.winend(j)],...
[6 6],'LineWidth',4.5,'Color',[204/255 0 102/255])
end
% changepoints avg_z
[cptsr.(cow).(day).ipoints, cptsr.(cow).(day).resid] = ...
findchangepts(data_day.avg_z_sm,...
'Statistic',cpts.CDstat,...
'MinDistance',cpts.MinDistance,...
'MinThreshold',1200);
cptsr.(cow).(day).ipoints = [[1;...
cptsr.(cow).(day).ipoints(1:end-1);...
cptsr.(cow).(day).ipoints(end)],...
[cptsr.(cow).(day).ipoints;height(data_day)]];
% calculate mean value + std
for j = 1:size(cptsr.(cow).(day).ipoints,1)
cptsr.(cow).(day).ipoints(j,3) = ...
mean(data_day.avg_z_sm(...
cptsr.(cow).(day).ipoints(j,1):...
cptsr.(cow).(day).ipoints(j,2)));
cptsr.(cow).(day).ipoints(j,4) = ...
std(data_day.avg_z_sm(...
cptsr.(cow).(day).ipoints(j,1):...
cptsr.(cow).(day).ipoints(j,2)));
% plot changepoints avg_z_sm
subplot(2,1,1)
plot([data_day.numdate(cptsr.(cow).(day).ipoints(j,2)) ...
data_day.numdate(cptsr.(cow).(day).ipoints(j,2))],...
[0 3],'r','LineWidth',1)
plot([data_day.numdate(cptsr.(cow).(day).ipoints(j,1)) ...
data_day.numdate(cptsr.(cow).(day).ipoints(j,2))],...
[cptsr.(cow).(day).ipoints(j,3) ...
cptsr.(cow).(day).ipoints(j,3)],'r','LineWidth',2)
end
% changepoints cd
[cptsr.(cow).(day).ipointsCD, cptsr.(cow).(day).residCD] = ...
findchangepts(data_day.centerdist,...
'Statistic',cpts.CDstat,...
'MinDistance',cpts.MinDistance,...
'MinThreshold',7154);
cptsr.(cow).(day).ipointsCD = [[1;...
cptsr.(cow).(day).ipointsCD(1:end-1);...
cptsr.(cow).(day).ipointsCD(end)],...
[cptsr.(cow).(day).ipointsCD;height(data_day)]];
% calculate mean value + std
for j = 1:size(cptsr.(cow).(day).ipointsCD,1)
cptsr.(cow).(day).ipointsCD(j,3) = mean(data_day.centerdist(...
cptsr.(cow).(day).ipointsCD(j,1):...
cptsr.(cow).(day).ipointsCD(j,2)));
cptsr.(cow).(day).ipointsCD(j,4) = std(data_day.centerdist(...
cptsr.(cow).(day).ipointsCD(j,1):...
cptsr.(cow).(day).ipointsCD(j,2)));
% plot changepoints centerdist
subplot(2,1,2)
plot([data_day.numdate(cptsr.(cow).(day).ipointsCD(j,2)) ...
data_day.numdate(cptsr.(cow).(day).ipointsCD(j,2))],...
[0 15],'b','LineWidth',1)
plot([data_day.numdate(cptsr.(cow).(day).ipointsCD(j,1)) ...
data_day.numdate(cptsr.(cow).(day).ipointsCD(j,2))],...
[cptsr.(cow).(day).ipointsCD(j,3)
cptsr.(cow).(day).ipointsCD(j,3)],'b',...
'LineWidth',2)
end
saveas(k,[init_.resdir 'figa_cpts_' cow '_' day '.fig'])
saveas(k,[init_.resdir 'figa_cpts_' cow '_' day '.tif'])
end
end
end
% save results
save([init_.resdir 'D5a_cpts_res.mat'],'cptsr','cpts')
% clear variables
clear ans cow data data_day daat_ld day days downevents f j i h ind
clear x y upevents seglist data2 k
%% Postprocessing - filter and combine results
% >>> SENSITIVITY
% >>> SPECIFICITY
% - change = in segment
% - difference in mean value (~std of the signal)
% - CD: difference at changepoints in level > 0.25m
% - CD: difference at changepoints in std >
% load results
% load([init_.resdir 'D4_cpts_res.mat'])
% cows
cows = fieldnames(cptsr);
for i = 1:length(fields_) % all cows with results
T=1; % count cpts for cd
T2=1; % count cpts for z
% all days for this cow
days = fieldnames(cptsr.(cows{i}));
close all
% add to changepoints the information for postprocessing
for j = 1:length(days)
% ref day numeric
day = days{j}; day = str2double(day(5:end)); % num day (~refday)
% select and store/split data
data.(cows{i}).(days{j}) = sdata.(cows{i})(:,[2 3 7 13 17 18 15]);
data.(cows{i}).(days{j})(isnan(data.(cows{i}).(days{j}).avg_z_sm) | ...
isnan(data.(cows{i}).(days{j}).centerdist),:) = [];
data.(cows{i}).(days{j}).numdate = ...
datenum(data.(cows{i}).(days{j}).date)-refdate; % numeric date
data.(cows{i}).(days{j}) = data.(cows{i}).(days{j})(...
floor(data.(cows{i}).(days{j}).numdate) == day,:);
% load figure and plot
% k = openfig([init_.resdir 'fig_cpts_' cows{i} '_' days{j} '.fig']);
% results z
rzdata = array2table(cptsr.(cows{i}).(days{j}).ipoints,...
'VariableNames',{'istart','iend','avg','std'});
% difference in mean for CD
rzdata.avgdif(1) = NaN;
rzdata.avgdif(2:end) = abs(diff(rzdata.avg));
% idx = find(rzdata.avgdif <= 0.25); % less than 0.25 meter diff
% if ~isempty(idx)
% rzdata.iend(idx-1) = rzdata.iend(idx);
% rzdata(idx,:) = [];
% end
% calculate mean value + std
for l = 1:height(rzdata)
rzdata.avgnew(l) = mean(...
data.(cows{i}).(days{j}).avg_z_sm(...
rzdata.istart(l):...
rzdata.iend(l)));
rzdata.stdnew(l) = std(...
data.(cows{i}).(days{j}).avg_z_sm(...
rzdata.istart(l):...
rzdata.iend(l)));
rzdata.numdatestart(l) = data.(cows{i}).(days{j}).numdate(...
rzdata.istart(l)); % start
rzdata.numdateend(l) = data.(cows{i}).(days{j}).numdate(...
rzdata.iend(l)); % end
% plot changepoints centerdist
% subplot(2,1,1)
% plot([data.(cows{i}).(days{j}).numdate(rzdata.iend(l)) ...
% data.(cows{i}).(days{j}).numdate(rzdata.iend(l))],...
% [0 2.5],'Color',[1 185/255 0],'LineWidth',1)
% plot([data.(cows{i}).(days{j}).numdate(rzdata.istart(l)) ...
% data.(cows{i}).(days{j}).numdate(rzdata.istart(l))],...
% [0 2.5],'Color',[1 185/255 0],'LineWidth',1)
% plot([data.(cows{i}).(days{j}).numdate(rzdata.istart(l)) ...
% data.(cows{i}).(days{j}).numdate(rzdata.iend(l))],...
% [rzdata.avgnew(l) rzdata.avgnew(l)],...
% 'Color',[1 185/255 0],'LineWidth',2.5)
%
end
% add that it's a changepoint to data
data.(cows{i}).(days{j}).ischangeZ(1) = T2;
data.(cows{i}).(days{j}).ischangeZ(rzdata.iend) = ...
(T2+1:T2+height(rzdata))';
% save new changepoints
fie = ['z_' days{j}];
data.(cows{i}).(fie) = rzdata;
T2 = T2+height(rzdata);
% results cd
rcpdata = array2table(cptsr.(cows{i}).(days{j}).ipointsCD,...
'VariableNames',{'istart','iend','avg','std'});
% difference in mean for CD
rcpdata.avgdif(1) = NaN;
rcpdata.avgdif(2:end) = abs(diff(rcpdata.avg));
% idx = find(rcpdata.avgdif <= 0.5);% less than 0.5 meter diff
% if ~isempty(idx)
% rcpdata.iend(idx-1) = rcpdata.iend(idx);
% rcpdata(idx,:) = [];
% end
% calculate mean value + std
for l = 1:height(rcpdata)
rcpdata.avgnew(l) = mean(...
data.(cows{i}).(days{j}).centerdist(...
rcpdata.istart(l):...
rcpdata.iend(l)));
rcpdata.stdnew(l) = std(...
data.(cows{i}).(days{j}).centerdist(...
rcpdata.istart(l):...
rcpdata.iend(l)));
rcpdata.numdatestart(l) = data.(cows{i}).(days{j}).numdate(...
rcpdata.istart(l)); % start
rcpdata.numdateend(l) = data.(cows{i}).(days{j}).numdate(...
rcpdata.iend(l)); % end
% plot changepoints centerdist
% subplot(2,1,2)
% plot([data.(cows{i}).(days{j}).numdate(rcpdata.iend(l)) ...
% data.(cows{i}).(days{j}).numdate(rcpdata.iend(l))],...
% [0 15],'Color',[1 185/255 0],'LineWidth',1)
% plot([data.(cows{i}).(days{j}).numdate(rcpdata.istart(l)) ...
% data.(cows{i}).(days{j}).numdate(rcpdata.istart(l))],...
% [0 15],'Color',[1 185/255 0],'LineWidth',1)
% plot([data.(cows{i}).(days{j}).numdate(rcpdata.istart(l)) ...
% data.(cows{i}).(days{j}).numdate(rcpdata.iend(l))],...
% [rcpdata.avgnew(l) rcpdata.avgnew(l)],...
% 'Color',[1 185/255 0],'LineWidth',2.5)
end
% add that it's a changepoint to data
data.(cows{i}).(days{j}).ischangeCD(1) = T;
data.(cows{i}).(days{j}).ischangeCD(rcpdata.iend) = ...
(T+1:T+height(rcpdata))';
T = T+height(rcpdata);
% save new changepoints
fie = ['cd_' days{j}];
data.(cows{i}).(fie) = rcpdata;
% save figure
% saveas(k, [init_.resdir 'nfig_cpts_' cows{i} '_' days{j} '.fig'])
% saveas(k, [init_.resdir 'nfig_cpts_' cows{i} '_' days{j} '.tif'])
end
end
save([init_.resdir 'D5a_dataday.mat'],'data', '-v7.3')
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment