source: orange/Orange/evaluation/scoring.py @ 11678:be71dd8a74e4

Revision 11678:be71dd8a74e4, 99.0 KB checked in by Ales Erjavec <ales.erjavec@…>, 8 months ago (diff)

Added test for 'split_by_classifiers' helper.

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