function crit = MM_CritBW(A, mmin, mmax, points, demarc)
% Written by Melanie Krause, this version 25.07.2016.

% For a nonparametric Kernel density estimation, this is the procedure to
% find the critical bandwidth as defined by Silverman (1981): the lowest
% bandwidth consistent with an m-modal density function (lower bandwidths
% lead to a density with at least m+1 modes).
% We want to find the mmax-element vector containing the critical
% bandwidths for densities with mmin < m < mmax modes. 
% We use a binary search procedure: The algorithm searches in the interval
% [Low, Up] with changing upper and lower bounds so that the interval halves 
% in each iteration until convergence.

% A is the vector containing the observations, mmax the number of modes up
% to which the critical bandwidth should be calculated and points is the
% number of points at which the density should be evaluated. 
% Horizontal demarcation line: If demarc = 0 (default) even tiny modes
% count, if it increases modes whose highest density value is below demarc
% will not count.

crit=zeros(mmax-mmin+1,1);
n=100; %iterations in the binary search procedure
Low = 0; % lower interval bound for critical bandwidth search
High = 2*max(A); % higher interval bound for critical bandwidth search
U=zeros(n,mmax-mmin+1); %matrix with the iterated bandwidths

for m=mmin:mmax %(mode-number whose critical value we look for)
U(1,m-mmin+1)=Low + (High-Low)/2; %initial value
for i=1:n %(iterations)
    [f,xi] = ksdensity(A, 'npoints', points, 'width', U(i,m-mmin+1)); % estimate Kernel density with i-th bandwidth 
   % Mode Counting: Method 1 - indicate mode at a point z if derivative is first > 0, then < 0:
   f_diff=diff(f'); % discrete derivative 
    mode_indicator=zeros(length(f_diff),1);
       for z=2:length(f_diff)
       if f_diff(z) < 0 & f_diff(z-1) >0 & max(f(z),f(z-1))> demarc
       mode_indicator(z)=1;  
       else mode_indicator(z)=0;
       end
       end
 % If we have m modes, see if we can decrease bandwidth further in next
  % iteration (set current bandwidth as upper bound); 
  % if we have already m+1 modes, increase bandwidth again 
  % (set current bandwidth as upper bound). 
       if sum(mode_indicator)<=m 
        High = U(i,m-mmin+1);
    else Low = U(i,m-mmin+1);
    end
  % Next bandwidth as midpoint of the new smaller interval
     U(i+1,m-mmin+1)=Low + (High-Low)/2; 
    %Break when we converge:
    if abs(U(i+1,m-mmin+1)-U(i,m-mmin+1))< 10^(-6)
        crit(m-mmin+1)=U(i,m-mmin+1); %store critical bandwidth
        High = 2*U(i,m-mmin+1); %critical values for densities with more modes must lower, so take this as starting value
        Low = 0; %reset lower bound to zero
        break
    end
end

end