# source:orange/Orange/evaluation/scoring.py@10367:058e9b1cc258

Revision 10367:058e9b1cc258, 98.2 KB checked in by Lan Zagar <lan.zagar@…>, 2 years ago (diff)

Some more code style and doc improvements.

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