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