1
2
3
4
5
6
7
8
9 """Utility class to compute the transfer error of classifiers."""
10
11 __docformat__ = 'restructuredtext'
12
13 import mvpa.support.copy as copy
14
15 import numpy as N
16
17 from sets import Set
18 from StringIO import StringIO
19 from math import log10, ceil
20
21 from mvpa.base import externals
22
23 from mvpa.misc.errorfx import meanPowerFx, rootMeanPowerFx, RMSErrorFx, \
24 CorrErrorFx, CorrErrorPFx, RelativeRMSErrorFx, MeanMismatchErrorFx, \
25 AUCErrorFx
26 from mvpa.base import warning
27 from mvpa.misc.state import StateVariable, ClassWithCollections
28 from mvpa.base.dochelpers import enhancedDocString, table2string
29 from mvpa.clfs.stats import autoNullDist
30
31 if __debug__:
32 from mvpa.base import debug
33
34 if externals.exists('scipy'):
35 from scipy.stats.stats import nanmean
36 else:
37 from mvpa.clfs.stats import nanmean
38
39 -def _p2(x, prec=2):
40 """Helper to print depending on the type nicely. For some
41 reason %.2g for 100 prints exponential form which is ugly
42 """
43 if isinstance(x, int):
44 return "%d" % x
45 elif isinstance(x, float):
46 s = ("%%.%df" % prec % x).rstrip('0').rstrip('.').lstrip()
47 if s == '':
48 s = '0'
49 return s
50 else:
51 return "%s" % x
52
56 """Basic class to collect targets/predictions and report summary statistics
57
58 It takes care about collecting the sets, which are just tuples
59 (targets, predictions, values). While 'computing' the matrix, all
60 sets are considered together. Children of the class are
61 responsible for computation and display.
62 """
63
64 _STATS_DESCRIPTION = (
65 ('# of sets',
66 'number of target/prediction sets which were provided',
67 None), )
68
69
70 - def __init__(self, targets=None, predictions=None, values=None, sets=None):
71 """Initialize SummaryStatistics
72
73 targets or predictions cannot be provided alone (ie targets
74 without predictions)
75
76 :Parameters:
77 targets
78 Optional set of targets
79 predictions
80 Optional set of predictions
81 values
82 Optional set of values (which served for prediction)
83 sets
84 Optional list of sets
85 """
86 self._computed = False
87 """Flag either it was computed for a given set of data"""
88
89 self.__sets = (sets, [])[int(sets is None)]
90 """Datasets (target, prediction) to compute confusion matrix on"""
91
92 if not targets is None or not predictions is None:
93 if not targets is None and not predictions is None:
94 self.add(targets=targets, predictions=predictions,
95 values=values)
96 else:
97 raise ValueError, \
98 "Please provide none or both targets and predictions"
99
100
101 - def add(self, targets, predictions, values=None):
102 """Add new results to the set of known results"""
103 if len(targets) != len(predictions):
104 raise ValueError, \
105 "Targets[%d] and predictions[%d]" % (len(targets),
106 len(predictions)) + \
107 " have different number of samples"
108
109 if values is not None and len(targets) != len(values):
110 raise ValueError, \
111 "Targets[%d] and values[%d]" % (len(targets),
112 len(values)) + \
113 " have different number of samples"
114
115
116
117
118 nonetype = type(None)
119 for i in xrange(len(targets)):
120 t1, t2 = type(targets[i]), type(predictions[i])
121
122
123 if t1 != t2 and t2 != nonetype:
124
125
126 if isinstance(predictions, tuple):
127 predictions = list(predictions)
128 predictions[i] = t1(predictions[i])
129
130 self.__sets.append( (targets, predictions, values) )
131 self._computed = False
132
133
134 - def asstring(self, short=False, header=True, summary=True,
135 description=False):
136 """'Pretty print' the matrix
137
138 :Parameters:
139 short : bool
140 if True, ignores the rest of the parameters and provides consise
141 1 line summary
142 header : bool
143 print header of the table
144 summary : bool
145 print summary (accuracy)
146 description : bool
147 print verbose description of presented statistics
148 """
149 raise NotImplementedError
150
151
153 """String summary over the `SummaryStatistics`
154
155 It would print description of the summary statistics if 'CM'
156 debug target is active
157 """
158 if __debug__:
159 description = ('CM' in debug.active)
160 else:
161 description = False
162 return self.asstring(short=False, header=True, summary=True,
163 description=description)
164
165
167 """Add the sets from `other` s `SummaryStatistics` to current one
168 """
169
170
171
172 othersets = copy.copy(other.__sets)
173 for set in othersets:
174 self.add(*set)
175 return self
176
177
179 """Add two `SummaryStatistics`s
180 """
181 result = copy.copy(self)
182 result += other
183 return result
184
185
187 """Actually compute the confusion matrix based on all the sets"""
188 if self._computed:
189 return
190
191 self._compute()
192 self._computed = True
193
194
196 self._stats = {'# of sets' : len(self.sets)}
197
198
199 @property
201 """Return a list of separate summaries per each stored set"""
202 return [ self.__class__(sets=[x]) for x in self.sets ]
203
204
205 @property
207 raise NotImplementedError
208
209
210 @property
212 self.compute()
213 return self._stats
214
215
217 """Cleans summary -- all data/sets are wiped out
218 """
219 self.__sets = []
220 self._computed = False
221
222
223 sets = property(lambda self:self.__sets)
224
227 """Generic class for ROC curve computation and plotting
228 """
229
231 """
232 :Parameters:
233 labels : list
234 labels which were used (in order of values if multiclass,
235 or 1 per class for binary problems (e.g. in SMLR))
236 sets : list of tuples
237 list of sets for the analysis
238 """
239 self._labels = labels
240 self._sets = sets
241 self.__computed = False
242
243
245 """Lazy computation if needed
246 """
247 if self.__computed:
248 return
249
250 labels = self._labels
251 Nlabels = len(labels)
252 sets = self._sets
253
254
255 def _checkValues(set_):
256 """Check if values are 'acceptable'"""
257 if len(set_)<3: return False
258 x = set_[2]
259
260 if (x is None) or len(x) == 0: return False
261 for v in x:
262 try:
263 if Nlabels <= 2 and N.isscalar(v):
264 continue
265 if (isinstance(v, dict) or
266 ((Nlabels>=2) and len(v)!=Nlabels)
267 ): return False
268 except Exception, e:
269
270
271
272 if __debug__:
273 debug('ROC', "Exception %s while checking "
274 "either %s are valid labels" % (str(e), x))
275 return False
276 return True
277
278 sets_wv = filter(_checkValues, sets)
279
280 Nsets_wv = len(sets_wv)
281 if Nsets_wv > 0 and len(sets) != Nsets_wv:
282 warning("Only %d sets have values assigned from %d sets" %
283 (Nsets_wv, len(sets)))
284
285
286
287
288 for iset,s in enumerate(sets_wv):
289
290 values = s[2]
291
292 if isinstance(values, N.ndarray) and len(values.shape)==1:
293 values = list(values)
294 rangev = None
295 for i in xrange(len(values)):
296 v = values[i]
297 if N.isscalar(v):
298 if Nlabels == 2:
299 def last_el(x):
300 """Helper function. Returns x if x is scalar, and
301 last element if x is not (ie list/tuple)"""
302 if N.isscalar(x): return x
303 else: return x[-1]
304 if rangev is None:
305
306
307 values_ = [last_el(x) for x in values]
308 rangev = N.min(values_) + N.max(values_)
309 values[i] = [rangev - v, v]
310 else:
311 raise ValueError, \
312 "Cannot have a single 'value' for multiclass" \
313 " classification. Got %s" % (v)
314 elif len(v) != Nlabels:
315 raise ValueError, \
316 "Got %d values whenever there is %d labels" % \
317 (len(v), Nlabels)
318
319 sets_wv[iset] = (s[0], s[1], N.asarray(values))
320
321
322
323
324
325 ROCs, aucs = [], []
326 for i,label in enumerate(labels):
327 aucs_pl = []
328 ROCs_pl = []
329 for s in sets_wv:
330 targets_pl = (s[0] == label).astype(int)
331
332 ROC = AUCErrorFx()
333 aucs_pl += [ROC([x[i] for x in s[2]], targets_pl)]
334 ROCs_pl.append(ROC)
335 if len(aucs_pl)>0:
336 ROCs += [ROCs_pl]
337 aucs += [nanmean(aucs_pl)]
338
339
340
341 self._ROCs = ROCs
342 self._aucs = aucs
343 self.__computed = True
344
345
346 @property
348 """Compute and return set of AUC values 1 per label
349 """
350 self._compute()
351 return self._aucs
352
353
354 @property
356 self._compute()
357 return self._ROCs
358
359
360 - def plot(self, label_index=0):
361 """
362
363 TODO: make it friendly to labels given by values?
364 should we also treat labels_map?
365 """
366 externals.exists("pylab", raiseException=True)
367 import pylab as P
368
369 self._compute()
370
371 labels = self._labels
372
373 ROCs = self.ROCs[label_index]
374
375 fig = P.gcf()
376 ax = P.gca()
377
378 P.plot([0, 1], [0, 1], 'k:')
379
380 for ROC in ROCs:
381 P.plot(ROC.fp, ROC.tp, linewidth=1)
382
383 P.axis((0.0, 1.0, 0.0, 1.0))
384 P.axis('scaled')
385 P.title('Label %s. Mean AUC=%.2f' % (label_index, self.aucs[label_index]))
386
387 P.xlabel('False positive rate')
388 P.ylabel('True positive rate')
389
392 """Class to contain information and display confusion matrix.
393
394 Implementation of the `SummaryStatistics` in the case of
395 classification problem. Actual computation of confusion matrix is
396 delayed until all data is acquired (to figure out complete set of
397 labels). If testing data doesn't have a complete set of labels,
398 but you like to include all labels, provide them as a parameter to
399 the constructor.
400
401 Confusion matrix provides a set of performance statistics (use
402 asstring(description=True) for the description of abbreviations),
403 as well ROC curve (http://en.wikipedia.org/wiki/ROC_curve)
404 plotting and analysis (AUC) in the limited set of problems:
405 binary, multiclass 1-vs-all.
406 """
407
408 _STATS_DESCRIPTION = (
409 ('TP', 'true positive (AKA hit)', None),
410 ('TN', 'true negative (AKA correct rejection)', None),
411 ('FP', 'false positive (AKA false alarm, Type I error)', None),
412 ('FN', 'false negative (AKA miss, Type II error)', None),
413 ('TPR', 'true positive rate (AKA hit rate, recall, sensitivity)',
414 'TPR = TP / P = TP / (TP + FN)'),
415 ('FPR', 'false positive rate (AKA false alarm rate, fall-out)',
416 'FPR = FP / N = FP / (FP + TN)'),
417 ('ACC', 'accuracy', 'ACC = (TP + TN) / (P + N)'),
418 ('SPC', 'specificity', 'SPC = TN / (FP + TN) = 1 - FPR'),
419 ('PPV', 'positive predictive value (AKA precision)',
420 'PPV = TP / (TP + FP)'),
421 ('NPV', 'negative predictive value', 'NPV = TN / (TN + FN)'),
422 ('FDR', 'false discovery rate', 'FDR = FP / (FP + TP)'),
423 ('MCC', "Matthews Correlation Coefficient",
424 "MCC = (TP*TN - FP*FN)/sqrt(P N P' N')"),
425 ('AUC', "Area under (AUC) curve", None),
426 ) + SummaryStatistics._STATS_DESCRIPTION
427
428
429 - def __init__(self, labels=None, labels_map=None, **kwargs):
430 """Initialize ConfusionMatrix with optional list of `labels`
431
432 :Parameters:
433 labels : list
434 Optional set of labels to include in the matrix
435 labels_map : None or dict
436 Dictionary from original dataset to show mapping into
437 numerical labels
438 targets
439 Optional set of targets
440 predictions
441 Optional set of predictions
442 """
443
444 SummaryStatistics.__init__(self, **kwargs)
445
446 if labels == None:
447 labels = []
448 self.__labels = labels
449 """List of known labels"""
450 self.__labels_map = labels_map
451 """Mapping from original into given labels"""
452 self.__matrix = None
453 """Resultant confusion matrix"""
454
455
456
457
458 @property
464
465
467 """Actually compute the confusion matrix based on all the sets"""
468
469 super(ConfusionMatrix, self)._compute()
470
471 if __debug__:
472 if not self.__matrix is None:
473 debug("LAZY",
474 "Have to recompute %s#%s" \
475 % (self.__class__.__name__, id(self)))
476
477
478
479
480 try:
481
482 labels = \
483 list(reduce(lambda x, y: x.union(Set(y[0]).union(Set(y[1]))),
484 self.sets,
485 Set(self.__labels)))
486 except:
487 labels = self.__labels
488
489
490 labels_map = self.__labels_map
491 if labels_map is not None:
492 labels_set = Set(labels)
493 map_labels_set = Set(labels_map.values())
494
495 if not map_labels_set.issuperset(labels_set):
496 warning("Provided labels_map %s is not coherent with labels "
497 "provided to ConfusionMatrix. No reverse mapping "
498 "will be provided" % labels_map)
499 labels_map = None
500
501
502 labels_map_rev = None
503 if labels_map is not None:
504 labels_map_rev = {}
505 for k,v in labels_map.iteritems():
506 v_mapping = labels_map_rev.get(v, [])
507 v_mapping.append(k)
508 labels_map_rev[v] = v_mapping
509 self.__labels_map_rev = labels_map_rev
510
511 labels.sort()
512 self.__labels = labels
513
514 Nlabels, Nsets = len(labels), len(self.sets)
515
516 if __debug__:
517 debug("CM", "Got labels %s" % labels)
518
519
520 mat_all = N.zeros( (Nsets, Nlabels, Nlabels), dtype=int )
521
522
523
524
525 counts_all = N.zeros( (Nsets, Nlabels) )
526
527
528 rev_map = dict([ (x[1], x[0]) for x in enumerate(labels)])
529 for iset, set_ in enumerate(self.sets):
530 for t,p in zip(*set_[:2]):
531 mat_all[iset, rev_map[p], rev_map[t]] += 1
532
533
534
535
536
537 self.__matrix = N.sum(mat_all, axis=0)
538 self.__Nsamples = N.sum(self.__matrix, axis=0)
539 self.__Ncorrect = sum(N.diag(self.__matrix))
540
541 TP = N.diag(self.__matrix)
542 offdiag = self.__matrix - N.diag(TP)
543 stats = {
544 '# of labels' : Nlabels,
545 'TP' : TP,
546 'FP' : N.sum(offdiag, axis=1),
547 'FN' : N.sum(offdiag, axis=0)}
548
549 stats['CORR'] = N.sum(TP)
550 stats['TN'] = stats['CORR'] - stats['TP']
551 stats['P'] = stats['TP'] + stats['FN']
552 stats['N'] = N.sum(stats['P']) - stats['P']
553 stats["P'"] = stats['TP'] + stats['FP']
554 stats["N'"] = stats['TN'] + stats['FN']
555 stats['TPR'] = stats['TP'] / (1.0*stats['P'])
556
557
558 stats['TPR'][stats['P'] == 0] = 0
559 stats['PPV'] = stats['TP'] / (1.0*stats["P'"])
560 stats['NPV'] = stats['TN'] / (1.0*stats["N'"])
561 stats['FDR'] = stats['FP'] / (1.0*stats["P'"])
562 stats['SPC'] = (stats['TN']) / (1.0*stats['FP'] + stats['TN'])
563
564 MCC_denom = N.sqrt(1.0*stats['P']*stats['N']*stats["P'"]*stats["N'"])
565 nz = MCC_denom!=0.0
566 stats['MCC'] = N.zeros(stats['TP'].shape)
567 stats['MCC'][nz] = \
568 (stats['TP'] * stats['TN'] - stats['FP'] * stats['FN'])[nz] \
569 / MCC_denom[nz]
570
571 stats['ACC'] = N.sum(TP)/(1.0*N.sum(stats['P']))
572 stats['ACC%'] = stats['ACC'] * 100.0
573
574
575
576 ROC = ROCCurve(labels=labels, sets=self.sets)
577 aucs = ROC.aucs
578 if len(aucs)>0:
579 stats['AUC'] = aucs
580 if len(aucs) != Nlabels:
581 raise RuntimeError, \
582 "We must got a AUC per label. Got %d instead of %d" % \
583 (len(aucs), Nlabels)
584 self.ROC = ROC
585 else:
586
587 stats['AUC'] = [N.nan] * Nlabels
588 self.ROC = None
589
590
591
592 for k,v in stats.items():
593 stats['mean(%s)' % k] = N.mean(v)
594
595 self._stats.update(stats)
596
597
598 - def asstring(self, short=False, header=True, summary=True,
599 description=False):
600 """'Pretty print' the matrix
601
602 :Parameters:
603 short : bool
604 if True, ignores the rest of the parameters and provides consise
605 1 line summary
606 header : bool
607 print header of the table
608 summary : bool
609 print summary (accuracy)
610 description : bool
611 print verbose description of presented statistics
612 """
613 if len(self.sets) == 0:
614 return "Empty"
615
616 self.compute()
617
618
619 labels = self.__labels
620 labels_map_rev = self.__labels_map_rev
621 matrix = self.__matrix
622
623 labels_rev = []
624 if labels_map_rev is not None:
625 labels_rev = [','.join([str(x) for x in labels_map_rev[l]])
626 for l in labels]
627
628 out = StringIO()
629
630 Nlabels = len(labels)
631 Nsamples = self.__Nsamples.astype(int)
632
633 stats = self._stats
634 if short:
635 return "%(# of sets)d sets %(# of labels)d labels " \
636 " ACC:%(ACC).2f" \
637 % stats
638
639 Ndigitsmax = int(ceil(log10(max(Nsamples))))
640 Nlabelsmax = max( [len(str(x)) for x in labels] )
641
642
643 L = max(Ndigitsmax+2, Nlabelsmax)
644 res = ""
645
646 stats_perpredict = ["P'", "N'", 'FP', 'FN', 'PPV', 'NPV', 'TPR',
647 'SPC', 'FDR', 'MCC']
648
649 if self.ROC is not None: stats_perpredict += [ 'AUC' ]
650 stats_pertarget = ['P', 'N', 'TP', 'TN']
651 stats_summary = ['ACC', 'ACC%', '# of sets']
652
653
654
655 prefixlen = Nlabelsmax + 1
656 pref = ' '*(prefixlen)
657
658 if matrix.shape != (Nlabels, Nlabels):
659 raise ValueError, \
660 "Number of labels %d doesn't correspond the size" + \
661 " of a confusion matrix %s" % (Nlabels, matrix.shape)
662
663
664 printed = []
665 underscores = [" %s" % ("-" * L)] * Nlabels
666 if header:
667
668 printed.append(['@l----------. '] + labels_rev)
669 printed.append(['@lpredictions\\targets'] + labels)
670
671 printed.append(['@l `------'] \
672 + underscores + stats_perpredict)
673
674
675 for i, line in enumerate(matrix):
676 l = labels[i]
677 if labels_rev != []:
678 l = '@r%10s / %s' % (labels_rev[i], l)
679 printed.append(
680 [l] +
681 [ str(x) for x in line ] +
682 [ _p2(stats[x][i]) for x in stats_perpredict])
683
684 if summary:
685
686
687
688
689
690
691 printed.append(['@lPer target:'] + underscores)
692 for stat in stats_pertarget:
693 printed.append([stat] + [
694 _p2(stats[stat][i]) for i in xrange(Nlabels)])
695
696
697
698
699 mean_stats = N.mean(N.array([stats[k] for k in stats_perpredict]),
700 axis=1)
701 printed.append(['@lSummary \ Means:'] + underscores
702 + [_p2(stats['mean(%s)' % x])
703 for x in stats_perpredict])
704
705 for stat in stats_summary:
706 printed.append([stat] + [_p2(stats[stat])])
707
708 table2string(printed, out)
709
710 if description:
711 out.write("\nStatistics computed in 1-vs-rest fashion per each " \
712 "target.\n")
713 out.write("Abbreviations (for details see " \
714 "http://en.wikipedia.org/wiki/ROC_curve):\n")
715 for d, val, eq in self._STATS_DESCRIPTION:
716 out.write(" %-3s: %s\n" % (d, val))
717 if eq is not None:
718 out.write(" " + eq + "\n")
719
720
721 result = out.getvalue()
722 out.close()
723 return result
724
725
726 - def plot(self, labels=None, numbers=False, origin='upper',
727 numbers_alpha=None, xlabels_vertical=True, numbers_kwargs={},
728 **kwargs):
729 """Provide presentation of confusion matrix in image
730
731 :Parameters:
732 labels : list of int or basestring
733 Optionally provided labels guarantee the order of
734 presentation. Also value of None places empty column/row,
735 thus provides visual groupping of labels (Thanks Ingo)
736 numbers : bool
737 Place values inside of confusion matrix elements
738 numbers_alpha : None or float
739 Controls textual output of numbers. If None -- all numbers
740 are plotted in the same intensity. If some float -- it controls
741 alpha level -- higher value would give higher contrast. (good
742 value is 2)
743 origin : basestring
744 Which left corner diagonal should start
745 xlabels_vertical : bool
746 Either to plot xlabels vertical (benefitial if number of labels
747 is large)
748 numbers_kwargs : dict
749 Additional keyword parameters to be added to numbers (if numbers
750 is True)
751 **kwargs
752 Additional arguments given to imshow (\eg me cmap)
753
754 :Returns:
755 (fig, im, cb) -- figure, imshow, colorbar
756 """
757
758 externals.exists("pylab", raiseException=True)
759 import pylab as P
760
761 self.compute()
762 labels_order = labels
763
764
765 labels = self.__labels
766 labels_map = self.__labels_map
767 labels_map_rev = self.__labels_map_rev
768 matrix = self.__matrix
769
770
771 labels_indexes = dict([(x,i) for i,x in enumerate(labels)])
772
773 labels_rev = []
774 if labels_map_rev is not None:
775 labels_rev = [','.join([str(x) for x in labels_map_rev[l]])
776 for l in labels]
777 labels_map_full = dict(zip(labels_rev, labels))
778
779 if labels_order is not None:
780 labels_order_filtered = filter(lambda x:x is not None, labels_order)
781 labels_order_filtered_set = Set(labels_order_filtered)
782
783 if Set(labels) == labels_order_filtered_set:
784
785 labels_plot = labels_order
786 elif len(labels_rev) \
787 and Set(labels_rev) == labels_order_filtered_set:
788
789
790 labels_plot = []
791 for l in labels_order:
792 v = None
793 if l is not None: v = labels_map_full[l]
794 labels_plot += [v]
795 else:
796 raise ValueError, \
797 "Provided labels %s do not match set of known " \
798 "original labels (%s) or mapped labels (%s)" % \
799 (labels_order, labels, labels_rev)
800 else:
801 labels_plot = labels
802
803
804 isempty = N.array([l is None for l in labels_plot])
805 non_empty = N.where(N.logical_not(isempty))[0]
806
807 NlabelsNN = len(non_empty)
808 Nlabels = len(labels_plot)
809
810 if matrix.shape != (NlabelsNN, NlabelsNN):
811 raise ValueError, \
812 "Number of labels %d doesn't correspond the size" + \
813 " of a confusion matrix %s" % (NlabelsNN, matrix.shape)
814
815 confusionmatrix = N.zeros((Nlabels, Nlabels))
816 mask = confusionmatrix.copy()
817 ticks = []
818 tick_labels = []
819
820 reordered_indexes = [labels_indexes[i] for i in labels_plot
821 if i is not None]
822 for i, l in enumerate(labels_plot):
823 if l is not None:
824 j = labels_indexes[l]
825 confusionmatrix[i, non_empty] = matrix[j, reordered_indexes]
826 confusionmatrix[non_empty, i] = matrix[reordered_indexes, j]
827 ticks += [i + 0.5]
828 if labels_map_rev is not None:
829 tick_labels += ['/'.join(labels_map_rev[l])]
830 else:
831 tick_labels += [str(l)]
832 else:
833 mask[i, :] = mask[:, i] = 1
834
835 confusionmatrix = N.ma.MaskedArray(confusionmatrix, mask=mask)
836
837
838 if P.matplotlib.get_backend() == 'TkAgg':
839 P.ioff()
840
841 fig = P.gcf()
842 ax = P.gca()
843 ax.axis('off')
844
845
846 xticks_position, yticks, ybottom = {
847 'upper': ('top', [Nlabels-x for x in ticks], 0.1),
848 'lower': ('bottom', ticks, 0.2)
849 }[origin]
850
851
852
853 axi = fig.add_axes([0.15, ybottom, 0.7, 0.7])
854 im = axi.imshow(confusionmatrix, interpolation="nearest", origin=origin,
855 aspect='equal', extent=(0, Nlabels, 0, Nlabels),
856 **kwargs)
857
858
859 if numbers:
860 numbers_kwargs_ = {'fontsize': 10,
861 'horizontalalignment': 'center',
862 'verticalalignment': 'center'}
863 maxv = float(N.max(confusionmatrix))
864 colors = [im.to_rgba(0), im.to_rgba(maxv)]
865 for i,j in zip(*N.logical_not(mask).nonzero()):
866 v = confusionmatrix[j, i]
867
868 if numbers_alpha is None:
869 alpha = 1.0
870 else:
871
872 alpha = 1 - N.array(1 - v / maxv) ** numbers_alpha
873 y = {'lower':j, 'upper':Nlabels-j-1}[origin]
874 numbers_kwargs_['color'] = colors[int(v<maxv/2)]
875 numbers_kwargs_.update(numbers_kwargs)
876 P.text(i+0.5, y+0.5, '%d' % v, alpha=alpha, **numbers_kwargs_)
877
878 maxv = N.max(confusionmatrix)
879 boundaries = N.linspace(0, maxv, N.min((maxv, 10)), True)
880
881
882 P.xlabel("targets")
883 P.ylabel("predictions")
884
885 P.setp(axi, xticks=ticks, yticks=yticks,
886 xticklabels=tick_labels, yticklabels=tick_labels)
887
888 axi.xaxis.set_ticks_position(xticks_position)
889 axi.xaxis.set_label_position(xticks_position)
890
891 if xlabels_vertical:
892 P.setp(P.getp(axi, 'xticklabels'), rotation='vertical')
893
894 axcb = fig.add_axes([0.8, ybottom, 0.02, 0.7])
895 cb = P.colorbar(im, cax=axcb, format='%d', ticks = boundaries)
896
897 if P.matplotlib.get_backend() == 'TkAgg':
898 P.ion()
899 P.draw()
900
901 self._plotted_confusionmatrix = confusionmatrix
902 return fig, im, cb
903
904
905 @property
907 self.compute()
908 return 1.0-self.__Ncorrect*1.0/sum(self.__Nsamples)
909
910
911 @property
913 self.compute()
914 return self.__labels
915
916
918 return self.__labels_map
919
920
922 if val is None or isinstance(val, dict):
923 self.__labels_map = val
924 else:
925 raise ValueError, "Cannot set labels_map to %s" % val
926
927 self.__labels_map_rev = None
928 self._computed = False
929
930
931 @property
933 self.compute()
934 return self.__matrix
935
936
937 @property
939 self.compute()
940 return 100.0*self.__Ncorrect/sum(self.__Nsamples)
941
942 labels_map = property(fget=getLabels_map, fset=setLabels_map)
943
946 """Class to contain information and display on regression results.
947
948 """
949
950 _STATS_DESCRIPTION = (
951 ('CCe', 'Error based on correlation coefficient',
952 '1 - corr_coef'),
953 ('CCp', 'Correlation coefficient (p-value)', None),
954 ('RMSE', 'Root mean squared error', None),
955 ('STD', 'Standard deviation', None),
956 ('RMP', 'Root mean power (compare to RMSE of results)',
957 'sqrt(mean( data**2 ))'),
958 ) + SummaryStatistics._STATS_DESCRIPTION
959
960
962 """Initialize RegressionStatistics
963
964 :Parameters:
965 targets
966 Optional set of targets
967 predictions
968 Optional set of predictions
969 """
970
971 SummaryStatistics.__init__(self, **kwargs)
972
973
975 """Actually compute the confusion matrix based on all the sets"""
976
977 super(RegressionStatistics, self)._compute()
978 sets = self.sets
979 Nsets = len(sets)
980
981 stats = {}
982
983 funcs = {
984 'RMP_t': lambda p,t:rootMeanPowerFx(t),
985 'STD_t': lambda p,t:N.std(t),
986 'RMP_p': lambda p,t:rootMeanPowerFx(p),
987 'STD_p': lambda p,t:N.std(p),
988 'CCe': CorrErrorFx(),
989 'CCp': CorrErrorPFx(),
990 'RMSE': RMSErrorFx(),
991 'RMSE/RMP_t': RelativeRMSErrorFx()
992 }
993
994 for funcname, func in funcs.iteritems():
995 funcname_all = funcname + '_all'
996 stats[funcname_all] = []
997 for i, (targets, predictions, values) in enumerate(sets):
998 stats[funcname_all] += [func(predictions, targets)]
999 stats[funcname_all] = N.array(stats[funcname_all])
1000 stats[funcname] = N.mean(stats[funcname_all])
1001 stats[funcname+'_std'] = N.std(stats[funcname_all])
1002 stats[funcname+'_max'] = N.max(stats[funcname_all])
1003 stats[funcname+'_min'] = N.min(stats[funcname_all])
1004
1005
1006
1007
1008 targets, predictions = [], []
1009 for i, (targets_, predictions_, values_) in enumerate(sets):
1010 targets += list(targets_)
1011 predictions += list(predictions_)
1012
1013 for funcname, func in funcs.iteritems():
1014 funcname_all = 'Summary ' + funcname
1015 stats[funcname_all] = func(predictions, targets)
1016
1017 self._stats.update(stats)
1018
1019
1020 - def plot(self,
1021 plot=True, plot_stats=True,
1022 splot=True
1023
1024
1025
1026
1027 ):
1028 """Provide presentation of regression performance in image
1029
1030 :Parameters:
1031 plot : bool
1032 Plot regular plot of values (targets/predictions)
1033 plot_stats : bool
1034 Print basic statistics in the title
1035 splot : bool
1036 Plot scatter plot
1037
1038 :Returns:
1039 (fig, im, cb) -- figure, imshow, colorbar
1040 """
1041 externals.exists("pylab", raiseException=True)
1042 import pylab as P
1043
1044 self.compute()
1045
1046 nplots = plot + splot
1047
1048
1049 if P.matplotlib.get_backend() == 'TkAgg':
1050 P.ioff()
1051
1052 fig = P.gcf()
1053 P.clf()
1054 sps = []
1055
1056 nplot = 0
1057 if plot:
1058 nplot += 1
1059 sps.append(P.subplot(nplots, 1, nplot))
1060 xstart = 0
1061 lines = []
1062 for s in self.sets:
1063 nsamples = len(s[0])
1064 xend = xstart+nsamples
1065 xs = xrange(xstart, xend)
1066 lines += [P.plot(xs, s[0], 'b')]
1067 lines += [P.plot(xs, s[1], 'r')]
1068
1069 P.plot([xend, xend], [N.min(s[0]), N.max(s[0])], 'k--')
1070 xstart = xend
1071 if len(lines)>1:
1072 P.legend(lines[:2], ('Target', 'Prediction'))
1073 if plot_stats:
1074 P.title(self.asstring(short='very'))
1075
1076 if splot:
1077 nplot += 1
1078 sps.append(P.subplot(nplots, 1, nplot))
1079 for s in self.sets:
1080 P.plot(s[0], s[1], 'o',
1081 markeredgewidth=0.2,
1082 markersize=2)
1083 P.gca().set_aspect('equal')
1084
1085 if P.matplotlib.get_backend() == 'TkAgg':
1086 P.ion()
1087 P.draw()
1088
1089 return fig, sps
1090
1091 - def asstring(self, short=False, header=True, summary=True,
1092 description=False):
1093 """'Pretty print' the statistics"""
1094
1095 if len(self.sets) == 0:
1096 return "Empty"
1097
1098 self.compute()
1099
1100 stats = self.stats
1101
1102 if short:
1103 if short == 'very':
1104
1105 return "%(# of sets)d sets CCe=%(CCe).2f p=%(CCp).2g" \
1106 " RMSE:%(RMSE).2f" \
1107 " Summary: " \
1108 "CCe=%(Summary CCe).2f p=%(Summary CCp).2g" \
1109 % stats
1110 else:
1111 return "%(# of sets)d sets CCe=%(CCe).2f+-%(CCe_std).3f" \
1112 " RMSE=%(RMSE).2f+-%(RMSE_std).3f" \
1113 " RMSE/RMP_t=%(RMSE/RMP_t).2f+-%(RMSE/RMP_t_std).3f" \
1114 % stats
1115
1116 stats_data = ['RMP_t', 'STD_t', 'RMP_p', 'STD_p']
1117
1118 stats_ = ['CCe', 'RMSE', 'RMSE/RMP_t']
1119 stats_summary = ['# of sets']
1120
1121 out = StringIO()
1122
1123 printed = []
1124 if header:
1125
1126 printed.append(['Statistics', 'Mean', 'Std', 'Min', 'Max'])
1127
1128 printed.append(['----------', '-----', '-----', '-----', '-----'])
1129
1130 def print_stats(printed, stats_):
1131
1132 for stat in stats_:
1133 s = [stat]
1134 for suffix in ['', '_std', '_min', '_max']:
1135 s += [ _p2(stats[stat+suffix], 3) ]
1136 printed.append(s)
1137
1138 printed.append(["Data: "])
1139 print_stats(printed, stats_data)
1140 printed.append(["Results: "])
1141 print_stats(printed, stats_)
1142 printed.append(["Summary: "])
1143 printed.append(["CCe", _p2(stats['Summary CCe']), "", "p=", '%g' % stats['Summary CCp']])
1144 printed.append(["RMSE", _p2(stats['Summary RMSE'])])
1145 printed.append(["RMSE/RMP_t", _p2(stats['Summary RMSE/RMP_t'])])
1146
1147 if summary:
1148 for stat in stats_summary:
1149 printed.append([stat] + [_p2(stats[stat])])
1150
1151 table2string(printed, out)
1152
1153 if description:
1154 out.write("\nDescription of printed statistics.\n"
1155 " Suffixes: _t - targets, _p - predictions\n")
1156
1157 for d, val, eq in self._STATS_DESCRIPTION:
1158 out.write(" %-3s: %s\n" % (d, val))
1159 if eq is not None:
1160 out.write(" " + eq + "\n")
1161
1162 result = out.getvalue()
1163 out.close()
1164 return result
1165
1166
1167 @property
1171
1175 """Compute (or return) some error of a (trained) classifier on a dataset.
1176 """
1177
1178 confusion = StateVariable(enabled=False)
1179 """TODO Think that labels might be also symbolic thus can't directly
1180 be indicies of the array
1181 """
1182
1183 training_confusion = StateVariable(enabled=False,
1184 doc="Proxy training_confusion from underlying classifier.")
1185
1186
1187 - def __init__(self, clf, labels=None, train=True, **kwargs):
1188 """Initialization.
1189
1190 :Parameters:
1191 clf : Classifier
1192 Either trained or untrained classifier
1193 labels : list
1194 if provided, should be a set of labels to add on top of the
1195 ones present in testdata
1196 train : bool
1197 unless train=False, classifier gets trained if
1198 trainingdata provided to __call__
1199 """
1200 ClassWithCollections.__init__(self, **kwargs)
1201 self.__clf = clf
1202
1203 self._labels = labels
1204 """Labels to add on top to existing in testing data"""
1205
1206 self.__train = train
1207 """Either to train classifier if trainingdata is provided"""
1208
1209
1210 __doc__ = enhancedDocString('ClassifierError', locals(), ClassWithCollections)
1211
1212
1218
1219
1220 - def _precall(self, testdataset, trainingdataset=None):
1221 """Generic part which trains the classifier if necessary
1222 """
1223 if not trainingdataset is None:
1224 if self.__train:
1225
1226
1227
1228
1229
1230
1231
1232
1233
1234
1235
1236 if self.states.isEnabled('training_confusion'):
1237 self.__clf.states._changeTemporarily(
1238 enable_states=['training_confusion'])
1239 self.__clf.train(trainingdataset)
1240 if self.states.isEnabled('training_confusion'):
1241 self.training_confusion = self.__clf.training_confusion
1242 self.__clf.states._resetEnabledTemporarily()
1243
1244 if self.__clf.states.isEnabled('trained_labels') and \
1245 not testdataset is None:
1246 newlabels = Set(testdataset.uniquelabels) \
1247 - Set(self.__clf.trained_labels)
1248 if len(newlabels)>0:
1249 warning("Classifier %s wasn't trained to classify labels %s" %
1250 (`self.__clf`, `newlabels`) +
1251 " present in testing dataset. Make sure that you have" +
1252 " not mixed order/names of the arguments anywhere")
1253
1254
1255
1256
1257
1258
1259
1260
1261
1262
1263
1264
1265 - def _call(self, testdataset, trainingdataset=None):
1266 raise NotImplementedError
1267
1268
1269 - def _postcall(self, testdataset, trainingdataset=None, error=None):
1271
1272
1273 - def __call__(self, testdataset, trainingdataset=None):
1274 """Compute the transfer error for a certain test dataset.
1275
1276 If `trainingdataset` is not `None` the classifier is trained using the
1277 provided dataset before computing the transfer error. Otherwise the
1278 classifier is used in it's current state to make the predictions on
1279 the test dataset.
1280
1281 Returns a scalar value of the transfer error.
1282 """
1283 self._precall(testdataset, trainingdataset)
1284 error = self._call(testdataset, trainingdataset)
1285 self._postcall(testdataset, trainingdataset, error)
1286 if __debug__:
1287 debug('CERR', 'Classifier error on %s: %.2f'
1288 % (testdataset, error))
1289 return error
1290
1291
1292 @property
1295
1296
1297 @property
1300
1304 """Compute the transfer error of a (trained) classifier on a dataset.
1305
1306 The actual error value is computed using a customizable error function.
1307 Optionally the classifier can be trained by passing an additional
1308 training dataset to the __call__() method.
1309 """
1310
1311 null_prob = StateVariable(enabled=True,
1312 doc="Stores the probability of an error result under "
1313 "the NULL hypothesis")
1314 samples_error = StateVariable(enabled=False,
1315 doc="Per sample errors computed by invoking the "
1316 "error function for each sample individually. "
1317 "Errors are available in a dictionary with each "
1318 "samples origid as key.")
1319
1322 """Initialization.
1323
1324 :Parameters:
1325 clf : Classifier
1326 Either trained or untrained classifier
1327 errorfx
1328 Functor that computes a scalar error value from the vectors of
1329 desired and predicted values (e.g. subclass of `ErrorFunction`)
1330 labels : list
1331 if provided, should be a set of labels to add on top of the
1332 ones present in testdata
1333 null_dist : instance of distribution estimator
1334 """
1335 ClassifierError.__init__(self, clf, labels, **kwargs)
1336 self.__errorfx = errorfx
1337 self.__null_dist = autoNullDist(null_dist)
1338
1339
1340 __doc__ = enhancedDocString('TransferError', locals(), ClassifierError)
1341
1342
1350
1351
1352 - def _call(self, testdataset, trainingdataset=None):
1353 """Compute the transfer error for a certain test dataset.
1354
1355 If `trainingdataset` is not `None` the classifier is trained using the
1356 provided dataset before computing the transfer error. Otherwise the
1357 classifier is used in it's current state to make the predictions on
1358 the test dataset.
1359
1360 Returns a scalar value of the transfer error.
1361 """
1362
1363 clf = self.clf
1364 if testdataset is None:
1365
1366
1367 import traceback as tb
1368 filenames = [x[0] for x in tb.extract_stack(limit=100)]
1369 rfe_matches = [f for f in filenames if f.endswith('/rfe.py')]
1370 cv_matches = [f for f in filenames if
1371 f.endswith('cvtranserror.py')]
1372 msg = ""
1373 if len(rfe_matches) > 0 and len(cv_matches):
1374 msg = " It is possible that you used RFE with stopping " \
1375 "criterion based on the TransferError and directly" \
1376 " from CrossValidatedTransferError, such approach" \
1377 " would require exposing testing dataset " \
1378 " to the classifier which might heavily bias " \
1379 " generalization performance estimate. If you are " \
1380 " sure to use it that way, create CVTE with " \
1381 " parameter expose_testdataset=True"
1382 raise ValueError, "Transfer error call obtained None " \
1383 "as a dataset for testing.%s" % msg
1384 predictions = clf.predict(testdataset.samples)
1385
1386
1387
1388
1389
1390
1391
1392 states = self.states
1393 if states.isEnabled('confusion'):
1394 confusion = clf._summaryClass(
1395
1396 targets=testdataset.labels,
1397 predictions=predictions,
1398 values=clf.states.get('values', None))
1399 try:
1400 confusion.labels_map = testdataset.labels_map
1401 except:
1402 pass
1403 states.confusion = confusion
1404
1405 if states.isEnabled('samples_error'):
1406 samples_error = []
1407 for i, p in enumerate(predictions):
1408 samples_error.append(self.__errorfx(p, testdataset.labels[i]))
1409
1410 states.samples_error = dict(zip(testdataset.origids, samples_error))
1411
1412
1413 error = self.__errorfx(predictions, testdataset.labels)
1414
1415 return error
1416
1417
1418 - def _postcall(self, vdata, wdata=None, error=None):
1419 """
1420 """
1421
1422
1423 if not self.__null_dist is None and not wdata is None:
1424
1425
1426
1427 null_terr = copy.copy(self)
1428 null_terr.__null_dist = None
1429 self.__null_dist.fit(null_terr, wdata, vdata)
1430
1431
1432
1433 if not error is None and not self.__null_dist is None:
1434 self.null_prob = self.__null_dist.p(error)
1435
1436
1437 @property
1438 - def errorfx(self): return self.__errorfx
1439
1440 @property
1441 - def null_dist(self): return self.__null_dist
1442
1446 """For a given classifier report an error based on internally
1447 computed error measure (given by some `ConfusionMatrix` stored in
1448 some state variable of `Classifier`).
1449
1450 This way we can perform feature selection taking as the error
1451 criterion either learning error, or transfer to splits error in
1452 the case of SplitClassifier
1453 """
1454
1455 - def __init__(self, clf, labels=None, confusion_state="training_confusion",
1456 **kwargs):
1457 """Initialization.
1458
1459 :Parameters:
1460 clf : Classifier
1461 Either trained or untrained classifier
1462 confusion_state
1463 Id of the state variable which stores `ConfusionMatrix`
1464 labels : list
1465 if provided, should be a set of labels to add on top of the
1466 ones present in testdata
1467 """
1468 ClassifierError.__init__(self, clf, labels, **kwargs)
1469
1470 self.__confusion_state = confusion_state
1471 """What state to extract from"""
1472
1473 if not clf.states.isKnown(confusion_state):
1474 raise ValueError, \
1475 "State variable %s is not defined for classifier %s" % \
1476 (confusion_state, `clf`)
1477 if not clf.states.isEnabled(confusion_state):
1478 if __debug__:
1479 debug('CERR', "Forcing state %s to be enabled for %s" %
1480 (confusion_state, `clf`))
1481 clf.states.enable(confusion_state)
1482
1483
1484 __doc__ = enhancedDocString('ConfusionBasedError', locals(),
1485 ClassifierError)
1486
1487
1488 - def _call(self, testdata, trainingdata=None):
1494