function [is_rejected,p_out,p_hi,p_lo,scale,shape] = ...
    wbl_tail_test_func(samples,block_labels,left_cens_prctile_thr,p_test,niter,censorams)
% 
%   F. Marra - Oct 2022
% 
% H0: block maxima are samples from a parent distribution with Weibull tail
%     (tail defined by a given left censoring threhsold)
% 
% INPUT
%   samples: data sample (array with all independent ordinary events) 
%   block_label: labels with info on the blocks of the samples (array with
%                same size as samples, same value for each block)
%   left_cens_prctile_thr: left censoring threshold (percentile of samples)  
%   p_test: probability to be used for the test (test will count the block
%           maxima out of the Y = 1-p_out confidence interval)
%   niter: number of stochastic realizations
%   censorams: logical variable (use true to censor the block maxima)
% 
% OUTPUT
%   is_rejected: outcome of the test (true means that the assumption of
%                Weibull tails for the given left-censoring threshold is
%                rejected) 
%   p_out: fraction of block maxima outside of the Y = 1-p_out confidence
%          interval 
%   p_hi: fraction of block maxima higher of the Y = 1-p_out confidence
%         interval 
%   p_lo: fraction of block maxima lower of the Y = 1-p_out confidence
%         interval 
%   scale: scale parameter of the Weibull distribution describing the
%          non-censored samples
%   shape: shape parameter of the Weibull distribution describing the
%          non-censored samples% 

[samples,ib] = sort(samples); % sorts samples
block_labels = block_labels(ib); % sorts block labels accordingly

% identifies block maxima
is_block_maximum = false(size(samples));
for id = unique(block_labels)'
    block_maximum=0;
    for l = find(block_labels==id)'
        if samples(l)>block_maximum; block_maximum = samples(l); i=l; end
    end
    is_block_maximum(i)=true;
end

% estimates Weibull parameters censoring the samples that do not exceed the
% given left-censoring & the block maxima  
to_use = max(1,fix(numel(samples)*left_cens_prctile_thr)) : ceil(numel(samples)); % censors samples below threhsold
if censorams
    to_use = to_use(ismember(to_use,find(~is_block_maximum))); % censors block maxima (optional)
end

% parameter estimation (least square linear regression in
% Weibull-transformed coordinates) 
ECDF = (1:numel(samples))'/(numel(samples)+1);  % empirical probabilities computed with Weibull plotting positions
X = log(log(1./(1-ECDF(to_use))));              % Weibull-tranformation for probabilities
Y = log(samples(to_use));                       % Weibull-tranformation for samples
tmp=regress(Y(:),[ones(size(X(:))) X(:)]);      % linear regression
scale = exp(tmp(1));                            % Weibull scale parameter
shape = 1/tmp(2);                               % Weibull shape parameter

istest = is_block_maximum; % values over which the hypothesis is tested, i.e. block maxima

randy = sort( wblinv(rand(niter,numel(samples)),scale,shape), 2);           % weibull-distributed stochastic samples
p_lo = nanmean( samples(istest)<quantile(randy(:,istest),p_test/2,1)' );    % fraction of block maxima below the (1-p) CI
p_hi = nanmean( samples(istest)>quantile(randy(:,istest),1-p_test/2,1)' );  % fraction of block maxima above the (1-p) CI

p_out = p_hi + p_lo; % fraction of block maxima out of the (1-p) CI
 
is_rejected = p_out>p_test; % outcome of the test

end