%% Function to compute scores for the 2018 PhysioNet/CinC Challenge.
%%
%% This function aggregates the outputs of a proposed classification
%% algorithm, and calculates the overall (gross) area under the
%% precision-recall curve, as well as the area under the receiver
%% operating characteristic curve.
%%
%% 'truth' is a vector of arousal values: zero for non-arousal
%% regions, positive for target arousal regions, and negative for
%% unscored regions.
%%
%% 'predictions' is a vector of probabilities produced by the
%% classification algorithm being tested. This vector must be the
%% same length as 'truth', and each value must be between 0 and 1.
%%
%% Call this function once for each record being tested; it will
%% return the gross AUPRC and AUROC (of all records so far), as well
%% as the per-record AUPRC and AUROC (which will be NaN if the record
%% contains no arousals.) Call the function with no arguments to
%% re-initialize.
function [auprc_g, auroc_g, auprc_r, auroc_r] = score2018(truth, predictions)
vscale = 1000;
persistent all_pos_values;
persistent all_neg_values;
if nargin == 0 || isempty(all_pos_values)
all_pos_values = zeros(1, vscale + 1);
all_neg_values = zeros(1, vscale + 1);
end
if nargin == 0
return;
end
%% Check if length is correct
if length(predictions) ~= length(truth)
error('length of predictions does not match truth');
end
%% Compute the histogram of all input probabilities
bins = [-1/2 : vscale + 1/2] / vscale;
all_values = xhistogram(predictions, bins);
%% Check if input contains any out-of-bounds or NaN values
if sum(all_values) ~= length(predictions)
error('invalid values in predictions');
end
%% Compute the histogram of probabilities within arousal regions
pos_samples = find(truth > 0);
pos_values = xhistogram(predictions(pos_samples), bins);
%% Compute the histogram of probabilities within unscored regions
ign_samples = find(truth < 0);
ign_values = xhistogram(predictions(ign_samples), bins);
%% Compute the histogram of probabilities in non-arousal regions,
%% given the above
neg_values = all_values - pos_values - ign_values;
all_pos_values = all_pos_values + pos_values;
all_neg_values = all_neg_values + neg_values;
[auprc_r, auroc_r] = xauc(pos_values, neg_values);
[auprc_g, auroc_g] = xauc(all_pos_values, all_neg_values);
end
function h = xhistogram(predictions, bins)
if isempty(predictions)
h = zeros(1, length(bins) - 1);
else
h = histc(predictions, bins);
h = h(1 : length(bins) - 1);
if ~isrow(h)
h = h';
end
end
end
function [auprc, auroc] = xauc(pos_values, neg_values)
%% Calculate areas under the ROC and PR curves by iterating over
%% the possible threshold values.
%% At the minimum threshold value, all samples are classified as
%% positive, and thus TPR = 1 and TNR = 0.
tp = sum(pos_values);
fp = sum(neg_values);
tn = 0;
fn = 0;
if tp == 0 || fp == 0
%% If either class is empty, scores are undefined.
auprc = nan;
auroc = nan;
return;
end
tpr = 1;
tnr = 0;
ppv = tp / (tp + fp);
auprc = 0;
auroc = 0;
%% As the threshold increases, TP decreases (and FN increases)
%% by pos_values(j), while TN increases (and FP decreases) by
%% neg_values(j).
for j = 1 : length(pos_values)
tp = tp - pos_values(j);
fn = fn + pos_values(j);
fp = fp - neg_values(j);
tn = tn + neg_values(j);
tpr_prev = tpr;
tnr_prev = tnr;
ppv_prev = ppv;
tpr = tp / (tp + fn);
tnr = tn / (tn + fp);
if tp + fp > 0
ppv = tp / (tp + fp);
end
%% ROC is piecewise linear; integrate by adding the areas
%% of trapezoids
auroc = auroc + (tpr_prev - tpr) * (tnr + tnr_prev) * 0.5;
%% PRC is piecewise constant; integrate by adding the
%% areas of rectangles
auprc = auprc + (tpr_prev - tpr) * ppv_prev;
end
end