matthiasm@8: function [FPrate, TPrate, AUC, thresholds] = computeROC(confidence, testClass)
matthiasm@8: % function [FPrate, TPrate, AUC, thresholds] = computeROC(confidence, testClass)
matthiasm@8: %
matthiasm@8: % computeROC computes the data for an ROC curve based on a classifier's confidence output.
matthiasm@8: % It returns the false positive rate and the true positive rate along with
matthiasm@8: % the area under the ROC curve, and the list of thresholds.
matthiasm@8: %
matthiasm@8: % Inputs: 
matthiasm@8: %           - confidence(i) is proportional to the probability that
matthiasm@8: %             testClass(i) is positive
matthiasm@8: %
matthiasm@8: % testClass = 0 => target absent
matthiasm@8: % testClass = 1 => target present
matthiasm@8: %
matthiasm@8: % Based on algorithms 2 and 4 from Tom Fawcett's paper "ROC Graphs: Notes and
matthiasm@8: % Practical Considerations for Data Mining Researchers" (2003)
matthiasm@8: % http://www.hpl.hp.com/techreports/2003/HPL-2003-4.pdf"
matthiasm@8: %
matthiasm@8: % Vlad Magdin, 21 Feb 2005
matthiasm@8: 
matthiasm@8: % break ties in scores
matthiasm@8: S = rand('state');
matthiasm@8: rand('state',0);
matthiasm@8: confidence = confidence + rand(size(confidence))*10^(-10);
matthiasm@8: rand('state',S)
matthiasm@8: [thresholds order] = sort(confidence, 'descend');
matthiasm@8: testClass = testClass(order);
matthiasm@8: 
matthiasm@8: %%% -- calculate TP/FP rates and totals -- %%%
matthiasm@8: AUC = 0;
matthiasm@8: faCnt = 0;
matthiasm@8: tpCnt = 0;
matthiasm@8: falseAlarms = zeros(1,size(thresholds,2));
matthiasm@8: detections = zeros(1,size(thresholds,2));
matthiasm@8: fPrev = -inf;
matthiasm@8: faPrev = 0;
matthiasm@8: tpPrev = 0;
matthiasm@8: 
matthiasm@8: P = max(size(find(testClass==1)));
matthiasm@8: N = max(size(find(testClass==0)));
matthiasm@8: 
matthiasm@8: for i=1:length(thresholds)
matthiasm@8:     if thresholds(i) ~= fPrev
matthiasm@8:         falseAlarms(i) = faCnt;
matthiasm@8:         detections(i) = tpCnt;
matthiasm@8: 
matthiasm@8:         AUC = AUC + polyarea([faPrev faPrev faCnt/N faCnt/N],[0 tpPrev tpCnt/P 0]);
matthiasm@8: 
matthiasm@8:         fPrev = thresholds(i);
matthiasm@8:         faPrev = faCnt/N;
matthiasm@8:         tpPrev = tpCnt/P;
matthiasm@8:     end
matthiasm@8:     
matthiasm@8:     if testClass(i) == 1
matthiasm@8:         tpCnt = tpCnt + 1;
matthiasm@8:     else
matthiasm@8:         faCnt = faCnt + 1;
matthiasm@8:     end
matthiasm@8: end
matthiasm@8: 
matthiasm@8: AUC = AUC + polyarea([faPrev faPrev 1 1],[0 tpPrev 1 0]);
matthiasm@8: 
matthiasm@8: FPrate = falseAlarms/N;
matthiasm@8: TPrate = detections/P;