%ROC Receiver-Operator Curve
%
% E = ROC(B,C,N)
%
% INPUT
% A Dataset
% W Trained classifier, or
% B Classification result, B = A*W*CLASSC
% C Index of desired class (default: C = 1)
% N Number of points on the Receiver-Operator Curve (default: 100)
%
% OUTPUT
% E Structure containing the error of the two classes
%
% DESCRIPTION
% Computes N points on the receiver-operator curve of the classifier W for
% class C in the labeled dataset B, which is typically the result of
% B = A*W; or for the dataset A labelled by applying the (cell array of)
% trained classifiers W.
%
% Note that a Receiver-Operator Curve is related to a specific class (class C)
% for which the errors are plotted horizontally. The total error on all other
% classes is plotted vertically. The class index C refers to its position in
% the label list of the dataset (A or B). It can be found by GETCLASSI.
%
% The curve is computed for N thresholds of the posteriori probabilities
% stored in B. The resulting error frequencies for the two classes are
% stored in the structure E. E.XVALUES contains the errors in the first
% class, E.ERROR contains the errors in the second class. In multi-class
% problems these are the mean values in a single class, respectively the
% mean values in all other classes. This may not be very useful, but not
% much more can be done as for multi-class cases the ROC is equivalent to a
% multi-dimensional surface.
%
% Use PLOTE(E) for plotting the result. In the plot the two types of error
% are annotated as 'Error I' (error of the first kind) and 'Error II' (error
% of the second kind). All error estimates are weighted according the class
% prior probabilities. Remove the priors in A or B (by setprior(A,[])) to
% produce a vanilla ROC.
%
% EXAMPLES
% Train set A and test set T:
% B = T*NMC(A); E = ROC(T,50); PLOTE(E); % Plots a single curve
% E = ROC(T,A*{NMC,UDC,QDC}); PLOTE(E); % Plots 3 curves
%
function [roc, thr] = roc_detection(a,n)
% Depending on the call, CLAS may the third or second argument.
% and N the third or the fourth.
if nargin < 2 || isempty(n), n = 100; end
datname = getname(a);
lablist = getlablist(a,'string');
required_labs = {'FP' 'TP'};
aux_lablist = intersect(cellstr(lablist), required_labs );
cant_required_labs = length(required_labs);
if( length(aux_lablist) ~= cant_required_labs )
fprintf(2, ['Esta funcion esta pensada para usarse en datasets de deteccion con labels:\n' colvec([char(required_labs) repmat('\n', cant_required_labs, 1)]')' ] );
error();
end
clas = find( strcmpi(cellstr(lablist), 'TP') );
clasname = lablist(clas,:);
%DXD: also check the class sizes:
cs = classsizes(a);
if any(cs == 0)
error('Ambas clases deben contener ejemplos');
end
% Set up the ROC structure.
% If a cell array of classifiers was given, apply each one here.
a = a*normm; % make sure we have a normalised classification matrix
[m,c] = size(a);
nlab = getnlab(a);
d = sort(a(:));
% Attempt to compute a good threshold range on the first classified
% dataset.
% thr = [max(0, min(d)-eps) rowvec(d(round(linspace(2,length(d)-1, n-1)))) 1];
thr = linspace(0,1,n+1);
% NLAB_OUT will be one where B is larger than THR.
I = matchlablist(getlablist(a),getfeatlab(a)); % Correct possible changes class orders
nlab_out = (repmat(+a(:,I(clas)),1,n+1) >= repmat(thr,m,1));
% aciertos will be 1 where the numeric label is unequal to NLAB_OUT
% (i.e., where errors occur).
bTP = nlab == clas;
aciertos = (repmat(bTP,1,n+1) == nlab_out);
sensitivity = mean(aciertos(bTP,:),1); % S
pospred = sum(aciertos(bTP,:)) ./ sum(nlab_out); % +P
roc = [ colvec(sensitivity) colvec(pospred) ];
mod = sqrt(sum(roc.^2,2));
[max_mod max_mod_idx] = max(mod);
figure(100);
h = subplot(1,2,1);
plot(h(1), sensitivity(:), pospred(:), 'bo-' )
hold(h(1), 'on');
plot(h(1), sensitivity(max_mod_idx), pospred(max_mod_idx), 'rx:', 'MarkerSize',11)
plot(h(1), sensitivity(max_mod_idx), pospred(max_mod_idx), 'mo:', 'MarkerSize',11)
hold(h(1), 'off');
axis(h(1), 'square');
box(h(1), 'off')
xlabel(h(1), 'Sensitivity')
ylabel(h(1), 'Positive predictivity')
% x_lim = xlim();
% y_lim = ylim();
% xlim([min([x_lim ylim]) 1.1]);
% xlim([min([x_lim ylim]) 1.1]);
h(2) = subplot(1,2,2);
plot(h(2), thr, mod, 'bo-' )
hold(h(2), 'on');
plot(h(2), thr(max_mod_idx), mod(max_mod_idx), 'rx', 'MarkerSize',11)
plot(h(2), thr(max_mod_idx), mod(max_mod_idx), 'mo', 'MarkerSize',11)
plot(h(2), [0.5 0.5], ylim(), 'k--')
hold(h(2), 'off');
box(h(2), 'off')
ylabel(h(2), 'mod(S,P+)')
xlabel(h(2), 'Operating point')
ConfusionMat = [ sum(aciertos(~bTP,max_mod_idx)) sum(~aciertos(~bTP,max_mod_idx)); sum(~aciertos(bTP,max_mod_idx)) sum(aciertos(bTP,max_mod_idx)) ];
DisplayResults('dsResult', ConfusionMat, 'SupportDataset', a);