%function outliers = outlier_bt500(scores_raw)

% Analyzes the scores and screens the observers finding the outliers as
% defined by ITU-R BT.500.
% input:
%  - scores_raw: cols = subjects, rows = presentations
% output:
%  - scores_screened: outlier scores set to NaN
%  - nb_outliers: number of outliers
%  - outliers: outliers

% number of subjects and test conditions
load('Exp3.mat')
scores_raw = Ann;

[L,N] = size(scores_raw);

u = scores_raw;

% mean and standard S
u_ = mean(scores_raw,2);
S = std(scores_raw,0,2);

% compute kurtosis coefficients
beta2 = kurtosis(scores_raw,1,2);

outliers = [];
for i = 1 : N
    Pi = 0;
    Qi = 0;
    for l = 1 : L
        if beta2(l) >= 2 && beta2(l) <= 4
            if u(l,i) >= u_(l) + 2 * S(l)
                Pi = Pi + 1;
            end
            if u(l,i) <= u_(l) - 2 * S(l)
                Qi = Qi + 1;
            end
        else
            if u(l,i) >= u_(l) + sqrt(20) * S(l)
                Pi = Pi + 1;
            end
            if u(l,i) <= u_(l) - sqrt(20) * S(l);
                Qi = Qi + 1;
            end
        end
    end
    
    if (Pi + Qi) / L > 0.05 && abs((Pi - Qi)/(Pi + Qi)) < 0.3
        outliers = [outliers ; i];
    end
end

 nb_outliers = length(outliers);
 scores_screened = scores_raw;
 scores_screened(:,outliers) = [];