Persistence Barcodes Versus Kolmogorov Signatures: Detecting Modes of One-Dimensional Signals

Abstract

We investigate the problem of estimating the number of modes (i.e., local maxima)-a well known question in statistical inference-and we show how to do so without presmoothing the data. To this end, we modify the ideas of persistence barcodes by first relating persistence values in dimension one to distances (with respect to the supremum norm) to the sets of functions with a given number of modes, and subsequently working with norms different from the supremum norm. As a particular case we investigate the Kolmogorov norm. We argue that this modification has certain statistical advantages. We offer confidence bands for the attendant Kolmogorov signatures, thereby allowing for the selection of relevant signatures with a statistically controllable error. As a result of independent interest, we show that taut strings minimize the number of critical points for a very general class of functions. We illustrate our results by several numerical examples.

Topological density estimation

“…The statistics toolbox is required for this script. %% Loop over PDFs for j = 1:size(X,1) %% Optimal kernel density estimate w/ knowledge of PDF opt = optimalkernelbandwidth(x,f(j,:),n(aa),kernelflag); if kernelflag > 0 pd_opt = fitdist(X(j,:)',kernelargstr,'BandWidth',opt.h); elseif kernelflag < 0 pd_opt = fitdist(X(j,:)','Kernel','Kernel',kernelargstr,... 'BandWidth',opt.h); end kde_opt = pdf(pd_opt,x); %% Numerically ISE minimizing bandwidth temp = vbise(x,f(j,:),X(j,:),kernelflag); ind = find(temp.ISE==min(temp.ISE),1,'first'); h_opt = temp.h(ind); if kernelflag > 0 pd_numopt = fitdist(X(j,:)',kernelargstr,'BandWidth',h_opt); elseif kernelflag < 0 pd_numopt = fitdist(X(j,:)','Kernel','Kernel',... kernelargstr,'BandWidth',h_opt); end kde_numopt = pdf(pd_numopt,x); %% Cross-validation KDE cv = cv1d(X(j,:)',kernelflag); if kernelflag > 0 % % For MATLAB estimator % pd = fitdist(X(j,:)',kernelargstr); pd = fitdist(X(j,:)',kernelargstr,'BandWidth',cv.h); elseif kernelflag < 0 % % For MATLAB estimator % pd = fitdist(X(j,:)','Kernel','Kernel',kernelargstr); pd = fitdist(X(j,:)','Kernel','Kernel',kernelargstr,'BandWidth',cv.h); end kde = pdf(pd,x); h_hat = pd.BandWidth; pre_c1_CV(ii,j,aa) = h_hat-h_opt; pre_c23_CV(ii,j,aa) = ise(x,f(j,:),X(j,:),h_hat,kernelflag); pre_c45_CV(ii,j,aa) = ... pre_c23_CV(ii,j,aa)-ise(x,f(j,:),X(j,:),h_opt,kernelflag); %% Topological KDE tde = tde1d(X(j,:),kernelflag); f_tde = tde.y/(sum(tde.y)*mean(diff(tde.x))); h_hat = tde.h; pre_c1_top(ii,j,aa) = h_hat-h_opt; pre_c23_top(ii,j,aa) = ise(x,f(j,:),X(j,:),h_hat,kernelflag); pre_c45_top(ii,j,aa) = ... pre_c23_top(ii,j,aa)-ise(x,f(j,:),X(j,:),h_opt,kernelflag); %% Unimodal category and number of local maxima ucat_CV(ii,j,aa) = nnz(sum(unidec(kde,0),2)>sqrt(eps)); lmax_CV(ii,j,aa) = nnz(diff(sign(diff([0,kde,0])))<0); ucat_top(ii,j,aa) = tde.mfuc; lmax_top(ii,j,aa) = nnz(diff(sign(diff([0,f_tde,0])))<0); end end end %% c1 for j = 1:6 % Normalized histogram lo = min([min(min(pre_c1_top(:,j,:))),min(min(pre_c1_CV(:,j,:)))]); hi = max([max(max(pre_c1_top(:,j,:))),max(max(pre_c1_CV(:,j,:)))]); L = linspace(lo,hi,50); for k = 1:5 % n = [25,50,100,200,500] A(:,k) = histc(squeeze(pre_c1_top(:,j,k)),L,1); B(:,k) = histc(squeeze(pre_c1_CV(:,j,k)),L,1); end eval(normmacro); % Plot subplot(3,6,j); eval(plotmacro); title(['$f_',num2str(j),'$'],'Interpreter','latex'); end subplot(3,6,1); ylabel('$\hat h -h_{opt}$','Interpreter','latex'); %% c2, c3 for j = 1:6 % Normalized histogram lo = 0; hi = max([max(max(pre_c23_top(:,j,:))),max(max(pre_c23_CV(:,j,:)))]); L = linspace(lo,hi,50); for k = 1:5 % n = [25,50,100,200,500] A(:,k) = histc(squeeze(pre_c23_top(:,j,k)),L,1); B(:,k) = histc(squeeze(pre_c23_CV(:,j,k)),L,1); end eval(normmacro); % Plot subplot (3,6,j+6); eval(plotmacro); end subplot (3,6,7); ylabel('ISE$(\hat h)$','Interpreter','latex'); %% c4, c5 for j = 1:6 % Form histograms (note logspace). NB.…”

“…NB. The only difference between L1 % and L2 versions is rescaling lo = -4; hi = 0.25; L = linspace(lo,hi,50); % abs turns out to not have any effect, but included for correctness for k = 1:5 % n = [25,50,100,200,500] A(:,k) = histc(squeeze(log10(abs(pre_c45_top(:,j,k)))),L,1); B(:,k) = histc(squeeze(log10(abs(pre_c45_CV(:,j,k)))),L,1); end eval(normmacro); % Plot subplot (3,6,j+12); eval(plotmacro); xlabel('$n$','Interpreter','latex'); end subplot (3,6,13); ylabel('$c_{45}$','Interpreter','latex');…”

“…% Script for plotting performance data for evaluating CV and TDE on the % densities $f_j$ ($1 \le j \le 6$). %% Plot macro % "box on" gives inadequate results but retained for tickmarks; manually % set axis and put in a box since "box on" doesn't give adequate results plotmacro = ['layertrans2(0:5,[L,Inf],A,B);box on;',... 'set(gca,''XTick'',1:5,''XTickLabel'',',... '{''25'',''50'',''100'',''200'',''500''});',... 'hold on;axis([0.5,5.5,lo,hi]);ax = axis;',... 'line ([ax(1),ax(1)]+1e-6,[ax(3),ax(4)],''Color'',''k'');',... 'line ([ax(2),ax (2)]-1e-6,[ax(3),ax(4)],''Color'',''k'');',... 'line ([ax(1),ax (2)],[ax(3),ax(3)]+1e-6,''Color'',''k'');',... 'line ([ax(1),ax (2)],[ax(4),ax(4)]-1e-6,''Color'',''k'');']; % and L2 versions is rescaling % The following is roughly equivalent in practice to % temp = cat(3,pre_c45_top(:,:,j),pre_c45_CV(:,:,j)); % qlog = .05; % quantile margin for log cutoff % lo = log10(quantile(temp(:),qlog)); % hi = max([0,log10(max(temp(:)))]); % but has the distinct advantage of being constant lo = -4; hi = 0.25; L = linspace(lo,hi,50); % abs turns out to not have any effect, but included for correctness A = histc(log10(abs(pre_c45_top(:,:,j))),L,1); B = histc(log10(abs(pre_c45_CV(:,:,j))),L,1); eval(normmacro); % Transparency plots subplot (3,3,3); ind = 7:16; eval(plotmacro); title('$c_{45}$','Interpreter','latex'); subplot (3,3,6); ind = 17:26; eval(plotmacro); subplot (3,3,9); ind = 27:36; eval(plotmacro);…”

“…Kernel density estimation (KDE) is such a ubiquitous and fundamental technique in statistics that our claim in this paper of an interesting and useful new contribution to the enormous body of literature (see, e.g., [6,[15][16][17]) almost inevitably entails some degree of hubris. Even the idea of using KDE to determine the multimodal structure of a probability density function (henceforth simply called "density") has by now a long history [8,10,14] that has very recently been explicitly coupled with ideas of topological persistence [3].…”

Topological data analysis for true step detection in periodic piecewise constant signals

“…Recent work has come in the form of classification and quantification of periodicity and quasi-periodicity [34][35][36][37][38]. Other applications using TSP include wheeze detection [39], computer performance [40], market prices [41,42], sonar [43], image processing [44], mode hunting [45], ice core analysis [46], machining dynamics [47][48][49][50] and gene expression [51,52].…”