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