# source:orange/Orange/evaluation/scoring.py@11725:983a29afcb85

Revision 11725:983a29afcb85, 99.1 KB checked in by Ales Erjavec <ales.erjavec@…>, 7 months ago (diff)

Fixed use of deprecated attribute names.

Line
1import math
2import functools
3from operator import add
4from collections import Iterable
5
6import numpy
7
8import Orange
9from Orange import statc, corn
10from Orange.utils import deprecated_keywords, deprecated_function_name, \
11    deprecation_warning, environ
12from Orange.evaluation import testing
13
14try:
15    import matplotlib
16    HAS_MATPLOTLIB = True
17except ImportError:
18    matplotlib = None
19    HAS_MATPLOTLIB = False
20
21#### Private stuff
22
23def log2(x):
24    """Calculate logarithm in base 2."""
25    return math.log(x) / math.log(2)
26
27def check_non_zero(x):
28    """Throw Value Error when x = 0."""
29    if x == 0.:
30        raise ValueError, "Cannot compute the score: no examples or sum of weights is 0."
31
32def gettotweight(res):
33    """Sum all the weights."""
34    totweight = reduce(lambda x, y: x + y.weight, res.results, 0)
35    if totweight == 0.:
36        raise ValueError, "Cannot compute the score: sum of weights is 0."
37    return totweight
38
39def gettotsize(res):
40    """Get number of result instances."""
41    if len(res.results):
42        return len(res.results)
43    else:
44        raise ValueError, "Cannot compute the score: no examples."
45
46# Backward compatibility
47def replace_use_weights(fun):
48    if environ.orange_no_deprecated_members:
49        return fun
50
51    @functools.wraps(fun)
52    def wrapped(*args, **kwargs):
53        use_weights = kwargs.pop("useWeights", None)
54        if use_weights is not None:
55            deprecation_warning("useWeights", "ignore_weights")
56            kwargs["ignore_weights"] = not use_weights
57        return fun(*args, **kwargs)
58    return wrapped
59
60def replace_discrete_probabilities_with_list(method):
61    if environ.orange_no_deprecated_members:
62        return lambda fun: fun
63
64    def decorator(fun):
65        @functools.wraps(fun)
66        def wrapped(*args, **kwargs):
67            res = args[method] if len(args)>method else kwargs.get("res", kwargs.get("test_results", None))
68            convert = res is not None
69
70            if convert:
71                old_probs = []
72                for r in res.results:
73                    old_probs.append(r.probabilities)
74                    r.probabilities = [list(p) if type(p) is Orange.statistics.distribution.Discrete
75                                       else p for p in r.probabilities]
76            result = fun(*args, **kwargs)
77            if convert:
78                for r, old in zip(res.results, old_probs):
79                    r.probabilities = old
80            return result
81        return wrapped
82    return decorator
83
84def split_by_iterations(res):
85    """Split ExperimentResults of a multiple iteratation test into a list
86    of ExperimentResults, one for each iteration.
87    """
88    if res.number_of_iterations < 2:
89        return [res]
90
91    ress = [Orange.evaluation.testing.ExperimentResults(
92                1, res.classifier_names, res.class_values,
93                res.weights, classifiers=res.classifiers,
95            for _ in range(res.number_of_iterations)]
96    for te in res.results:
97        ress[te.iteration_number].results.append(te)
98    return ress
99
100
101def split_by_classifiers(res):
102    """Split an instance of :obj:ExperimentResults into a list of
103    :obj:ExperimentResults, one for each classifier.
104    """
105    split_res = []
106    for i in range(len(res.classifier_names)):
107        r = Orange.evaluation.testing.ExperimentResults(
108                    res.number_of_iterations,
109                    [res.classifier_names[i]], res.class_values,
110                    weights=res.weights, base_class=res.base_class,
111                    classifiers=[res.classifiers[i]] if res.classifiers else [],
112                    test_type=res.test_type, labels=res.labels)
113        r.results = []
114        for te in res.results:
115            r.results.append(Orange.evaluation.testing.TestedExample(
116                te.iteration_number, te.actual_class, n=1, weight=te.weight))
117            r.results[-1].classes = [te.classes[i]]
118            r.results[-1].probabilities = [te.probabilities[i]]
119
120        if hasattr(res, "examples"):
121            r.examples = res.examples
122
123        split_res.append(r)
124
125    return split_res
126
127
128def class_probabilities_from_res(res, **argkw):
129    """Calculate class probabilities."""
130    probs = [0.] * len(res.class_values)
131    if argkw.get("unweighted", 0) or not res.weights:
132        for tex in res.results:
133            probs[int(tex.actual_class)] += 1.
134        totweight = gettotsize(res)
135    else:
136        totweight = 0.
137        for tex in res.results:
138            probs[tex.actual_class] += tex.weight
139            totweight += tex.weight
140        check_non_zero(totweight)
141    return [prob / totweight for prob in probs]
142
143
144@deprecated_keywords({
145    "foldN": "fold_n",
146    "reportSE": "report_se",
147    "iterationIsOuter": "iteration_is_outer"})
148def statistics_by_folds(stats, fold_n, report_se, iteration_is_outer):
149    # remove empty folds, turn the matrix so that learner is outer
150    if iteration_is_outer:
151        if not stats:
152            raise ValueError, "Cannot compute the score: no examples or sum of weights is 0."
153        number_of_learners = len(stats[0])
154        stats = filter(lambda (x, fN): fN > 0, zip(stats, fold_n))
155        stats = [[x[lrn] / fN for x, fN in stats]
156                 for lrn in range(number_of_learners)]
157    else:
158        stats = [[x / Fn for x, Fn in filter(lambda (x, Fn): Fn > 0,
159                 zip(lrnD, fold_n))] for lrnD in stats]
160
161    if not stats:
162        raise ValueError, "Cannot compute the score: no classifiers"
163    if not stats[0]:
164        raise ValueError, "Cannot compute the score: no examples or sum of weights is 0."
165
166    if report_se:
167        return [(statc.mean(x), statc.sterr(x)) for x in stats]
168    else:
169        return [statc.mean(x) for x in stats]
170
171def ME(res, **argkw):
172    MEs = [0.] * res.number_of_learners
173
174    if argkw.get("unweighted", 0) or not res.weights:
175        for tex in res.results:
176            MEs = map(lambda res, cls, ac=float(tex.actual_class):
177                      res + abs(float(cls) - ac), MEs, tex.classes)
178        totweight = gettotsize(res)
179    else:
180        for tex in res.results:
181            MEs = map(lambda res, cls, ac=float(tex.actual_class), tw=tex.weight:
182                       res + tw * abs(float(cls) - ac), MEs, tex.classes)
183        totweight = gettotweight(res)
184
185    return [x / totweight for x in MEs]
186
187MAE = ME
188
189
190class ConfusionMatrix:
191    """
192    Classification result summary.
193    """
194    #: True Positive predictions
195    TP = 0.
196    #:True Negative predictions
197    TN = 0.
198    #:False Positive predictions
199    FP = 0.
200    #: False Negative predictions
201    FN = 0.
202
203    @deprecated_keywords({"predictedPositive": "predicted_positive",
204                          "isPositive": "is_positive"})
205    def addTFPosNeg(self, predicted_positive, is_positive, weight=1.0):
206        """
207        Update confusion matrix with result of a single classification.
208
209        :param predicted_positive: positive class value was predicted
210        :param is_positive: correct class value is positive
211        :param weight: weight of the selected instance
212         """
213        if predicted_positive:
214            if is_positive:
215                self.TP += weight
216            else:
217                self.FP += weight
218        else:
219            if is_positive:
220                self.FN += weight
221            else:
222                self.TN += weight
223
224
225
226#########################################################################
227# PERFORMANCE MEASURES:
228# Scores for evaluation of numeric predictions
229
230def check_argkw(dct, lst):
231    """Return True if any item from lst has a non-zero value in dct."""
232    return reduce(lambda x, y: x or y, [dct.get(k, 0) for k in lst])
233
234def regression_error(res, **argkw):
235    """Return the regression error (default: MSE)."""
236    if argkw.get("SE", 0) and res.number_of_iterations > 1:
237        # computes the scores for each iteration, then averages
238        scores = [[0.] * res.number_of_iterations
239                  for _ in range(res.number_of_learners)]
240        norm = None
241        if argkw.get("norm-abs", 0) or argkw.get("norm-sqr", 0):
242            norm = [0.] * res.number_of_iterations
243
244        # counts examples in each iteration
245        nIter = [0] * res.number_of_iterations
246        # average class in each iteration
247        a = [0] * res.number_of_iterations
248        for tex in res.results:
249            nIter[tex.iteration_number] += 1
250            a[tex.iteration_number] += float(tex.actual_class)
251        a = [a[i] / nIter[i] for i in range(res.number_of_iterations)]
252
253        if argkw.get("unweighted", 0) or not res.weights:
254            # iterate accross test cases
255            for tex in res.results:
256                ai = float(tex.actual_class)
257                nIter[tex.iteration_number] += 1
258
259                # compute normalization, if required
260                if argkw.get("norm-abs", 0):
261                    norm[tex.iteration_number] += abs(ai - a[tex.iteration_number])
262                elif argkw.get("norm-sqr", 0):
263                    norm[tex.iteration_number] += (ai - a[tex.iteration_number]) ** 2
264
265                # iterate accross results of different regressors
266                for i, cls in enumerate(tex.classes):
267                    if argkw.get("abs", 0):
268                        scores[i][tex.iteration_number] += abs(float(cls) - ai)
269                    else:
270                        scores[i][tex.iteration_number] += (float(cls) - ai) ** 2
271        else: # unweighted != 0
272            raise NotImplementedError, "weighted error scores with SE not implemented yet"
273
274        if argkw.get("norm-abs") or argkw.get("norm-sqr"):
275            scores = [[x / n for x, n in zip(y, norm)] for y in scores]
276        else:
277            scores = [[x / ni for x, ni in zip(y, nIter)] for y in scores]
278
279        if argkw.get("R2"):
280            scores = [[1.0 - x for x in y] for y in scores]
281
282        if argkw.get("sqrt", 0):
283            scores = [[math.sqrt(x) for x in y] for y in scores]
284
285        return [(statc.mean(x), statc.std(x)) for x in scores]
286
287    else: # single iteration (testing on a single test set)
288        scores = [0.] * res.number_of_learners
289        norm = 0.
290
291        if argkw.get("unweighted", 0) or not res.weights:
292            a = sum([tex.actual_class for tex in res.results]) \
293                / len(res.results)
294            for tex in res.results:
295                if argkw.get("abs", 0):
296                    scores = map(lambda res, cls, ac=float(tex.actual_class):
297                                 res + abs(float(cls) - ac), scores, tex.classes)
298                else:
299                    scores = map(lambda res, cls, ac=float(tex.actual_class):
300                                 res + (float(cls) - ac) ** 2, scores, tex.classes)
301
302                if argkw.get("norm-abs", 0):
303                    norm += abs(tex.actual_class - a)
304                elif argkw.get("norm-sqr", 0):
305                    norm += (tex.actual_class - a) ** 2
306            totweight = gettotsize(res)
307        else:
308            # UNFINISHED
309            MSEs = [0.] * res.number_of_learners
310            for tex in res.results:
311                MSEs = map(lambda res, cls, ac=float(tex.actual_class),
312                           tw=tex.weight:
313                           res + tw * (float(cls) - ac) ** 2, MSEs, tex.classes)
314            totweight = gettotweight(res)
315
316        if argkw.get("norm-abs", 0) or argkw.get("norm-sqr", 0):
317            scores = [s / norm for s in scores]
318        else: # normalize by number of instances (or sum of weights)
319            scores = [s / totweight for s in scores]
320
321        if argkw.get("R2"):
322            scores = [1. - s for s in scores]
323
324        if argkw.get("sqrt", 0):
325            scores = [math.sqrt(x) for x in scores]
326
327        return scores
328
329def MSE(res, **argkw):
330    """Compute mean-squared error."""
331    return regression_error(res, **argkw)
332
333def RMSE(res, **argkw):
334    """Compute root mean-squared error."""
335    argkw.setdefault("sqrt", True)
336    return regression_error(res, **argkw)
337
338def MAE(res, **argkw):
339    """Compute mean absolute error."""
340    argkw.setdefault("abs", True)
341    return regression_error(res, **argkw)
342
343def RSE(res, **argkw):
344    """Compute relative squared error."""
345    argkw.setdefault("norm-sqr", True)
346    return regression_error(res, **argkw)
347
348def RRSE(res, **argkw):
349    """Compute relative squared error."""
350    argkw.setdefault("norm-sqr", True)
351    argkw.setdefault("sqrt", True)
352    return regression_error(res, **argkw)
353
354def RAE(res, **argkw):
355    """Compute relative absolute error."""
356    argkw.setdefault("abs", True)
357    argkw.setdefault("norm-abs", True)
358    return regression_error(res, **argkw)
359
360def R2(res, **argkw):
361    """Compute the coefficient of determination, R-squared."""
362    argkw.setdefault("norm-sqr", True)
363    argkw.setdefault("R2", True)
364    return regression_error(res, **argkw)
365
366def MSE_old(res, **argkw):
367    """Compute mean-squared error."""
368    if argkw.get("SE", 0) and res.number_of_iterations > 1:
369        MSEs = [[0.] * res.number_of_iterations
370                for _ in range(res.number_of_learners)]
371        nIter = [0] * res.number_of_iterations
372        if argkw.get("unweighted", 0) or not res.weights:
373            for tex in res.results:
374                ac = float(tex.actual_class)
375                nIter[tex.iteration_number] += 1
376                for i, cls in enumerate(tex.classes):
377                    MSEs[i][tex.iteration_number] += (float(cls) - ac) ** 2
378        else:
379            raise ValueError, "weighted RMSE with SE not implemented yet"
380        MSEs = [[x / ni for x, ni in zip(y, nIter)] for y in MSEs]
381        if argkw.get("sqrt", 0):
382            MSEs = [[math.sqrt(x) for x in y] for y in MSEs]
383        return [(statc.mean(x), statc.std(x)) for x in MSEs]
384
385    else:
386        MSEs = [0.] * res.number_of_learners
387        if argkw.get("unweighted", 0) or not res.weights:
388            for tex in res.results:
389                MSEs = map(lambda res, cls, ac=float(tex.actual_class):
390                           res + (float(cls) - ac) ** 2, MSEs, tex.classes)
391            totweight = gettotsize(res)
392        else:
393            for tex in res.results:
394                MSEs = map(lambda res, cls, ac=float(tex.actual_class),
395                           tw=tex.weight: res + tw * (float(cls) - ac) ** 2,
396                           MSEs, tex.classes)
397            totweight = gettotweight(res)
398
399        if argkw.get("sqrt", 0):
400            MSEs = [math.sqrt(x) for x in MSEs]
401        return [x / totweight for x in MSEs]
402
403def RMSE_old(res, **argkw):
404    """Compute root mean-squared error."""
405    argkw.setdefault("sqrt", 1)
406    return MSE_old(res, **argkw)
407
408#########################################################################
409# PERFORMANCE MEASURES:
410# Scores for evaluation of classifiers
411
412class CA(list):
413    """
414    Compute percentage of matches between predicted and actual class.
415
416    :param test_results: :obj:~Orange.evaluation.testing.ExperimentResults
417                         or list of :obj:ConfusionMatrix.
418    :param report_se: include standard error in result.
419    :param ignore_weights: ignore instance weights.
420    :rtype: list of scores, one for each learner.
421
422    Standard errors are estimated from deviation of CAs across folds (if
423    test_results were produced by cross_validation) or approximated under
424    the assumption of normal distribution otherwise.
425    """
426    CONFUSION_MATRIX = 0
427    CONFUSION_MATRIX_LIST = 1
428    CLASSIFICATION = 2
429    CROSS_VALIDATION = 3
430
431    @deprecated_keywords({"reportSE": "report_se",
432                          "unweighted": "ignore_weights"})
433    def __init__(self, test_results, report_se=False, ignore_weights=False):
434        super(CA, self).__init__()
435        self.report_se = report_se
436        self.ignore_weights = ignore_weights
437
438        input_type = self.get_input_type(test_results)
439        if input_type == self.CONFUSION_MATRIX:
440            self[:] = [self.from_confusion_matrix(test_results)]
441        elif input_type == self.CONFUSION_MATRIX_LIST:
442            self[:] = self.from_confusion_matrix_list(test_results)
443        elif input_type == self.CLASSIFICATION:
444            self[:] = self.from_classification_results(test_results)
445        elif input_type == self.CROSS_VALIDATION:
446            self[:] = self.from_crossvalidation_results(test_results)
447
448    def from_confusion_matrix(self, cm):
449        all_predictions = 0.
450        correct_predictions = 0.
451        if isinstance(cm, ConfusionMatrix):
452            all_predictions += cm.TP + cm.FN + cm.FP + cm.TN
453            correct_predictions += cm.TP + cm.TN
454        else:
455            for r, row in enumerate(cm):
456                for c, column in enumerate(row):
457                    if r == c:
458                        correct_predictions += column
459                    all_predictions += column
460
461        check_non_zero(all_predictions)
462        ca = correct_predictions / all_predictions
463
464        if self.report_se:
465            return ca, ca * (1 - ca) / math.sqrt(all_predictions)
466        else:
467            return ca
468
469    def from_confusion_matrix_list(self, confusion_matrices):
470        return [self.from_confusion_matrix(cm) for cm in confusion_matrices]
471
472    def from_classification_results(self, test_results):
473        CAs = [0.] * test_results.number_of_learners
474        totweight = 0.
475        for tex in test_results.results:
476            w = 1. if self.ignore_weights else tex.weight
477            CAs = map(lambda res, cls: res + (cls == tex.actual_class and w),
478                      CAs, tex.classes)
479            totweight += w
480        check_non_zero(totweight)
481        ca = [x / totweight for x in CAs]
482
483        if self.report_se:
484            return [(x, x * (1 - x) / math.sqrt(totweight)) for x in ca]
485        else:
486            return ca
487
488    def from_crossvalidation_results(self, test_results):
489        CAsByFold = [[0.] * test_results.number_of_iterations
490                     for _ in range(test_results.number_of_learners)]
491        foldN = [0.] * test_results.number_of_iterations
492
493        for tex in test_results.results:
494            w = 1. if self.ignore_weights else tex.weight
495            for lrn in range(test_results.number_of_learners):
496                CAsByFold[lrn][tex.iteration_number] += (tex.classes[lrn] ==
497                    tex.actual_class) and w
498            foldN[tex.iteration_number] += w
499
500        return statistics_by_folds(CAsByFold, foldN, self.report_se, False)
501
502    def get_input_type(self, test_results):
503        if isinstance(test_results, ConfusionMatrix):
504            return self.CONFUSION_MATRIX
505        elif isinstance(test_results, testing.ExperimentResults):
506            if test_results.number_of_iterations == 1:
507                return self.CLASSIFICATION
508            else:
509                return self.CROSS_VALIDATION
510        elif isinstance(test_results, list):
511            return self.CONFUSION_MATRIX_LIST
512
513
514@deprecated_keywords({"reportSE": "report_se",
515                      "unweighted": "ignore_weights"})
516def AP(res, report_se=False, ignore_weights=False, **argkw):
517    """Compute the average probability assigned to the correct class."""
518    if res.number_of_iterations == 1:
519        APs = [0.] * res.number_of_learners
520        if ignore_weights or not res.weights:
521            for tex in res.results:
522                APs = map(lambda res, probs: res + probs[tex.actual_class],
523                          APs, tex.probabilities)
524            totweight = gettotsize(res)
525        else:
526            totweight = 0.
527            for tex in res.results:
528                APs = map(lambda res, probs: res + probs[tex.actual_class] *
529                          tex.weight, APs, tex.probabilities)
530                totweight += tex.weight
531        check_non_zero(totweight)
532        return [AP / totweight for AP in APs]
533
534    APsByFold = [[0.] * res.number_of_learners
535                 for _ in range(res.number_of_iterations)]
536    foldN = [0.] * res.number_of_iterations
537    if ignore_weights or not res.weights:
538        for tex in res.results:
539            APsByFold[tex.iteration_number] = map(lambda res, probs:
540                res + probs[tex.actual_class],
541                APsByFold[tex.iteration_number], tex.probabilities)
542            foldN[tex.iteration_number] += 1
543    else:
544        for tex in res.results:
545            APsByFold[tex.iteration_number] = map(lambda res, probs:
546                res + probs[tex.actual_class] * tex.weight,
547                APsByFold[tex.iteration_number], tex.probabilities)
548            foldN[tex.iteration_number] += tex.weight
549
550    return statistics_by_folds(APsByFold, foldN, report_se, True)
551
552
553@deprecated_keywords({"reportSE": "report_se",
554                      "unweighted": "ignore_weights"})
555def Brier_score(res, report_se=False, ignore_weights=False, **argkw):
556    """Compute the Brier score, defined as the average (over test instances)
557    of :math:\sum_x(t(x) - p(x))^2, where x is a class value, t(x) is 1 for
558    the actual class value and 0 otherwise, and p(x) is the predicted
559    probability of x.
560    """
561    # Computes an average (over examples) of sum_x(t(x) - p(x))^2, where
562    #    x is class,
563    #    t(x) is 0 for 'wrong' and 1 for 'correct' class
564    #    p(x) is predicted probabilty.
565    # There's a trick: since t(x) is zero for all classes but the
566    # correct one (c), we compute the sum as sum_x(p(x)^2) - 2*p(c) + 1
567    # Since +1 is there for each example, it adds 1 to the average
568    # We skip the +1 inside the sum and add it just at the end of the function
569    # We take max(result, 0) to avoid -0.0000x due to rounding errors
570
571    if res.number_of_iterations == 1:
572        MSEs = [0.] * res.number_of_learners
573        if ignore_weights or not res.weights:
574            totweight = 0.
575            for tex in res.results:
576                MSEs = map(lambda res, probs: res + reduce(
577                    lambda s, pi: s + pi ** 2, probs, 0) -
578                    2 * probs[tex.actual_class], MSEs, tex.probabilities)
579                totweight += tex.weight
580        else:
581            for tex in res.results:
582                MSEs = map(lambda res, probs: res + tex.weight * reduce(
583                    lambda s, pi: s + pi ** 2, probs, 0) -
584                    2 * probs[tex.actual_class], MSEs, tex.probabilities)
585            totweight = gettotweight(res)
586        check_non_zero(totweight)
587        if report_se:
588            ## change this, not zero!!!
589            return [(max(x / totweight + 1., 0), 0) for x in MSEs]
590        else:
591            return [max(x / totweight + 1., 0) for x in MSEs]
592
593    BSs = [[0.] * res.number_of_learners
594           for _ in range(res.number_of_iterations)]
595    foldN = [0.] * res.number_of_iterations
596
597    if ignore_weights or not res.weights:
598        for tex in res.results:
599            BSs[tex.iteration_number] = map(lambda rr, probs: rr + reduce(
600                lambda s, pi: s + pi ** 2, probs, 0) -
601                2 * probs[tex.actual_class], BSs[tex.iteration_number],
602                tex.probabilities)
603            foldN[tex.iteration_number] += 1
604    else:
605        for tex in res.results:
606            BSs[tex.iteration_number] = map(lambda res, probs:
607                res + tex.weight * reduce(lambda s, pi: s + pi ** 2, probs, 0) -
608                2 * probs[tex. actual_class], BSs[tex.iteration_number],
609                tex.probabilities)
610            foldN[tex.iteration_number] += tex.weight
611
612    stats = statistics_by_folds(BSs, foldN, report_se, True)
613    if report_se:
614        return [(x + 1., y) for x, y in stats]
615    else:
616        return [x + 1. for x in stats]
617
618def BSS(res, **argkw):
619    return [1 - x / 2 for x in apply(Brier_score, (res,), argkw)]
620
621def IS_ex(Pc, P):
622    """Pc aposterior probability, P aprior"""
623    if Pc >= P:
624        return -log2(P) + log2(Pc)
625    else:
626        return -(-log2(1 - P) + log2(1 - Pc))
627
628
629@deprecated_keywords({"reportSE": "report_se"})
630def IS(res, apriori=None, report_se=False, **argkw):
631    """Compute the information score as defined by
632    Kononenko and Bratko (1991) \
634    Argument :obj:apriori gives the apriori class
635    distribution; if it is omitted, the class distribution is computed from
636    the actual classes of examples in :obj:res.
637    """
638    if not apriori:
639        apriori = class_probabilities_from_res(res)
640
641    if res.number_of_iterations == 1:
642        ISs = [0.] * res.number_of_learners
643        if argkw.get("unweighted", 0) or not res.weights:
644            for tex in res.results:
645              for i in range(len(tex.probabilities)):
646                    cls = tex.actual_class
647                    ISs[i] += IS_ex(tex.probabilities[i][cls], apriori[cls])
648            totweight = gettotsize(res)
649        else:
650            for tex in res.results:
651              for i in range(len(tex.probabilities)):
652                    cls = tex.actual_class
653                    ISs[i] += (IS_ex(tex.probabilities[i][cls], apriori[cls]) *
654                               tex.weight)
655            totweight = gettotweight(res)
656        if report_se:
657            return [(IS / totweight, 0) for IS in ISs]
658        else:
659            return [IS / totweight for IS in ISs]
660
661
662    ISs = [[0.] * res.number_of_iterations
663           for _ in range(res.number_of_learners)]
664    foldN = [0.] * res.number_of_iterations
665
666    # compute info scores for each fold
667    if argkw.get("unweighted", 0) or not res.weights:
668        for tex in res.results:
669            for i in range(len(tex.probabilities)):
670                cls = tex.actual_class
671                ISs[i][tex.iteration_number] += IS_ex(tex.probabilities[i][cls],
672                                                apriori[cls])
673            foldN[tex.iteration_number] += 1
674    else:
675        for tex in res.results:
676            for i in range(len(tex.probabilities)):
677                cls = tex.actual_class
678                ISs[i][tex.iteration_number] += IS_ex(tex.probabilities[i][cls],
679                                                apriori[cls]) * tex.weight
680            foldN[tex.iteration_number] += tex.weight
681
682    return statistics_by_folds(ISs, foldN, report_se, False)
683
684
685def Wilcoxon(res, statistics, **argkw):
686    res1, res2 = [], []
687    for ri in split_by_iterations(res):
688        stats = statistics(ri, **argkw)
689        if len(stats) != 2:
690            raise TypeError, "Wilcoxon compares two classifiers, no more, no less"
691        res1.append(stats[0])
692        res2.append(stats[1])
693    return statc.wilcoxont(res1, res2)
694
695def rank_difference(res, statistics, **argkw):
696    if not res.results:
697        raise TypeError, "no experiments"
698
699    k = len(res.results[0].classes)
700    if k < 2:
701        raise TypeError, "nothing to compare (less than two classifiers given)"
702    if k == 2:
703        return apply(Wilcoxon, (res, statistics), argkw)
704    else:
705        return apply(Friedman, (res, statistics), argkw)
706
707
708@deprecated_keywords({"res": "test_results",
709                      "classIndex": "class_index",
710                      "unweighted": "ignore_weights"})
711def confusion_matrices(test_results, class_index= -1,
712                       ignore_weights=False, cutoff=.5):
713    """
714    Return confusion matrices for test_results.
715
716    :param test_results: test results
717    :param class_index: index of class value for which the confusion matrices
718        are to be computed (by default unspecified - see note below).
719    :param ignore_weights: ignore instance weights.
720    :param cutoff: cutoff for probability
721
722    :rtype: list of :obj:ConfusionMatrix or list-of-list-of-lists (see
723        note below)
724
725    .. note:: If class_index is not specified and test_results
726        contain predictions for multi-class problem, then the return
727        value is a list of 2D tables (list-of-lists) of all class
728        value pairwise misclassifications.
729
730    """
731    tfpns = [ConfusionMatrix() for _ in range(test_results.number_of_learners)]
732
733    if class_index < 0:
734        numberOfClasses = len(test_results.class_values)
735        if class_index < -1 or numberOfClasses > 2:
736            cm = [[[0.] * numberOfClasses for _ in range(numberOfClasses)]
737                  for _ in range(test_results.number_of_learners)]
738            if ignore_weights or not test_results.weights:
739                for tex in test_results.results:
740                    trueClass = int(tex.actual_class)
741                    for li, pred in enumerate(tex.classes):
742                        predClass = int(pred)
743                        if predClass < numberOfClasses:
744                            cm[li][trueClass][predClass] += 1
745            else:
746                for tex in test_results.results:
747                    trueClass = int(tex.actual_class)
748                    for li, pred in tex.classes:
749                        predClass = int(pred)
750                        if predClass < numberOfClasses:
751                            cm[li][trueClass][predClass] += tex.weight
752            return cm
753
754        elif test_results.base_class >= 0:
755            class_index = test_results.base_class
756        else:
757            class_index = 1
758
759    if cutoff != .5:
760        if ignore_weights or not test_results.weights:
761            for lr in test_results.results:
762                isPositive = (lr.actual_class == class_index)
763                for i in range(test_results.number_of_learners):
765                                         cutoff, isPositive)
766        else:
767            for lr in test_results.results:
768                isPositive = (lr.actual_class == class_index)
769                for i in range(test_results.number_of_learners):
771                                         cutoff, isPositive, lr.weight)
772    else:
773        if ignore_weights or not test_results.weights:
774            for lr in test_results.results:
775                isPositive = (lr.actual_class == class_index)
776                for i in range(test_results.number_of_learners):
777                    tfpns[i].addTFPosNeg(lr.classes[i] == class_index,
778                                         isPositive)
779        else:
780            for lr in test_results.results:
781                isPositive = (lr.actual_class == class_index)
782                for i in range(test_results.number_of_learners):
783                    tfpns[i].addTFPosNeg(lr.classes[i] == class_index,
784                                         isPositive, lr.weight)
785    return tfpns
786
787
788# obsolete (renamed)
789compute_confusion_matrices = confusion_matrices
790
791
792@deprecated_keywords({"confusionMatrix": "confusion_matrix"})
793def confusion_chi_square(confusion_matrix):
794    """
795    Return chi square statistic of the confusion matrix
796    (higher value indicates that prediction is not by chance).
797    """
798    if isinstance(confusion_matrix, ConfusionMatrix) or \
799       not isinstance(confusion_matrix[1], list):
800        return _confusion_chi_square(confusion_matrix)
801    else:
802        return map(_confusion_chi_square, confusion_matrix)
803
804def _confusion_chi_square(confusion_matrix):
805    if isinstance(confusion_matrix, ConfusionMatrix):
806        c = confusion_matrix
807        confusion_matrix = [[c.TP, c.FN], [c.FP, c.TN]]
808    dim = len(confusion_matrix)
809    rowPriors = [sum(r) for r in confusion_matrix]
810    colPriors = [sum(r[i] for r in confusion_matrix) for i in range(dim)]
811    total = sum(rowPriors)
812    rowPriors = [r / total for r in rowPriors]
813    colPriors = [r / total for r in colPriors]
814    ss = 0
815    for ri, row in enumerate(confusion_matrix):
816        for ci, o in enumerate(row):
817            e = total * rowPriors[ri] * colPriors[ci]
818            if not e:
819                return -1, -1, -1
820            ss += (o - e) ** 2 / e
821    df = (dim - 1) ** 2
822    return ss, df, statc.chisqprob(ss, df)
823
824class CMScore(list):
825    """
826    :param test_results: :obj:~Orange.evaluation.testing.ExperimentResults
827                         or list of :obj:ConfusionMatrix.
828    :rtype: list of scores, one for each learner."""
829    def __new__(cls, test_results, **kwargs):
830        self = list.__new__(cls)
831        if isinstance(test_results, ConfusionMatrix):
832            self.__init__(test_results, **kwargs)
833            return self[0]
834        return self
835
836
837    @deprecated_keywords({"confm": "test_results"})
838    def __init__(self, test_results=None):
839        super(CMScore, self).__init__()
840
841        if test_results is not None:
842            self[:] = self.__call__(test_results)
843
844    def __call__(self, test_results):
845        if isinstance(test_results, testing.ExperimentResults):
846            test_results = confusion_matrices(test_results, class_index=1)
847        if isinstance(test_results, ConfusionMatrix):
848            test_results = [test_results]
849
850        return map(self.compute, test_results)
851
852
853
854class Sensitivity(CMScore):
855    __doc__ = """Compute sensitivity
856    <http://en.wikipedia.org/wiki/Sensitivity_and_specificity>_ (proportion
857    of actual positives which are correctly identified as such).
858    """ + CMScore.__doc__
859    @classmethod
860    def compute(self, confusion_matrix):
861        tot = confusion_matrix.TP + confusion_matrix.FN
862        if tot < 1e-6:
863            import warnings
864            warnings.warn("Can't compute sensitivity: one or both classes have no instances")
865            return None
866
867        return confusion_matrix.TP / tot
868
869
870class Recall(Sensitivity):
871    __doc__ = """ Compute recall
872    <http://en.wikipedia.org/wiki/Precision_and_recall>_
873    (fraction of relevant instances that are retrieved).
874    """ + CMScore.__doc__
875    pass # Recall == Sensitivity
876
877
878class Specificity(CMScore):
879    __doc__ = """Compute specificity
880    <http://en.wikipedia.org/wiki/Sensitivity_and_specificity>_
881    (proportion of negatives which are correctly identified).
882    """ + CMScore.__doc__
883    @classmethod
884    def compute(self, confusion_matrix):
885        tot = confusion_matrix.FP + confusion_matrix.TN
886        if tot < 1e-6:
887            import warnings
888            warnings.warn("Can't compute specificity: one or both classes have no instances")
889            return None
890        return confusion_matrix.TN / tot
891
892
893class PPV(CMScore):
894    __doc__ = """Compute positive predictive value
895    <http://en.wikipedia.org/wiki/Positive_predictive_value>_ (proportion of
896    subjects with positive test results who are correctly diagnosed).
897    """ + CMScore.__doc__
898    @classmethod
899    def compute(self, confusion_matrix):
900        tot = confusion_matrix.TP + confusion_matrix.FP
901        if tot < 1e-6:
902            import warnings
903            warnings.warn("Can't compute PPV: one or both classes have no instances")
904            return None
905        return confusion_matrix.TP / tot
906
907
908class Precision(PPV):
909    __doc__ = """Compute precision <http://en.wikipedia.org/wiki/Precision_and_recall>_
910    (retrieved instances that are relevant).
911    """ + CMScore.__doc__
912    pass # Precision == PPV
913
914
915class NPV(CMScore):
916    __doc__ = """Compute negative predictive value
917    <http://en.wikipedia.org/wiki/Negative_predictive_value>_ (proportion of
918    subjects with a negative test result who are correctly diagnosed).
919     """ + CMScore.__doc__
920    @classmethod
921    def compute(self, confusion_matrix):
922        tot = confusion_matrix.FN + confusion_matrix.TN
923        if tot < 1e-6:
924            import warnings
925            warnings.warn("Can't compute NPV: one or both classes have no instances")
926            return None
927        return confusion_matrix.TN / tot
928
929
930class F1(CMScore):
931    __doc__ = """Return F1 score
932    <http://en.wikipedia.org/wiki/F1_score>_
933    (harmonic mean of precision and recall).
934    """ + CMScore.__doc__
935    @classmethod
936    def compute(self, confusion_matrix):
937        p = Precision.compute(confusion_matrix)
938        r = Recall.compute(confusion_matrix)
939        if p is not None and r is not None and (p + r) != 0:
940            return 2. * p * r / (p + r)
941        else:
942            import warnings
943            warnings.warn("Can't compute F1: P + R is zero or not defined")
944            return None
945
946
947class Falpha(CMScore):
948    __doc__ = """Compute the alpha-mean of precision and recall over the given confusion
949    matrix.
950    """ + CMScore.__doc__
951
952    def __init__(self, test_results, alpha=1.):
953        self.alpha = alpha
954        super(Falpha, self).__init__(test_results)
955
956    def compute(self, confusion_matrix):
957        p = Precision.compute(confusion_matrix)
958        r = Recall.compute(confusion_matrix)
959        return (1. + self.alpha) * p * r / (self.alpha * p + r)
960
961
962class MCC(CMScore):
963    __doc__ = """Compute Matthew correlation coefficient
964    <http://en.wikipedia.org/wiki/Matthews_correlation_coefficient>_
965    (correlation coefficient between the observed and predicted binary
966    classifications).
967    """ + CMScore.__doc__
968    @classmethod
969    def compute(self, cm):
970        # code by Boris Gorelik
971        TP, TN, FP, FN = cm.TP, cm.TN, cm.FP, cm.FN
972
973        try:
974            return (TP * TN - FP * FN) / \
975                 math.sqrt((TP + FP) * (TP + FN) * (TN + FP) * (TN + FN))
976        except ZeroDivisionError:
977            # Zero division occurs when there is either no true positives
978            # or no true negatives i.e. the problem contains only one
979            # type of classes.
980            import warnings
981            warnings.warn("Can't compute MCC: TP or TN is zero or not defined")
982
983
984@deprecated_keywords({"bIsListOfMatrices": "b_is_list_of_matrices"})
985def scotts_pi(confusion_matrix, b_is_list_of_matrices=True):
986   """Compute Scott's Pi for measuring inter-rater agreement for nominal data
987
988   http://en.wikipedia.org/wiki/Scott%27s_Pi
989   Scott's Pi is a statistic for measuring inter-rater reliability for nominal
990   raters.
991
992   @param confusion_matrix: confusion matrix, or list of confusion matrices.
993                            To obtain non-binary confusion matrix, call
994                            Orange.evaluation.scoring.confusion_matrices
995                            and set the class_index parameter to -2.
996   @param b_is_list_of_matrices: specifies whether confm is list of matrices.
997                            This function needs to operate on non-binary
998                            confusion matrices, which are represented by python
999                            lists, therefore one needs a way to distinguish
1000                            between a single matrix and list of matrices
1001   """
1002
1003   if b_is_list_of_matrices:
1004       try:
1005           return [scotts_pi(cm, b_is_list_of_matrices=False)
1006                   for cm in confusion_matrix]
1007       except TypeError:
1008           # Nevermind the parameter, maybe this is a "conventional" binary
1009           # confusion matrix and bIsListOfMatrices was specified by mistake
1010           return scottsPiSingle(confusion_matrix, bIsListOfMatrices=False)
1011   else:
1012       if isinstance(confusion_matrix, ConfusionMatrix):
1013           confusion_matrix = numpy.array([[confusion_matrix.TP,
1014               confusion_matrix.FN], [confusion_matrix.FP,
1015               confusion_matrix.TN]], dtype=float)
1016       else:
1017           confusion_matrix = numpy.array(confusion_matrix, dtype=float)
1018
1019       marginalSumOfRows = numpy.sum(confusion_matrix, axis=0)
1020       marginalSumOfColumns = numpy.sum(confusion_matrix, axis=1)
1021       jointProportion = (marginalSumOfColumns + marginalSumOfRows) / \
1022                           (2. * numpy.sum(confusion_matrix))
1023       # In the eq. above, 2. is what the Wikipedia page calls
1024       # the number of annotators. Here we have two annotators:
1025       # the observed (true) labels (annotations) and the predicted by
1026       # the learners.
1027
1028       prExpected = numpy.sum(jointProportion ** 2)
1029       prActual = numpy.sum(numpy.diag(confusion_matrix)) / \
1030                  numpy.sum(confusion_matrix)
1031
1032       ret = (prActual - prExpected) / (1.0 - prExpected)
1033       return ret
1034
1035# Backward compatibility
1036sens = Sensitivity
1037spec = Specificity
1038precision = Precision
1039recall = Recall
1040
1041
1042
1043@deprecated_keywords({"classIndex": "class_index",
1044                      "unweighted": "ignore_weights"})
1045@replace_discrete_probabilities_with_list(False)
1046def AUCWilcoxon(res, class_index= -1, ignore_weights=False, **argkw):
1047    """Compute the area under ROC (AUC) and its standard error using
1048    Wilcoxon's approach proposed by Hanley and McNeal (1982). If
1049    :obj:class_index is not specified, the first class is used as
1050    "the positive" and others are negative. The result is a list of
1051    tuples (aROC, standard error).
1052
1053    If test results consist of multiple folds, you need to split them using
1054    :obj:split_by_iterations and perform this test on each fold separately.
1055    """
1056    useweights = res.weights and not ignore_weights
1057    problists, tots = corn.computeROCCumulative(res, class_index, useweights)
1058
1059    results = []
1060
1061    totPos, totNeg = tots[1], tots[0]
1062    N = totPos + totNeg
1063    for plist in problists:
1064        highPos, lowNeg = totPos, 0.
1065        W, Q1, Q2 = 0., 0., 0.
1066        for prob in plist:
1067            thisPos, thisNeg = prob[1][1], prob[1][0]
1068            highPos -= thisPos
1069            W += thisNeg * (highPos + thisPos / 2.)
1070            Q2 += thisPos * (lowNeg ** 2 + lowNeg * thisNeg + thisNeg ** 2 / 3.)
1071            Q1 += thisNeg * (highPos ** 2 + highPos * thisPos + thisPos ** 2 / 3.)
1072
1073            lowNeg += thisNeg
1074
1075        W /= (totPos * totNeg)
1076        Q1 /= (totNeg * totPos ** 2)
1077        Q2 /= (totPos * totNeg ** 2)
1078
1079        SE = math.sqrt((W * (1 - W) + (totPos - 1) * (Q1 - W ** 2) +
1080                       (totNeg - 1) * (Q2 - W ** 2)) / (totPos * totNeg))
1081        results.append((W, SE))
1082    return results
1083
1084AROC = AUCWilcoxon # for backward compatibility, AROC is obsolete
1085
1086
1087@deprecated_keywords({"classIndex": "class_index",
1088                      "unweighted": "ignore_weights"})
1089@replace_discrete_probabilities_with_list(False)
1090def compare_2_AUCs(res, lrn1, lrn2, class_index= -1,
1091                   ignore_weights=False, **argkw):
1092    return corn.compare2ROCs(res, lrn1, lrn2, class_index,
1093                             res.weights and not ignore_weights)
1094
1095# for backward compatibility, compare_2_AROCs is obsolete
1096compare_2_AROCs = compare_2_AUCs
1097
1098
1099@deprecated_keywords({"classIndex": "class_index"})
1100@replace_discrete_probabilities_with_list(False)
1101def compute_ROC(res, class_index= -1):
1102    """Compute a ROC curve as a list of (x, y) tuples, where x is
1103    1-specificity and y is sensitivity.
1104    """
1105    problists, tots = corn.computeROCCumulative(res, class_index)
1106
1107    results = []
1108    totPos, totNeg = tots[1], tots[0]
1109
1110    for plist in problists:
1111        curve = [(1., 1.)]
1112        TP, TN = totPos, 0.
1113        FN, FP = 0., totNeg
1114        for prob in plist:
1115            thisPos, thisNeg = prob[1][1], prob[1][0]
1116            # thisPos go from TP to FN
1117            TP -= thisPos
1118            FN += thisPos
1119            # thisNeg go from FP to TN
1120            TN += thisNeg
1121            FP -= thisNeg
1122
1123            sens = TP / (TP + FN)
1124            spec = TN / (FP + TN)
1125            curve.append((1 - spec, sens))
1126        results.append(curve)
1127
1128    return results
1129
1130## TC's implementation of algorithms, taken from:
1131## T Fawcett: ROC Graphs: Notes and Practical Considerations for
1132## Data Mining Researchers, submitted to KDD Journal.
1133def ROC_slope((P1x, P1y, P1fscore), (P2x, P2y, P2fscore)):
1134    if P1x == P2x:
1135        return 1e300
1136    return (P1y - P2y) / (P1x - P2x)
1137
1138
1139@deprecated_keywords({"keepConcavities": "keep_concavities"})
1140def ROC_add_point(P, R, keep_concavities=1):
1141    if keep_concavities:
1142        R.append(P)
1143    else:
1144        while True:
1145            if len(R) < 2:
1146                R.append(P)
1147                return R
1148            else:
1149                T = R.pop()
1150                T2 = R[-1]
1151                if ROC_slope(T2, T) > ROC_slope(T, P):
1152                    R.append(T)
1153                    R.append(P)
1154                    return R
1155    return R
1156
1157
1158@deprecated_keywords({"classIndex": "class_index",
1159                      "keepConcavities": "keep_concavities"})
1160@replace_discrete_probabilities_with_list(False)
1161def TC_compute_ROC(res, class_index= -1, keep_concavities=1):
1162    problists, tots = corn.computeROCCumulative(res, class_index)
1163
1164    results = []
1165    P, N = tots[1], tots[0]
1166
1167    for plist in problists:
1168        ## corn gives an increasing by scores list, we need a decreasing
1169        plist.reverse()
1170        TP = 0.
1171        FP = 0.
1172        curve = []
1173        fPrev = 10e300 # "infinity" score at 0., 0.
1174        for prob in plist:
1175            f = prob[0]
1176            if f != fPrev:
1177                if P:
1178                    tpr = TP / P
1179                else:
1180                    tpr = 0.
1181                if N:
1182                    fpr = FP / N
1183                else:
1184                    fpr = 0.
1185                curve = ROC_add_point((fpr, tpr, fPrev), curve,
1186                                      keep_concavities)
1187                fPrev = f
1188            thisPos, thisNeg = prob[1][1], prob[1][0]
1189            TP += thisPos
1190            FP += thisNeg
1191        if P:
1192            tpr = TP / P
1193        else:
1194            tpr = 0.
1195        if N:
1196            fpr = FP / N
1197        else:
1198            fpr = 0.
1199        curve = ROC_add_point((fpr, tpr, f), curve, keep_concavities) ## ugly
1200        results.append(curve)
1201
1202    return results
1203
1204## returns a list of points at the intersection of the tangential
1205## iso-performance line and the given ROC curve
1206## for given values of FPcost, FNcost and pval
1207def TC_best_thresholds_on_ROC_curve(FPcost, FNcost, pval, curve):
1208    m = (FPcost * (1. - pval)) / (FNcost * pval)
1209
1210    ## put the iso-performance line in point (0., 1.)
1211    x0, y0 = (0., 1.)
1212    x1, y1 = (1., 1. + m)
1213    d01 = math.sqrt((x1 - x0) * (x1 - x0) + (y1 - y0) * (y1 - y0))
1214
1215    ## calculate and find the closest point to the line
1216    firstp = 1
1217    mind = 0.
1218    a = x0 * y1 - x1 * y0
1219    closestPoints = []
1220    for (x, y, fscore) in curve:
1221        d = ((y0 - y1) * x + (x1 - x0) * y + a) / d01
1222        d = abs(d)
1223        if firstp or d < mind:
1224            mind, firstp = d, 0
1225            closestPoints = [(x, y, fscore)]
1226        else:
1227            if abs(d - mind) <= 0.0001: ## close enough
1228                closestPoints.append((x, y, fscore))
1229    return closestPoints
1230
1231def frange(start, end=None, inc=None):
1232    """A range function, that does accept float increments..."""
1233
1234    if end is None:
1235        end = start + 0.
1236        start = 0.
1237
1238    if inc is None or inc == 0:
1239        inc = 1.
1240
1241    L = [start]
1242    while 1:
1243        next = start + len(L) * inc
1244        if inc > 0 and next >= end:
1245            L.append(end)
1246            break
1247        elif inc < 0 and next <= end:
1248            L.append(end)
1249            break
1250        L.append(next)
1251
1252    return L
1253
1254## input ROCcurves are of form [ROCcurves1, ROCcurves2, ... ROCcurvesN],
1255## where ROCcurvesX is a set of ROC curves,
1256## where a (one) ROC curve is a set of (FP, TP) points
1257##
1258## for each (sub)set of input ROC curves
1259## returns the average ROC curve and an array of (vertical) standard deviations
1260@deprecated_keywords({"ROCcurves": "roc_curves"})
1261def TC_vertical_average_ROC(roc_curves, samples=10):
1262    def INTERPOLATE((P1x, P1y, P1fscore), (P2x, P2y, P2fscore), X):
1263        if (P1x == P2x) or P1x < X > P2x or P1x > X < P2x:
1264            raise ValueError, "assumptions for interpolation are not met: P1 = %f,%f P2 = %f,%f X = %f" % (P1x, P1y, P2x, P2y, X)
1265        dx = float(P2x) - float(P1x)
1266        dy = float(P2y) - float(P1y)
1267        m = dy / dx
1268        return P1y + m * (X - P1x)
1269
1270    def TP_FOR_FP(FPsample, ROC, npts):
1271        i = 0
1272        while i < npts - 1:
1273            (fp, _, _) = ROC[i + 1]
1274            if fp <= FPsample:
1275                i += 1
1276            else:
1277                break
1278        (fp, tp, _) = ROC[i]
1279        if fp == FPsample:
1280            return tp
1281        elif fp < FPsample and i + 1 < len(ROC):
1282            return INTERPOLATE(ROC[i], ROC[i + 1], FPsample)
1283        elif fp < FPsample and i + 1 == len(ROC): # return the last
1284            return ROC[i][1]
1285        raise ValueError, "cannot compute: TP_FOR_FP in TC_vertical_average_ROC"
1286        #return 0.
1287
1288    average = []
1289    stdev = []
1290    for ROCS in roc_curves:
1291        npts = []
1292        for c in ROCS:
1293            npts.append(len(c))
1294        nrocs = len(ROCS)
1295
1296        TPavg = []
1297        TPstd = []
1298        for FPsample in frange(0., 1., 1. / samples):
1299            TPsum = []
1300            for i in range(nrocs):
1301                ##TPsum = TPsum + TP_FOR_FP(FPsample, ROCS[i], npts[i])
1302                TPsum.append(TP_FOR_FP(FPsample, ROCS[i], npts[i]))
1303            TPavg.append((FPsample, statc.mean(TPsum)))
1304            if len(TPsum) > 1:
1305                stdv = statc.std(TPsum)
1306            else:
1307                stdv = 0.
1308            TPstd.append(stdv)
1309
1310        average.append(TPavg)
1311        stdev.append(TPstd)
1312
1313    return average, stdev
1314
1315## input ROCcurves are of form [ROCcurves1, ROCcurves2, ... ROCcurvesN],
1316## where ROCcurvesX is a set of ROC curves,
1317## where a (one) ROC curve is a set of (FP, TP) points
1318##
1319## for each (sub)set of input ROC curves
1320## returns the average ROC curve, an array of vertical standard deviations and an array of horizontal standard deviations
1321@deprecated_keywords({"ROCcurves": "roc_curves"})
1322def TC_threshold_average_ROC(roc_curves, samples=10):
1323    def POINT_AT_THRESH(ROC, npts, thresh):
1324        i = 0
1325        while i < npts - 1:
1326            (px, py, pfscore) = ROC[i]
1327            if pfscore > thresh:
1328                i += 1
1329            else:
1330                break
1331        return ROC[i]
1332
1333    average = []
1334    stdevV = []
1335    stdevH = []
1336    for ROCS in roc_curves:
1337        npts = []
1338        for c in ROCS:
1339            npts.append(len(c))
1340        nrocs = len(ROCS)
1341
1342        T = []
1343        for c in ROCS:
1344            for (px, py, pfscore) in c:
1345##                try:
1346##                    T.index(pfscore)
1347##                except:
1348                T.append(pfscore)
1349        T.sort()
1350        T.reverse() ## ugly
1351
1352        TPavg = []
1353        TPstdV = []
1354        TPstdH = []
1355        for tidx in frange(0, (len(T) - 1.), float(len(T)) / samples):
1356            FPsum = []
1357            TPsum = []
1358            for i in range(nrocs):
1359                (fp, tp, _) = POINT_AT_THRESH(ROCS[i], npts[i], T[int(tidx)])
1360                FPsum.append(fp)
1361                TPsum.append(tp)
1362            TPavg.append((statc.mean(FPsum), statc.mean(TPsum)))
1363            ## vertical standard deviation
1364            if len(TPsum) > 1:
1365                stdv = statc.std(TPsum)
1366            else:
1367                stdv = 0.
1368            TPstdV.append(stdv)
1369            ## horizontal standard deviation
1370            if len(FPsum) > 1:
1371                stdh = statc.std(FPsum)
1372            else:
1373                stdh = 0.
1374            TPstdH.append(stdh)
1375
1376        average.append(TPavg)
1377        stdevV.append(TPstdV)
1378        stdevH.append(TPstdH)
1379
1380    return average, stdevV, stdevH
1381
1382## Calibration Curve
1383## returns an array of (curve, yesClassPredictions, noClassPredictions)
1384## elements, where:
1385##  - curve is an array of points (x, y) on the calibration curve
1386##  - yesClassRugPoints is an array of (x, 1) points
1387##  - noClassRugPoints is an array of (x, 0) points
1388@deprecated_keywords({"classIndex": "class_index"})
1389@replace_discrete_probabilities_with_list(False)
1390def compute_calibration_curve(res, class_index= -1):
1391    ## merge multiple iterations into one
1392    mres = Orange.evaluation.testing.ExperimentResults(1, res.classifier_names,
1393        res.class_values, res.weights, classifiers=res.classifiers,
1395    for te in res.results:
1396        mres.results.append(te)
1397
1398    problists, tots = corn.computeROCCumulative(mres, class_index)
1399
1400    results = []
1401
1402    bins = 10 ## divide interval between 0. and 1. into N bins
1403
1404    for plist in problists:
1405        yesClassRugPoints = []
1406        noClassRugPoints = []
1407
1408        yesBinsVals = [0] * bins
1409        noBinsVals = [0] * bins
1410        for (f, (thisNeg, thisPos)) in plist:
1411            yesClassRugPoints.append((f, thisPos)) # 1.
1412            noClassRugPoints.append((f, thisNeg)) # 1.
1413
1414            index = int(f * bins)
1415            index = min(index, bins - 1) ## just in case for value 1.
1416            yesBinsVals[index] += thisPos
1417            noBinsVals[index] += thisNeg
1418
1419        curve = []
1420        for cn in range(bins):
1421            f = float(cn * 1. / bins) + (1. / 2. / bins)
1422            yesVal = yesBinsVals[cn]
1423            noVal = noBinsVals[cn]
1424            allVal = yesVal + noVal
1425            if allVal == 0.: continue
1426            y = float(yesVal) / float(allVal)
1427            curve.append((f, y))
1428
1429        ## smooth the curve
1430        maxnPoints = 100
1431        if len(curve) >= 3:
1432#            loessCurve = statc.loess(curve, -3, 0.6)
1433            loessCurve = statc.loess(curve, maxnPoints, 0.5, 3)
1434        else:
1435            loessCurve = curve
1436        clen = len(loessCurve)
1437        if clen > maxnPoints:
1438            df = clen / maxnPoints
1439            if df < 1: df = 1
1440            curve = [loessCurve[i]  for i in range(0, clen, df)]
1441        else:
1442            curve = loessCurve
1443        ## remove the third value (variance of epsilon?) that suddenly
1444        ## appeared in the output of the statc.loess function
1445        curve = [c[:2] for c in curve]
1446        results.append((curve, yesClassRugPoints, noClassRugPoints))
1447
1448    return results
1449
1450
1451## Lift Curve
1452## returns an array of curve elements, where:
1453##  - curve is an array of points ((TP + FP) / (P + N), TP / P, (th, FP / N))
1454##    on the Lift Curve
1455@deprecated_keywords({"classIndex": "class_index"})
1456@replace_discrete_probabilities_with_list(False)
1457def compute_lift_curve(res, class_index= -1):
1458    ## merge multiple iterations into one
1459    mres = Orange.evaluation.testing.ExperimentResults(1, res.classifier_names,
1460        res.class_values, res.weights, classifiers=res.classifiers,
1462    for te in res.results:
1463        mres.results.append(te)
1464
1465    problists, tots = corn.computeROCCumulative(mres, class_index)
1466
1467    results = []
1468    P, N = tots[1], tots[0]
1469    for plist in problists:
1470        ## corn gives an increasing by scores list, we need a decreasing
1471        plist.reverse()
1472        TP = 0.
1473        FP = 0.
1474        curve = [(0., 0., (10e300, 0.))]
1475        for (f, (thisNeg, thisPos)) in plist:
1476            TP += thisPos
1477            FP += thisNeg
1478            curve.append(((TP + FP) / (P + N), TP, (f, FP / (N or 1))))
1479        results.append(curve)
1480
1481    return P, N, results
1482
1483
1484class CDT:
1485  """Store the number of concordant (C), discordant (D) and tied (T) pairs."""
1486  def __init__(self, C=0., D=0., T=0.):
1487    self.C, self.D, self.T = C, D, T
1488
1489def is_CDT_empty(cdt):
1490    return cdt.C + cdt.D + cdt.T < 1e-20
1491
1492
1493@deprecated_keywords({"classIndex": "class_index",
1494                      "unweighted": "ignore_weights"})
1495@replace_discrete_probabilities_with_list(False)
1496def compute_CDT(res, class_index= -1, ignore_weights=False, **argkw):
1497    """Obsolete, don't use."""
1498    if class_index < 0:
1499        if res.base_class >= 0:
1500            class_index = res.base_class
1501        else:
1502            class_index = 1
1503
1504    useweights = res.weights and not ignore_weights
1505    weightByClasses = argkw.get("weightByClasses", True)
1506
1507    if res.number_of_iterations > 1:
1508        CDTs = [CDT() for _ in range(res.number_of_learners)]
1509        iterationExperiments = split_by_iterations(res)
1510        for exp in iterationExperiments:
1511            expCDTs = corn.computeCDT(exp, class_index, useweights)
1512            for i in range(len(CDTs)):
1513                CDTs[i].C += expCDTs[i].C
1514                CDTs[i].D += expCDTs[i].D
1515                CDTs[i].T += expCDTs[i].T
1516        for i in range(res.number_of_learners):
1517            if is_CDT_empty(CDTs[0]):
1518                return corn.computeCDT(res, class_index, useweights)
1519
1520        return CDTs
1521    else:
1522        return corn.computeCDT(res, class_index, useweights)
1523
1524class AUC(list):
1525    """
1526    Compute the area under ROC curve given a set of experimental results.
1527    If testing consisted of multiple folds, each fold is scored and the
1528    average score is returned. If a fold contains only instances with the
1529    same class value, folds will be merged.
1530
1531    :param test_results: test results to score
1532    :param ignore_weights: ignore instance weights when calculating score
1533    :param multiclass: tells what kind of averaging to perform if the target
1534                       class has more than 2 values.
1535    """
1536
1537    #!Compute AUC for each pair of classes (ignoring instances of all other
1538    #!classes) and average the results, weighting them by the number of
1539    #!pairs of instances from these two classes (e.g. by the product of
1540    #!probabilities of the two classes). AUC computed in this way still
1541    #!behaves as the concordance index, e.g., gives the probability that two
1542    #!randomly chosen instances from different classes will be correctly
1543    #!recognized (if the classifier knows from which two classes the
1544    #!instances came).
1545    ByWeightedPairs = 0
1546
1547    #!Similar to ByWeightedPairs, except that the average over class pairs
1548    #!is not weighted. This AUC is, like the binary version, independent of
1549    #!class distributions, but it is not related to the concordance index
1550    #!any more.
1551    ByPairs = 1
1552
1553    #!For each class, it computes AUC for this class against all others (that
1554    #!is, treating other classes as one class). The AUCs are then averaged by
1555    #!the class probabilities. This is related to the concordance index in
1556    #!which we test the classifier's (average) capability of distinguishing
1557    #!the instances from a specified class from those that come from other
1558    #!classes.
1559    #!Unlike the binary AUC, the measure is not independent of class
1560    #!distributions.
1561    WeightedOneAgainstAll = 2
1562
1563    #!Similar to weighted_one_against_all, except that the average
1564    #!is not weighted.
1565    OneAgainstAll = 3
1566
1567    @replace_use_weights
1568    @deprecated_keywords({"method": "multiclass"})
1569    def __init__(self, test_results=None, multiclass=ByWeightedPairs, ignore_weights=False):
1570
1571        super(AUC, self).__init__()
1572
1573        self.ignore_weights = ignore_weights
1574        self.method = multiclass
1575
1576        if test_results is not None:
1577            r = self.__call__(test_results)
1578            if r == False: #when the test is invalid it can return a single False
1579                r = [ False ] * test_results.number_of_learners
1580            self[:] = r
1581
1582    @replace_discrete_probabilities_with_list(method=True)
1583    def __call__(self, test_results):
1584        if len(test_results.class_values) < 2:
1585            raise ValueError("Cannot compute AUC on a single-class problem")
1586        elif len(test_results.class_values) == 2:
1587            return self._compute_for_binary_class(test_results)
1588        else:
1589            return self._compute_for_multi_value_class(test_results, self.method)
1590
1591    def _compute_for_binary_class(self, res):
1592        """AUC for binary classification problems"""
1593        if res.number_of_iterations > 1:
1594            return self._compute_for_multiple_folds(
1595                self._compute_one_class_against_all,
1596                split_by_iterations(res),
1597                (-1, res, res.number_of_iterations))
1598        else:
1599            return self._compute_one_class_against_all(res, -1)[0]
1600
1601    def _compute_for_multi_value_class(self, res, method=0):
1602        """AUC for multiclass classification problems"""
1603        numberOfClasses = len(res.class_values)
1604
1605        if res.number_of_iterations > 1:
1606            iterations = split_by_iterations(res)
1607            all_ite = res
1608        else:
1609            iterations = [res]
1610            all_ite = None
1611
1612        # by pairs
1613        sum_aucs = [0.] * res.number_of_learners
1614        usefulClassPairs = 0.
1615
1616        prob = None
1617        if method in [self.ByWeightedPairs, self.WeightedOneAgainstAll]:
1618            prob = class_probabilities_from_res(res)
1619
1620        if method in [self.ByWeightedPairs, self.ByPairs]:
1621            for classIndex1 in range(numberOfClasses):
1622                for classIndex2 in range(classIndex1):
1623                    subsum_aucs = self._compute_for_multiple_folds(
1624                        self._compute_one_class_against_another, iterations,
1625                        (classIndex1, classIndex2, all_ite,
1626                        res.number_of_iterations))
1627                    if subsum_aucs:
1628                        if method == self.ByWeightedPairs:
1629                            p_ij = prob[classIndex1] * prob[classIndex2]
1630                            subsum_aucs = [x * p_ij  for x in subsum_aucs]
1631                            usefulClassPairs += p_ij
1632                        else:
1633                            usefulClassPairs += 1
1634                        sum_aucs = map(add, sum_aucs, subsum_aucs)
1635        else:
1636            for classIndex in range(numberOfClasses):
1637                subsum_aucs = self._compute_for_multiple_folds(
1638                    self._compute_one_class_against_all,
1639                    iterations, (classIndex, all_ite,
1640                                 res.number_of_iterations))
1641                if subsum_aucs:
1642                    if method == self.ByWeightedPairs:
1643                        p_i = prob[classIndex]
1644                        subsum_aucs = [x * p_i  for x in subsum_aucs]
1645                        usefulClassPairs += p_i
1646                    else:
1647                        usefulClassPairs += 1
1648                    sum_aucs = map(add, sum_aucs, subsum_aucs)
1649
1650        if usefulClassPairs > 0:
1651            sum_aucs = [x / usefulClassPairs for x in sum_aucs]
1652
1653        return sum_aucs
1654
1655    # computes the average AUC over folds using "AUCcomputer" (AUC_i or AUC_ij)
1656    # it returns the sum of what is returned by the computer,
1657    # unless at a certain fold the computer has to resort to computing
1658    # over all folds or even this failed;
1659    # in these cases the result is returned immediately
1660    def _compute_for_multiple_folds(self, auc_computer, iterations,
1661                                 computer_args):
1662        """Compute the average AUC over folds using :obj:auc_computer."""
1663        subsum_aucs = [0.] * iterations[0].number_of_learners
1664        for ite in iterations:
1665            aucs, foldsUsed = auc_computer(*(ite,) + computer_args)
1666            if not aucs:
1667                import warnings
1668                warnings.warn("AUC cannot be computed (all instances belong to the same class).")
1669                return
1670            if not foldsUsed:
1671                self[:] = aucs
1672                return aucs
1673            subsum_aucs = map(add, subsum_aucs, aucs)
1674        return subsum_aucs
1675
1676    # Computes AUC
1677    # in multivalued class problem, AUC is computed as one against all
1678    # results over folds are averages
1679    # if some folds examples from one class only, the folds are merged
1680    def _compute_for_single_class(self, res, class_index):
1681        if res.number_of_iterations > 1:
1682            return self._compute_for_multiple_folds(
1683                self._compute_one_class_against_all, split_by_iterations(res),
1684                (class_index, res, res.number_of_iterations))
1685        else:
1686            return self._compute_one_class_against_all(res, class_index)[0]
1687
1688    # Computes AUC for a pair of classes (as if there were no other classes)
1689    # results over folds are averages
1690    # if some folds have examples from one class only, the folds are merged
1691    def _compute_for_pair_of_classes(self, res, class_index1, class_index2):
1692        if res.number_of_iterations > 1:
1693            return self._compute_for_multiple_folds(
1694                self._compute_one_class_against_another,
1695                split_by_iterations(res),
1696                (class_index1, class_index2, res, res.number_of_iterations))
1697        else:
1698            return self._compute_one_class_against_another(res, class_index1,
1699                                                    class_index2)
1700
1701    # computes AUC between class i and the other classes
1702    # (treating them as the same class)
1703    @deprecated_keywords({"classIndex": "class_index",
1704                          "divideByIfIte": "divide_by_if_ite"})
1705    def _compute_one_class_against_all(self, ite, class_index, all_ite=None,
1706                                      divide_by_if_ite=1.0):
1707        """Compute AUC between class i and all the other classes)"""
1708        return self._compute_auc(corn.computeCDT, ite, all_ite,
1709            divide_by_if_ite, (class_index, not self.ignore_weights))
1710
1711
1712    # computes AUC between classes i and j as if there are no other classes
1713    def _compute_one_class_against_another(
1714        self, ite, class_index1, class_index2,
1715        all_ite=None, divide_by_if_ite=1.):
1716        """
1717        Compute AUC between classes i and j as if there are no other classes.
1718        """
1719        return self._compute_auc(corn.computeCDTPair, ite,
1720            all_ite, divide_by_if_ite,
1721            (class_index1, class_index2, not self.ignore_weights))
1722
1723    # computes AUC using a specified 'cdtComputer' function
1724    # It tries to compute AUCs from 'ite' (examples from a single iteration)
1725    # and, if C+D+T=0, from 'all_ite' (entire test set). In the former case,
1726    # the AUCs are divided by 'divideByIfIte'.
1727    # Additional flag is returned which is True in the former case,
1728    # or False in the latter.
1729    @deprecated_keywords({"cdt_computer": "cdtComputer",
1730                          "divideByIfIte": "divide_by_if_ite",
1731                          "computerArgs": "computer_args"})
1732    def _compute_auc(self, cdt_computer, ite, all_ite, divide_by_if_ite,
1733                     computer_args):
1734        """
1735        Compute AUC using a :obj:cdt_computer.
1736        """
1737        cdts = cdt_computer(*(ite,) + computer_args)
1738        if not is_CDT_empty(cdts[0]):
1739            return [(cdt.C + cdt.T / 2) / (cdt.C + cdt.D + cdt.T) /
1740                    divide_by_if_ite for cdt in cdts], True
1741
1742        if all_ite:
1743            cdts = cdt_computer(*(all_ite,) + computer_args)
1744            if not is_CDT_empty(cdts[0]):
1745                return [(cdt.C + cdt.T / 2) / (cdt.C + cdt.D + cdt.T)
1746                        for cdt in cdts], False
1747
1748        return False, False
1749
1750class AUC_for_single_class(AUC):
1751    """
1752    Compute AUC where the class with the given class_index is singled
1753    out and all other classes are treated as a single class.
1754    """
1755    def __init__(self, test_results=None, class_index= -1, ignore_weights=False):
1756        if class_index < 0:
1757            if test_results and test_results.base_class >= 0:
1758                self.class_index = test_results.base_class
1759            else:
1760                self.class_index = 1
1761        else:
1762            self.class_index = class_index
1763
1764        super(AUC_for_single_class, self).__init__(test_results, ignore_weights=ignore_weights)
1765
1766    @replace_discrete_probabilities_with_list(method=True)
1767    def __call__(self, test_results):
1768        return self._compute_for_single_class(test_results, self.class_index)
1769
1770
1771class AUC_for_pair_of_classes(AUC):
1772    """
1773    Computes AUC between a pair of classes, ignoring instances from all
1774    other classes.
1775    """
1776    def __init__(self, test_results, class_index1, class_index2, ignore_weights=False):
1777        self.class_index1 = class_index1
1778        self.class_index2 = class_index2
1779
1780        super(AUC_for_pair_of_classes, self).__init__(test_results, ignore_weights=ignore_weights)
1781
1782    @replace_discrete_probabilities_with_list(method=True)
1783    def __call__(self, test_results):
1784        return self._compute_for_pair_of_classes(test_results, self.class_index1, self.class_index2)
1785
1786
1787class AUC_matrix(AUC):
1788    """
1789    Compute a (lower diagonal) matrix with AUCs for all pairs of classes.
1790    If there are empty classes, the corresponding elements in the matrix
1791    are -1.
1792    """
1793
1794    @replace_discrete_probabilities_with_list(method=True)
1795    def __call__(self, test_results):
1796        numberOfClasses = len(test_results.class_values)
1797        number_of_learners = test_results.number_of_learners
1798        if test_results.number_of_iterations > 1:
1799            iterations, all_ite = split_by_iterations(test_results), test_results
1800        else:
1801            iterations, all_ite = [test_results], None
1802        aucs = [[[] for _ in range(numberOfClasses)]
1803        for _ in range(number_of_learners)]
1804        for classIndex1 in range(numberOfClasses):
1805            for classIndex2 in range(classIndex1):
1806                pair_aucs = self._compute_for_multiple_folds(
1807                    self._compute_one_class_against_another, iterations,
1808                    (classIndex1, classIndex2, all_ite,
1809                     test_results.number_of_iterations))
1810                if pair_aucs:
1811                    for lrn in range(number_of_learners):
1812                        aucs[lrn][classIndex1].append(pair_aucs[lrn])
1813                else:
1814                    for lrn in range(number_of_learners):
1815                        aucs[lrn][classIndex1].append(-1)
1816        return aucs
1817
1818#Backward compatibility
1819@replace_use_weights
1820@replace_discrete_probabilities_with_list(False)
1821def AUC_binary(res, ignore_weights=False):
1822    auc = deprecated_function_name(AUC)(ignore_weights=ignore_weights)
1823    auc._compute_for_binary_class(res)
1824    return auc
1825
1826@replace_use_weights
1827@replace_discrete_probabilities_with_list(False)
1828def AUC_multi(res, ignore_weights=False, method=0):
1829    auc = deprecated_function_name(AUC)(ignore_weights=ignore_weights,
1830        method=method)
1831    auc._compute_for_multi_value_class(res)
1832    return auc
1833
1834
1835@deprecated_keywords({"AUCcomputer": "auc_computer",
1836                      "computerArgs": "computer_args"})
1837def AUC_iterations(auc_computer, iterations, computer_args):
1838    auc = deprecated_function_name(AUC)()
1839    auc._compute_for_multiple_folds(auc_computer, iterations, computer_args)
1840    return auc
1841
1842def AUC_x(cdtComputer, ite, all_ite, divide_by_if_ite, computer_args):
1843    auc = deprecated_function_name(AUC)()
1844    result = auc._compute_auc(cdtComputer, ite, all_ite, divide_by_if_ite,
1845                              computer_args)
1846    return result
1847
1848@replace_use_weights
1849def AUC_i(ite, class_index, ignore_weights=False, all_ite=None,
1850          divide_by_if_ite=1.):
1851    auc = deprecated_function_name(AUC)(ignore_weights=ignore_weights)
1852    result = auc._compute_one_class_against_another(ite, class_index,
1853        all_ite=all_ite, divide_by_if_ite=divide_by_if_ite)
1854    return result
1855
1856
1857@replace_use_weights
1858def AUC_ij(ite, class_index1, class_index2, ignore_weights=False,
1859           all_ite=None, divide_by_if_ite=1.):
1860    auc = deprecated_function_name(AUC)(ignore_weights=ignore_weights)
1861    result = auc._compute_one_class_against_another(
1862        ite, class_index1, class_index2, all_ite=all_ite, divide_by_if_ite=divide_by_if_ite)
1863    return result
1864
1865AUC_single = replace_use_weights(
1866             deprecated_keywords({"classIndex": "class_index"})(
1867             deprecated_function_name(AUC_for_single_class)))
1868AUC_pair = replace_use_weights(
1869           deprecated_keywords({"classIndex1": "class_index1",
1870                                "classIndex2": "class_index2"})(
1871           deprecated_function_name(AUC_for_pair_of_classes)))
1872AUC_matrix = replace_use_weights(AUC_matrix)
1873
1874
1875@deprecated_keywords({"unweighted": "ignore_weights"})
1876@replace_discrete_probabilities_with_list(False)
1877def McNemar(res, ignore_weights=False, **argkw):
1878    """
1879    Compute a triangular matrix with McNemar statistics for each pair of
1880    classifiers. The statistics is distributed by chi-square distribution with
1881    one degree of freedom; critical value for 5% significance is around 3.84.
1882    """
1883    nLearners = res.number_of_learners
1884    mcm = []
1885    for i in range(nLearners):
1886       mcm.append([0.] * res.number_of_learners)
1887
1888    if not res.weights or ignore_weights:
1889        for i in res.results:
1890            actual = i.actual_class
1891            classes = i.classes
1892            for l1 in range(nLearners):
1893                for l2 in range(l1, nLearners):
1894                    if classes[l1] == actual:
1895                        if classes[l2] != actual:
1896                            mcm[l1][l2] += 1
1897                    elif classes[l2] == actual:
1898                        mcm[l2][l1] += 1
1899    else:
1900        for i in res.results:
1901            actual = i.actual_class
1902            classes = i.classes
1903            for l1 in range(nLearners):
1904                for l2 in range(l1, nLearners):
1905                    if classes[l1] == actual:
1906                        if classes[l2] != actual:
1907                            mcm[l1][l2] += i.weight
1908                    elif classes[l2] == actual:
1909                        mcm[l2][l1] += i.weight
1910
1911    for l1 in range(nLearners):
1912        for l2 in range(l1, nLearners):
1913            su = mcm[l1][l2] + mcm[l2][l1]
1914            if su:
1915                mcm[l2][l1] = (abs(mcm[l1][l2] - mcm[l2][l1]) - 1) ** 2 / su
1916            else:
1917                mcm[l2][l1] = 0
1918
1919    for l1 in range(nLearners):
1920        mcm[l1] = mcm[l1][:l1]
1921
1922    return mcm
1923
1924@replace_discrete_probabilities_with_list(False)
1925def McNemar_of_two(res, lrn1, lrn2, ignore_weights=False):
1926    """
1927    McNemar_of_two computes a McNemar statistics for a pair of classifier,
1928    specified by indices learner1 and learner2.
1929    """
1930    tf = ft = 0.
1931    if not res.weights or ignore_weights:
1932        for i in res.results:
1933            actual = i.actual_class
1934            if i.classes[lrn1] == actual:
1935                if i.classes[lrn2] != actual:
1936                    tf += i.weight
1937            elif i.classes[lrn2] == actual:
1938                    ft += i.weight
1939    else:
1940        for i in res.results:
1941            actual = i.actual_class
1942            if i.classes[lrn1] == actual:
1943                if i.classes[lrn2] != actual:
1944                    tf += 1.
1945            elif i.classes[lrn2] == actual:
1946                    ft += 1.
1947
1948    su = tf + ft
1949    if su:
1950        return (abs(tf - ft) - 1) ** 2 / su
1951    else:
1952        return 0
1953
1954@replace_discrete_probabilities_with_list(False)
1955def Friedman(res, stat=CA):
1956    """
1957    Compare classifiers with Friedman test, treating folds as different examles.
1958    Returns F, p and average ranks.
1959    """
1960    res_split = split_by_iterations(res)
1961    res = [stat(r) for r in res_split]
1962
1963    N = len(res)
1964    k = len(res[0])
1965    sums = [0.] * k
1966    for r in res:
1967        ranks = [k - x + 1 for x in statc.rankdata(r)]
1968        if stat == Brier_score: # reverse ranks for Brier_score (lower better)
1969            ranks = [k + 1 - x for x in ranks]
1970        sums = map(add, ranks, sums)
1971
1972    T = sum(x * x for x in sums)
1973    sums = [x / N for x in sums]
1974
1975    F = 12. / (N * k * (k + 1)) * T - 3 * N * (k + 1)
1976
1977    return F, statc.chisqprob(F, k - 1), sums
1978
1979@replace_discrete_probabilities_with_list(False)
1980def Wilcoxon_pairs(res, avgranks, stat=CA):
1981    """
1982    Return a triangular matrix, where element[i][j] stores significance of
1983    difference between the i-th and the j-th classifier, as computed by the
1984    Wilcoxon test. The element is positive if the i-th is better than the j-th,
1985    negative if it is worse, and 1 if they are equal.
1986    Arguments are ExperimentResults, average ranks (as returned by Friedman)
1987    and, optionally, a statistics; greater values should mean better results.
1988    """
1989    res_split = split_by_iterations(res)
1990    res = [stat(r) for r in res_split]
1991
1992    k = len(res[0])
1993    bt = []
1994    for m1 in range(k):
1995        nl = []
1996        for m2 in range(m1 + 1, k):
1997            t, p = statc.wilcoxont([r[m1] for r in res], [r[m2] for r in res])
1998            if avgranks[m1] < avgranks[m2]:
1999                nl.append(p)
2000            elif avgranks[m2] < avgranks[m1]:
2001                nl.append(-p)
2002            else:
2003                nl.append(1)
2004        bt.append(nl)
2005    return bt
2006
2007
2008@deprecated_keywords({"allResults": "all_results",
2009                      "noConfidence": "no_confidence"})
2010def plot_learning_curve_learners(file, all_results, proportions, learners,
2011                                 no_confidence=0):
2012    plot_learning_curve(file, all_results, proportions,
2013        [Orange.misc.getobjectname(learners[i], "Learner %i" % i)
2014        for i in range(len(learners))], no_confidence)
2015
2016
2017@deprecated_keywords({"allResults": "all_results",
2018                      "noConfidence": "no_confidence"})
2019def plot_learning_curve(file, all_results, proportions, legend,
2020                        no_confidence=0):
2021    import types
2022    fopened = 0
2023    if type(file) == types.StringType:
2024        file = open(file, "wt")
2025        fopened = 1
2026
2027    file.write("set yrange [0:1]\n")
2028    file.write("set xrange [%f:%f]\n" % (proportions[0], proportions[-1]))
2029    file.write("set multiplot\n\n")
2030    CAs = [CA(x, report_se=True) for x in all_results]
2031
2032    file.write("plot \\\n")
2033    for i in range(len(legend) - 1):
2034        if not no_confidence:
2035            file.write("'-' title '' with yerrorbars pointtype %i,\\\n" % (i + 1))
2036        file.write("'-' title '%s' with linespoints pointtype %i,\\\n" % (legend[i], i + 1))
2037    if not no_confidence:
2038        file.write("'-' title '' with yerrorbars pointtype %i,\\\n" % (len(legend)))
2039    file.write("'-' title '%s' with linespoints pointtype %i\n" % (legend[-1], len(legend)))
2040
2041    for i in range(len(legend)):
2042        if not no_confidence:
2043            for p in range(len(proportions)):
2044                file.write("%f\t%f\t%f\n" % (proportions[p], CAs[p][i][0], 1.96 * CAs[p][i][1]))
2045            file.write("e\n\n")
2046
2047        for p in range(len(proportions)):
2048            file.write("%f\t%f\n" % (proportions[p], CAs[p][i][0]))
2049        file.write("e\n\n")
2050
2051    if fopened:
2052        file.close()
2053
2054
2055def print_single_ROC_curve_coordinates(file, curve):
2056    import types
2057    fopened = 0
2058    if type(file) == types.StringType:
2059        file = open(file, "wt")
2060        fopened = 1
2061
2062    for coord in curve:
2063        file.write("%5.3f\t%5.3f\n" % tuple(coord))
2064
2065    if fopened:
2066        file.close()
2067
2068
2069def plot_ROC_learners(file, curves, learners):
2070    plot_ROC(file, curves, [Orange.misc.getobjectname(learners[i], "Learner %i" % i) for i in range(len(learners))])
2071
2072def plot_ROC(file, curves, legend):
2073    import types
2074    fopened = 0
2075    if type(file) == types.StringType:
2076        file = open(file, "wt")
2077        fopened = 1
2078
2079    file.write("set yrange [0:1]\n")
2080    file.write("set xrange [0:1]\n")
2081    file.write("set multiplot\n\n")
2082
2083    file.write("plot \\\n")
2084    for leg in legend:
2085        file.write("'-' title '%s' with lines,\\\n" % leg)
2086    file.write("'-' title '' with lines\n")
2087
2088    for curve in curves:
2089        for coord in curve:
2090            file.write("%5.3f\t%5.3f\n" % tuple(coord))
2091        file.write("e\n\n")
2092
2093    file.write("1.0\t1.0\n0.0\t0.0e\n\n")
2094
2095    if fopened:
2096        file.close()
2097
2098
2099@deprecated_keywords({"allResults": "all_results"})
2100def plot_McNemar_curve_learners(file, all_results, proportions, learners, reference= -1):
2101    plot_McNemar_curve(file, all_results, proportions, [Orange.misc.getobjectname(learners[i], "Learner %i" % i) for i in range(len(learners))], reference)
2102
2103
2104@deprecated_keywords({"allResults": "all_results"})
2105def plot_McNemar_curve(file, all_results, proportions, legend, reference= -1):
2106    if reference < 0:
2107        reference = len(legend) - 1
2108
2109    import types
2110    fopened = 0
2111    if type(file) == types.StringType:
2112        file = open(file, "wt")
2113        fopened = 1
2114
2115    #file.write("set yrange [0:1]\n")
2116    #file.write("set xrange [%f:%f]\n" % (proportions[0], proportions[-1]))
2117    file.write("set multiplot\n\n")
2118    file.write("plot \\\n")
2119    tmap = range(reference) + range(reference + 1, len(legend))
2120    for i in tmap[:-1]:
2121        file.write("'-' title '%s' with linespoints pointtype %i,\\\n" % (legend[i], i + 1))
2122    file.write("'-' title '%s' with linespoints pointtype %i\n" % (legend[tmap[-1]], tmap[-1]))
2123    file.write("\n")
2124
2125    for i in tmap:
2126        for p in range(len(proportions)):
2127            file.write("%f\t%f\n" % (proportions[p], McNemar_of_two(all_results[p], i, reference)))
2128        file.write("e\n\n")
2129
2130    if fopened:
2131        file.close()
2132
2133default_point_types = ("{$\\circ$}", "{$\\diamond$}", "{$+$}", "{$\\times$}", "{$|$}") + tuple([chr(x) for x in range(97, 122)])
2134default_line_types = ("\\setsolid", "\\setdashpattern <4pt, 2pt>", "\\setdashpattern <8pt, 2pt>", "\\setdashes", "\\setdots")
2135
2136@deprecated_keywords({"allResults": "all_results"})
2137def learning_curve_learners_to_PiCTeX(file, all_results, proportions, **options):
2138    return apply(learning_curve_to_PiCTeX, (file, all_results, proportions), options)
2139
2140
2141@deprecated_keywords({"allResults": "all_results"})
2142def learning_curve_to_PiCTeX(file, all_results, proportions, **options):
2143    import types
2144    fopened = 0
2145    if type(file) == types.StringType:
2146        file = open(file, "wt")
2147        fopened = 1
2148
2149    nexamples = len(all_results[0].results)
2150    CAs = [CA(x, report_se=True) for x in all_results]
2151
2152    graphsize = float(options.get("graphsize", 10.0)) #cm
2153    difprop = proportions[-1] - proportions[0]
2154    ntestexamples = nexamples * proportions[-1]
2155    xunit = graphsize / ntestexamples
2156
2157    yshift = float(options.get("yshift", -ntestexamples / 20.))
2158
2159    pointtypes = options.get("pointtypes", default_point_types)
2160    linetypes = options.get("linetypes", default_line_types)
2161
2162    if options.has_key("numberedx"):
2163        numberedx = options["numberedx"]
2164        if type(numberedx) == types.IntType:
2165            if numberedx > 0:
2166                numberedx = [nexamples * proportions[int(i / float(numberedx) * len(proportions))] for i in range(numberedx)] + [proportions[-1] * nexamples]
2167            elif numberedx < 0:
2168                numberedx = -numberedx
2169                newn = []
2170                for i in range(numberedx + 1):
2171                    wanted = proportions[0] + float(i) / numberedx * difprop
2172                    best = (10, 0)
2173                    for t in proportions:
2174                        td = abs(wanted - t)
2175                        if td < best[0]:
2176                            best = (td, t)
2177                    if not best[1] in newn:
2178                        newn.append(best[1])
2179                newn.sort()
2180                numberedx = [nexamples * x for x in newn]
2181        elif type(numberedx[0]) == types.FloatType:
2182            numberedx = [nexamples * x for x in numberedx]
2183    else:
2184        numberedx = [nexamples * x for x in proportions]
2185
2186    file.write("\\mbox{\n")
2187    file.write(\\beginpicture\n")
2188    file.write(\\setcoordinatesystem units <%10.8fcm, %5.3fcm>\n\n" % (xunit, graphsize))
2189    file.write(\\setplotarea x from %5.3f to %5.3f, y from 0 to 1\n" % (0, ntestexamples))
2190    file.write(\\axis bottom invisible\n")# label {#examples}\n")
2191    file.write("      ticks short at %s /\n" % reduce(lambda x, y:x + " " + y, ["%i" % (x * nexamples + 0.5) for x in proportions]))
2192    if numberedx:
2193        file.write("            long numbered at %s /\n" % reduce(lambda x, y:x + y, ["%i " % int(x + 0.5) for x in numberedx]))
2194    file.write("  /\n")
2195    file.write(\\axis left invisible\n")# label {classification accuracy}\n")
2196    file.write("      shiftedto y=%5.3f\n" % yshift)
2197    file.write("      ticks short from 0.0 to 1.0 by 0.05\n")
2198    file.write("            long numbered from 0.0 to 1.0 by 0.25\n")
2199    file.write("  /\n")
2200    if options.has_key("default"):
2201        file.write(\\setdashpattern<1pt, 1pt>\n")
2202        file.write(\\plot %5.3f %5.3f %5.3f %5.3f /\n" % (0., options["default"], ntestexamples, options["default"]))
2203
2204    for i in range(len(CAs[0])):
2205        coordinates = reduce(lambda x, y:x + " " + y, ["%i %5.3f" % (proportions[p] * nexamples, CAs[p][i][0]) for p in range(len(proportions))])
2206        if linetypes:
2207            file.write(%s\n" % linetypes[i])
2208            file.write(\\plot %s /\n" % coordinates)
2209        if pointtypes:
2210            file.write(\\multiput %s at %s /\n" % (pointtypes[i], coordinates))
2211
2212    file.write(\\endpicture\n")
2213    file.write("}\n")
2214    if fopened:
2215        file.close()
2216    file.close()
2217    del file
2218
2219def legend_learners_to_PiCTeX(file, learners, **options):
2220  return apply(legend_to_PiCTeX, (file, [Orange.misc.getobjectname(learners[i], "Learner %i" % i) for i in range(len(learners))]), options)
2221
2222def legend_to_PiCTeX(file, legend, **options):
2223    import types
2224    fopened = 0
2225    if type(file) == types.StringType:
2226        file = open(file, "wt")
2227        fopened = 1
2228
2229    pointtypes = options.get("pointtypes", default_point_types)
2230    linetypes = options.get("linetypes", default_line_types)
2231
2232    file.write("\\mbox{\n")
2233    file.write(\\beginpicture\n")
2234    file.write(\\setcoordinatesystem units <5cm, 1pt>\n\n")
2235    file.write(\\setplotarea x from 0.000 to %5.3f, y from 0 to 12\n" % len(legend))
2236
2237    for i in range(len(legend)):
2238        if linetypes:
2239            file.write(%s\n" % linetypes[i])
2240            file.write(\\plot %5.3f 6 %5.3f 6 /\n" % (i, i + 0.2))
2241        if pointtypes:
2242            file.write(\\put {%s} at %5.3f 6\n" % (pointtypes[i], i + 0.1))
2243        file.write(\\put {%s} [lb] at %5.3f 0\n" % (legend[i], i + 0.25))
2244
2245    file.write(\\endpicture\n")
2246    file.write("}\n")
2247    if fopened:
2248        file.close()
2249    file.close()
2250    del file
2251
2252
2253def compute_friedman(avranks, N):
2254    """ Returns a tuple composed of (friedman statistic, degrees of freedom)
2255    and (Iman statistic - F-distribution, degrees of freedoma) given average
2256    ranks and a number of tested data sets N.
2257    """
2258
2259    k = len(avranks)
2260
2261    def friedman(N, k, ranks):
2262        return 12 * N * (sum([rank ** 2.0 for rank in ranks]) - (k * (k + 1) * (k + 1) / 4.0)) / (k * (k + 1))
2263
2264    def iman(fried, N, k):
2265        return (N - 1) * fried / (N * (k - 1) - fried)
2266
2267    f = friedman(N, k, avranks)
2268    im = iman(f, N, k)
2269    fdistdof = (k - 1, (k - 1) * (N - 1))
2270
2271    return (f, k - 1), (im, fdistdof)
2272
2273def compute_CD(avranks, N, alpha="0.05", type="nemenyi"):
2274    """ Returns critical difference for Nemenyi or Bonferroni-Dunn test
2275    according to given alpha (either alpha="0.05" or alpha="0.1") for average
2276    ranks and number of tested data sets N. Type can be either "nemenyi" for
2277    for Nemenyi two tailed test or "bonferroni-dunn" for Bonferroni-Dunn test.
2278    """
2279
2280    k = len(avranks)
2281
2282    d = {("nemenyi", "0.05"): [0, 0, 1.959964, 2.343701, 2.569032, 2.727774,
2283                               2.849705, 2.94832, 3.030879, 3.101730, 3.163684,
2284                               3.218654, 3.268004, 3.312739, 3.353618, 3.39123,
2285                               3.426041, 3.458425, 3.488685, 3.517073, 3.543799]
2286        , ("nemenyi", "0.1"): [0, 0, 1.644854, 2.052293, 2.291341, 2.459516,
2287                               2.588521, 2.692732, 2.779884, 2.854606, 2.919889,
2288                               2.977768, 3.029694, 3.076733, 3.119693, 3.159199,
2289                               3.195743, 3.229723, 3.261461, 3.291224, 3.319233]
2290        , ("bonferroni-dunn", "0.05"): [0, 0, 1.960, 2.241, 2.394, 2.498, 2.576,
2291                                        2.638, 2.690, 2.724, 2.773],
2292         ("bonferroni-dunn", "0.1"): [0, 0, 1.645, 1.960, 2.128, 2.241, 2.326,
2293                                      2.394, 2.450, 2.498, 2.539]}
2294
2295    #can be computed in R as qtukey(0.95, n, Inf)**0.5
2296    #for (x in c(2:20)) print(qtukey(0.95, x, Inf)/(2**0.5)
2297
2298    q = d[(type, alpha)]
2299
2300    cd = q[k] * (k * (k + 1) / (6.0 * N)) ** 0.5
2301
2302    return cd
2303
2304
2305def graph_ranks(filename, avranks, names, cd=None, cdmethod=None, lowv=None, highv=None, width=6, textspace=1, reverse=False, **kwargs):
2306    """
2307    Draws a CD graph, which is used to display  the differences in methods'
2308    performance.
2309    See Janez Demsar, Statistical Comparisons of Classifiers over
2310    Multiple Data Sets, 7(Jan):1--30, 2006.
2311
2312    Needs matplotlib to work.
2313
2314    :param filename: Output file name (with extension). Formats supported
2315                     by matplotlib can be used.
2316    :param avranks: List of average methods' ranks.
2317    :param names: List of methods' names.
2318
2319    :param cd: Critical difference. Used for marking methods that whose
2320               difference is not statistically significant.
2321    :param lowv: The lowest shown rank, if None, use 1.
2322    :param highv: The highest shown rank, if None, use len(avranks).
2323    :param width: Width of the drawn figure in inches, default 6 in.
2324    :param textspace: Space on figure sides left for the description
2325                      of methods, default 1 in.
2326    :param reverse:  If True, the lowest rank is on the right. Default\: False.
2327    :param cdmethod: None by default. It can be an index of element in avranks
2328                     or or names which specifies the method which should be
2329                     marked with an interval.
2330    """
2331    if not HAS_MATPLOTLIB:
2332        import sys
2333        print >> sys.stderr, "Function requires matplotlib. Please install it."
2334        return
2335
2336    width = float(width)
2337    textspace = float(textspace)
2338
2339    def nth(l, n):
2340        """
2341        Returns only nth elemnt in a list.
2342        """
2343        n = lloc(l, n)
2344        return [ a[n] for a in l ]
2345
2346    def lloc(l, n):
2347        """
2348        List location in list of list structure.
2349        Enable the use of negative locations:
2350        -1 is the last element, -2 second last...
2351        """
2352        if n < 0:
2353            return len(l[0]) + n
2354        else:
2355            return n
2356
2357    def mxrange(lr):
2358        """
2359        Multiple xranges. Can be used to traverse matrices.
2360        This function is very slow due to unknown number of
2361        parameters.
2362
2363        >>> mxrange([3,5])
2364        [(0, 0), (0, 1), (0, 2), (1, 0), (1, 1), (1, 2)]
2365
2366        >>> mxrange([[3,5,1],[9,0,-3]])
2367        [(3, 9), (3, 6), (3, 3), (4, 9), (4, 6), (4, 3)]
2368
2369        """
2370        if not len(lr):
2371            yield ()
2372        else:
2373            #it can work with single numbers
2374            index = lr[0]
2375            if type(1) == type(index):
2376                index = [ index ]
2377            for a in range(*index):
2378                for b in mxrange(lr[1:]):
2379                    yield tuple([a] + list(b))
2380
2381
2382    from matplotlib.figure import Figure
2383    from matplotlib.backends.backend_agg import FigureCanvasAgg
2384
2385    def print_figure(fig, *args, **kwargs):
2386        canvas = FigureCanvasAgg(fig)
2387        canvas.print_figure(*args, **kwargs)
2388
2389    sums = avranks
2390
2391    tempsort = sorted([ (a, i) for i, a in  enumerate(sums) ], reverse=reverse)
2392    ssums = nth(tempsort, 0)
2393    sortidx = nth(tempsort, 1)
2394    nnames = [ names[x] for x in sortidx ]
2395
2396    if lowv is None:
2397        lowv = min(1, int(math.floor(min(ssums))))
2398    if highv is None:
2399        highv = max(len(avranks), int(math.ceil(max(ssums))))
2400
2401    cline = 0.4
2402
2403    k = len(sums)
2404
2405    lines = None
2406
2407    linesblank = 0
2408    scalewidth = width - 2 * textspace
2409
2410    def rankpos(rank):
2411        if not reverse:
2412            a = rank - lowv
2413        else:
2414            a = highv - rank
2415        return textspace + scalewidth / (highv - lowv) * a
2416
2417    distanceh = 0.25
2418
2419    if cd and cdmethod is None:
2420
2421        #get pairs of non significant methods
2422
2423        def get_lines(sums, hsd):
2424
2425            #get all pairs
2426            lsums = len(sums)
2427            allpairs = [ (i, j) for i, j in mxrange([[lsums], [lsums]]) if j > i ]
2428            #remove not significant
2429            notSig = [ (i, j) for i, j in allpairs if abs(sums[i] - sums[j]) <= hsd ]
2430            #keep only longest
2431
2432            def no_longer((i, j), notSig):
2433                for i1, j1 in notSig:
2434                    if (i1 <= i and j1 > j) or (i1 < i and j1 >= j):
2435                        return False
2436                return True
2437
2438            longest = [ (i, j) for i, j in notSig if no_longer((i, j), notSig) ]
2439
2440            return longest
2441
2442        lines = get_lines(ssums, cd)
2443        linesblank = 0.2 + 0.2 + (len(lines) - 1) * 0.1
2444
2446        distanceh = 0.25
2447        cline += distanceh
2448
2449    #calculate height needed height of an image
2450    minnotsignificant = max(2 * 0.2, linesblank)
2451    height = cline + ((k + 1) / 2) * 0.2 + minnotsignificant
2452
2453    fig = Figure(figsize=(width, height))
2454    ax = fig.add_axes([0, 0, 1, 1]) #reverse y axis
2455    ax.set_axis_off()
2456
2457    hf = 1. / height # height factor
2458    wf = 1. / width
2459
2460    def hfl(l):
2461        return [ a * hf for a in l ]
2462
2463    def wfl(l):
2464        return [ a * wf for a in l ]
2465
2466
2467    # Upper left corner is (0,0).
2468
2469    ax.plot([0, 1], [0, 1], c="w")
2470    ax.set_xlim(0, 1)
2471    ax.set_ylim(1, 0)
2472
2473    def line(l, color='k', **kwargs):
2474        """
2475        Input is a list of pairs of points.
2476        """
2477        ax.plot(wfl(nth(l, 0)), hfl(nth(l, 1)), color=color, **kwargs)
2478
2479    def text(x, y, s, *args, **kwargs):
2480        ax.text(wf * x, hf * y, s, *args, **kwargs)
2481
2482    line([(textspace, cline), (width - textspace, cline)], linewidth=0.7)
2483
2484    bigtick = 0.1
2485    smalltick = 0.05
2486
2487
2488    import numpy
2489    tick = None
2490    for a in list(numpy.arange(lowv, highv, 0.5)) + [highv]:
2491        tick = smalltick
2492        if a == int(a): tick = bigtick
2493        line([(rankpos(a), cline - tick / 2), (rankpos(a), cline)], linewidth=0.7)
2494
2495    for a in range(lowv, highv + 1):
2496        text(rankpos(a), cline - tick / 2 - 0.05, str(a), ha="center", va="bottom")
2497
2498    k = len(ssums)
2499
2500    for i in range((k + 1) / 2):
2501        chei = cline + minnotsignificant + i * 0.2
2502        line([(rankpos(ssums[i]), cline), (rankpos(ssums[i]), chei), (textspace - 0.1, chei)], linewidth=0.7)
2503        text(textspace - 0.2, chei, nnames[i], ha="right", va="center")
2504
2505    for i in range((k + 1) / 2, k):
2506        chei = cline + minnotsignificant + (k - i - 1) * 0.2
2507        line([(rankpos(ssums[i]), cline), (rankpos(ssums[i]), chei), (textspace + scalewidth + 0.1, chei)], linewidth=0.7)
2508        text(textspace + scalewidth + 0.2, chei, nnames[i], ha="left", va="center")
2509
2510    if cd and cdmethod is None:
2511
2512        #upper scale
2513        if not reverse:
2514            begin, end = rankpos(lowv), rankpos(lowv + cd)
2515        else:
2516            begin, end = rankpos(highv), rankpos(highv - cd)
2517
2518        line([(begin, distanceh), (end, distanceh)], linewidth=0.7)
2519        line([(begin, distanceh + bigtick / 2), (begin, distanceh - bigtick / 2)], linewidth=0.7)
2520        line([(end, distanceh + bigtick / 2), (end, distanceh - bigtick / 2)], linewidth=0.7)
2521        text((begin + end) / 2, distanceh - 0.05, "CD", ha="center", va="bottom")
2522
2523        #non significance lines
2524        def draw_lines(lines, side=0.05, height=0.1):
2525            start = cline + 0.2
2526            for l, r in lines:
2527                line([(rankpos(ssums[l]) - side, start), (rankpos(ssums[r]) + side, start)], linewidth=2.5)
2528                start += height
2529
2530        draw_lines(lines)
2531
2532    elif cd:
2533        begin = rankpos(avranks[cdmethod] - cd)
2534        end = rankpos(avranks[cdmethod] + cd)
2535        line([(begin, cline), (end, cline)], linewidth=2.5)
2536        line([(begin, cline + bigtick / 2), (begin, cline - bigtick / 2)], linewidth=2.5)
2537        line([(end, cline + bigtick / 2), (end, cline - bigtick / 2)], linewidth=2.5)
2538
2539    print_figure(fig, filename, **kwargs)
2540
2541def mlc_hamming_loss(res):
2542    """
2543    Schapire and Singer (2000) presented Hamming Loss, which id defined as:
2544
2545    :math:HammingLoss(H,D)=\\frac{1}{|D|} \\sum_{i=1}^{|D|} \\frac{Y_i \\vartriangle Z_i}{|L|}
2546    """
2547    losses = [0] * res.number_of_learners
2548    label_num = len(res.labels)
2549    example_num = gettotsize(res)
2550
2551    for e in res.results:
2552        aclass = e.actual_class
2553        for i, labels in enumerate(e.classes):
2554            labels = map(int, labels)
2555            if len(labels) <> len(aclass):
2556                raise ValueError, "The dimensions of the classified output and the actual class array do not match."
2557            for j in range(label_num):
2558                if labels[j] != aclass[j]:
2559                    losses[i] += 1
2560
2561    return [float(x) / (label_num * example_num) for x in losses]
2562
2563def mlc_accuracy(res, forgiveness_rate=1.0):
2564    """
2565    Godbole & Sarawagi, 2004 uses the metrics accuracy, precision, recall as follows:
2566
2567    :math:Accuracy(H,D)=\\frac{1}{|D|} \\sum_{i=1}^{|D|} \\frac{|Y_i \\cap Z_i|}{|Y_i \\cup Z_i|}
2568
2569    Boutell et al. (2004) give a more generalized version using a parameter :math:\\alpha \\ge 0,
2570    called forgiveness rate:
2571
2572    :math:Accuracy(H,D)=\\frac{1}{|D|} \\sum_{i=1}^{|D|} (\\frac{|Y_i \\cap Z_i|}{|Y_i \\cup Z_i|})^{\\alpha}
2573    """
2574    accuracies = [0.0] * res.number_of_learners
2575    example_num = gettotsize(res)
2576
2577    for e in res.results:
2578        aclass = e.actual_class
2579        for i, labels in enumerate(e.classes):
2580            labels = map(int, labels)
2581            if len(labels) <> len(aclass):
2582                raise ValueError, "The dimensions of the classified output and the actual class array do not match."
2583
2584            intersection = 0.0
2585            union = 0.0
2586            for real, pred in zip(labels, aclass):
2587                if real and pred:
2588                    intersection += 1
2589                if real or pred:
2590                    union += 1
2591
2592            if union:
2593                accuracies[i] += intersection / union
2594
2595    return [math.pow(x / example_num, forgiveness_rate) for x in accuracies]
2596
2597def mlc_precision(res):
2598    """
2599    :math:Precision(H,D)=\\frac{1}{|D|} \\sum_{i=1}^{|D|} \\frac{|Y_i \\cap Z_i|}{|Z_i|}
2600    """
2601    precisions = [0.0] * res.number_of_learners
2602    example_num = gettotsize(res)
2603
2604    for e in res.results:
2605        aclass = e.actual_class
2606        for i, labels in enumerate(e.classes):
2607            labels = map(int, labels)
2608            if len(labels) <> len(aclass):
2609                raise ValueError, "The dimensions of the classified output and the actual class array do not match."
2610
2611            intersection = 0.0
2612            predicted = 0.0
2613            for real, pred in zip(labels, aclass):
2614                if real and pred:
2615                    intersection += 1
2616                if real:
2617                    predicted += 1
2618            if predicted:
2619                precisions[i] += intersection / predicted
2620
2621    return [x / example_num for x in precisions]
2622
2623def mlc_recall(res):
2624    """
2625    :math:Recall(H,D)=\\frac{1}{|D|} \\sum_{i=1}^{|D|} \\frac{|Y_i \\cap Z_i|}{|Y_i|}
2626    """
2627    recalls = [0.0] * res.number_of_learners
2628    example_num = gettotsize(res)
2629
2630    for e in res.results:
2631        aclass = e.actual_class
2632        for i, labels in enumerate(e.classes):
2633            labels = map(int, labels)
2634            if len(labels) <> len(aclass):
2635                raise ValueError, "The dimensions of the classified output and the actual class array do not match."
2636
2637            intersection = 0.0
2638            actual = 0.0
2639            for real, pred in zip(labels, aclass):
2640                if real and pred:
2641                    intersection += 1
2642                if pred:
2643                    actual += 1
2644            if actual:
2645                recalls[i] += intersection / actual
2646    return [x / example_num for x in recalls]
2647
2648#def mlc_ranking_loss(res):
2649#    pass
2650#
2651#def mlc_average_precision(res):
2652#    pass
2653#
2654#def mlc_hierarchical_loss(res):
2655#    pass
2656
2657def logloss(res):
2658    """
2659    Calculates LogLoss, n is the number of all test results and :math:p_{i} is the probability
2660     withw hich the classifier predicted the actual class.
2661     :math:LogLoss = \\frac{1}{n}\\sum_{i = 1}^{n} -max(log(p_{i}), log \\frac{1}{n}) \\newline
2662    """
2663    results = []
2664    n_results = len(res.results)
2665    min_log = math.log(1.0/n_results)
2666    for l in xrange(res.number_of_learners):
2667        temp = 0.0
2668        for r in res.results:
2669            if not r.probabilities[l]:
2670                raise ValueError, "Probabilities are needed to compute logloss"
2671            temp-=max(math.log(max(r.probabilities[l][int(r.actual_class)],1e-20)),min_log)
2672
2673        results.append(temp/n_results)
2674    return results
2675
2676
2677def mlc_F1_micro(res):
2678    """
2679    F1_{micro} = 2 * \frac{\overline{precision}  * \overline{recall}}{\overline{precision} + \overline{recall}}
2680    """
2681
2682    precision = mlc_precision(res)
2683    recall = mlc_recall(res)
2684    return [0.0 if p == 0 and r == 0 else 2 * p * r / (p + r) for p,r in zip(precision, recall)]
2685
2686
2687def mlc_F1_macro(res):
2688    """
2689    F1_{macro} = \frac{1}{d}\sum_{j=0}^{d} 2 * \frac{precision_j * recall_j}{precision_j + recall_j}
2690    """
2691
2692    results = []
2693    n_results = gettotsize(res)
2694    n_classes =  len(res.results[0].actual_class)
2695
2696    for l in xrange(res.number_of_learners):
2697        true_positive = [0.0000001] * n_classes
2698        sum_fptp = [0.0000001] * n_classes
2699        sum_fntp = [0.0000001] * n_classes
2700        for r in res.results:
2701            aclass = r.actual_class
2702
2703            for i, cls_val in enumerate(r.classes[l]):
2704                if aclass[i]==1 and int(cls_val)==1:
2705                    true_positive[i] += 1
2706                if int(cls_val)==1:
2707                    sum_fptp[i] += 1
2708                if aclass[i]==1:
2709                    sum_fntp[i] += 1
2710
2711        results.append(sum([ 2*(tp/fptp * tp/fntp)/(tp/fptp + tp/fntp) for tp, fptp, fntp in \
2712            zip(true_positive, sum_fptp, sum_fntp)] ) / n_classes)
2713    return results
2714
2715
2716
2717################################################################################
2718if __name__ == "__main__":
2719    avranks = [3.143, 2.000, 2.893, 1.964]
2720    names = ["prva", "druga", "tretja", "cetrta" ]
2721    cd = compute_CD(avranks, 14)
2722    #cd = compute_CD(avranks, 10, type="bonferroni-dunn")
2723    print cd
2724
2725    print compute_friedman(avranks, 14)
2726
2727    #graph_ranks("test.eps", avranks, names, cd=cd, cdmethod=0, width=6, textspace=1.5)
Note: See TracBrowser for help on using the repository browser.