function [Mhat, Phat, Nhat, lhat] = func_empiricalEstimators_Histo(m, Y)
% Computes the empirical estimators of the spectral method for the histogram basis.

disp('Computing empirical estimators.');
tic;

temp = size(Y);
n = temp(2);

indicator=zeros(m,m,m);
for a=1:m
    indicA = (Y(:,1:(n-2))>(a-1)/m) .* (Y(:,1:(n-2))<=a/m);
    for b=1:m
        indicB = (Y(:,2:(n-1))>(b-1)/m) .* (Y(:,2:(n-1))<=b/m);
        for c=1:m
            indicator(a,b,c) = mean( indicA .*...
                indicB .*...
                (Y(:,3:n)>(c-1)/m) .* (Y(:,3:n)<=c/m ) );
        end
    end
end
Mhat = (m^(3/2)) * indicator;
Phat = (m) * squeeze(sum(indicator,2));
Nhat = (m) * squeeze(sum(indicator,3));
lhat = (m^(1/2)) * squeeze(sum(sum(indicator,3),2));

elapsedTime = toc;
disp([' ... finished in ' num2str(elapsedTime) ' seconds.']);
