source: orange/Orange/evaluation/scoring.py @ 10362:539f49e3586d

Revision 10362:539f49e3586d, 97.9 KB checked in by Lan Zagar <lan.zagar@…>, 2 years ago (diff)

Small (but numerous) code style improvements in the first third of scoring.py

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