source: orange/Orange/evaluation/scoring.py @ 10425:40e3b077c358

Revision 10425:40e3b077c358, 98.4 KB checked in by anzeh <anze.staric@…>, 2 years ago (diff)

Scoring functions ((based on confusion matrices) are now classes, but documented as functions.

Line 
1import math
2import functools
3from operator import add
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,
56                loaded=res.loaded, test_type=res.test_type, labels=res.labels)
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) \
588    <http://www.springerlink.com/content/g5p7473160476612/>`_.
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):
712                    tfpns[i].addTFPosNeg(lr.probabilities[i][class_index] >
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):
718                    tfpns[i].addTFPosNeg(lr.probabilities[i][class_index] > 
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):
725                    tfpns[i].addTFPosNeg(lr.classes[i] == class_index,
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):
731                    tfpns[i].addTFPosNeg(lr.classes[i] == class_index,
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
772class CMScore(list):
773    """
774    :param test_results: :obj:`~Orange.evaluation.testing.ExperimentResults`
775                         or list of :obj:`ConfusionMatrix`.
776    :rtype: list of scores, one for each learner."""
777    @deprecated_keywords({"confm": "test_results"})
778    def __init__(self, test_results=None):
779        super(CMScore, self).__init__()
780
781        if test_results is not None:
782            self[:] = self.__call__(test_results)
783
784    def __call__(self, test_results):
785        if isinstance(test_results, testing.ExperimentResults):
786            test_results = confusion_matrices(test_results, class_index=1)
787        if isinstance(test_results, ConfusionMatrix):
788            test_results = [test_results]
789
790        return map(self.compute, test_results)
791
792
793
794class Sensitivity(CMScore):
795    __doc__ = """Compute `sensitivity
796    <http://en.wikipedia.org/wiki/Sensitivity_and_specificity>`_ (proportion
797    of actual positives which are correctly identified as such).
798    """ + CMScore.__doc__
799    @classmethod
800    def compute(self, confusion_matrix):
801        tot = confusion_matrix.TP+confusion_matrix.FN
802        if tot < 1e-6:
803            import warnings
804            warnings.warn("Can't compute sensitivity: one or both classes have no instances")
805            return None
806
807        return confusion_matrix.TP / tot
808
809
810class Recall(Sensitivity):
811    __doc__ = """ Compute `recall
812    <http://en.wikipedia.org/wiki/Precision_and_recall>`_
813    (fraction of relevant instances that are retrieved).
814    """ + CMScore.__doc__
815    pass # Recall == Sensitivity
816
817
818class Specificity(CMScore):
819    __doc__ = """Compute `specificity
820    <http://en.wikipedia.org/wiki/Sensitivity_and_specificity>`_
821    (proportion of negatives which are correctly identified).
822    """ + CMScore.__doc__
823    @classmethod
824    def compute(self, confusion_matrix):
825        tot = confusion_matrix.FP+confusion_matrix.TN
826        if tot < 1e-6:
827            import warnings
828            warnings.warn("Can't compute specificity: one or both classes have no instances")
829            return None
830        return confusion_matrix.TN / tot
831
832
833class PPV(CMScore):
834    __doc__ = """Compute `positive predictive value
835    <http://en.wikipedia.org/wiki/Positive_predictive_value>`_ (proportion of
836    subjects with positive test results who are correctly diagnosed).
837    """ + CMScore.__doc__
838    @classmethod
839    def compute(self, confusion_matrix):
840        tot = confusion_matrix.TP + confusion_matrix.FP
841        if tot < 1e-6:
842            import warnings
843            warnings.warn("Can't compute PPV: one or both classes have no instances")
844            return None
845        return confusion_matrix.TP/tot
846
847
848class Precision(PPV):
849    __doc__ = """Compute `precision <http://en.wikipedia.org/wiki/Precision_and_recall>`_
850    (retrieved instances that are relevant).
851    """ + CMScore.__doc__
852    pass # Precision == PPV
853
854
855class NPV(CMScore):
856    __doc__ = """Compute `negative predictive value
857    <http://en.wikipedia.org/wiki/Negative_predictive_value>`_ (proportion of
858    subjects with a negative test result who are correctly diagnosed).
859     """ + CMScore.__doc__
860    @classmethod
861    def compute(self, confusion_matrix):
862        tot = confusion_matrix.FN + confusion_matrix.TN
863        if tot < 1e-6:
864            import warnings
865            warnings.warn("Can't compute NPV: one or both classes have no instances")
866            return None
867        return confusion_matrix.TN / tot
868
869
870class F1(CMScore):
871    __doc__ = """Return `F1 score
872    <http://en.wikipedia.org/wiki/F1_score>`_
873    (harmonic mean of precision and recall).
874    """ + CMScore.__doc__
875    @classmethod
876    def compute(self, confusion_matrix):
877        p = Precision.compute(confusion_matrix)
878        r = Recall.compute(confusion_matrix)
879        if p is not None and r is not None and (p+r) != 0:
880            return 2. * p * r / (p + r)
881        else:
882            import warnings
883            warnings.warn("Can't compute F1: P + R is zero or not defined")
884            return None
885
886
887class Falpha(CMScore):
888    __doc__ = """Compute the alpha-mean of precision and recall over the given confusion
889    matrix.
890    """ + CMScore.__doc__
891
892    def __init__(self, test_results, alpha=1.):
893        self.alpha = alpha
894        super(Falpha, self).__init__(test_results)
895
896    def compute(self, confusion_matrix):
897        p = Precision.compute(confusion_matrix)
898        r = Recall.compute(confusion_matrix)
899        return (1. + self.alpha) * p * r / (self.alpha * p + r)
900
901
902class MCC(CMScore):
903    __doc__ = """Compute `Matthew correlation coefficient
904    <http://en.wikipedia.org/wiki/Matthews_correlation_coefficient>`_
905    (correlation coefficient between the observed and predicted binary
906    classifications).
907    """ + CMScore.__doc__
908    @classmethod
909    def compute(self, cm):
910        # code by Boris Gorelik
911        TP, TN, FP, FN = cm.TP, cm.TN, cm.FP, cm.FN
912         
913        try:
914            return (TP*TN - FP*FN) /\
915                 math.sqrt((TP+FP) * (TP+FN) * (TN+ FP) * (TN+FN))
916        except ZeroDivisionError:
917            # Zero division occurs when there is either no true positives
918            # or no true negatives i.e. the problem contains only one
919            # type of classes.
920            import warnings
921            warnings.warn("Can't compute MCC: TP or TN is zero or not defined")
922
923
924@deprecated_keywords({"bIsListOfMatrices": "b_is_list_of_matrices"})
925def scotts_pi(confusion_matrix, b_is_list_of_matrices=True):
926   """Compute Scott's Pi for measuring inter-rater agreement for nominal data
927
928   http://en.wikipedia.org/wiki/Scott%27s_Pi
929   Scott's Pi is a statistic for measuring inter-rater reliability for nominal
930   raters.
931
932   @param confusion_matrix: confusion matrix, or list of confusion matrices.
933                            To obtain non-binary confusion matrix, call
934                            Orange.evaluation.scoring.confusion_matrices
935                            and set the class_index parameter to -2.
936   @param b_is_list_of_matrices: specifies whether confm is list of matrices.
937                            This function needs to operate on non-binary
938                            confusion matrices, which are represented by python
939                            lists, therefore one needs a way to distinguish
940                            between a single matrix and list of matrices
941   """
942
943   if b_is_list_of_matrices:
944       try:
945           return [scotts_pi(cm, b_is_list_of_matrices=False)
946                   for cm in confusion_matrix]
947       except TypeError:
948           # Nevermind the parameter, maybe this is a "conventional" binary
949           # confusion matrix and bIsListOfMatrices was specified by mistake
950           return scottsPiSingle(confusion_matrix, bIsListOfMatrices=False)
951   else:
952       if isinstance(confusion_matrix, ConfusionMatrix):
953           confusion_matrix = numpy.array([[confusion_matrix.TP,
954               confusion_matrix.FN], [confusion_matrix.FP,
955               confusion_matrix.TN]], dtype=float)
956       else:
957           confusion_matrix = numpy.array(confusion_matrix, dtype=float)
958
959       marginalSumOfRows = numpy.sum(confusion_matrix, axis=0)
960       marginalSumOfColumns = numpy.sum(confusion_matrix, axis=1)
961       jointProportion = (marginalSumOfColumns + marginalSumOfRows) / \
962                           (2. * numpy.sum(confusion_matrix))
963       # In the eq. above, 2. is what the Wikipedia page calls
964       # the number of annotators. Here we have two annotators:
965       # the observed (true) labels (annotations) and the predicted by
966       # the learners.
967
968       prExpected = numpy.sum(jointProportion ** 2)
969       prActual = numpy.sum(numpy.diag(confusion_matrix)) / \
970                  numpy.sum(confusion_matrix)
971
972       ret = (prActual - prExpected) / (1.0 - prExpected)
973       return ret
974
975# Backward compatibility
976sens = Sensitivity
977spec = Specificity
978precision = Precision
979recall = Recall
980
981
982
983@deprecated_keywords({"classIndex": "class_index",
984                      "unweighted": "ignore_weights"})
985def AUCWilcoxon(res, class_index=-1, ignore_weights=False, **argkw):
986    """Compute the area under ROC (AUC) and its standard error using
987    Wilcoxon's approach proposed by Hanley and McNeal (1982). If
988    :obj:`class_index` is not specified, the first class is used as
989    "the positive" and others are negative. The result is a list of
990    tuples (aROC, standard error).
991
992    If test results consist of multiple folds, you need to split them using
993    :obj:`split_by_iterations` and perform this test on each fold separately.
994    """
995    useweights = res.weights and not ignore_weights
996    problists, tots = corn.computeROCCumulative(res, class_index, useweights)
997
998    results = []
999
1000    totPos, totNeg = tots[1], tots[0]
1001    N = totPos + totNeg
1002    for plist in problists:
1003        highPos, lowNeg = totPos, 0.
1004        W, Q1, Q2 = 0., 0., 0.
1005        for prob in plist:
1006            thisPos, thisNeg = prob[1][1], prob[1][0]
1007            highPos -= thisPos
1008            W += thisNeg * (highPos + thisPos / 2.)
1009            Q2 += thisPos * (lowNeg**2  + lowNeg*thisNeg  + thisNeg**2 / 3.)
1010            Q1 += thisNeg * (highPos**2 + highPos*thisPos + thisPos**2 / 3.)
1011
1012            lowNeg += thisNeg
1013
1014        W  /= (totPos*totNeg)
1015        Q1 /= (totNeg*totPos**2)
1016        Q2 /= (totPos*totNeg**2)
1017
1018        SE = math.sqrt((W * (1 - W) + (totPos - 1) * (Q1 - W**2) +
1019                       (totNeg - 1) * (Q2 - W**2)) / (totPos * totNeg))
1020        results.append((W, SE))
1021    return results
1022
1023AROC = AUCWilcoxon # for backward compatibility, AROC is obsolete
1024
1025
1026@deprecated_keywords({"classIndex": "class_index",
1027                      "unweighted": "ignore_weights"})
1028def compare_2_AUCs(res, lrn1, lrn2, class_index=-1,
1029                   ignore_weights=False, **argkw):
1030    return corn.compare2ROCs(res, lrn1, lrn2, class_index,
1031                             res.weights and not ignore_weights)
1032
1033# for backward compatibility, compare_2_AROCs is obsolete
1034compare_2_AROCs = compare_2_AUCs
1035
1036
1037@deprecated_keywords({"classIndex": "class_index"})
1038def compute_ROC(res, class_index=-1):
1039    """Compute a ROC curve as a list of (x, y) tuples, where x is
1040    1-specificity and y is sensitivity.
1041    """
1042    problists, tots = corn.computeROCCumulative(res, class_index)
1043
1044    results = []
1045    totPos, totNeg = tots[1], tots[0]
1046
1047    for plist in problists:
1048        curve=[(1., 1.)]
1049        TP, TN = totPos, 0.
1050        FN, FP = 0., totNeg
1051        for prob in plist:
1052            thisPos, thisNeg = prob[1][1], prob[1][0]
1053            # thisPos go from TP to FN
1054            TP -= thisPos
1055            FN += thisPos
1056            # thisNeg go from FP to TN
1057            TN += thisNeg
1058            FP -= thisNeg
1059
1060            sens = TP / (TP + FN)
1061            spec = TN / (FP + TN)
1062            curve.append((1 - spec, sens))
1063        results.append(curve)
1064
1065    return results   
1066
1067## TC's implementation of algorithms, taken from:
1068## T Fawcett: ROC Graphs: Notes and Practical Considerations for
1069## Data Mining Researchers, submitted to KDD Journal.
1070def ROC_slope((P1x, P1y, P1fscore), (P2x, P2y, P2fscore)):
1071    if P1x == P2x:
1072        return 1e300
1073    return (P1y - P2y) / (P1x - P2x)
1074
1075
1076@deprecated_keywords({"keepConcavities": "keep_concavities"})
1077def ROC_add_point(P, R, keep_concavities=1):
1078    if keep_concavities:
1079        R.append(P)
1080    else:
1081        while True:
1082            if len(R) < 2:
1083                R.append(P)
1084                return R
1085            else:
1086                T = R.pop()
1087                T2 = R[-1]
1088                if ROC_slope(T2, T) > ROC_slope(T, P):
1089                    R.append(T)
1090                    R.append(P)
1091                    return R
1092    return R
1093
1094
1095@deprecated_keywords({"classIndex": "class_index",
1096                      "keepConcavities": "keep_concavities"})
1097def TC_compute_ROC(res, class_index=-1, keep_concavities=1):
1098    problists, tots = corn.computeROCCumulative(res, class_index)
1099
1100    results = []
1101    P, N = tots[1], tots[0]
1102
1103    for plist in problists:
1104        ## corn gives an increasing by scores list, we need a decreasing
1105        plist.reverse()
1106        TP = 0.
1107        FP = 0.
1108        curve=[]
1109        fPrev = 10e300 # "infinity" score at 0., 0.
1110        for prob in plist:
1111            f = prob[0]
1112            if f != fPrev:
1113                if P:
1114                    tpr = TP / P
1115                else:
1116                    tpr = 0.
1117                if N:
1118                    fpr = FP / N
1119                else:
1120                    fpr = 0.
1121                curve = ROC_add_point((fpr, tpr, fPrev), curve,
1122                                      keep_concavities)
1123                fPrev = f
1124            thisPos, thisNeg = prob[1][1], prob[1][0]
1125            TP += thisPos
1126            FP += thisNeg
1127        if P:
1128            tpr = TP / P
1129        else:
1130            tpr = 0.
1131        if N:
1132            fpr = FP / N
1133        else:
1134            fpr = 0.
1135        curve = ROC_add_point((fpr, tpr, f), curve, keep_concavities) ## ugly
1136        results.append(curve)
1137
1138    return results
1139
1140## returns a list of points at the intersection of the tangential
1141## iso-performance line and the given ROC curve
1142## for given values of FPcost, FNcost and pval
1143def TC_best_thresholds_on_ROC_curve(FPcost, FNcost, pval, curve):
1144    m = (FPcost * (1. - pval)) / (FNcost * pval)
1145
1146    ## put the iso-performance line in point (0., 1.)
1147    x0, y0 = (0., 1.)
1148    x1, y1 = (1., 1. + m)
1149    d01 = math.sqrt((x1 - x0) * (x1 - x0) + (y1 - y0) * (y1 - y0))
1150
1151    ## calculate and find the closest point to the line
1152    firstp = 1
1153    mind = 0.
1154    a = x0 * y1 - x1 * y0
1155    closestPoints = []
1156    for (x, y, fscore) in curve:
1157        d = ((y0 - y1) * x + (x1 - x0) * y + a) / d01
1158        d = abs(d)
1159        if firstp or d < mind:
1160            mind, firstp = d, 0
1161            closestPoints = [(x, y, fscore)]
1162        else:
1163            if abs(d - mind) <= 0.0001: ## close enough
1164                closestPoints.append( (x, y, fscore) )
1165    return closestPoints         
1166
1167def frange(start, end=None, inc=None):
1168    """A range function, that does accept float increments..."""
1169
1170    if end is None:
1171        end = start + 0.
1172        start = 0.
1173
1174    if inc is None or inc == 0:
1175        inc = 1.
1176
1177    L = [start]
1178    while 1:
1179        next = start + len(L) * inc
1180        if inc > 0 and next >= end:
1181            L.append(end)
1182            break
1183        elif inc < 0 and next <= end:
1184            L.append(end)
1185            break
1186        L.append(next)
1187       
1188    return L
1189
1190## input ROCcurves are of form [ROCcurves1, ROCcurves2, ... ROCcurvesN],
1191## where ROCcurvesX is a set of ROC curves,
1192## where a (one) ROC curve is a set of (FP, TP) points
1193##
1194## for each (sub)set of input ROC curves
1195## returns the average ROC curve and an array of (vertical) standard deviations
1196@deprecated_keywords({"ROCcurves": "roc_curves"})
1197def TC_vertical_average_ROC(roc_curves, samples = 10):
1198    def INTERPOLATE((P1x, P1y, P1fscore), (P2x, P2y, P2fscore), X):
1199        if (P1x == P2x) or P1x < X > P2x or P1x > X < P2x:
1200            raise ValueError, "assumptions for interpolation are not met: P1 = %f,%f P2 = %f,%f X = %f" % (P1x, P1y, P2x, P2y, X)
1201        dx = float(P2x) - float(P1x)
1202        dy = float(P2y) - float(P1y)
1203        m = dy / dx
1204        return P1y + m * (X - P1x)
1205
1206    def TP_FOR_FP(FPsample, ROC, npts):
1207        i = 0
1208        while i < npts - 1:
1209            (fp, _, _) = ROC[i + 1]
1210            if fp <= FPsample:
1211                i += 1
1212            else:
1213                break
1214        (fp, tp, _) = ROC[i]
1215        if fp == FPsample:
1216            return tp
1217        elif fp < FPsample and i + 1 < len(ROC):
1218            return INTERPOLATE(ROC[i], ROC[i + 1], FPsample)
1219        elif fp < FPsample and i + 1 == len(ROC): # return the last
1220            return ROC[i][1]
1221        raise ValueError, "cannot compute: TP_FOR_FP in TC_vertical_average_ROC"
1222        #return 0.
1223
1224    average = []
1225    stdev = []
1226    for ROCS in roc_curves:
1227        npts = []
1228        for c in ROCS:
1229            npts.append(len(c))
1230        nrocs = len(ROCS)
1231
1232        TPavg = []
1233        TPstd = []
1234        for FPsample in frange(0., 1., 1. / samples):
1235            TPsum = []
1236            for i in range(nrocs):
1237                ##TPsum = TPsum + TP_FOR_FP(FPsample, ROCS[i], npts[i])
1238                TPsum.append(TP_FOR_FP(FPsample, ROCS[i], npts[i]))
1239            TPavg.append((FPsample, statc.mean(TPsum)))
1240            if len(TPsum) > 1:
1241                stdv = statc.std(TPsum)
1242            else:
1243                stdv = 0.
1244            TPstd.append(stdv)
1245
1246        average.append(TPavg)
1247        stdev.append(TPstd)
1248
1249    return average, stdev
1250
1251## input ROCcurves are of form [ROCcurves1, ROCcurves2, ... ROCcurvesN],
1252## where ROCcurvesX is a set of ROC curves,
1253## where a (one) ROC curve is a set of (FP, TP) points
1254##
1255## for each (sub)set of input ROC curves
1256## returns the average ROC curve, an array of vertical standard deviations and an array of horizontal standard deviations
1257@deprecated_keywords({"ROCcurves": "roc_curves"})
1258def TC_threshold_average_ROC(roc_curves, samples=10):
1259    def POINT_AT_THRESH(ROC, npts, thresh):
1260        i = 0
1261        while i < npts - 1:
1262            (px, py, pfscore) = ROC[i]
1263            if pfscore > thresh:
1264                i += 1
1265            else:
1266                break
1267        return ROC[i]
1268
1269    average = []
1270    stdevV = []
1271    stdevH = []
1272    for ROCS in roc_curves:
1273        npts = []
1274        for c in ROCS:
1275            npts.append(len(c))
1276        nrocs = len(ROCS)
1277
1278        T = []
1279        for c in ROCS:
1280            for (px, py, pfscore) in c:
1281##                try:
1282##                    T.index(pfscore)
1283##                except:
1284                T.append(pfscore)
1285        T.sort()
1286        T.reverse() ## ugly
1287
1288        TPavg = []
1289        TPstdV = []
1290        TPstdH = []
1291        for tidx in frange(0, (len(T) - 1.), float(len(T)) / samples):
1292            FPsum = []
1293            TPsum = []
1294            for i in range(nrocs):
1295                (fp, tp, _) = POINT_AT_THRESH(ROCS[i], npts[i], T[int(tidx)])
1296                FPsum.append(fp)
1297                TPsum.append(tp)
1298            TPavg.append((statc.mean(FPsum), statc.mean(TPsum)))
1299            ## vertical standard deviation
1300            if len(TPsum) > 1:
1301                stdv = statc.std(TPsum)
1302            else:
1303                stdv = 0.
1304            TPstdV.append( stdv )
1305            ## horizontal standard deviation
1306            if len(FPsum) > 1:
1307                stdh = statc.std(FPsum)
1308            else:
1309                stdh = 0.
1310            TPstdH.append( stdh )
1311
1312        average.append(TPavg)
1313        stdevV.append(TPstdV)
1314        stdevH.append(TPstdH)
1315
1316    return average, stdevV, stdevH
1317
1318## Calibration Curve
1319## returns an array of (curve, yesClassPredictions, noClassPredictions)
1320## elements, where:
1321##  - curve is an array of points (x, y) on the calibration curve
1322##  - yesClassRugPoints is an array of (x, 1) points
1323##  - noClassRugPoints is an array of (x, 0) points
1324@deprecated_keywords({"classIndex": "class_index"})
1325def compute_calibration_curve(res, class_index=-1):
1326    ## merge multiple iterations into one
1327    mres = Orange.evaluation.testing.ExperimentResults(1, res.classifier_names,
1328        res.class_values, res.weights, classifiers=res.classifiers,
1329        loaded=res.loaded, test_type=res.test_type, labels=res.labels)
1330    for te in res.results:
1331        mres.results.append(te)
1332
1333    problists, tots = corn.computeROCCumulative(mres, class_index)
1334
1335    results = []
1336
1337    bins = 10 ## divide interval between 0. and 1. into N bins
1338
1339    for plist in problists:
1340        yesClassRugPoints = [] 
1341        noClassRugPoints = []
1342
1343        yesBinsVals = [0] * bins
1344        noBinsVals = [0] * bins
1345        for (f, (thisNeg, thisPos)) in plist:
1346            yesClassRugPoints.append((f, thisPos)) # 1.
1347            noClassRugPoints.append((f, thisNeg)) # 1.
1348
1349            index = int(f * bins )
1350            index = min(index, bins - 1) ## just in case for value 1.
1351            yesBinsVals[index] += thisPos
1352            noBinsVals[index] += thisNeg
1353
1354        curve = []
1355        for cn in range(bins):
1356            f = float(cn * 1. / bins) + (1. / 2. / bins)
1357            yesVal = yesBinsVals[cn]
1358            noVal = noBinsVals[cn]
1359            allVal = yesVal + noVal
1360            if allVal == 0.: continue
1361            y = float(yesVal) / float(allVal)
1362            curve.append((f,  y))
1363
1364        ## smooth the curve
1365        maxnPoints = 100
1366        if len(curve) >= 3:
1367#            loessCurve = statc.loess(curve, -3, 0.6)
1368            loessCurve = statc.loess(curve, maxnPoints, 0.5, 3)
1369        else:
1370            loessCurve = curve
1371        clen = len(loessCurve)
1372        if clen > maxnPoints:
1373            df = clen / maxnPoints
1374            if df < 1: df = 1
1375            curve = [loessCurve[i]  for i in range(0, clen, df)]
1376        else:
1377            curve = loessCurve
1378        ## remove the third value (variance of epsilon?) that suddenly
1379        ## appeared in the output of the statc.loess function
1380        curve = [c[:2] for c in curve]
1381        results.append((curve, yesClassRugPoints, noClassRugPoints))
1382
1383    return results
1384
1385
1386## Lift Curve
1387## returns an array of curve elements, where:
1388##  - curve is an array of points ((TP + FP) / (P + N), TP / P, (th, FP / N))
1389##    on the Lift Curve
1390@deprecated_keywords({"classIndex": "class_index"})
1391def compute_lift_curve(res, class_index=-1):
1392    ## merge multiple iterations into one
1393    mres = Orange.evaluation.testing.ExperimentResults(1, res.classifier_names,
1394        res.class_values, res.weights, classifiers=res.classifiers,
1395        loaded=res.loaded, test_type=res.test_type, labels=res.labels)
1396    for te in res.results:
1397        mres.results.append(te)
1398
1399    problists, tots = corn.computeROCCumulative(mres, class_index)
1400
1401    results = []
1402    P, N = tots[1], tots[0]
1403    for plist in problists:
1404        ## corn gives an increasing by scores list, we need a decreasing
1405        plist.reverse()
1406        TP = 0.
1407        FP = 0.
1408        curve = [(0., 0., (10e300, 0.))]
1409        for (f, (thisNeg, thisPos)) in plist:
1410            TP += thisPos
1411            FP += thisNeg
1412            curve.append(((TP + FP) / (P + N), TP, (f, FP / (N or 1))))
1413        results.append(curve)
1414
1415    return P, N, results
1416
1417
1418class CDT:
1419  """Store the number of concordant (C), discordant (D) and tied (T) pairs."""
1420  def __init__(self, C=0., D=0., T=0.):
1421    self.C, self.D, self.T = C, D, T
1422   
1423def is_CDT_empty(cdt):
1424    return cdt.C + cdt.D + cdt.T < 1e-20
1425
1426
1427@deprecated_keywords({"classIndex": "class_index",
1428                      "unweighted": "ignore_weights"})
1429def compute_CDT(res, class_index=-1, ignore_weights=False, **argkw):
1430    """Obsolete, don't use."""
1431    if class_index < 0:
1432        if res.baseClass >= 0:
1433            class_index = res.baseClass
1434        else:
1435            class_index = 1
1436           
1437    useweights = res.weights and not ignore_weights
1438    weightByClasses = argkw.get("weightByClasses", True)
1439
1440    if res.number_of_iterations>1:
1441        CDTs = [CDT() for _ in range(res.number_of_learners)]
1442        iterationExperiments = split_by_iterations(res)
1443        for exp in iterationExperiments:
1444            expCDTs = corn.computeCDT(exp, class_index, useweights)
1445            for i in range(len(CDTs)):
1446                CDTs[i].C += expCDTs[i].C
1447                CDTs[i].D += expCDTs[i].D
1448                CDTs[i].T += expCDTs[i].T
1449        for i in range(res.number_of_learners):
1450            if is_CDT_empty(CDTs[0]):
1451                return corn.computeCDT(res, class_index, useweights)
1452       
1453        return CDTs
1454    else:
1455        return corn.computeCDT(res, class_index, useweights)
1456
1457# Backward compatibility
1458def replace_use_weights(fun):
1459    if environ.orange_no_deprecated_members:
1460        return fun
1461
1462    @functools.wraps(fun)
1463    def wrapped(*args, **kwargs):
1464        use_weights = kwargs.pop("useWeights", None)
1465        if use_weights is not None:
1466            deprecation_warning("useWeights", "ignore_weights")
1467            kwargs["ignore_weights"] = not use_weights
1468        return fun(*args, **kwargs)
1469    return wrapped
1470
1471class AUC(list):
1472    """
1473    Compute the area under ROC curve given a set of experimental results.
1474    For multivalued class problems, return the result of
1475    :obj:`by_weighted_pairs`.
1476    If testing consisted of multiple folds, each fold is scored and the
1477    average score is returned. If a fold contains only instances with the
1478    same class value, folds will be merged.
1479
1480    :param test_results: test results to score
1481    :param ignore_weights: ignore instance weights when calculating score
1482    :param method: DEPRECATED, call the appropriate method directly.
1483    """
1484
1485    @replace_use_weights
1486    def __init__(self, test_results=None, method=0, ignore_weights=False):
1487
1488        super(AUC, self).__init__()
1489
1490        self.ignore_weights=ignore_weights
1491        self.method=method
1492
1493        if test_results is not None:
1494            self.__call__(test_results)
1495
1496    def __call__(self, test_results):
1497        if len(test_results.class_values) < 2:
1498            raise ValueError("Cannot compute AUC on a single-class problem")
1499        elif len(test_results.class_values) == 2:
1500            self._compute_for_binary_class(test_results)
1501        else:
1502            self._compute_for_multi_value_class(test_results, self.method)
1503
1504    @classmethod
1505    def by_weighted_pairs(cls, res, ignore_weights=False):
1506        """
1507        Compute AUC for each pair of classes (ignoring instances of all other
1508        classes) and average the results, weighting them by the number of
1509        pairs of instances from these two classes (e.g. by the product of
1510        probabilities of the two classes). AUC computed in this way still
1511        behaves as the concordance index, e.g., gives the probability that two
1512        randomly chosen instances from different classes will be correctly
1513        recognized (if the classifier knows from which two classes the
1514        instances came).
1515        """
1516        auc = AUC(ignore_weights=ignore_weights)
1517        auc._compute_for_multi_value_class(res, method=cls.ByWeightedPairs)
1518        return auc
1519
1520    @classmethod
1521    def by_pairs(cls, res, ignore_weights=False):
1522        """
1523        Similar to by_weighted_pairs, except that the average over class pairs
1524        is not weighted. This AUC is, like the binary version, independent of
1525        class distributions, but it is not related to the concordance index
1526        any more.
1527        """
1528        auc = AUC(ignore_weights=ignore_weights)
1529        auc._compute_for_multi_value_class(res, method=cls.ByPairs)
1530        return auc
1531
1532    @classmethod
1533    def weighted_one_against_all(cls, res, ignore_weights=False):
1534        """
1535        For each class, it computes AUC for this class against all others (that
1536        is, treating other classes as one class). The AUCs are then averaged by
1537        the class probabilities. This is related to the concordance index in
1538        which we test the classifier's (average) capability of distinguishing
1539        the instances from a specified class from those that come from other
1540        classes.
1541        Unlike the binary AUC, the measure is not independent of class
1542        distributions.
1543        """
1544        auc = AUC(ignore_weights=ignore_weights)
1545        auc._compute_for_multi_value_class(res,
1546            method=cls.WeightedOneAgainstAll)
1547        return auc
1548
1549    @classmethod
1550    def one_against_all(cls, res, ignore_weights=False):
1551        """
1552        Similar to weighted_one_against_all, except that the average
1553        is not weighted.
1554        """
1555        auc = AUC(ignore_weights=ignore_weights)
1556        auc._compute_for_multi_value_class(res, method=cls.OneAgainstAll)
1557        return auc
1558
1559    @classmethod
1560    def single_class(cls, res, class_index=-1, ignore_weights=False):
1561        """
1562        Compute AUC where the class with the given class_index is singled
1563        out and all other classes are treated as a single class.
1564        """
1565        if class_index < 0:
1566            if res.base_class >= 0:
1567                class_index = res.base_class
1568            else:
1569                class_index = 1
1570
1571        auc = AUC(ignore_weights=ignore_weights)
1572        auc._compute_for_single_class(res, class_index)
1573        return auc
1574
1575    @classmethod
1576    def pair(cls, res, class_index1, class_index2, ignore_weights=False):
1577        """
1578        Computes AUC between a pair of classes, ignoring instances from all
1579        other classes.
1580        """
1581        auc = AUC(ignore_weights=ignore_weights)
1582        auc._compute_for_pair_of_classes(res, class_index1, class_index2)
1583        return auc
1584
1585    @classmethod
1586    def matrix(cls, res, ignore_weights=False):
1587        """
1588        Compute a (lower diagonal) matrix with AUCs for all pairs of classes.
1589        If there are empty classes, the corresponding elements in the matrix
1590        are -1.
1591        """
1592        auc = AUC(ignore_weights=ignore_weights)
1593        auc._compute_matrix(res)
1594        return auc
1595
1596    def _compute_for_binary_class(self, res):
1597        """AUC for binary classification problems"""
1598        if res.number_of_iterations > 1:
1599            return self._compute_for_multiple_folds(
1600                self._compute_one_class_against_all,
1601                split_by_iterations(res),
1602                (-1, res, res.number_of_iterations))
1603        else:
1604            auc, _ = self._compute_one_class_against_all(res, -1)
1605            self[:] = auc
1606            return self
1607
1608    def _compute_for_multi_value_class(self, res, method=0):
1609        """AUC for multiclass classification problems"""
1610        numberOfClasses = len(res.class_values)
1611
1612        if res.number_of_iterations > 1:
1613            iterations = split_by_iterations(res)
1614            all_ite = res
1615        else:
1616            iterations = [res]
1617            all_ite = None
1618
1619        # by pairs
1620        sum_aucs = [0.] * res.number_of_learners
1621        usefulClassPairs = 0.
1622
1623        prob = None
1624        if method in [self.ByWeightedPairs, self.WeightedOneAgainstAll]:
1625            prob = class_probabilities_from_res(res)
1626
1627        if method in [self.ByWeightedPairs, self.ByPairs]:
1628            for classIndex1 in range(numberOfClasses):
1629                for classIndex2 in range(classIndex1):
1630                    subsum_aucs = self._compute_for_multiple_folds(
1631                        self._compute_one_class_against_another, iterations,
1632                        (classIndex1, classIndex2, all_ite,
1633                        res.number_of_iterations))
1634                    if subsum_aucs:
1635                        if method == self.ByWeightedPairs:
1636                            p_ij = prob[classIndex1] * prob[classIndex2]
1637                            subsum_aucs = [x * p_ij  for x in subsum_aucs]
1638                            usefulClassPairs += p_ij
1639                        else:
1640                            usefulClassPairs += 1
1641                        sum_aucs = map(add, sum_aucs, subsum_aucs)
1642        else:
1643            for classIndex in range(numberOfClasses):
1644                subsum_aucs = self._compute_for_multiple_folds(
1645                    self._compute_one_class_against_all,
1646                    iterations, (classIndex, all_ite,
1647                                 res.number_of_iterations))
1648                if subsum_aucs:
1649                    if method == self.ByWeightedPairs:
1650                        p_i = prob[classIndex]
1651                        subsum_aucs = [x * p_i  for x in subsum_aucs]
1652                        usefulClassPairs += p_i
1653                    else:
1654                        usefulClassPairs += 1
1655                    sum_aucs = map(add, sum_aucs, subsum_aucs)
1656
1657        if usefulClassPairs > 0:
1658            sum_aucs = [x/usefulClassPairs for x in sum_aucs]
1659
1660        self[:] = sum_aucs
1661        return self
1662
1663    # computes the average AUC over folds using "AUCcomputer" (AUC_i or AUC_ij)
1664    # it returns the sum of what is returned by the computer,
1665    # unless at a certain fold the computer has to resort to computing
1666    # over all folds or even this failed;
1667    # in these cases the result is returned immediately
1668    def _compute_for_multiple_folds(self, auc_computer, iterations,
1669                                 computer_args):
1670        """Compute the average AUC over folds using :obj:`auc_computer`."""
1671        subsum_aucs = [0.] * iterations[0].number_of_learners
1672        for ite in iterations:
1673            aucs, foldsUsed = auc_computer(*(ite, ) + computer_args)
1674            if not aucs:
1675                import warnings
1676                warnings.warn("AUC cannot be computed (all instances belong to the same class).")
1677                return
1678            if not foldsUsed:
1679                self[:] = aucs
1680                return aucs
1681            subsum_aucs = map(add, subsum_aucs, aucs)
1682        self[:] = subsum_aucs
1683        return self
1684
1685    # Computes AUC
1686    # in multivalued class problem, AUC is computed as one against all
1687    # results over folds are averages
1688    # if some folds examples from one class only, the folds are merged
1689    def _compute_for_single_class(self, res, class_index):
1690        if res.number_of_iterations > 1:
1691            self._compute_for_multiple_folds(
1692                self._compute_one_class_against_all, split_by_iterations(res),
1693                (class_index, res, res.number_of_iterations))
1694        else:
1695            self._compute_one_class_against_all(res, class_index)
1696
1697    # Computes AUC for a pair of classes (as if there were no other classes)
1698    # results over folds are averages
1699    # if some folds have examples from one class only, the folds are merged
1700    def _compute_for_pair_of_classes(self, res, class_index1, class_index2):
1701        if res.number_of_iterations > 1:
1702            self._compute_for_multiple_folds(
1703                self._compute_one_class_against_another,
1704                split_by_iterations(res),
1705                (class_index1, class_index2, res, res.number_of_iterations))
1706        else:
1707            self._compute_one_class_against_another(res, class_index1,
1708                                                    class_index2)
1709
1710    # computes AUC between class i and the other classes
1711    # (treating them as the same class)
1712    @deprecated_keywords({"classIndex": "class_index",
1713                          "divideByIfIte": "divide_by_if_ite"})
1714    def _compute_one_class_against_all(self, ite, class_index, all_ite=None,
1715                                      divide_by_if_ite=1.0):
1716        """Compute AUC between class i and all the other classes)"""
1717        return self._compute_auc(corn.computeCDT, ite, all_ite,
1718            divide_by_if_ite, (class_index, not self.ignore_weights))
1719
1720
1721    # computes AUC between classes i and j as if there are no other classes
1722    def _compute_one_class_against_another(
1723        self, ite, class_index1, class_index2,
1724        all_ite=None, divide_by_if_ite=1.):
1725        """
1726        Compute AUC between classes i and j as if there are no other classes.
1727        """
1728        return self._compute_auc(corn.computeCDTPair, ite,
1729            all_ite, divide_by_if_ite,
1730            (class_index1, class_index2, not self.ignore_weights))
1731
1732    # computes AUC using a specified 'cdtComputer' function
1733    # It tries to compute AUCs from 'ite' (examples from a single iteration)
1734    # and, if C+D+T=0, from 'all_ite' (entire test set). In the former case,
1735    # the AUCs are divided by 'divideByIfIte'.
1736    # Additional flag is returned which is True in the former case,
1737    # or False in the latter.
1738    @deprecated_keywords({"cdt_computer": "cdtComputer",
1739                          "divideByIfIte": "divide_by_if_ite",
1740                          "computerArgs": "computer_args"})
1741    def _compute_auc(self, cdt_computer, ite, all_ite, divide_by_if_ite,
1742                     computer_args):
1743        """
1744        Compute AUC using a :obj:`cdt_computer`.
1745        """
1746        cdts = cdt_computer(*(ite, ) + computer_args)
1747        if not is_CDT_empty(cdts[0]):
1748            return [(cdt.C + cdt.T / 2) / (cdt.C + cdt.D + cdt.T) /
1749                    divide_by_if_ite for cdt in cdts], True
1750
1751        if all_ite:
1752            cdts = cdt_computer(*(all_ite, ) + computer_args)
1753            if not is_CDT_empty(cdts[0]):
1754                return [(cdt.C + cdt.T / 2) / (cdt.C + cdt.D + cdt.T)
1755                        for cdt in cdts], False
1756
1757        return False, False
1758
1759    def _compute_matrix(self, res):
1760        numberOfClasses = len(res.class_values)
1761        number_of_learners = res.number_of_learners
1762        if res.number_of_iterations > 1:
1763            iterations, all_ite = split_by_iterations(res), res
1764        else:
1765            iterations, all_ite = [res], None
1766        aucs = [[[] for _ in range(numberOfClasses)]
1767                for _ in range(number_of_learners)]
1768        for classIndex1 in range(numberOfClasses):
1769            for classIndex2 in range(classIndex1):
1770                pair_aucs = self._compute_for_multiple_folds(
1771                    self._compute_one_class_against_another, iterations,
1772                    (classIndex1, classIndex2, all_ite,
1773                     res.number_of_iterations))
1774                if pair_aucs:
1775                    for lrn in range(number_of_learners):
1776                        aucs[lrn][classIndex1].append(pair_aucs[lrn])
1777                else:
1778                    for lrn in range(number_of_learners):
1779                        aucs[lrn][classIndex1].append(-1)
1780        self[:] = aucs
1781        return aucs
1782
1783#Backward compatibility
1784AUC.ByWeightedPairs = 0
1785AUC.ByPairs = 1
1786AUC.WeightedOneAgainstAll = 2
1787AUC.OneAgainstAll = 3
1788
1789@replace_use_weights
1790def AUC_binary(res, ignore_weights=False):
1791    auc = deprecated_function_name(AUC)(ignore_weights=ignore_weights)
1792    auc._compute_for_binary_class(res)
1793    return auc
1794
1795@replace_use_weights
1796def AUC_multi(res, ignore_weights=False, method=0):
1797    auc = deprecated_function_name(AUC)(ignore_weights=ignore_weights,
1798        method=method)
1799    auc._compute_for_multi_value_class(res)
1800    return auc
1801
1802
1803@deprecated_keywords({"AUCcomputer": "auc_computer",
1804                      "computerArgs": "computer_args"})
1805def AUC_iterations(auc_computer, iterations, computer_args):
1806    auc = deprecated_function_name(AUC)()
1807    auc._compute_for_multiple_folds(auc_computer, iterations, computer_args)
1808    return auc
1809
1810def AUC_x(cdtComputer, ite, all_ite, divide_by_if_ite, computer_args):
1811    auc = deprecated_function_name(AUC)()
1812    result = auc._compute_auc(cdtComputer, ite, all_ite, divide_by_if_ite,
1813                              computer_args)
1814    return result
1815
1816@replace_use_weights
1817def AUC_i(ite, class_index, ignore_weights=False, all_ite=None,
1818          divide_by_if_ite=1.):
1819    auc = deprecated_function_name(AUC)()
1820    result = auc._compute_one_class_against_another(ite, class_index,
1821        all_ite=None, divide_by_if_ite=1.)
1822    return result
1823
1824
1825@replace_use_weights
1826def AUC_ij(ite, class_index1, class_index2, ignore_weights=False,
1827           all_ite=None, divide_by_if_ite=1.):
1828    auc = deprecated_function_name(AUC)(ignore_weights=ignore_weights)
1829    result = auc._compute_one_class_against_another(
1830        ite, class_index1, class_index2, all_ite=None, divide_by_if_ite=1.)
1831    return result
1832
1833
1834
1835
1836#AUC_binary = replace_use_weights(deprecated_function_name(AUC()._compute_for_binary_class))
1837#AUC_multi = replace_use_weights(deprecated_function_name(AUC._compute_for_multi_value_class))
1838#AUC_iterations = replace_use_weights(deprecated_function_name(AUC._compute_for_multiple_folds))
1839#AUC_x = replace_use_weights(deprecated_function_name(AUC._compute_auc))
1840#AUC_i = replace_use_weights(deprecated_function_name(AUC._compute_one_class_against_all))
1841#AUC_ij = replace_use_weights(deprecated_function_name(AUC._compute_one_class_against_another))
1842
1843AUC_single = replace_use_weights(
1844             deprecated_keywords({"classIndex": "class_index"})(
1845             deprecated_function_name(AUC.single_class)))
1846AUC_pair = replace_use_weights(
1847           deprecated_keywords({"classIndex1": "class_index1",
1848                                "classIndex2": "class_index2"})(
1849           deprecated_function_name(AUC.pair)))
1850AUC_matrix = replace_use_weights(deprecated_function_name(AUC.matrix))
1851
1852
1853@deprecated_keywords({"unweighted": "ignore_weights"})
1854def McNemar(res, ignore_weights=False, **argkw):
1855    """
1856    Compute a triangular matrix with McNemar statistics for each pair of
1857    classifiers. The statistics is distributed by chi-square distribution with
1858    one degree of freedom; critical value for 5% significance is around 3.84.
1859    """
1860    nLearners = res.number_of_learners
1861    mcm = []
1862    for i in range(nLearners):
1863       mcm.append([0.] * res.number_of_learners)
1864
1865    if not res.weights or ignore_weights:
1866        for i in res.results:
1867            actual = i.actual_class
1868            classes = i.classes
1869            for l1 in range(nLearners):
1870                for l2 in range(l1, nLearners):
1871                    if classes[l1] == actual:
1872                        if classes[l2] != actual:
1873                            mcm[l1][l2] += 1
1874                    elif classes[l2] == actual:
1875                        mcm[l2][l1] += 1
1876    else:
1877        for i in res.results:
1878            actual = i.actual_class
1879            classes = i.classes
1880            for l1 in range(nLearners):
1881                for l2 in range(l1, nLearners):
1882                    if classes[l1] == actual:
1883                        if classes[l2] != actual:
1884                            mcm[l1][l2] += i.weight
1885                    elif classes[l2] == actual:
1886                        mcm[l2][l1] += i.weight
1887
1888    for l1 in range(nLearners):
1889        for l2 in range(l1, nLearners):
1890            su = mcm[l1][l2] + mcm[l2][l1]
1891            if su:
1892                mcm[l2][l1] = (abs(mcm[l1][l2] - mcm[l2][l1]) - 1)**2 / su
1893            else:
1894                mcm[l2][l1] = 0
1895
1896    for l1 in range(nLearners):
1897        mcm[l1]=mcm[l1][:l1]
1898
1899    return mcm
1900
1901
1902def McNemar_of_two(res, lrn1, lrn2, ignore_weights=False):
1903    """
1904    McNemar_of_two computes a McNemar statistics for a pair of classifier,
1905    specified by indices learner1 and learner2.
1906    """
1907    tf = ft = 0.
1908    if not res.weights or ignore_weights:
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 += i.weight
1914            elif i.classes[lrn2] == actual:
1915                    ft += i.weight
1916    else:
1917        for i in res.results:
1918            actual = i.actual_class
1919            if i.classes[lrn1] == actual:
1920                if i.classes[lrn2] != actual:
1921                    tf += 1.
1922            elif i.classes[lrn2] == actual:
1923                    ft += 1.
1924
1925    su = tf + ft
1926    if su:
1927        return (abs(tf - ft) - 1)**2 / su
1928    else:
1929        return 0
1930
1931
1932def Friedman(res, stat=CA):
1933    """
1934    Compare classifiers with Friedman test, treating folds as different examles.
1935    Returns F, p and average ranks.
1936    """
1937    res_split = split_by_iterations(res)
1938    res = [stat(r) for r in res_split]
1939   
1940    N = len(res)
1941    k = len(res[0])
1942    sums = [0.] * k
1943    for r in res:
1944        ranks = [k - x + 1 for x in statc.rankdata(r)]
1945        if stat == Brier_score: # reverse ranks for Brier_score (lower better)
1946            ranks = [k + 1 - x for x in ranks]
1947        sums = map(add, ranks, sums)
1948
1949    T = sum(x * x for x in sums)
1950    sums = [x / N for x in sums]
1951
1952    F = 12. / (N * k * (k + 1)) * T  - 3 * N * (k + 1)
1953
1954    return F, statc.chisqprob(F, k - 1), sums
1955
1956
1957def Wilcoxon_pairs(res, avgranks, stat=CA):
1958    """
1959    Return a triangular matrix, where element[i][j] stores significance of
1960    difference between the i-th and the j-th classifier, as computed by the
1961    Wilcoxon test. The element is positive if the i-th is better than the j-th,
1962    negative if it is worse, and 1 if they are equal.
1963    Arguments are ExperimentResults, average ranks (as returned by Friedman)
1964    and, optionally, a statistics; greater values should mean better results.
1965    """
1966    res_split = split_by_iterations(res)
1967    res = [stat(r) for r in res_split]
1968
1969    k = len(res[0])
1970    bt = []
1971    for m1 in range(k):
1972        nl = []
1973        for m2 in range(m1 + 1, k):
1974            t, p = statc.wilcoxont([r[m1] for r in res], [r[m2] for r in res])
1975            if avgranks[m1]<avgranks[m2]:
1976                nl.append(p)
1977            elif avgranks[m2]<avgranks[m1]:
1978                nl.append(-p)
1979            else:
1980                nl.append(1)
1981        bt.append(nl)
1982    return bt
1983
1984
1985@deprecated_keywords({"allResults": "all_results",
1986                      "noConfidence": "no_confidence"})
1987def plot_learning_curve_learners(file, all_results, proportions, learners,
1988                                 no_confidence=0):
1989    plot_learning_curve(file, all_results, proportions,
1990        [Orange.misc.getobjectname(learners[i], "Learner %i" % i)
1991        for i in range(len(learners))], no_confidence)
1992
1993
1994@deprecated_keywords({"allResults": "all_results",
1995                      "noConfidence": "no_confidence"})
1996def plot_learning_curve(file, all_results, proportions, legend,
1997                        no_confidence=0):
1998    import types
1999    fopened = 0
2000    if type(file) == types.StringType:
2001        file = open(file, "wt")
2002        fopened = 1
2003       
2004    file.write("set yrange [0:1]\n")
2005    file.write("set xrange [%f:%f]\n" % (proportions[0], proportions[-1]))
2006    file.write("set multiplot\n\n")
2007    CAs = [CA(x, report_se=True) for x in all_results]
2008
2009    file.write("plot \\\n")
2010    for i in range(len(legend) - 1):
2011        if not no_confidence:
2012            file.write("'-' title '' with yerrorbars pointtype %i,\\\n" % (i+1))
2013        file.write("'-' title '%s' with linespoints pointtype %i,\\\n" % (legend[i], i+1))
2014    if not no_confidence:
2015        file.write("'-' title '' with yerrorbars pointtype %i,\\\n" % (len(legend)))
2016    file.write("'-' title '%s' with linespoints pointtype %i\n" % (legend[-1], len(legend)))
2017
2018    for i in range(len(legend)):
2019        if not no_confidence:
2020            for p in range(len(proportions)):
2021                file.write("%f\t%f\t%f\n" % (proportions[p], CAs[p][i][0], 1.96*CAs[p][i][1]))
2022            file.write("e\n\n")
2023
2024        for p in range(len(proportions)):
2025            file.write("%f\t%f\n" % (proportions[p], CAs[p][i][0]))
2026        file.write("e\n\n")
2027
2028    if fopened:
2029        file.close()
2030
2031
2032def print_single_ROC_curve_coordinates(file, curve):
2033    import types
2034    fopened=0
2035    if type(file)==types.StringType:
2036        file=open(file, "wt")
2037        fopened=1
2038
2039    for coord in curve:
2040        file.write("%5.3f\t%5.3f\n" % tuple(coord))
2041
2042    if fopened:
2043        file.close()
2044
2045
2046def plot_ROC_learners(file, curves, learners):
2047    plot_ROC(file, curves, [Orange.misc.getobjectname(learners[i], "Learner %i" % i) for i in range(len(learners))])
2048   
2049def plot_ROC(file, curves, legend):
2050    import types
2051    fopened=0
2052    if type(file)==types.StringType:
2053        file=open(file, "wt")
2054        fopened=1
2055
2056    file.write("set yrange [0:1]\n")
2057    file.write("set xrange [0:1]\n")
2058    file.write("set multiplot\n\n")
2059
2060    file.write("plot \\\n")
2061    for leg in legend:
2062        file.write("'-' title '%s' with lines,\\\n" % leg)
2063    file.write("'-' title '' with lines\n")
2064
2065    for curve in curves:
2066        for coord in curve:
2067            file.write("%5.3f\t%5.3f\n" % tuple(coord))
2068        file.write("e\n\n")
2069
2070    file.write("1.0\t1.0\n0.0\t0.0e\n\n")         
2071
2072    if fopened:
2073        file.close()
2074
2075
2076@deprecated_keywords({"allResults": "all_results"})
2077def plot_McNemar_curve_learners(file, all_results, proportions, learners, reference=-1):
2078    plot_McNemar_curve(file, all_results, proportions, [Orange.misc.getobjectname(learners[i], "Learner %i" % i) for i in range(len(learners))], reference)
2079
2080
2081@deprecated_keywords({"allResults": "all_results"})
2082def plot_McNemar_curve(file, all_results, proportions, legend, reference=-1):
2083    if reference<0:
2084        reference=len(legend)-1
2085       
2086    import types
2087    fopened=0
2088    if type(file)==types.StringType:
2089        file=open(file, "wt")
2090        fopened=1
2091       
2092    #file.write("set yrange [0:1]\n")
2093    #file.write("set xrange [%f:%f]\n" % (proportions[0], proportions[-1]))
2094    file.write("set multiplot\n\n")
2095    file.write("plot \\\n")
2096    tmap=range(reference)+range(reference+1, len(legend))
2097    for i in tmap[:-1]:
2098        file.write("'-' title '%s' with linespoints pointtype %i,\\\n" % (legend[i], i+1))
2099    file.write("'-' title '%s' with linespoints pointtype %i\n" % (legend[tmap[-1]], tmap[-1]))
2100    file.write("\n")
2101
2102    for i in tmap:
2103        for p in range(len(proportions)):
2104            file.write("%f\t%f\n" % (proportions[p], McNemar_of_two(all_results[p], i, reference)))
2105        file.write("e\n\n")
2106
2107    if fopened:
2108        file.close()
2109
2110default_point_types=("{$\\circ$}", "{$\\diamond$}", "{$+$}", "{$\\times$}", "{$|$}")+tuple([chr(x) for x in range(97, 122)])
2111default_line_types=("\\setsolid", "\\setdashpattern <4pt, 2pt>", "\\setdashpattern <8pt, 2pt>", "\\setdashes", "\\setdots")
2112
2113@deprecated_keywords({"allResults": "all_results"})
2114def learning_curve_learners_to_PiCTeX(file, all_results, proportions, **options):
2115    return apply(learning_curve_to_PiCTeX, (file, all_results, proportions), options)
2116
2117
2118@deprecated_keywords({"allResults": "all_results"})
2119def learning_curve_to_PiCTeX(file, all_results, proportions, **options):
2120    import types
2121    fopened=0
2122    if type(file)==types.StringType:
2123        file=open(file, "wt")
2124        fopened=1
2125
2126    nexamples=len(all_results[0].results)
2127    CAs = [CA(x, report_se=True) for x in all_results]
2128
2129    graphsize=float(options.get("graphsize", 10.0)) #cm
2130    difprop=proportions[-1]-proportions[0]
2131    ntestexamples=nexamples*proportions[-1]
2132    xunit=graphsize/ntestexamples
2133
2134    yshift=float(options.get("yshift", -ntestexamples/20.))
2135   
2136    pointtypes=options.get("pointtypes", default_point_types)
2137    linetypes=options.get("linetypes", default_line_types)
2138
2139    if options.has_key("numberedx"):
2140        numberedx=options["numberedx"]
2141        if type(numberedx)==types.IntType:
2142            if numberedx>0:
2143                numberedx=[nexamples*proportions[int(i/float(numberedx)*len(proportions))] for i in range(numberedx)]+[proportions[-1]*nexamples]
2144            elif numberedx<0:
2145                numberedx = -numberedx
2146                newn=[]
2147                for i in range(numberedx+1):
2148                    wanted=proportions[0]+float(i)/numberedx*difprop
2149                    best=(10, 0)
2150                    for t in proportions:
2151                        td=abs(wanted-t)
2152                        if td<best[0]:
2153                            best=(td, t)
2154                    if not best[1] in newn:
2155                        newn.append(best[1])
2156                newn.sort()
2157                numberedx=[nexamples*x for x in newn]
2158        elif type(numberedx[0])==types.FloatType:
2159            numberedx=[nexamples*x for x in numberedx]
2160    else:
2161        numberedx=[nexamples*x for x in proportions]
2162
2163    file.write("\\mbox{\n")
2164    file.write(\\beginpicture\n")
2165    file.write(\\setcoordinatesystem units <%10.8fcm, %5.3fcm>\n\n" % (xunit, graphsize))   
2166    file.write(\\setplotarea x from %5.3f to %5.3f, y from 0 to 1\n" % (0, ntestexamples))   
2167    file.write(\\axis bottom invisible\n")# label {#examples}\n")
2168    file.write("      ticks short at %s /\n" % reduce(lambda x,y:x+" "+y, ["%i"%(x*nexamples+0.5) for x in proportions]))
2169    if numberedx:
2170        file.write("            long numbered at %s /\n" % reduce(lambda x,y:x+y, ["%i " % int(x+0.5) for x in numberedx]))
2171    file.write("  /\n")
2172    file.write(\\axis left invisible\n")# label {classification accuracy}\n")
2173    file.write("      shiftedto y=%5.3f\n" % yshift)
2174    file.write("      ticks short from 0.0 to 1.0 by 0.05\n")
2175    file.write("            long numbered from 0.0 to 1.0 by 0.25\n")
2176    file.write("  /\n")
2177    if options.has_key("default"):
2178        file.write(\\setdashpattern<1pt, 1pt>\n")
2179        file.write(\\plot %5.3f %5.3f %5.3f %5.3f /\n" % (0., options["default"], ntestexamples, options["default"]))
2180   
2181    for i in range(len(CAs[0])):
2182        coordinates=reduce(lambda x,y:x+" "+y, ["%i %5.3f" % (proportions[p]*nexamples, CAs[p][i][0]) for p in range(len(proportions))])
2183        if linetypes:
2184            file.write(%s\n" % linetypes[i])
2185            file.write(\\plot %s /\n" % coordinates)
2186        if pointtypes:
2187            file.write(\\multiput %s at %s /\n" % (pointtypes[i], coordinates))
2188
2189    file.write(\\endpicture\n")
2190    file.write("}\n")
2191    if fopened:
2192        file.close()
2193    file.close()
2194    del file
2195
2196def legend_learners_to_PiCTeX(file, learners, **options):
2197  return apply(legend_to_PiCTeX, (file, [Orange.misc.getobjectname(learners[i], "Learner %i" % i) for i in range(len(learners))]), options)
2198   
2199def legend_to_PiCTeX(file, legend, **options):
2200    import types
2201    fopened=0
2202    if type(file)==types.StringType:
2203        file=open(file, "wt")
2204        fopened=1
2205
2206    pointtypes=options.get("pointtypes", default_point_types)
2207    linetypes=options.get("linetypes", default_line_types)
2208
2209    file.write("\\mbox{\n")
2210    file.write(\\beginpicture\n")
2211    file.write(\\setcoordinatesystem units <5cm, 1pt>\n\n")
2212    file.write(\\setplotarea x from 0.000 to %5.3f, y from 0 to 12\n" % len(legend))
2213
2214    for i in range(len(legend)):
2215        if linetypes:
2216            file.write(%s\n" % linetypes[i])
2217            file.write(\\plot %5.3f 6 %5.3f 6 /\n" % (i, i+0.2))
2218        if pointtypes:
2219            file.write(\\put {%s} at %5.3f 6\n" % (pointtypes[i], i+0.1))
2220        file.write(\\put {%s} [lb] at %5.3f 0\n" % (legend[i], i+0.25))
2221
2222    file.write(\\endpicture\n")
2223    file.write("}\n")
2224    if fopened:
2225        file.close()
2226    file.close()
2227    del file
2228
2229
2230def compute_friedman(avranks, N):
2231    """ Returns a tuple composed of (friedman statistic, degrees of freedom)
2232    and (Iman statistic - F-distribution, degrees of freedoma) given average
2233    ranks and a number of tested data sets N.
2234    """
2235
2236    k = len(avranks)
2237
2238    def friedman(N, k, ranks):
2239        return 12*N*(sum([rank**2.0 for rank in ranks]) - (k*(k+1)*(k+1)/4.0) )/(k*(k+1))
2240
2241    def iman(fried, N, k):
2242        return (N-1)*fried/(N*(k-1) - fried)
2243
2244    f = friedman(N, k, avranks)
2245    im = iman(f, N, k)
2246    fdistdof = (k-1, (k-1)*(N-1))
2247
2248    return (f, k-1), (im, fdistdof)
2249
2250def compute_CD(avranks, N, alpha="0.05", type="nemenyi"):
2251    """ Returns critical difference for Nemenyi or Bonferroni-Dunn test
2252    according to given alpha (either alpha="0.05" or alpha="0.1") for average
2253    ranks and number of tested data sets N. Type can be either "nemenyi" for
2254    for Nemenyi two tailed test or "bonferroni-dunn" for Bonferroni-Dunn test.
2255    """
2256
2257    k = len(avranks)
2258   
2259    d = {("nemenyi", "0.05"): [0, 0, 1.959964, 2.343701, 2.569032, 2.727774,
2260                               2.849705, 2.94832, 3.030879, 3.101730, 3.163684,
2261                               3.218654, 3.268004, 3.312739, 3.353618, 3.39123,
2262                               3.426041, 3.458425, 3.488685, 3.517073, 3.543799]
2263        , ("nemenyi", "0.1"): [0, 0, 1.644854, 2.052293, 2.291341, 2.459516,
2264                               2.588521, 2.692732, 2.779884, 2.854606, 2.919889,
2265                               2.977768, 3.029694, 3.076733, 3.119693, 3.159199,
2266                               3.195743, 3.229723, 3.261461, 3.291224, 3.319233]
2267        , ("bonferroni-dunn", "0.05"): [0, 0, 1.960, 2.241, 2.394, 2.498, 2.576,
2268                                        2.638, 2.690, 2.724, 2.773],
2269         ("bonferroni-dunn", "0.1"): [0, 0, 1.645, 1.960, 2.128, 2.241, 2.326,
2270                                      2.394, 2.450, 2.498, 2.539]}
2271
2272    #can be computed in R as qtukey(0.95, n, Inf)**0.5
2273    #for (x in c(2:20)) print(qtukey(0.95, x, Inf)/(2**0.5)
2274
2275    q = d[(type, alpha)]
2276
2277    cd = q[k]*(k*(k+1)/(6.0*N))**0.5
2278
2279    return cd
2280 
2281
2282def graph_ranks(filename, avranks, names, cd=None, cdmethod=None, lowv=None, highv=None, width=6, textspace=1, reverse=False, **kwargs):
2283    """
2284    Draws a CD graph, which is used to display  the differences in methods'
2285    performance.
2286    See Janez Demsar, Statistical Comparisons of Classifiers over
2287    Multiple Data Sets, 7(Jan):1--30, 2006.
2288
2289    Needs matplotlib to work.
2290
2291    :param filename: Output file name (with extension). Formats supported
2292                     by matplotlib can be used.
2293    :param avranks: List of average methods' ranks.
2294    :param names: List of methods' names.
2295
2296    :param cd: Critical difference. Used for marking methods that whose
2297               difference is not statistically significant.
2298    :param lowv: The lowest shown rank, if None, use 1.
2299    :param highv: The highest shown rank, if None, use len(avranks).
2300    :param width: Width of the drawn figure in inches, default 6 in.
2301    :param textspace: Space on figure sides left for the description
2302                      of methods, default 1 in.
2303    :param reverse:  If True, the lowest rank is on the right. Default\: False.
2304    :param cdmethod: None by default. It can be an index of element in avranks
2305                     or or names which specifies the method which should be
2306                     marked with an interval.
2307    """
2308    if not HAS_MATPLOTLIB:
2309        import sys
2310        print >> sys.stderr, "Function requires matplotlib. Please install it."
2311        return
2312
2313    width = float(width)
2314    textspace = float(textspace)
2315
2316    def nth(l,n):
2317        """
2318        Returns only nth elemnt in a list.
2319        """
2320        n = lloc(l,n)
2321        return [ a[n] for a in l ]
2322
2323    def lloc(l,n):
2324        """
2325        List location in list of list structure.
2326        Enable the use of negative locations:
2327        -1 is the last element, -2 second last...
2328        """
2329        if n < 0:
2330            return len(l[0])+n
2331        else:
2332            return n
2333
2334    def mxrange(lr):
2335        """
2336        Multiple xranges. Can be used to traverse matrices.
2337        This function is very slow due to unknown number of
2338        parameters.
2339
2340        >>> mxrange([3,5])
2341        [(0, 0), (0, 1), (0, 2), (1, 0), (1, 1), (1, 2)]
2342
2343        >>> mxrange([[3,5,1],[9,0,-3]])
2344        [(3, 9), (3, 6), (3, 3), (4, 9), (4, 6), (4, 3)]
2345
2346        """
2347        if not len(lr):
2348            yield ()
2349        else:
2350            #it can work with single numbers
2351            index = lr[0]
2352            if type(1) == type(index):
2353                index = [ index ]
2354            for a in range(*index):
2355                for b in mxrange(lr[1:]):
2356                    yield tuple([a] + list(b))
2357
2358
2359    from matplotlib.figure import Figure
2360    from matplotlib.backends.backend_agg import FigureCanvasAgg
2361
2362    def print_figure(fig, *args, **kwargs):
2363        canvas = FigureCanvasAgg(fig)
2364        canvas.print_figure(*args, **kwargs)
2365
2366    sums = avranks
2367
2368    tempsort =  sorted([ (a,i) for i,a in  enumerate(sums) ], reverse=reverse)
2369    ssums = nth(tempsort, 0)
2370    sortidx = nth(tempsort, 1)
2371    nnames = [ names[x] for x in sortidx ]
2372   
2373    if lowv is None:
2374        lowv = min(1, int(math.floor(min(ssums))))
2375    if highv is None:
2376        highv = max(len(avranks), int(math.ceil(max(ssums))))
2377
2378    cline = 0.4
2379
2380    k = len(sums)
2381
2382    lines = None
2383
2384    linesblank = 0
2385    scalewidth = width - 2*textspace
2386
2387    def rankpos(rank):
2388        if not reverse:
2389            a = rank - lowv
2390        else:
2391            a = highv - rank
2392        return textspace+scalewidth/(highv-lowv)*a
2393
2394    distanceh = 0.25
2395
2396    if cd and cdmethod is None:
2397   
2398        #get pairs of non significant methods
2399
2400        def get_lines(sums, hsd):
2401
2402            #get all pairs
2403            lsums = len(sums)
2404            allpairs = [ (i,j) for i,j in mxrange([[lsums], [lsums]]) if j > i ]
2405            #remove not significant
2406            notSig = [ (i,j) for i,j in allpairs if abs(sums[i]-sums[j]) <= hsd ]
2407            #keep only longest
2408           
2409            def no_longer((i,j), notSig):
2410                for i1,j1 in notSig:
2411                    if (i1 <= i and j1 > j) or (i1 < i and j1 >= j):
2412                        return False
2413                return True
2414
2415            longest = [ (i,j) for i,j in notSig if no_longer((i,j),notSig) ]
2416           
2417            return longest
2418
2419        lines = get_lines(ssums, cd)
2420        linesblank = 0.2 + 0.2 + (len(lines)-1)*0.1
2421
2422        #add scale
2423        distanceh = 0.25
2424        cline += distanceh
2425
2426    #calculate height needed height of an image
2427    minnotsignificant = max(2*0.2, linesblank)
2428    height = cline + ((k+1)/2)*0.2 + minnotsignificant
2429
2430    fig = Figure(figsize=(width, height))
2431    ax = fig.add_axes([0,0,1,1]) #reverse y axis
2432    ax.set_axis_off()
2433
2434    hf = 1./height # height factor
2435    wf = 1./width
2436
2437    def hfl(l): 
2438        return [ a*hf for a in l ]
2439
2440    def wfl(l): 
2441        return [ a*wf for a in l ]
2442
2443
2444    # Upper left corner is (0,0).
2445
2446    ax.plot([0,1], [0,1], c="w")
2447    ax.set_xlim(0, 1)
2448    ax.set_ylim(1, 0)
2449
2450    def line(l, color='k', **kwargs):
2451        """
2452        Input is a list of pairs of points.
2453        """
2454        ax.plot(wfl(nth(l,0)), hfl(nth(l,1)), color=color, **kwargs)
2455
2456    def text(x, y, s, *args, **kwargs):
2457        ax.text(wf*x, hf*y, s, *args, **kwargs)
2458
2459    line([(textspace, cline), (width-textspace, cline)], linewidth=0.7)
2460   
2461    bigtick = 0.1
2462    smalltick = 0.05
2463
2464
2465    import numpy
2466    tick = None
2467    for a in list(numpy.arange(lowv, highv, 0.5)) + [highv]:
2468        tick = smalltick
2469        if a == int(a): tick = bigtick
2470        line([(rankpos(a), cline-tick/2),(rankpos(a), cline)], linewidth=0.7)
2471
2472    for a in range(lowv, highv+1):
2473        text(rankpos(a), cline-tick/2-0.05, str(a), ha="center", va="bottom")
2474
2475    k = len(ssums)
2476
2477    for i in range((k+1)/2):
2478        chei = cline+ minnotsignificant + i *0.2
2479        line([(rankpos(ssums[i]), cline), (rankpos(ssums[i]), chei), (textspace-0.1, chei)], linewidth=0.7)
2480        text(textspace-0.2, chei, nnames[i], ha="right", va="center")
2481
2482    for i in range((k+1)/2, k):
2483        chei = cline + minnotsignificant + (k-i-1)*0.2
2484        line([(rankpos(ssums[i]), cline), (rankpos(ssums[i]), chei), (textspace+scalewidth+0.1, chei)], linewidth=0.7)
2485        text(textspace+scalewidth+0.2, chei, nnames[i], ha="left", va="center")
2486
2487    if cd and cdmethod is None:
2488
2489        #upper scale
2490        if not reverse:
2491            begin, end = rankpos(lowv), rankpos(lowv+cd)
2492        else:
2493            begin, end = rankpos(highv), rankpos(highv - cd)
2494           
2495        line([(begin, distanceh), (end, distanceh)], linewidth=0.7)
2496        line([(begin, distanceh + bigtick/2), (begin, distanceh - bigtick/2)], linewidth=0.7)
2497        line([(end, distanceh + bigtick/2), (end, distanceh - bigtick/2)], linewidth=0.7)
2498        text((begin+end)/2, distanceh - 0.05, "CD", ha="center", va="bottom")
2499
2500        #non significance lines   
2501        def draw_lines(lines, side=0.05, height=0.1):
2502            start = cline + 0.2
2503            for l,r in lines: 
2504                line([(rankpos(ssums[l])-side, start), (rankpos(ssums[r])+side, start)], linewidth=2.5) 
2505                start += height
2506
2507        draw_lines(lines)
2508
2509    elif cd:
2510        begin = rankpos(avranks[cdmethod]-cd)
2511        end = rankpos(avranks[cdmethod]+cd)
2512        line([(begin, cline), (end, cline)], linewidth=2.5) 
2513        line([(begin, cline + bigtick/2), (begin, cline - bigtick/2)], linewidth=2.5)
2514        line([(end, cline + bigtick/2), (end, cline - bigtick/2)], linewidth=2.5)
2515 
2516    print_figure(fig, filename, **kwargs)
2517
2518def mlc_hamming_loss(res):
2519    """
2520    Schapire and Singer (2000) presented Hamming Loss, which id defined as:
2521   
2522    :math:`HammingLoss(H,D)=\\frac{1}{|D|} \\sum_{i=1}^{|D|} \\frac{Y_i \\vartriangle Z_i}{|L|}`
2523    """
2524    losses = [0]*res.number_of_learners
2525    label_num = len(res.labels)
2526    example_num = gettotsize(res)
2527
2528    for e in res.results:
2529        aclass = e.actual_class
2530        for i, labels in enumerate(e.classes):
2531            labels = map(int, labels)
2532            if len(labels) <> len(aclass):
2533                raise ValueError, "The dimensions of the classified output and the actual class array do not match."
2534            for j in range(label_num):
2535                if labels[j] != aclass[j]:
2536                    losses[i] += 1
2537           
2538    return [float(x)/(label_num*example_num) for x in losses]
2539
2540def mlc_accuracy(res, forgiveness_rate = 1.0):
2541    """
2542    Godbole & Sarawagi, 2004 uses the metrics accuracy, precision, recall as follows:
2543     
2544    :math:`Accuracy(H,D)=\\frac{1}{|D|} \\sum_{i=1}^{|D|} \\frac{|Y_i \\cap Z_i|}{|Y_i \\cup Z_i|}`
2545   
2546    Boutell et al. (2004) give a more generalized version using a parameter :math:`\\alpha \\ge 0`,
2547    called forgiveness rate:
2548   
2549    :math:`Accuracy(H,D)=\\frac{1}{|D|} \\sum_{i=1}^{|D|} (\\frac{|Y_i \\cap Z_i|}{|Y_i \\cup Z_i|})^{\\alpha}`
2550    """
2551    accuracies = [0.0]*res.number_of_learners
2552    example_num = gettotsize(res)
2553   
2554    for e in res.results:
2555        aclass = e.actual_class
2556        for i, labels in enumerate(e.classes):
2557            labels = map(int, labels)
2558            if len(labels) <> len(aclass):
2559                raise ValueError, "The dimensions of the classified output and the actual class array do not match."
2560           
2561            intersection = 0.0
2562            union = 0.0
2563            for real, pred in zip(labels, aclass):
2564                if real and pred:
2565                    intersection += 1
2566                if real or pred:
2567                    union += 1
2568
2569            if union:
2570                accuracies[i] += intersection / union
2571           
2572    return [math.pow(x/example_num,forgiveness_rate) for x in accuracies]
2573
2574def mlc_precision(res):
2575    """
2576    :math:`Precision(H,D)=\\frac{1}{|D|} \\sum_{i=1}^{|D|} \\frac{|Y_i \\cap Z_i|}{|Z_i|}`
2577    """
2578    precisions = [0.0]*res.number_of_learners
2579    example_num = gettotsize(res)
2580   
2581    for e in res.results:
2582        aclass = e.actual_class
2583        for i, labels in enumerate(e.classes):
2584            labels = map(int, labels)
2585            if len(labels) <> len(aclass):
2586                raise ValueError, "The dimensions of the classified output and the actual class array do not match."
2587           
2588            intersection = 0.0
2589            predicted = 0.0
2590            for real, pred in zip(labels, aclass):
2591                if real and pred:
2592                    intersection += 1
2593                if real:
2594                    predicted += 1
2595            if predicted:
2596                precisions[i] += intersection / predicted
2597           
2598    return [x/example_num for x in precisions]
2599
2600def mlc_recall(res):
2601    """
2602    :math:`Recall(H,D)=\\frac{1}{|D|} \\sum_{i=1}^{|D|} \\frac{|Y_i \\cap Z_i|}{|Y_i|}`
2603    """
2604    recalls = [0.0]*res.number_of_learners
2605    example_num = gettotsize(res)
2606   
2607    for e in res.results:
2608        aclass = e.actual_class
2609        for i, labels in enumerate(e.classes):
2610            labels = map(int, labels)
2611            if len(labels) <> len(aclass):
2612                raise ValueError, "The dimensions of the classified output and the actual class array do not match."
2613           
2614            intersection = 0.0
2615            actual = 0.0
2616            for real, pred in zip(labels, aclass):
2617                if real and pred:
2618                    intersection += 1
2619                if pred:
2620                    actual += 1
2621            if actual:
2622                recalls[i] += intersection / actual
2623           
2624    return [x/example_num for x in recalls]
2625
2626#def mlc_ranking_loss(res):
2627#    pass
2628#
2629#def mlc_average_precision(res):
2630#    pass
2631#
2632#def mlc_hierarchical_loss(res):
2633#    pass
2634
2635
2636def mt_average_score(res, score, weights=None):
2637    """
2638    Compute individual scores for each target and return the (weighted) average.
2639
2640    One method can be used to compute scores for all targets or a list of
2641    scoring methods can be passed to use different methods for different
2642    targets. In the latter case, care has to be taken if the ranges of scoring
2643    methods differ.
2644    For example, when the first target is scored from -1 to 1 (1 best) and the
2645    second from 0 to 1 (0 best), using `weights=[0.5,-1]` would scale both
2646    to a span of 1, and invert the second so that higher scores are better.
2647
2648    :param score: Single-target scoring method or a list of such methods
2649                  (one for each target).
2650    :param weights: List of real weights, one for each target,
2651                    for a weighted average.
2652
2653    """
2654    if not len(res.results):
2655        raise ValueError, "Cannot compute the score: no examples."
2656    if res.number_of_learners < 1:
2657        return []
2658    n_classes = len(res.results[0].actual_class)
2659    if weights is None:
2660        weights = [1.] * n_classes
2661    if not hasattr(score, '__len__'):
2662        score = [score] * n_classes
2663    elif len(score) != n_classes:
2664        raise ValueError, "Number of scoring methods and targets do not match."
2665    # save original classes
2666    clsss = [te.classes for te in res.results]
2667    aclsss = [te.actual_class for te in res.results]
2668    # compute single target scores
2669    single_scores = []
2670    for i in range(n_classes):
2671        for te, clss, aclss in zip(res.results, clsss, aclsss):
2672            te.classes = [cls[i] for cls in clss]
2673            te.actual_class = aclss[i]
2674        single_scores.append(score[i](res))
2675    # restore original classes
2676    for te, clss, aclss in zip(res.results, clsss, aclsss):
2677        te.classes = clss
2678        te.actual_class = aclss
2679    return [sum(w * s for w, s in zip(weights, scores)) / sum(weights)
2680        for scores in zip(*single_scores)]
2681
2682def mt_flattened_score(res, score):
2683    """
2684    Flatten (concatenate into a single list) the predictions of multiple
2685    targets and compute a single-target score.
2686   
2687    :param score: Single-target scoring method.
2688    """
2689    res2 = Orange.evaluation.testing.ExperimentResults(res.number_of_iterations,
2690        res.classifier_names, class_values=res.class_values,
2691        weights=res.weights, classifiers=res.classifiers, loaded=res.loaded,
2692        test_type=Orange.evaluation.testing.TEST_TYPE_SINGLE, labels=res.labels)
2693    for te in res.results:
2694        for i, ac in enumerate(te.actual_class):
2695            te2 = Orange.evaluation.testing.TestedExample(
2696                iteration_number=te.iteration_number, actual_class=ac)
2697            for c, p in zip(te.classes, te.probabilities):
2698                te2.add_result(c[i], p[i])
2699            res2.results.append(te2)
2700    return score(res2)
2701
2702
2703################################################################################
2704if __name__ == "__main__":
2705    avranks =  [3.143, 2.000, 2.893, 1.964]
2706    names = ["prva", "druga", "tretja", "cetrta" ]
2707    cd = compute_CD(avranks, 14)
2708    #cd = compute_CD(avranks, 10, type="bonferroni-dunn")
2709    print cd
2710
2711    print compute_friedman(avranks, 14)
2712
2713    #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.