source: orange/Orange/evaluation/scoring.py @ 10426:ce57e8dbcc18

Revision 10426:ce57e8dbcc18, 97.2 KB checked in by anzeh <anze.staric@…>, 2 years ago (diff)

Refactored AUC.

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    If testing consisted of multiple folds, each fold is scored and the
1475    average score is returned. If a fold contains only instances with the
1476    same class value, folds will be merged.
1477
1478    :param test_results: test results to score
1479    :param ignore_weights: ignore instance weights when calculating score
1480    :param multiclass: tells what kind of averaging to perform if the target
1481                       class has more than 2 values.
1482    """
1483
1484    #!Compute AUC for each pair of classes (ignoring instances of all other
1485    #!classes) and average the results, weighting them by the number of
1486    #!pairs of instances from these two classes (e.g. by the product of
1487    #!probabilities of the two classes). AUC computed in this way still
1488    #!behaves as the concordance index, e.g., gives the probability that two
1489    #!randomly chosen instances from different classes will be correctly
1490    #!recognized (if the classifier knows from which two classes the
1491    #!instances came).
1492    ByWeightedPairs = 0
1493
1494    #!Similar to ByWeightedPairs, except that the average over class pairs
1495    #!is not weighted. This AUC is, like the binary version, independent of
1496    #!class distributions, but it is not related to the concordance index
1497    #!any more.
1498    ByPairs = 1
1499
1500    #!For each class, it computes AUC for this class against all others (that
1501    #!is, treating other classes as one class). The AUCs are then averaged by
1502    #!the class probabilities. This is related to the concordance index in
1503    #!which we test the classifier's (average) capability of distinguishing
1504    #!the instances from a specified class from those that come from other
1505    #!classes.
1506    #!Unlike the binary AUC, the measure is not independent of class
1507    #!distributions.
1508    WeightedOneAgainstAll = 2
1509
1510    #!Similar to weighted_one_against_all, except that the average
1511    #!is not weighted.
1512    OneAgainstAll = 3
1513
1514    @replace_use_weights
1515    @deprecated_keywords({"method": "multiclass"})
1516    def __init__(self, test_results=None, multiclass=ByWeightedPairs, ignore_weights=False):
1517
1518        super(AUC, self).__init__()
1519
1520        self.ignore_weights=ignore_weights
1521        self.method=multiclass
1522
1523        if test_results is not None:
1524            self[:] = self.__call__(test_results)
1525
1526    def __call__(self, test_results):
1527        if len(test_results.class_values) < 2:
1528            raise ValueError("Cannot compute AUC on a single-class problem")
1529        elif len(test_results.class_values) == 2:
1530            return self._compute_for_binary_class(test_results)
1531        else:
1532            return self._compute_for_multi_value_class(test_results, self.method)
1533
1534    def _compute_for_binary_class(self, res):
1535        """AUC for binary classification problems"""
1536        if res.number_of_iterations > 1:
1537            return self._compute_for_multiple_folds(
1538                self._compute_one_class_against_all,
1539                split_by_iterations(res),
1540                (-1, res, res.number_of_iterations))
1541        else:
1542            return self._compute_one_class_against_all(res, -1)[0]
1543
1544    def _compute_for_multi_value_class(self, res, method=0):
1545        """AUC for multiclass classification problems"""
1546        numberOfClasses = len(res.class_values)
1547
1548        if res.number_of_iterations > 1:
1549            iterations = split_by_iterations(res)
1550            all_ite = res
1551        else:
1552            iterations = [res]
1553            all_ite = None
1554
1555        # by pairs
1556        sum_aucs = [0.] * res.number_of_learners
1557        usefulClassPairs = 0.
1558
1559        prob = None
1560        if method in [self.ByWeightedPairs, self.WeightedOneAgainstAll]:
1561            prob = class_probabilities_from_res(res)
1562
1563        if method in [self.ByWeightedPairs, self.ByPairs]:
1564            for classIndex1 in range(numberOfClasses):
1565                for classIndex2 in range(classIndex1):
1566                    subsum_aucs = self._compute_for_multiple_folds(
1567                        self._compute_one_class_against_another, iterations,
1568                        (classIndex1, classIndex2, all_ite,
1569                        res.number_of_iterations))
1570                    if subsum_aucs:
1571                        if method == self.ByWeightedPairs:
1572                            p_ij = prob[classIndex1] * prob[classIndex2]
1573                            subsum_aucs = [x * p_ij  for x in subsum_aucs]
1574                            usefulClassPairs += p_ij
1575                        else:
1576                            usefulClassPairs += 1
1577                        sum_aucs = map(add, sum_aucs, subsum_aucs)
1578        else:
1579            for classIndex in range(numberOfClasses):
1580                subsum_aucs = self._compute_for_multiple_folds(
1581                    self._compute_one_class_against_all,
1582                    iterations, (classIndex, all_ite,
1583                                 res.number_of_iterations))
1584                if subsum_aucs:
1585                    if method == self.ByWeightedPairs:
1586                        p_i = prob[classIndex]
1587                        subsum_aucs = [x * p_i  for x in subsum_aucs]
1588                        usefulClassPairs += p_i
1589                    else:
1590                        usefulClassPairs += 1
1591                    sum_aucs = map(add, sum_aucs, subsum_aucs)
1592
1593        if usefulClassPairs > 0:
1594            sum_aucs = [x/usefulClassPairs for x in sum_aucs]
1595
1596        return sum_aucs
1597
1598    # computes the average AUC over folds using "AUCcomputer" (AUC_i or AUC_ij)
1599    # it returns the sum of what is returned by the computer,
1600    # unless at a certain fold the computer has to resort to computing
1601    # over all folds or even this failed;
1602    # in these cases the result is returned immediately
1603    def _compute_for_multiple_folds(self, auc_computer, iterations,
1604                                 computer_args):
1605        """Compute the average AUC over folds using :obj:`auc_computer`."""
1606        subsum_aucs = [0.] * iterations[0].number_of_learners
1607        for ite in iterations:
1608            aucs, foldsUsed = auc_computer(*(ite, ) + computer_args)
1609            if not aucs:
1610                import warnings
1611                warnings.warn("AUC cannot be computed (all instances belong to the same class).")
1612                return
1613            if not foldsUsed:
1614                self[:] = aucs
1615                return aucs
1616            subsum_aucs = map(add, subsum_aucs, aucs)
1617        return subsum_aucs
1618
1619    # Computes AUC
1620    # in multivalued class problem, AUC is computed as one against all
1621    # results over folds are averages
1622    # if some folds examples from one class only, the folds are merged
1623    def _compute_for_single_class(self, res, class_index):
1624        if res.number_of_iterations > 1:
1625            return self._compute_for_multiple_folds(
1626                self._compute_one_class_against_all, split_by_iterations(res),
1627                (class_index, res, res.number_of_iterations))
1628        else:
1629            return self._compute_one_class_against_all(res, class_index)
1630
1631    # Computes AUC for a pair of classes (as if there were no other classes)
1632    # results over folds are averages
1633    # if some folds have examples from one class only, the folds are merged
1634    def _compute_for_pair_of_classes(self, res, class_index1, class_index2):
1635        if res.number_of_iterations > 1:
1636            return self._compute_for_multiple_folds(
1637                self._compute_one_class_against_another,
1638                split_by_iterations(res),
1639                (class_index1, class_index2, res, res.number_of_iterations))
1640        else:
1641            return self._compute_one_class_against_another(res, class_index1,
1642                                                    class_index2)
1643
1644    # computes AUC between class i and the other classes
1645    # (treating them as the same class)
1646    @deprecated_keywords({"classIndex": "class_index",
1647                          "divideByIfIte": "divide_by_if_ite"})
1648    def _compute_one_class_against_all(self, ite, class_index, all_ite=None,
1649                                      divide_by_if_ite=1.0):
1650        """Compute AUC between class i and all the other classes)"""
1651        return self._compute_auc(corn.computeCDT, ite, all_ite,
1652            divide_by_if_ite, (class_index, not self.ignore_weights))
1653
1654
1655    # computes AUC between classes i and j as if there are no other classes
1656    def _compute_one_class_against_another(
1657        self, ite, class_index1, class_index2,
1658        all_ite=None, divide_by_if_ite=1.):
1659        """
1660        Compute AUC between classes i and j as if there are no other classes.
1661        """
1662        return self._compute_auc(corn.computeCDTPair, ite,
1663            all_ite, divide_by_if_ite,
1664            (class_index1, class_index2, not self.ignore_weights))
1665
1666    # computes AUC using a specified 'cdtComputer' function
1667    # It tries to compute AUCs from 'ite' (examples from a single iteration)
1668    # and, if C+D+T=0, from 'all_ite' (entire test set). In the former case,
1669    # the AUCs are divided by 'divideByIfIte'.
1670    # Additional flag is returned which is True in the former case,
1671    # or False in the latter.
1672    @deprecated_keywords({"cdt_computer": "cdtComputer",
1673                          "divideByIfIte": "divide_by_if_ite",
1674                          "computerArgs": "computer_args"})
1675    def _compute_auc(self, cdt_computer, ite, all_ite, divide_by_if_ite,
1676                     computer_args):
1677        """
1678        Compute AUC using a :obj:`cdt_computer`.
1679        """
1680        cdts = cdt_computer(*(ite, ) + computer_args)
1681        if not is_CDT_empty(cdts[0]):
1682            return [(cdt.C + cdt.T / 2) / (cdt.C + cdt.D + cdt.T) /
1683                    divide_by_if_ite for cdt in cdts], True
1684
1685        if all_ite:
1686            cdts = cdt_computer(*(all_ite, ) + computer_args)
1687            if not is_CDT_empty(cdts[0]):
1688                return [(cdt.C + cdt.T / 2) / (cdt.C + cdt.D + cdt.T)
1689                        for cdt in cdts], False
1690
1691        return False, False
1692
1693class AUC_for_single_class(AUC):
1694    """
1695    Compute AUC where the class with the given class_index is singled
1696    out and all other classes are treated as a single class.
1697    """
1698    def __init__(self, test_results=None, class_index=-1, ignore_weights=False):
1699        if class_index < 0:
1700            if test_results and test_results.base_class >= 0:
1701                self.class_index = test_results.base_class
1702            else:
1703                self.class_index = 1
1704        else:
1705            self.class_index = class_index
1706
1707        super(AUC_for_single_class, self).__init__(test_results, ignore_weights=ignore_weights)
1708
1709    def __call__(self, test_results):
1710        return self._compute_for_single_class(test_results, self.class_index)
1711
1712
1713class AUC_for_pair_of_classes(AUC):
1714    """
1715    Computes AUC between a pair of classes, ignoring instances from all
1716    other classes.
1717    """
1718    def __init__(self, test_results, class_index1, class_index2, ignore_weights=False):
1719        self.class_index1 = class_index1
1720        self.class_index2 = class_index2
1721
1722        super(AUC_for_pair_of_classes, self).__init__(test_results, ignore_weights=ignore_weights)
1723
1724    def __call__(self, test_results):
1725        return self._compute_for_pair_of_classes(test_results, self.class_index1, self.class_index2)
1726
1727
1728class AUC_matrix(AUC):
1729    """
1730    Compute a (lower diagonal) matrix with AUCs for all pairs of classes.
1731    If there are empty classes, the corresponding elements in the matrix
1732    are -1.
1733    """
1734
1735    def __call__(self, test_results):
1736        numberOfClasses = len(test_results.class_values)
1737        number_of_learners = test_results.number_of_learners
1738        if test_results.number_of_iterations > 1:
1739            iterations, all_ite = split_by_iterations(test_results), test_results
1740        else:
1741            iterations, all_ite = [test_results], None
1742        aucs = [[[] for _ in range(numberOfClasses)]
1743        for _ in range(number_of_learners)]
1744        for classIndex1 in range(numberOfClasses):
1745            for classIndex2 in range(classIndex1):
1746                pair_aucs = self._compute_for_multiple_folds(
1747                    self._compute_one_class_against_another, iterations,
1748                    (classIndex1, classIndex2, all_ite,
1749                     test_results.number_of_iterations))
1750                if pair_aucs:
1751                    for lrn in range(number_of_learners):
1752                        aucs[lrn][classIndex1].append(pair_aucs[lrn])
1753                else:
1754                    for lrn in range(number_of_learners):
1755                        aucs[lrn][classIndex1].append(-1)
1756        return aucs
1757
1758#Backward compatibility
1759@replace_use_weights
1760def AUC_binary(res, ignore_weights=False):
1761    auc = deprecated_function_name(AUC)(ignore_weights=ignore_weights)
1762    auc._compute_for_binary_class(res)
1763    return auc
1764
1765@replace_use_weights
1766def AUC_multi(res, ignore_weights=False, method=0):
1767    auc = deprecated_function_name(AUC)(ignore_weights=ignore_weights,
1768        method=method)
1769    auc._compute_for_multi_value_class(res)
1770    return auc
1771
1772
1773@deprecated_keywords({"AUCcomputer": "auc_computer",
1774                      "computerArgs": "computer_args"})
1775def AUC_iterations(auc_computer, iterations, computer_args):
1776    auc = deprecated_function_name(AUC)()
1777    auc._compute_for_multiple_folds(auc_computer, iterations, computer_args)
1778    return auc
1779
1780def AUC_x(cdtComputer, ite, all_ite, divide_by_if_ite, computer_args):
1781    auc = deprecated_function_name(AUC)()
1782    result = auc._compute_auc(cdtComputer, ite, all_ite, divide_by_if_ite,
1783                              computer_args)
1784    return result
1785
1786@replace_use_weights
1787def AUC_i(ite, class_index, ignore_weights=False, all_ite=None,
1788          divide_by_if_ite=1.):
1789    auc = deprecated_function_name(AUC)(ignore_weights=ignore_weights)
1790    result = auc._compute_one_class_against_another(ite, class_index,
1791        all_ite=all_ite, divide_by_if_ite=divide_by_if_ite)
1792    return result
1793
1794
1795@replace_use_weights
1796def AUC_ij(ite, class_index1, class_index2, ignore_weights=False,
1797           all_ite=None, divide_by_if_ite=1.):
1798    auc = deprecated_function_name(AUC)(ignore_weights=ignore_weights)
1799    result = auc._compute_one_class_against_another(
1800        ite, class_index1, class_index2, all_ite=all_ite, divide_by_if_ite=divide_by_if_ite)
1801    return result
1802
1803AUC_single = replace_use_weights(
1804             deprecated_keywords({"classIndex": "class_index"})(
1805             deprecated_function_name(AUC_for_single_class)))
1806AUC_pair = replace_use_weights(
1807           deprecated_keywords({"classIndex1": "class_index1",
1808                                "classIndex2": "class_index2"})(
1809           deprecated_function_name(AUC_for_pair_of_classes)))
1810AUC_matrix = replace_use_weights(deprecated_function_name(AUC_matrix))
1811
1812
1813@deprecated_keywords({"unweighted": "ignore_weights"})
1814def McNemar(res, ignore_weights=False, **argkw):
1815    """
1816    Compute a triangular matrix with McNemar statistics for each pair of
1817    classifiers. The statistics is distributed by chi-square distribution with
1818    one degree of freedom; critical value for 5% significance is around 3.84.
1819    """
1820    nLearners = res.number_of_learners
1821    mcm = []
1822    for i in range(nLearners):
1823       mcm.append([0.] * res.number_of_learners)
1824
1825    if not res.weights or ignore_weights:
1826        for i in res.results:
1827            actual = i.actual_class
1828            classes = i.classes
1829            for l1 in range(nLearners):
1830                for l2 in range(l1, nLearners):
1831                    if classes[l1] == actual:
1832                        if classes[l2] != actual:
1833                            mcm[l1][l2] += 1
1834                    elif classes[l2] == actual:
1835                        mcm[l2][l1] += 1
1836    else:
1837        for i in res.results:
1838            actual = i.actual_class
1839            classes = i.classes
1840            for l1 in range(nLearners):
1841                for l2 in range(l1, nLearners):
1842                    if classes[l1] == actual:
1843                        if classes[l2] != actual:
1844                            mcm[l1][l2] += i.weight
1845                    elif classes[l2] == actual:
1846                        mcm[l2][l1] += i.weight
1847
1848    for l1 in range(nLearners):
1849        for l2 in range(l1, nLearners):
1850            su = mcm[l1][l2] + mcm[l2][l1]
1851            if su:
1852                mcm[l2][l1] = (abs(mcm[l1][l2] - mcm[l2][l1]) - 1)**2 / su
1853            else:
1854                mcm[l2][l1] = 0
1855
1856    for l1 in range(nLearners):
1857        mcm[l1]=mcm[l1][:l1]
1858
1859    return mcm
1860
1861
1862def McNemar_of_two(res, lrn1, lrn2, ignore_weights=False):
1863    """
1864    McNemar_of_two computes a McNemar statistics for a pair of classifier,
1865    specified by indices learner1 and learner2.
1866    """
1867    tf = ft = 0.
1868    if not res.weights or ignore_weights:
1869        for i in res.results:
1870            actual = i.actual_class
1871            if i.classes[lrn1] == actual:
1872                if i.classes[lrn2] != actual:
1873                    tf += i.weight
1874            elif i.classes[lrn2] == actual:
1875                    ft += i.weight
1876    else:
1877        for i in res.results:
1878            actual = i.actual_class
1879            if i.classes[lrn1] == actual:
1880                if i.classes[lrn2] != actual:
1881                    tf += 1.
1882            elif i.classes[lrn2] == actual:
1883                    ft += 1.
1884
1885    su = tf + ft
1886    if su:
1887        return (abs(tf - ft) - 1)**2 / su
1888    else:
1889        return 0
1890
1891
1892def Friedman(res, stat=CA):
1893    """
1894    Compare classifiers with Friedman test, treating folds as different examles.
1895    Returns F, p and average ranks.
1896    """
1897    res_split = split_by_iterations(res)
1898    res = [stat(r) for r in res_split]
1899   
1900    N = len(res)
1901    k = len(res[0])
1902    sums = [0.] * k
1903    for r in res:
1904        ranks = [k - x + 1 for x in statc.rankdata(r)]
1905        if stat == Brier_score: # reverse ranks for Brier_score (lower better)
1906            ranks = [k + 1 - x for x in ranks]
1907        sums = map(add, ranks, sums)
1908
1909    T = sum(x * x for x in sums)
1910    sums = [x / N for x in sums]
1911
1912    F = 12. / (N * k * (k + 1)) * T  - 3 * N * (k + 1)
1913
1914    return F, statc.chisqprob(F, k - 1), sums
1915
1916
1917def Wilcoxon_pairs(res, avgranks, stat=CA):
1918    """
1919    Return a triangular matrix, where element[i][j] stores significance of
1920    difference between the i-th and the j-th classifier, as computed by the
1921    Wilcoxon test. The element is positive if the i-th is better than the j-th,
1922    negative if it is worse, and 1 if they are equal.
1923    Arguments are ExperimentResults, average ranks (as returned by Friedman)
1924    and, optionally, a statistics; greater values should mean better results.
1925    """
1926    res_split = split_by_iterations(res)
1927    res = [stat(r) for r in res_split]
1928
1929    k = len(res[0])
1930    bt = []
1931    for m1 in range(k):
1932        nl = []
1933        for m2 in range(m1 + 1, k):
1934            t, p = statc.wilcoxont([r[m1] for r in res], [r[m2] for r in res])
1935            if avgranks[m1]<avgranks[m2]:
1936                nl.append(p)
1937            elif avgranks[m2]<avgranks[m1]:
1938                nl.append(-p)
1939            else:
1940                nl.append(1)
1941        bt.append(nl)
1942    return bt
1943
1944
1945@deprecated_keywords({"allResults": "all_results",
1946                      "noConfidence": "no_confidence"})
1947def plot_learning_curve_learners(file, all_results, proportions, learners,
1948                                 no_confidence=0):
1949    plot_learning_curve(file, all_results, proportions,
1950        [Orange.misc.getobjectname(learners[i], "Learner %i" % i)
1951        for i in range(len(learners))], no_confidence)
1952
1953
1954@deprecated_keywords({"allResults": "all_results",
1955                      "noConfidence": "no_confidence"})
1956def plot_learning_curve(file, all_results, proportions, legend,
1957                        no_confidence=0):
1958    import types
1959    fopened = 0
1960    if type(file) == types.StringType:
1961        file = open(file, "wt")
1962        fopened = 1
1963       
1964    file.write("set yrange [0:1]\n")
1965    file.write("set xrange [%f:%f]\n" % (proportions[0], proportions[-1]))
1966    file.write("set multiplot\n\n")
1967    CAs = [CA(x, report_se=True) for x in all_results]
1968
1969    file.write("plot \\\n")
1970    for i in range(len(legend) - 1):
1971        if not no_confidence:
1972            file.write("'-' title '' with yerrorbars pointtype %i,\\\n" % (i+1))
1973        file.write("'-' title '%s' with linespoints pointtype %i,\\\n" % (legend[i], i+1))
1974    if not no_confidence:
1975        file.write("'-' title '' with yerrorbars pointtype %i,\\\n" % (len(legend)))
1976    file.write("'-' title '%s' with linespoints pointtype %i\n" % (legend[-1], len(legend)))
1977
1978    for i in range(len(legend)):
1979        if not no_confidence:
1980            for p in range(len(proportions)):
1981                file.write("%f\t%f\t%f\n" % (proportions[p], CAs[p][i][0], 1.96*CAs[p][i][1]))
1982            file.write("e\n\n")
1983
1984        for p in range(len(proportions)):
1985            file.write("%f\t%f\n" % (proportions[p], CAs[p][i][0]))
1986        file.write("e\n\n")
1987
1988    if fopened:
1989        file.close()
1990
1991
1992def print_single_ROC_curve_coordinates(file, curve):
1993    import types
1994    fopened=0
1995    if type(file)==types.StringType:
1996        file=open(file, "wt")
1997        fopened=1
1998
1999    for coord in curve:
2000        file.write("%5.3f\t%5.3f\n" % tuple(coord))
2001
2002    if fopened:
2003        file.close()
2004
2005
2006def plot_ROC_learners(file, curves, learners):
2007    plot_ROC(file, curves, [Orange.misc.getobjectname(learners[i], "Learner %i" % i) for i in range(len(learners))])
2008   
2009def plot_ROC(file, curves, legend):
2010    import types
2011    fopened=0
2012    if type(file)==types.StringType:
2013        file=open(file, "wt")
2014        fopened=1
2015
2016    file.write("set yrange [0:1]\n")
2017    file.write("set xrange [0:1]\n")
2018    file.write("set multiplot\n\n")
2019
2020    file.write("plot \\\n")
2021    for leg in legend:
2022        file.write("'-' title '%s' with lines,\\\n" % leg)
2023    file.write("'-' title '' with lines\n")
2024
2025    for curve in curves:
2026        for coord in curve:
2027            file.write("%5.3f\t%5.3f\n" % tuple(coord))
2028        file.write("e\n\n")
2029
2030    file.write("1.0\t1.0\n0.0\t0.0e\n\n")         
2031
2032    if fopened:
2033        file.close()
2034
2035
2036@deprecated_keywords({"allResults": "all_results"})
2037def plot_McNemar_curve_learners(file, all_results, proportions, learners, reference=-1):
2038    plot_McNemar_curve(file, all_results, proportions, [Orange.misc.getobjectname(learners[i], "Learner %i" % i) for i in range(len(learners))], reference)
2039
2040
2041@deprecated_keywords({"allResults": "all_results"})
2042def plot_McNemar_curve(file, all_results, proportions, legend, reference=-1):
2043    if reference<0:
2044        reference=len(legend)-1
2045       
2046    import types
2047    fopened=0
2048    if type(file)==types.StringType:
2049        file=open(file, "wt")
2050        fopened=1
2051       
2052    #file.write("set yrange [0:1]\n")
2053    #file.write("set xrange [%f:%f]\n" % (proportions[0], proportions[-1]))
2054    file.write("set multiplot\n\n")
2055    file.write("plot \\\n")
2056    tmap=range(reference)+range(reference+1, len(legend))
2057    for i in tmap[:-1]:
2058        file.write("'-' title '%s' with linespoints pointtype %i,\\\n" % (legend[i], i+1))
2059    file.write("'-' title '%s' with linespoints pointtype %i\n" % (legend[tmap[-1]], tmap[-1]))
2060    file.write("\n")
2061
2062    for i in tmap:
2063        for p in range(len(proportions)):
2064            file.write("%f\t%f\n" % (proportions[p], McNemar_of_two(all_results[p], i, reference)))
2065        file.write("e\n\n")
2066
2067    if fopened:
2068        file.close()
2069
2070default_point_types=("{$\\circ$}", "{$\\diamond$}", "{$+$}", "{$\\times$}", "{$|$}")+tuple([chr(x) for x in range(97, 122)])
2071default_line_types=("\\setsolid", "\\setdashpattern <4pt, 2pt>", "\\setdashpattern <8pt, 2pt>", "\\setdashes", "\\setdots")
2072
2073@deprecated_keywords({"allResults": "all_results"})
2074def learning_curve_learners_to_PiCTeX(file, all_results, proportions, **options):
2075    return apply(learning_curve_to_PiCTeX, (file, all_results, proportions), options)
2076
2077
2078@deprecated_keywords({"allResults": "all_results"})
2079def learning_curve_to_PiCTeX(file, all_results, proportions, **options):
2080    import types
2081    fopened=0
2082    if type(file)==types.StringType:
2083        file=open(file, "wt")
2084        fopened=1
2085
2086    nexamples=len(all_results[0].results)
2087    CAs = [CA(x, report_se=True) for x in all_results]
2088
2089    graphsize=float(options.get("graphsize", 10.0)) #cm
2090    difprop=proportions[-1]-proportions[0]
2091    ntestexamples=nexamples*proportions[-1]
2092    xunit=graphsize/ntestexamples
2093
2094    yshift=float(options.get("yshift", -ntestexamples/20.))
2095   
2096    pointtypes=options.get("pointtypes", default_point_types)
2097    linetypes=options.get("linetypes", default_line_types)
2098
2099    if options.has_key("numberedx"):
2100        numberedx=options["numberedx"]
2101        if type(numberedx)==types.IntType:
2102            if numberedx>0:
2103                numberedx=[nexamples*proportions[int(i/float(numberedx)*len(proportions))] for i in range(numberedx)]+[proportions[-1]*nexamples]
2104            elif numberedx<0:
2105                numberedx = -numberedx
2106                newn=[]
2107                for i in range(numberedx+1):
2108                    wanted=proportions[0]+float(i)/numberedx*difprop
2109                    best=(10, 0)
2110                    for t in proportions:
2111                        td=abs(wanted-t)
2112                        if td<best[0]:
2113                            best=(td, t)
2114                    if not best[1] in newn:
2115                        newn.append(best[1])
2116                newn.sort()
2117                numberedx=[nexamples*x for x in newn]
2118        elif type(numberedx[0])==types.FloatType:
2119            numberedx=[nexamples*x for x in numberedx]
2120    else:
2121        numberedx=[nexamples*x for x in proportions]
2122
2123    file.write("\\mbox{\n")
2124    file.write(\\beginpicture\n")
2125    file.write(\\setcoordinatesystem units <%10.8fcm, %5.3fcm>\n\n" % (xunit, graphsize))   
2126    file.write(\\setplotarea x from %5.3f to %5.3f, y from 0 to 1\n" % (0, ntestexamples))   
2127    file.write(\\axis bottom invisible\n")# label {#examples}\n")
2128    file.write("      ticks short at %s /\n" % reduce(lambda x,y:x+" "+y, ["%i"%(x*nexamples+0.5) for x in proportions]))
2129    if numberedx:
2130        file.write("            long numbered at %s /\n" % reduce(lambda x,y:x+y, ["%i " % int(x+0.5) for x in numberedx]))
2131    file.write("  /\n")
2132    file.write(\\axis left invisible\n")# label {classification accuracy}\n")
2133    file.write("      shiftedto y=%5.3f\n" % yshift)
2134    file.write("      ticks short from 0.0 to 1.0 by 0.05\n")
2135    file.write("            long numbered from 0.0 to 1.0 by 0.25\n")
2136    file.write("  /\n")
2137    if options.has_key("default"):
2138        file.write(\\setdashpattern<1pt, 1pt>\n")
2139        file.write(\\plot %5.3f %5.3f %5.3f %5.3f /\n" % (0., options["default"], ntestexamples, options["default"]))
2140   
2141    for i in range(len(CAs[0])):
2142        coordinates=reduce(lambda x,y:x+" "+y, ["%i %5.3f" % (proportions[p]*nexamples, CAs[p][i][0]) for p in range(len(proportions))])
2143        if linetypes:
2144            file.write(%s\n" % linetypes[i])
2145            file.write(\\plot %s /\n" % coordinates)
2146        if pointtypes:
2147            file.write(\\multiput %s at %s /\n" % (pointtypes[i], coordinates))
2148
2149    file.write(\\endpicture\n")
2150    file.write("}\n")
2151    if fopened:
2152        file.close()
2153    file.close()
2154    del file
2155
2156def legend_learners_to_PiCTeX(file, learners, **options):
2157  return apply(legend_to_PiCTeX, (file, [Orange.misc.getobjectname(learners[i], "Learner %i" % i) for i in range(len(learners))]), options)
2158   
2159def legend_to_PiCTeX(file, legend, **options):
2160    import types
2161    fopened=0
2162    if type(file)==types.StringType:
2163        file=open(file, "wt")
2164        fopened=1
2165
2166    pointtypes=options.get("pointtypes", default_point_types)
2167    linetypes=options.get("linetypes", default_line_types)
2168
2169    file.write("\\mbox{\n")
2170    file.write(\\beginpicture\n")
2171    file.write(\\setcoordinatesystem units <5cm, 1pt>\n\n")
2172    file.write(\\setplotarea x from 0.000 to %5.3f, y from 0 to 12\n" % len(legend))
2173
2174    for i in range(len(legend)):
2175        if linetypes:
2176            file.write(%s\n" % linetypes[i])
2177            file.write(\\plot %5.3f 6 %5.3f 6 /\n" % (i, i+0.2))
2178        if pointtypes:
2179            file.write(\\put {%s} at %5.3f 6\n" % (pointtypes[i], i+0.1))
2180        file.write(\\put {%s} [lb] at %5.3f 0\n" % (legend[i], i+0.25))
2181
2182    file.write(\\endpicture\n")
2183    file.write("}\n")
2184    if fopened:
2185        file.close()
2186    file.close()
2187    del file
2188
2189
2190def compute_friedman(avranks, N):
2191    """ Returns a tuple composed of (friedman statistic, degrees of freedom)
2192    and (Iman statistic - F-distribution, degrees of freedoma) given average
2193    ranks and a number of tested data sets N.
2194    """
2195
2196    k = len(avranks)
2197
2198    def friedman(N, k, ranks):
2199        return 12*N*(sum([rank**2.0 for rank in ranks]) - (k*(k+1)*(k+1)/4.0) )/(k*(k+1))
2200
2201    def iman(fried, N, k):
2202        return (N-1)*fried/(N*(k-1) - fried)
2203
2204    f = friedman(N, k, avranks)
2205    im = iman(f, N, k)
2206    fdistdof = (k-1, (k-1)*(N-1))
2207
2208    return (f, k-1), (im, fdistdof)
2209
2210def compute_CD(avranks, N, alpha="0.05", type="nemenyi"):
2211    """ Returns critical difference for Nemenyi or Bonferroni-Dunn test
2212    according to given alpha (either alpha="0.05" or alpha="0.1") for average
2213    ranks and number of tested data sets N. Type can be either "nemenyi" for
2214    for Nemenyi two tailed test or "bonferroni-dunn" for Bonferroni-Dunn test.
2215    """
2216
2217    k = len(avranks)
2218   
2219    d = {("nemenyi", "0.05"): [0, 0, 1.959964, 2.343701, 2.569032, 2.727774,
2220                               2.849705, 2.94832, 3.030879, 3.101730, 3.163684,
2221                               3.218654, 3.268004, 3.312739, 3.353618, 3.39123,
2222                               3.426041, 3.458425, 3.488685, 3.517073, 3.543799]
2223        , ("nemenyi", "0.1"): [0, 0, 1.644854, 2.052293, 2.291341, 2.459516,
2224                               2.588521, 2.692732, 2.779884, 2.854606, 2.919889,
2225                               2.977768, 3.029694, 3.076733, 3.119693, 3.159199,
2226                               3.195743, 3.229723, 3.261461, 3.291224, 3.319233]
2227        , ("bonferroni-dunn", "0.05"): [0, 0, 1.960, 2.241, 2.394, 2.498, 2.576,
2228                                        2.638, 2.690, 2.724, 2.773],
2229         ("bonferroni-dunn", "0.1"): [0, 0, 1.645, 1.960, 2.128, 2.241, 2.326,
2230                                      2.394, 2.450, 2.498, 2.539]}
2231
2232    #can be computed in R as qtukey(0.95, n, Inf)**0.5
2233    #for (x in c(2:20)) print(qtukey(0.95, x, Inf)/(2**0.5)
2234
2235    q = d[(type, alpha)]
2236
2237    cd = q[k]*(k*(k+1)/(6.0*N))**0.5
2238
2239    return cd
2240 
2241
2242def graph_ranks(filename, avranks, names, cd=None, cdmethod=None, lowv=None, highv=None, width=6, textspace=1, reverse=False, **kwargs):
2243    """
2244    Draws a CD graph, which is used to display  the differences in methods'
2245    performance.
2246    See Janez Demsar, Statistical Comparisons of Classifiers over
2247    Multiple Data Sets, 7(Jan):1--30, 2006.
2248
2249    Needs matplotlib to work.
2250
2251    :param filename: Output file name (with extension). Formats supported
2252                     by matplotlib can be used.
2253    :param avranks: List of average methods' ranks.
2254    :param names: List of methods' names.
2255
2256    :param cd: Critical difference. Used for marking methods that whose
2257               difference is not statistically significant.
2258    :param lowv: The lowest shown rank, if None, use 1.
2259    :param highv: The highest shown rank, if None, use len(avranks).
2260    :param width: Width of the drawn figure in inches, default 6 in.
2261    :param textspace: Space on figure sides left for the description
2262                      of methods, default 1 in.
2263    :param reverse:  If True, the lowest rank is on the right. Default\: False.
2264    :param cdmethod: None by default. It can be an index of element in avranks
2265                     or or names which specifies the method which should be
2266                     marked with an interval.
2267    """
2268    if not HAS_MATPLOTLIB:
2269        import sys
2270        print >> sys.stderr, "Function requires matplotlib. Please install it."
2271        return
2272
2273    width = float(width)
2274    textspace = float(textspace)
2275
2276    def nth(l,n):
2277        """
2278        Returns only nth elemnt in a list.
2279        """
2280        n = lloc(l,n)
2281        return [ a[n] for a in l ]
2282
2283    def lloc(l,n):
2284        """
2285        List location in list of list structure.
2286        Enable the use of negative locations:
2287        -1 is the last element, -2 second last...
2288        """
2289        if n < 0:
2290            return len(l[0])+n
2291        else:
2292            return n
2293
2294    def mxrange(lr):
2295        """
2296        Multiple xranges. Can be used to traverse matrices.
2297        This function is very slow due to unknown number of
2298        parameters.
2299
2300        >>> mxrange([3,5])
2301        [(0, 0), (0, 1), (0, 2), (1, 0), (1, 1), (1, 2)]
2302
2303        >>> mxrange([[3,5,1],[9,0,-3]])
2304        [(3, 9), (3, 6), (3, 3), (4, 9), (4, 6), (4, 3)]
2305
2306        """
2307        if not len(lr):
2308            yield ()
2309        else:
2310            #it can work with single numbers
2311            index = lr[0]
2312            if type(1) == type(index):
2313                index = [ index ]
2314            for a in range(*index):
2315                for b in mxrange(lr[1:]):
2316                    yield tuple([a] + list(b))
2317
2318
2319    from matplotlib.figure import Figure
2320    from matplotlib.backends.backend_agg import FigureCanvasAgg
2321
2322    def print_figure(fig, *args, **kwargs):
2323        canvas = FigureCanvasAgg(fig)
2324        canvas.print_figure(*args, **kwargs)
2325
2326    sums = avranks
2327
2328    tempsort =  sorted([ (a,i) for i,a in  enumerate(sums) ], reverse=reverse)
2329    ssums = nth(tempsort, 0)
2330    sortidx = nth(tempsort, 1)
2331    nnames = [ names[x] for x in sortidx ]
2332   
2333    if lowv is None:
2334        lowv = min(1, int(math.floor(min(ssums))))
2335    if highv is None:
2336        highv = max(len(avranks), int(math.ceil(max(ssums))))
2337
2338    cline = 0.4
2339
2340    k = len(sums)
2341
2342    lines = None
2343
2344    linesblank = 0
2345    scalewidth = width - 2*textspace
2346
2347    def rankpos(rank):
2348        if not reverse:
2349            a = rank - lowv
2350        else:
2351            a = highv - rank
2352        return textspace+scalewidth/(highv-lowv)*a
2353
2354    distanceh = 0.25
2355
2356    if cd and cdmethod is None:
2357   
2358        #get pairs of non significant methods
2359
2360        def get_lines(sums, hsd):
2361
2362            #get all pairs
2363            lsums = len(sums)
2364            allpairs = [ (i,j) for i,j in mxrange([[lsums], [lsums]]) if j > i ]
2365            #remove not significant
2366            notSig = [ (i,j) for i,j in allpairs if abs(sums[i]-sums[j]) <= hsd ]
2367            #keep only longest
2368           
2369            def no_longer((i,j), notSig):
2370                for i1,j1 in notSig:
2371                    if (i1 <= i and j1 > j) or (i1 < i and j1 >= j):
2372                        return False
2373                return True
2374
2375            longest = [ (i,j) for i,j in notSig if no_longer((i,j),notSig) ]
2376           
2377            return longest
2378
2379        lines = get_lines(ssums, cd)
2380        linesblank = 0.2 + 0.2 + (len(lines)-1)*0.1
2381
2382        #add scale
2383        distanceh = 0.25
2384        cline += distanceh
2385
2386    #calculate height needed height of an image
2387    minnotsignificant = max(2*0.2, linesblank)
2388    height = cline + ((k+1)/2)*0.2 + minnotsignificant
2389
2390    fig = Figure(figsize=(width, height))
2391    ax = fig.add_axes([0,0,1,1]) #reverse y axis
2392    ax.set_axis_off()
2393
2394    hf = 1./height # height factor
2395    wf = 1./width
2396
2397    def hfl(l): 
2398        return [ a*hf for a in l ]
2399
2400    def wfl(l): 
2401        return [ a*wf for a in l ]
2402
2403
2404    # Upper left corner is (0,0).
2405
2406    ax.plot([0,1], [0,1], c="w")
2407    ax.set_xlim(0, 1)
2408    ax.set_ylim(1, 0)
2409
2410    def line(l, color='k', **kwargs):
2411        """
2412        Input is a list of pairs of points.
2413        """
2414        ax.plot(wfl(nth(l,0)), hfl(nth(l,1)), color=color, **kwargs)
2415
2416    def text(x, y, s, *args, **kwargs):
2417        ax.text(wf*x, hf*y, s, *args, **kwargs)
2418
2419    line([(textspace, cline), (width-textspace, cline)], linewidth=0.7)
2420   
2421    bigtick = 0.1
2422    smalltick = 0.05
2423
2424
2425    import numpy
2426    tick = None
2427    for a in list(numpy.arange(lowv, highv, 0.5)) + [highv]:
2428        tick = smalltick
2429        if a == int(a): tick = bigtick
2430        line([(rankpos(a), cline-tick/2),(rankpos(a), cline)], linewidth=0.7)
2431
2432    for a in range(lowv, highv+1):
2433        text(rankpos(a), cline-tick/2-0.05, str(a), ha="center", va="bottom")
2434
2435    k = len(ssums)
2436
2437    for i in range((k+1)/2):
2438        chei = cline+ minnotsignificant + i *0.2
2439        line([(rankpos(ssums[i]), cline), (rankpos(ssums[i]), chei), (textspace-0.1, chei)], linewidth=0.7)
2440        text(textspace-0.2, chei, nnames[i], ha="right", va="center")
2441
2442    for i in range((k+1)/2, k):
2443        chei = cline + minnotsignificant + (k-i-1)*0.2
2444        line([(rankpos(ssums[i]), cline), (rankpos(ssums[i]), chei), (textspace+scalewidth+0.1, chei)], linewidth=0.7)
2445        text(textspace+scalewidth+0.2, chei, nnames[i], ha="left", va="center")
2446
2447    if cd and cdmethod is None:
2448
2449        #upper scale
2450        if not reverse:
2451            begin, end = rankpos(lowv), rankpos(lowv+cd)
2452        else:
2453            begin, end = rankpos(highv), rankpos(highv - cd)
2454           
2455        line([(begin, distanceh), (end, distanceh)], linewidth=0.7)
2456        line([(begin, distanceh + bigtick/2), (begin, distanceh - bigtick/2)], linewidth=0.7)
2457        line([(end, distanceh + bigtick/2), (end, distanceh - bigtick/2)], linewidth=0.7)
2458        text((begin+end)/2, distanceh - 0.05, "CD", ha="center", va="bottom")
2459
2460        #non significance lines   
2461        def draw_lines(lines, side=0.05, height=0.1):
2462            start = cline + 0.2
2463            for l,r in lines: 
2464                line([(rankpos(ssums[l])-side, start), (rankpos(ssums[r])+side, start)], linewidth=2.5) 
2465                start += height
2466
2467        draw_lines(lines)
2468
2469    elif cd:
2470        begin = rankpos(avranks[cdmethod]-cd)
2471        end = rankpos(avranks[cdmethod]+cd)
2472        line([(begin, cline), (end, cline)], linewidth=2.5) 
2473        line([(begin, cline + bigtick/2), (begin, cline - bigtick/2)], linewidth=2.5)
2474        line([(end, cline + bigtick/2), (end, cline - bigtick/2)], linewidth=2.5)
2475 
2476    print_figure(fig, filename, **kwargs)
2477
2478def mlc_hamming_loss(res):
2479    """
2480    Schapire and Singer (2000) presented Hamming Loss, which id defined as:
2481   
2482    :math:`HammingLoss(H,D)=\\frac{1}{|D|} \\sum_{i=1}^{|D|} \\frac{Y_i \\vartriangle Z_i}{|L|}`
2483    """
2484    losses = [0]*res.number_of_learners
2485    label_num = len(res.labels)
2486    example_num = gettotsize(res)
2487
2488    for e in res.results:
2489        aclass = e.actual_class
2490        for i, labels in enumerate(e.classes):
2491            labels = map(int, labels)
2492            if len(labels) <> len(aclass):
2493                raise ValueError, "The dimensions of the classified output and the actual class array do not match."
2494            for j in range(label_num):
2495                if labels[j] != aclass[j]:
2496                    losses[i] += 1
2497           
2498    return [float(x)/(label_num*example_num) for x in losses]
2499
2500def mlc_accuracy(res, forgiveness_rate = 1.0):
2501    """
2502    Godbole & Sarawagi, 2004 uses the metrics accuracy, precision, recall as follows:
2503     
2504    :math:`Accuracy(H,D)=\\frac{1}{|D|} \\sum_{i=1}^{|D|} \\frac{|Y_i \\cap Z_i|}{|Y_i \\cup Z_i|}`
2505   
2506    Boutell et al. (2004) give a more generalized version using a parameter :math:`\\alpha \\ge 0`,
2507    called forgiveness rate:
2508   
2509    :math:`Accuracy(H,D)=\\frac{1}{|D|} \\sum_{i=1}^{|D|} (\\frac{|Y_i \\cap Z_i|}{|Y_i \\cup Z_i|})^{\\alpha}`
2510    """
2511    accuracies = [0.0]*res.number_of_learners
2512    example_num = gettotsize(res)
2513   
2514    for e in res.results:
2515        aclass = e.actual_class
2516        for i, labels in enumerate(e.classes):
2517            labels = map(int, labels)
2518            if len(labels) <> len(aclass):
2519                raise ValueError, "The dimensions of the classified output and the actual class array do not match."
2520           
2521            intersection = 0.0
2522            union = 0.0
2523            for real, pred in zip(labels, aclass):
2524                if real and pred:
2525                    intersection += 1
2526                if real or pred:
2527                    union += 1
2528
2529            if union:
2530                accuracies[i] += intersection / union
2531           
2532    return [math.pow(x/example_num,forgiveness_rate) for x in accuracies]
2533
2534def mlc_precision(res):
2535    """
2536    :math:`Precision(H,D)=\\frac{1}{|D|} \\sum_{i=1}^{|D|} \\frac{|Y_i \\cap Z_i|}{|Z_i|}`
2537    """
2538    precisions = [0.0]*res.number_of_learners
2539    example_num = gettotsize(res)
2540   
2541    for e in res.results:
2542        aclass = e.actual_class
2543        for i, labels in enumerate(e.classes):
2544            labels = map(int, labels)
2545            if len(labels) <> len(aclass):
2546                raise ValueError, "The dimensions of the classified output and the actual class array do not match."
2547           
2548            intersection = 0.0
2549            predicted = 0.0
2550            for real, pred in zip(labels, aclass):
2551                if real and pred:
2552                    intersection += 1
2553                if real:
2554                    predicted += 1
2555            if predicted:
2556                precisions[i] += intersection / predicted
2557           
2558    return [x/example_num for x in precisions]
2559
2560def mlc_recall(res):
2561    """
2562    :math:`Recall(H,D)=\\frac{1}{|D|} \\sum_{i=1}^{|D|} \\frac{|Y_i \\cap Z_i|}{|Y_i|}`
2563    """
2564    recalls = [0.0]*res.number_of_learners
2565    example_num = gettotsize(res)
2566   
2567    for e in res.results:
2568        aclass = e.actual_class
2569        for i, labels in enumerate(e.classes):
2570            labels = map(int, labels)
2571            if len(labels) <> len(aclass):
2572                raise ValueError, "The dimensions of the classified output and the actual class array do not match."
2573           
2574            intersection = 0.0
2575            actual = 0.0
2576            for real, pred in zip(labels, aclass):
2577                if real and pred:
2578                    intersection += 1
2579                if pred:
2580                    actual += 1
2581            if actual:
2582                recalls[i] += intersection / actual
2583           
2584    return [x/example_num for x in recalls]
2585
2586#def mlc_ranking_loss(res):
2587#    pass
2588#
2589#def mlc_average_precision(res):
2590#    pass
2591#
2592#def mlc_hierarchical_loss(res):
2593#    pass
2594
2595
2596def mt_average_score(res, score, weights=None):
2597    """
2598    Compute individual scores for each target and return the (weighted) average.
2599
2600    One method can be used to compute scores for all targets or a list of
2601    scoring methods can be passed to use different methods for different
2602    targets. In the latter case, care has to be taken if the ranges of scoring
2603    methods differ.
2604    For example, when the first target is scored from -1 to 1 (1 best) and the
2605    second from 0 to 1 (0 best), using `weights=[0.5,-1]` would scale both
2606    to a span of 1, and invert the second so that higher scores are better.
2607
2608    :param score: Single-target scoring method or a list of such methods
2609                  (one for each target).
2610    :param weights: List of real weights, one for each target,
2611                    for a weighted average.
2612
2613    """
2614    if not len(res.results):
2615        raise ValueError, "Cannot compute the score: no examples."
2616    if res.number_of_learners < 1:
2617        return []
2618    n_classes = len(res.results[0].actual_class)
2619    if weights is None:
2620        weights = [1.] * n_classes
2621    if not hasattr(score, '__len__'):
2622        score = [score] * n_classes
2623    elif len(score) != n_classes:
2624        raise ValueError, "Number of scoring methods and targets do not match."
2625    # save original classes
2626    clsss = [te.classes for te in res.results]
2627    aclsss = [te.actual_class for te in res.results]
2628    # compute single target scores
2629    single_scores = []
2630    for i in range(n_classes):
2631        for te, clss, aclss in zip(res.results, clsss, aclsss):
2632            te.classes = [cls[i] for cls in clss]
2633            te.actual_class = aclss[i]
2634        single_scores.append(score[i](res))
2635    # restore original classes
2636    for te, clss, aclss in zip(res.results, clsss, aclsss):
2637        te.classes = clss
2638        te.actual_class = aclss
2639    return [sum(w * s for w, s in zip(weights, scores)) / sum(weights)
2640        for scores in zip(*single_scores)]
2641
2642def mt_flattened_score(res, score):
2643    """
2644    Flatten (concatenate into a single list) the predictions of multiple
2645    targets and compute a single-target score.
2646   
2647    :param score: Single-target scoring method.
2648    """
2649    res2 = Orange.evaluation.testing.ExperimentResults(res.number_of_iterations,
2650        res.classifier_names, class_values=res.class_values,
2651        weights=res.weights, classifiers=res.classifiers, loaded=res.loaded,
2652        test_type=Orange.evaluation.testing.TEST_TYPE_SINGLE, labels=res.labels)
2653    for te in res.results:
2654        for i, ac in enumerate(te.actual_class):
2655            te2 = Orange.evaluation.testing.TestedExample(
2656                iteration_number=te.iteration_number, actual_class=ac)
2657            for c, p in zip(te.classes, te.probabilities):
2658                te2.add_result(c[i], p[i])
2659            res2.results.append(te2)
2660    return score(res2)
2661
2662
2663################################################################################
2664if __name__ == "__main__":
2665    avranks =  [3.143, 2.000, 2.893, 1.964]
2666    names = ["prva", "druga", "tretja", "cetrta" ]
2667    cd = compute_CD(avranks, 14)
2668    #cd = compute_CD(avranks, 10, type="bonferroni-dunn")
2669    print cd
2670
2671    print compute_friedman(avranks, 14)
2672
2673    #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.