Package mvpa :: Package clfs :: Module transerror
[hide private]
[frames] | no frames]

Source Code for Module mvpa.clfs.transerror

   1  #emacs: -*- mode: python-mode; py-indent-offset: 4; indent-tabs-mode: nil -*- 
   2  #ex: set sts=4 ts=4 sw=4 et: 
   3  ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ## 
   4  # 
   5  #   See COPYING file distributed along with the PyMVPA package for the 
   6  #   copyright and license terms. 
   7  # 
   8  ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ## 
   9  """Utility class to compute the transfer error of classifiers.""" 
  10   
  11  __docformat__ = 'restructuredtext' 
  12   
  13  import mvpa.misc.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, Stateful 
  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
53 54 55 -class SummaryStatistics(object):
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 # enforce labels in predictions to be of the same datatype as in 116 # targets, since otherwise we are getting doubles for unknown at a 117 # given moment labels 118 nonetype = type(None) 119 for i in xrange(len(targets)): 120 t1, t2 = type(targets[i]), type(predictions[i]) 121 # if there were no prediction made - leave None, otherwise 122 # convert to appropriate type 123 if t1 != t2 and t2 != nonetype: 124 #warning("Obtained target %s and prediction %s are of " % 125 # (t1, t2) + "different datatypes.") 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
152 - def __str__(self):
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
166 - def __iadd__(self, other):
167 """Add the sets from `other` s `SummaryStatistics` to current one 168 """ 169 #print "adding ", other, " to ", self 170 # need to do shallow copy, or otherwise smth like "cm += cm" 171 # would loop forever and exhaust memory eventually 172 othersets = copy.copy(other.__sets) 173 for set in othersets: 174 self.add(*set)#[0], set[1]) 175 return self
176 177
178 - def __add__(self, other):
179 """Add two `SummaryStatistics`s 180 """ 181 result = copy.copy(self) 182 result += other 183 return result
184 185
186 - def compute(self):
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
195 - def _compute(self):
196 self._stats = {'# of sets' : len(self.sets)}
197 198 199 @property
200 - def summaries(self):
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
206 - def error(self):
207 raise NotImplementedError
208 209 210 @property
211 - def stats(self):
212 self.compute() 213 return self._stats
214 215
216 - def reset(self):
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
225 226 -class ROCCurve(object):
227 """Generic class for ROC curve computation and plotting 228 """ 229
230 - def __init__(self, labels, sets=None):
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
244 - def _compute(self):
245 """Lazy computation if needed 246 """ 247 if self.__computed: 248 return 249 # local bindings 250 labels = self._labels 251 Nlabels = len(labels) 252 sets = self._sets 253 254 # take sets which have values in the shape we can handle 255 def _checkValues(set_): 256 """Check if values are 'acceptable'""" 257 x = set_[2] 258 # TODO: OPT: need optimization 259 if (x is None) or len(x) == 0: return False # undefined 260 for v in x: 261 try: 262 if Nlabels <= 2 and N.isscalar(v): 263 continue 264 if (isinstance(v, dict) or # not dict for pairs 265 ((Nlabels>=2) and len(v)!=Nlabels) # 1 per each label for multiclass 266 ): return False 267 except Exception, e: 268 # Something else which is not supported, like 269 # in shogun interface we don't yet extract values per each label or 270 # in pairs in the case of built-in multiclass 271 if __debug__: 272 debug('ROC', "Exception %s while checking " 273 "either %s are valid labels" % (str(e), x)) 274 return False 275 return True
276 277 sets_wv = filter(_checkValues, sets) 278 # check if all had values, if not -- complain 279 Nsets_wv = len(sets_wv) 280 if Nsets_wv > 0 and len(sets) != Nsets_wv: 281 warning("Only %d sets have values assigned from %d sets" % 282 (Nsets_wv, len(sets))) 283 # bring all values to the same 'shape': 284 # 1 value per each label. In binary classifier, if only a single 285 # value is provided, add '0' for 0th label 'value'... it should 286 # work taking drunk Yarik logic ;-) 287 for iset,s in enumerate(sets_wv): 288 # we will do inplace modification, thus go by index 289 values = s[2] 290 # we would need it to be a list to reassign element with a list 291 if isinstance(values, N.ndarray) and len(values.shape)==1: 292 values = list(values) 293 rangev = None 294 for i in xrange(len(values)): 295 v = values[i] 296 if N.isscalar(v): 297 if Nlabels == 2: 298 def last_el(x): 299 """Helper function. Returns x if x is scalar, and 300 last element if x is not (ie list/tuple)""" 301 if N.isscalar(x): return x 302 else: return x[-1]
303 if rangev is None: 304 # we need to figure out min/max values 305 # to invert for the 0th label 306 values_ = [last_el(x) for x in values] 307 rangev = N.min(values_) + N.max(values_) 308 values[i] = [rangev - v, v] 309 else: 310 raise ValueError, \ 311 "Cannot have a single 'value' for multiclass" \ 312 " classification. Got %s" % (v) 313 elif len(v) != Nlabels: 314 raise ValueError, \ 315 "Got %d values whenever there is %d labels" % \ 316 (len(v), Nlabels) 317 # reassign possibly adjusted values 318 sets_wv[iset] = (s[0], s[1], N.asarray(values)) 319 320 321 # we need to estimate ROC per each label 322 # XXX order of labels might not correspond to the one among 'values' 323 # which were used to make a decision... check 324 ROCs, aucs = [], [] # 1 per label 325 for i,label in enumerate(labels): 326 aucs_pl = [] 327 ROCs_pl = [] 328 for s in sets_wv: 329 targets_pl = (s[0] == label).astype(int) 330 # XXX we might unify naming between AUC/ROC 331 ROC = AUCErrorFx() 332 aucs_pl += [ROC([x[i] for x in s[2]], targets_pl)] 333 ROCs_pl.append(ROC) 334 if len(aucs_pl)>0: 335 ROCs += [ROCs_pl] 336 aucs += [nanmean(aucs_pl)] 337 #aucs += [N.mean(aucs_pl)] 338 339 # store results within the object 340 self._ROCs = ROCs 341 self._aucs = aucs 342 self.__computed = True 343 344 345 @property
346 - def aucs(self):
347 """Compute and return set of AUC values 1 per label 348 """ 349 self._compute() 350 return self._aucs
351 352 353 @property
354 - def ROCs(self):
355 self._compute() 356 return self._ROCs
357 358
359 - def plot(self, label_index=0):
360 """ 361 362 TODO: make it friendly to labels given by values? 363 should we also treat labels_map? 364 """ 365 externals.exists("pylab", raiseException=True) 366 import pylab as P 367 368 self._compute() 369 370 labels = self._labels 371 # select only ROCs for the given label 372 ROCs = self.ROCs[label_index] 373 374 fig = P.gcf() 375 ax = P.gca() 376 377 P.plot([0, 1], [0, 1], 'k:') 378 379 for ROC in ROCs: 380 P.plot(ROC.fp, ROC.tp, linewidth=1) 381 382 P.axis((0.0, 1.0, 0.0, 1.0)) 383 P.axis('scaled') 384 P.title('Label %s. Mean AUC=%.2f' % (label_index, self.aucs[label_index])) 385 386 P.xlabel('False positive rate') 387 P.ylabel('True positive rate')
388
389 390 -class ConfusionMatrix(SummaryStatistics):
391 """Class to contain information and display confusion matrix. 392 393 Implementation of the `SummaryStatistics` in the case of 394 classification problem. Actual computation of confusion matrix is 395 delayed until all data is acquired (to figure out complete set of 396 labels). If testing data doesn't have a complete set of labels, 397 but you like to include all labels, provide them as a parameter to 398 the constructor. 399 400 Confusion matrix provides a set of performance statistics (use 401 asstring(description=True) for the description of abbreviations), 402 as well ROC curve (http://en.wikipedia.org/wiki/ROC_curve) 403 plotting and analysis (AUC) in the limited set of problems: 404 binary, multiclass 1-vs-all. 405 """ 406 407 _STATS_DESCRIPTION = ( 408 ('TP', 'true positive (AKA hit)', None), 409 ('TN', 'true negative (AKA correct rejection)', None), 410 ('FP', 'false positive (AKA false alarm, Type I error)', None), 411 ('FN', 'false negative (AKA miss, Type II error)', None), 412 ('TPR', 'true positive rate (AKA hit rate, recall, sensitivity)', 413 'TPR = TP / P = TP / (TP + FN)'), 414 ('FPR', 'false positive rate (AKA false alarm rate, fall-out)', 415 'FPR = FP / N = FP / (FP + TN)'), 416 ('ACC', 'accuracy', 'ACC = (TP + TN) / (P + N)'), 417 ('SPC', 'specificity', 'SPC = TN / (FP + TN) = 1 - FPR'), 418 ('PPV', 'positive predictive value (AKA precision)', 419 'PPV = TP / (TP + FP)'), 420 ('NPV', 'negative predictive value', 'NPV = TN / (TN + FN)'), 421 ('FDR', 'false discovery rate', 'FDR = FP / (FP + TP)'), 422 ('MCC', "Matthews Correlation Coefficient", 423 "MCC = (TP*TN - FP*FN)/sqrt(P N P' N')"), 424 ('AUC', "Area under (AUC) curve", None), 425 ) + SummaryStatistics._STATS_DESCRIPTION 426 427
428 - def __init__(self, labels=None, labels_map=None, **kwargs):
429 """Initialize ConfusionMatrix with optional list of `labels` 430 431 :Parameters: 432 labels : list 433 Optional set of labels to include in the matrix 434 labels_map : None or dict 435 Dictionary from original dataset to show mapping into 436 numerical labels 437 targets 438 Optional set of targets 439 predictions 440 Optional set of predictions 441 """ 442 443 SummaryStatistics.__init__(self, **kwargs) 444 445 if labels == None: 446 labels = [] 447 self.__labels = labels 448 """List of known labels""" 449 self.__labels_map = labels_map 450 """Mapping from original into given labels""" 451 self.__matrix = None 452 """Resultant confusion matrix"""
453 454 455 # XXX might want to remove since summaries does the same, just without 456 # supplying labels 457 @property
458 - def matrices(self):
459 """Return a list of separate confusion matrix per each stored set""" 460 return [ self.__class__(labels=self.labels, 461 labels_map=self.labels_map, 462 sets=[x]) for x in self.sets]
463 464
465 - def _compute(self):
466 """Actually compute the confusion matrix based on all the sets""" 467 468 super(ConfusionMatrix, self)._compute() 469 470 if __debug__: 471 if not self.__matrix is None: 472 debug("LAZY", 473 "Have to recompute %s#%s" \ 474 % (self.__class__.__name__, id(self))) 475 476 477 # TODO: BinaryClassifier might spit out a list of predictions for each 478 # value need to handle it... for now just keep original labels 479 try: 480 # figure out what labels we have 481 labels = \ 482 list(reduce(lambda x, y: x.union(Set(y[0]).union(Set(y[1]))), 483 self.sets, 484 Set(self.__labels))) 485 except: 486 labels = self.__labels 487 488 # Check labels_map if it was provided if it covers all the labels 489 labels_map = self.__labels_map 490 if labels_map is not None: 491 labels_set = Set(labels) 492 map_labels_set = Set(labels_map.values()) 493 494 if not map_labels_set.issuperset(labels_set): 495 warning("Provided labels_map %s is not coherent with labels " 496 "provided to ConfusionMatrix. No reverse mapping " 497 "will be provided" % labels_map) 498 labels_map = None 499 500 # Create reverse map 501 labels_map_rev = None 502 if labels_map is not None: 503 labels_map_rev = {} 504 for k,v in labels_map.iteritems(): 505 v_mapping = labels_map_rev.get(v, []) 506 v_mapping.append(k) 507 labels_map_rev[v] = v_mapping 508 self.__labels_map_rev = labels_map_rev 509 510 labels.sort() 511 self.__labels = labels # store the recomputed labels 512 513 Nlabels, Nsets = len(labels), len(self.sets) 514 515 if __debug__: 516 debug("CM", "Got labels %s" % labels) 517 518 # Create a matrix for all votes 519 mat_all = N.zeros( (Nsets, Nlabels, Nlabels), dtype=int ) 520 521 # create total number of samples of each label counts 522 # just for convinience I guess since it can always be 523 # computed from mat_all 524 counts_all = N.zeros( (Nsets, Nlabels) ) 525 526 # reverse mapping from label into index in the list of labels 527 rev_map = dict([ (x[1], x[0]) for x in enumerate(labels)]) 528 for iset, (targets, predictions, values) in enumerate(self.sets): 529 for t,p in zip(targets, predictions): 530 mat_all[iset, rev_map[p], rev_map[t]] += 1 531 532 533 # for now simply compute a sum of votes across different sets 534 # we might do something more sophisticated later on, and this setup 535 # should easily allow it 536 self.__matrix = N.sum(mat_all, axis=0) 537 self.__Nsamples = N.sum(self.__matrix, axis=0) 538 self.__Ncorrect = sum(N.diag(self.__matrix)) 539 540 TP = N.diag(self.__matrix) 541 offdiag = self.__matrix - N.diag(TP) 542 stats = { 543 '# of labels' : Nlabels, 544 'TP' : TP, 545 'FP' : N.sum(offdiag, axis=1), 546 'FN' : N.sum(offdiag, axis=0)} 547 548 stats['CORR'] = N.sum(TP) 549 stats['TN'] = stats['CORR'] - stats['TP'] 550 stats['P'] = stats['TP'] + stats['FN'] 551 stats['N'] = N.sum(stats['P']) - stats['P'] 552 stats["P'"] = stats['TP'] + stats['FP'] 553 stats["N'"] = stats['TN'] + stats['FN'] 554 stats['TPR'] = stats['TP'] / (1.0*stats['P']) 555 # reset nans in TPRs to 0s whenever there is no entries 556 # for those labels among the targets 557 stats['TPR'][stats['P'] == 0] = 0 558 stats['PPV'] = stats['TP'] / (1.0*stats["P'"]) 559 stats['NPV'] = stats['TN'] / (1.0*stats["N'"]) 560 stats['FDR'] = stats['FP'] / (1.0*stats["P'"]) 561 stats['SPC'] = (stats['TN']) / (1.0*stats['FP'] + stats['TN']) 562 563 MCC_denom = N.sqrt(1.0*stats['P']*stats['N']*stats["P'"]*stats["N'"]) 564 nz = MCC_denom!=0.0 565 stats['MCC'] = N.zeros(stats['TP'].shape) 566 stats['MCC'][nz] = \ 567 (stats['TP'] * stats['TN'] - stats['FP'] * stats['FN'])[nz] \ 568 / MCC_denom[nz] 569 570 stats['ACC'] = N.sum(TP)/(1.0*N.sum(stats['P'])) 571 stats['ACC%'] = stats['ACC'] * 100.0 572 573 # 574 # ROC computation if available 575 ROC = ROCCurve(labels=labels, sets=self.sets) 576 aucs = ROC.aucs 577 if len(aucs)>0: 578 stats['AUC'] = aucs 579 if len(aucs) != Nlabels: 580 raise RuntimeError, \ 581 "We must got a AUC per label. Got %d instead of %d" % \ 582 (len(aucs), Nlabels) 583 self.ROC = ROC 584 else: 585 # we don't want to provide ROC if it is bogus 586 stats['AUC'] = [N.nan] * Nlabels 587 self.ROC = None 588 589 590 # compute mean stats 591 for k,v in stats.items(): 592 stats['mean(%s)' % k] = N.mean(v) 593 594 self._stats.update(stats)
595 596
597 - def asstring(self, short=False, header=True, summary=True, 598 description=False):
599 """'Pretty print' the matrix 600 601 :Parameters: 602 short : bool 603 if True, ignores the rest of the parameters and provides consise 604 1 line summary 605 header : bool 606 print header of the table 607 summary : bool 608 print summary (accuracy) 609 description : bool 610 print verbose description of presented statistics 611 """ 612 if len(self.sets) == 0: 613 return "Empty" 614 615 self.compute() 616 617 # some shortcuts 618 labels = self.__labels 619 labels_map_rev = self.__labels_map_rev 620 matrix = self.__matrix 621 622 labels_rev = [] 623 if labels_map_rev is not None: 624 labels_rev = [','.join([str(x) for x in labels_map_rev[l]]) 625 for l in labels] 626 627 out = StringIO() 628 # numbers of different entries 629 Nlabels = len(labels) 630 Nsamples = self.__Nsamples.astype(int) 631 632 stats = self._stats 633 if short: 634 return "%(# of sets)d sets %(# of labels)d labels " \ 635 " ACC:%(ACC).2f" \ 636 % stats 637 638 Ndigitsmax = int(ceil(log10(max(Nsamples)))) 639 Nlabelsmax = max( [len(str(x)) for x in labels] ) 640 641 # length of a single label/value 642 L = max(Ndigitsmax+2, Nlabelsmax) #, len("100.00%")) 643 res = "" 644 645 stats_perpredict = ["P'", "N'", 'FP', 'FN', 'PPV', 'NPV', 'TPR', 646 'SPC', 'FDR', 'MCC'] 647 # print AUC only if ROC was computed 648 if self.ROC is not None: stats_perpredict += [ 'AUC' ] 649 stats_pertarget = ['P', 'N', 'TP', 'TN'] 650 stats_summary = ['ACC', 'ACC%', '# of sets'] 651 652 653 #prefixlen = Nlabelsmax + 2 + Ndigitsmax + 1 654 prefixlen = Nlabelsmax + 1 655 pref = ' '*(prefixlen) # empty prefix 656 657 if matrix.shape != (Nlabels, Nlabels): 658 raise ValueError, \ 659 "Number of labels %d doesn't correspond the size" + \ 660 " of a confusion matrix %s" % (Nlabels, matrix.shape) 661 662 # list of lists of what is printed 663 printed = [] 664 underscores = [" %s" % ("-" * L)] * Nlabels 665 if header: 666 # labels 667 printed.append(['@l----------. '] + labels_rev) 668 printed.append(['@lpredictions\\targets'] + labels) 669 # underscores 670 printed.append(['@l `------'] \ 671 + underscores + stats_perpredict) 672 673 # matrix itself 674 for i, line in enumerate(matrix): 675 l = labels[i] 676 if labels_rev != []: 677 l = '@r%10s / %s' % (labels_rev[i], l) 678 printed.append( 679 [l] + 680 [ str(x) for x in line ] + 681 [ _p2(stats[x][i]) for x in stats_perpredict]) 682 683 if summary: 684 ## Various alternative schemes ;-) 685 # printed.append([''] + underscores) 686 # printed.append(['@lPer target \ Means:'] + underscores + \ 687 # [_p2(x) for x in mean_stats]) 688 # printed.append(['Means:'] + [''] * len(labels) 689 # + [_p2(x) for x in mean_stats]) 690 printed.append(['@lPer target:'] + underscores) 691 for stat in stats_pertarget: 692 printed.append([stat] + [ 693 _p2(stats[stat][i]) for i in xrange(Nlabels)]) 694 695 # compute mean stats 696 # XXX refactor to expose them in stats as well, as 697 # mean(FCC) 698 mean_stats = N.mean(N.array([stats[k] for k in stats_perpredict]), 699 axis=1) 700 printed.append(['@lSummary \ Means:'] + underscores 701 + [_p2(stats['mean(%s)' % x]) 702 for x in stats_perpredict]) 703 704 for stat in stats_summary: 705 printed.append([stat] + [_p2(stats[stat])]) 706 707 table2string(printed, out) 708 709 if description: 710 out.write("\nStatistics computed in 1-vs-rest fashion per each " \ 711 "target.\n") 712 out.write("Abbreviations (for details see " \ 713 "http://en.wikipedia.org/wiki/ROC_curve):\n") 714 for d, val, eq in self._STATS_DESCRIPTION: 715 out.write(" %-3s: %s\n" % (d, val)) 716 if eq is not None: 717 out.write(" " + eq + "\n") 718 719 #out.write("%s" % printed) 720 result = out.getvalue() 721 out.close() 722 return result
723 724
725 - def plot(self, labels=None, numbers=False, origin='upper', 726 numbers_alpha=None, xlabels_vertical=True, numbers_kwargs={}, 727 **kwargs):
728 """Provide presentation of confusion matrix in image 729 730 :Parameters: 731 labels : list of int or basestring 732 Optionally provided labels guarantee the order of 733 presentation. Also value of None places empty column/row, 734 thus provides visual groupping of labels (Thanks Ingo) 735 numbers : bool 736 Place values inside of confusion matrix elements 737 numbers_alpha : None or float 738 Controls textual output of numbers. If None -- all numbers 739 are plotted in the same intensity. If some float -- it controls 740 alpha level -- higher value would give higher contrast. (good 741 value is 2) 742 origin : basestring 743 Which left corner diagonal should start 744 xlabels_vertical : bool 745 Either to plot xlabels vertical (benefitial if number of labels 746 is large) 747 numbers_kwargs : dict 748 Additional keyword parameters to be added to numbers (if numbers 749 is True) 750 **kwargs 751 Additional arguments given to imshow (\eg me cmap) 752 753 :Returns: 754 (fig, im, cb) -- figure, imshow, colorbar 755 """ 756 757 externals.exists("pylab", raiseException=True) 758 import pylab as P 759 760 self.compute() 761 labels_order = labels 762 763 # some shortcuts 764 labels = self.__labels 765 labels_map = self.__labels_map 766 labels_map_rev = self.__labels_map_rev 767 matrix = self.__matrix 768 769 # craft original mapping from a label into index in the matrix 770 labels_indexes = dict([(x,i) for i,x in enumerate(labels)]) 771 772 labels_rev = [] 773 if labels_map_rev is not None: 774 labels_rev = [','.join([str(x) for x in labels_map_rev[l]]) 775 for l in labels] 776 777 if labels_order is not None: 778 labels_order_filtered = filter(lambda x:x is not None, labels_order) 779 labels_order_filtered_set = Set(labels_order_filtered) 780 # Verify if all labels provided in labels 781 if Set(labels) == labels_order_filtered_set: 782 # We were provided numerical (most probably) set 783 labels_plot = labels_order 784 elif len(labels_rev) \ 785 and Set(labels_map.keys()) == labels_order_filtered_set: 786 # not clear if right whenever there were multiple labels 787 # mapped into the same 788 labels_plot = [] 789 for l in labels_order: 790 v = None 791 if l is not None: v = labels_map[l] 792 labels_plot += [v] 793 else: 794 raise ValueError, \ 795 "Provided labels %s do not match set of known " \ 796 "original labels (%s) or mapped labels (%s)" % \ 797 (labels_order, labels, labels_rev) 798 else: 799 labels_plot = labels 800 801 # where we have Nones? 802 isempty = N.array([l is None for l in labels_plot]) 803 non_empty = N.where(N.logical_not(isempty))[0] 804 # numbers of different entries 805 NlabelsNN = len(non_empty) 806 Nlabels = len(labels_plot) 807 808 if matrix.shape != (NlabelsNN, NlabelsNN): 809 raise ValueError, \ 810 "Number of labels %d doesn't correspond the size" + \ 811 " of a confusion matrix %s" % (NlabelsNN, matrix.shape) 812 813 confusionmatrix = N.zeros((Nlabels, Nlabels)) 814 mask = confusionmatrix.copy() 815 ticks = [] 816 tick_labels = [] 817 # populate in a silly way 818 819 for i, l in enumerate(labels_plot): 820 if l is not None: 821 j = labels_indexes[l] 822 confusionmatrix[i, non_empty] = matrix[j, :] 823 confusionmatrix[non_empty, i] = matrix[:, j] 824 ticks += [i + 0.5] 825 if labels_map_rev is not None: 826 tick_labels += ['/'.join(labels_map_rev[l])] 827 else: 828 tick_labels += [str(l)] 829 else: 830 mask[i, :] = mask[:, i] = 1 831 832 confusionmatrix = N.ma.MaskedArray(confusionmatrix, mask=mask) 833 834 # turn off automatic update if interactive 835 if P.matplotlib.get_backend() == 'TkAgg': 836 P.ioff() 837 838 fig = P.gcf() 839 ax = P.gca() 840 ax.axis('off') 841 842 # some customization depending on the origin 843 xticks_position, yticks, ybottom = { 844 'upper': ('top', ticks[::-1], 0.1), 845 'lower': ('bottom', ticks, 0.2) 846 }[origin] 847 848 849 # Plot 850 axi = fig.add_axes([0.15, ybottom, 0.7, 0.7]) 851 im = axi.imshow(confusionmatrix, interpolation="nearest", origin=origin, 852 aspect='equal', **kwargs) 853 854 # plot numbers 855 if numbers: 856 numbers_kwargs_ = {'fontsize': 10, 857 'horizontalalignment': 'center', 858 'verticalalignment': 'center'} 859 maxv = float(N.max(confusionmatrix)) 860 colors = [im.to_rgba(0), im.to_rgba(maxv)] 861 for i,j in zip(*N.logical_not(mask).nonzero()): 862 v = confusionmatrix[j, i] 863 # scale alpha non-linearly 864 if numbers_alpha is None: 865 alpha = 1.0 866 else: 867 # scale according to value 868 alpha = 1 - N.array(1 - v / maxv) ** numbers_alpha 869 y = {'lower':j, 'upper':Nlabels-j-1}[origin] 870 numbers_kwargs_['color'] = colors[int(v<maxv/2)] 871 numbers_kwargs_.update(numbers_kwargs) 872 P.text(i+0.5, y+0.5, '%d' % v, alpha=alpha, **numbers_kwargs_) 873 874 maxv = N.max(confusionmatrix) 875 boundaries = N.linspace(0, maxv, N.min(maxv, 10), True) 876 877 # Label axes 878 P.xlabel("targets") 879 P.ylabel("predictions") 880 881 P.setp(axi, xticks=ticks, yticks=yticks, 882 xticklabels=tick_labels, yticklabels=tick_labels) 883 884 axi.xaxis.set_ticks_position(xticks_position) 885 axi.xaxis.set_label_position(xticks_position) 886 887 if xlabels_vertical: 888 P.setp(P.getp(axi, 'xticklabels'), rotation='vertical') 889 890 axcb = fig.add_axes([0.8, ybottom, 0.02, 0.7]) 891 cb = P.colorbar(im, cax=axcb, format='%d', ticks = boundaries) 892 893 if P.matplotlib.get_backend() == 'TkAgg': 894 P.ion() 895 P.draw() 896 return fig, im, cb
897 898 899 @property
900 - def error(self):
901 self.compute() 902 return 1.0-self.__Ncorrect*1.0/sum(self.__Nsamples)
903 904 905 @property
906 - def labels(self):
907 self.compute() 908 return self.__labels
909 910
911 - def getLabels_map(self):
912 return self.__labels_map
913 914
915 - def setLabels_map(self, val):
916 if val is None or isinstance(val, dict): 917 self.__labels_map = val 918 else: 919 raise ValueError, "Cannot set labels_map to %s" % val 920 # reset it just in case 921 self.__labels_map_rev = None 922 self._computed = False
923 924 925 @property
926 - def matrix(self):
927 self.compute() 928 return self.__matrix
929 930 931 @property
932 - def percentCorrect(self):
933 self.compute() 934 return 100.0*self.__Ncorrect/sum(self.__Nsamples)
935 936 labels_map = property(fget=getLabels_map, fset=setLabels_map)
937
938 939 -class RegressionStatistics(SummaryStatistics):
940 """Class to contain information and display on regression results. 941 942 """ 943 944 _STATS_DESCRIPTION = ( 945 ('CCe', 'Error based on correlation coefficient', 946 '1 - corr_coef'), 947 ('CCp', 'Correlation coefficient (p-value)', None), 948 ('RMSE', 'Root mean squared error', None), 949 ('STD', 'Standard deviation', None), 950 ('RMP', 'Root mean power (compare to RMSE of results)', 951 'sqrt(mean( data**2 ))'), 952 ) + SummaryStatistics._STATS_DESCRIPTION 953 954
955 - def __init__(self, **kwargs):
956 """Initialize RegressionStatistics 957 958 :Parameters: 959 targets 960 Optional set of targets 961 predictions 962 Optional set of predictions 963 """ 964 965 SummaryStatistics.__init__(self, **kwargs)
966 967
968 - def _compute(self):
969 """Actually compute the confusion matrix based on all the sets""" 970 971 super(RegressionStatistics, self)._compute() 972 sets = self.sets 973 Nsets = len(sets) 974 975 stats = {} 976 977 funcs = { 978 'RMP_t': lambda p,t:rootMeanPowerFx(t), 979 'STD_t': lambda p,t:N.std(t), 980 'RMP_p': lambda p,t:rootMeanPowerFx(p), 981 'STD_p': lambda p,t:N.std(p), 982 'CCe': CorrErrorFx(), 983 'CCp': CorrErrorPFx(), 984 'RMSE': RMSErrorFx(), 985 'RMSE/RMP_t': RelativeRMSErrorFx() 986 } 987 988 for funcname, func in funcs.iteritems(): 989 funcname_all = funcname + '_all' 990 stats[funcname_all] = [] 991 for i, (targets, predictions, values) in enumerate(sets): 992 stats[funcname_all] += [func(predictions, targets)] 993 stats[funcname_all] = N.array(stats[funcname_all]) 994 stats[funcname] = N.mean(stats[funcname_all]) 995 stats[funcname+'_std'] = N.std(stats[funcname_all]) 996 stats[funcname+'_max'] = N.max(stats[funcname_all]) 997 stats[funcname+'_min'] = N.min(stats[funcname_all]) 998 999 self._stats.update(stats)
1000 1001
1002 - def asstring(self, short=False, header=True, summary=True, 1003 description=False):
1004 """'Pretty print' the statistics""" 1005 1006 if len(self.sets) == 0: 1007 return "Empty" 1008 1009 self.compute() 1010 1011 stats = self.stats 1012 1013 if short: 1014 if short == 'very': 1015 return "%(# of sets)d sets CCe:%(CCe).2f CCp:%(CCp).2g" \ 1016 " RMSE:%(RMSE).2f" \ 1017 " RMSE/RMP_t:%(RMSE/RMP_t).2f" \ 1018 % stats 1019 else: 1020 return "%(# of sets)d sets CCe:%(CCe).2f+-%(CCe_std).3f" \ 1021 " RMSE:%(RMSE).2f+-%(RMSE_std).3f" \ 1022 " RMSE/RMP_t:%(RMSE/RMP_t).2f+-%(RMSE/RMP_t_std).3f" \ 1023 % stats 1024 1025 stats_data = ['RMP_t', 'STD_t', 'RMP_p', 'STD_p'] 1026 # CCp needs tune up of format so excluded 1027 stats_ = ['CCe', 'RMSE', 'RMSE/RMP_t'] 1028 stats_summary = ['# of sets'] 1029 1030 out = StringIO() 1031 1032 printed = [] 1033 if header: 1034 # labels 1035 printed.append(['Statistics', 'Mean', 'Std', 'Min', 'Max']) 1036 # underscores 1037 printed.append(['----------', '-----', '-----', '-----', '-----']) 1038 1039 def print_stats(printed, stats_): 1040 # Statistics itself 1041 for stat in stats_: 1042 s = [stat] 1043 for suffix in ['', '_std', '_min', '_max']: 1044 s += [ _p2(stats[stat+suffix], 3) ] 1045 printed.append(s)
1046 1047 printed.append(["Data: "]) 1048 print_stats(printed, stats_data) 1049 printed.append(["Results: "]) 1050 print_stats(printed, stats_) 1051 1052 if summary: 1053 for stat in stats_summary: 1054 printed.append([stat] + [_p2(stats[stat])]) 1055 1056 table2string(printed, out) 1057 1058 if description: 1059 out.write("\nDescription of printed statistics.\n" 1060 " Suffixes: _t - targets, _p - predictions\n") 1061 1062 for d, val, eq in self._STATS_DESCRIPTION: 1063 out.write(" %-3s: %s\n" % (d, val)) 1064 if eq is not None: 1065 out.write(" " + eq + "\n") 1066 1067 result = out.getvalue() 1068 out.close() 1069 return result
1070 1071 1072 @property
1073 - def error(self):
1074 self.compute() 1075 return self.stats['RMSE']
1076
1077 1078 1079 -class ClassifierError(Stateful):
1080 """Compute (or return) some error of a (trained) classifier on a dataset. 1081 """ 1082 1083 confusion = StateVariable(enabled=False) 1084 """TODO Think that labels might be also symbolic thus can't directly 1085 be indicies of the array 1086 """ 1087 1088 training_confusion = StateVariable(enabled=False, 1089 doc="Proxy training_confusion from underlying classifier.") 1090 1091
1092 - def __init__(self, clf, labels=None, train=True, **kwargs):
1093 """Initialization. 1094 1095 :Parameters: 1096 clf : Classifier 1097 Either trained or untrained classifier 1098 labels : list 1099 if provided, should be a set of labels to add on top of the 1100 ones present in testdata 1101 train : bool 1102 unless train=False, classifier gets trained if 1103 trainingdata provided to __call__ 1104 """ 1105 Stateful.__init__(self, **kwargs) 1106 self.__clf = clf 1107 1108 self._labels = labels 1109 """Labels to add on top to existing in testing data""" 1110 1111 self.__train = train 1112 """Either to train classifier if trainingdata is provided"""
1113 1114 1115 __doc__ = enhancedDocString('ClassifierError', locals(), Stateful) 1116 1117
1118 - def __copy__(self):
1119 """TODO: think... may be we need to copy self.clf""" 1120 out = ClassifierError.__new__(TransferError) 1121 ClassifierError.__init__(out, self.clf) 1122 return out
1123 1124
1125 - def _precall(self, testdataset, trainingdataset=None):
1126 """Generic part which trains the classifier if necessary 1127 """ 1128 if not trainingdataset is None: 1129 if self.__train: 1130 # XXX can be pretty annoying if triggered inside an algorithm 1131 # where it cannot be switched of, but retraining might be 1132 # intended or at least not avoidable. 1133 # Additonally isTrained docs say: 1134 # MUST BE USED WITH CARE IF EVER 1135 # 1136 # switching it off for now 1137 #if self.__clf.isTrained(trainingdataset): 1138 # warning('It seems that classifier %s was already trained' % 1139 # self.__clf + ' on dataset %s. Please inspect' \ 1140 # % trainingdataset) 1141 if self.states.isEnabled('training_confusion'): 1142 self.__clf.states._changeTemporarily( 1143 enable_states=['training_confusion']) 1144 self.__clf.train(trainingdataset) 1145 if self.states.isEnabled('training_confusion'): 1146 self.training_confusion = self.__clf.training_confusion 1147 self.__clf.states._resetEnabledTemporarily() 1148 1149 if self.__clf.states.isEnabled('trained_labels') and \ 1150 not testdataset is None: 1151 newlabels = Set(testdataset.uniquelabels) \ 1152 - Set(self.__clf.trained_labels) 1153 if len(newlabels)>0: 1154 warning("Classifier %s wasn't trained to classify labels %s" % 1155 (`self.__clf`, `newlabels`) + 1156 " present in testing dataset. Make sure that you have" + 1157 " not mixed order/names of the arguments anywhere")
1158 1159 ### Here checking for if it was trained... might be a cause of trouble 1160 # XXX disabled since it is unreliable.. just rely on explicit 1161 # self.__train 1162 # if not self.__clf.isTrained(trainingdataset): 1163 # self.__clf.train(trainingdataset) 1164 # elif __debug__: 1165 # debug('CERR', 1166 # 'Not training classifier %s since it was ' % `self.__clf` 1167 # + ' already trained on dataset %s' % `trainingdataset`) 1168 1169
1170 - def _call(self, testdataset, trainingdataset=None):
1171 raise NotImplementedError
1172 1173
1174 - def _postcall(self, testdataset, trainingdataset=None, error=None):
1175 pass
1176 1177
1178 - def __call__(self, testdataset, trainingdataset=None):
1179 """Compute the transfer error for a certain test dataset. 1180 1181 If `trainingdataset` is not `None` the classifier is trained using the 1182 provided dataset before computing the transfer error. Otherwise the 1183 classifier is used in it's current state to make the predictions on 1184 the test dataset. 1185 1186 Returns a scalar value of the transfer error. 1187 """ 1188 self._precall(testdataset, trainingdataset) 1189 error = self._call(testdataset, trainingdataset) 1190 self._postcall(testdataset, trainingdataset, error) 1191 return error
1192 1193 1194 @property
1195 - def clf(self):
1196 return self.__clf
1197 1198 1199 @property
1200 - def labels(self):
1201 return self._labels
1202
1203 1204 1205 -class TransferError(ClassifierError):
1206 """Compute the transfer error of a (trained) classifier on a dataset. 1207 1208 The actual error value is computed using a customizable error function. 1209 Optionally the classifier can be trained by passing an additional 1210 training dataset to the __call__() method. 1211 """ 1212 1213 null_prob = StateVariable(enabled=True, 1214 doc="Stores the probability of an error result under " 1215 "the NULL hypothesis") 1216 samples_error = StateVariable(enabled=False, 1217 doc="Per sample errors computed by invoking the " 1218 "error function for each sample individually. " 1219 "Errors are available in a dictionary with each " 1220 "samples origid as key.") 1221
1222 - def __init__(self, clf, errorfx=MeanMismatchErrorFx(), labels=None, 1223 null_dist=None, **kwargs):
1224 """Initialization. 1225 1226 :Parameters: 1227 clf : Classifier 1228 Either trained or untrained classifier 1229 errorfx 1230 Functor that computes a scalar error value from the vectors of 1231 desired and predicted values (e.g. subclass of `ErrorFunction`) 1232 labels : list 1233 if provided, should be a set of labels to add on top of the 1234 ones present in testdata 1235 null_dist : instance of distribution estimator 1236 """ 1237 ClassifierError.__init__(self, clf, labels, **kwargs) 1238 self.__errorfx = errorfx 1239 self.__null_dist = autoNullDist(null_dist)
1240 1241 1242 __doc__ = enhancedDocString('TransferError', locals(), ClassifierError) 1243 1244
1245 - def __copy__(self):
1246 """Performs deepcopying of the classifier.""" 1247 # TODO -- use ClassifierError.__copy__ 1248 from mvpa.clfs.base import _deepcopyclf 1249 out = TransferError.__new__(TransferError) 1250 TransferError.__init__(out, _deepcopyclf(self.clf), self.errorfx, self._labels) 1251 1252 return out
1253 1254
1255 - def _call(self, testdataset, trainingdataset=None):
1256 """Compute the transfer error for a certain test dataset. 1257 1258 If `trainingdataset` is not `None` the classifier is trained using the 1259 provided dataset before computing the transfer error. Otherwise the 1260 classifier is used in it's current state to make the predictions on 1261 the test dataset. 1262 1263 Returns a scalar value of the transfer error. 1264 """ 1265 # OPT: local binding 1266 clf = self.clf 1267 if testdataset is None: 1268 # We cannot do anythin, but we can try to figure out WTF and 1269 # warn the user accordingly in some usecases 1270 import traceback as tb 1271 filenames = [x[0] for x in tb.extract_stack(limit=100)] 1272 rfe_matches = [f for f in filenames if f.endswith('/rfe.py')] 1273 cv_matches = [f for f in filenames if 1274 f.endswith('cvtranserror.py')] 1275 msg = "" 1276 if len(rfe_matches) > 0 and len(cv_matches): 1277 msg = " It is possible that you used RFE with stopping " \ 1278 "criterion based on the TransferError and directly" \ 1279 " from CrossValidatedTransferError, such approach" \ 1280 " would require exposing testing dataset " \ 1281 " to the classifier which might heavily bias " \ 1282 " generalization performance estimate. If you are " \ 1283 " sure to use it that way, create CVTE with " \ 1284 " parameter expose_testdataset=True" 1285 raise ValueError, "Transfer error call obtained None " \ 1286 "as a dataset for testing.%s" % msg 1287 predictions = clf.predict(testdataset.samples) 1288 1289 # compute confusion matrix 1290 # Should it migrate into ClassifierError.__postcall? 1291 # -> Probably not because other childs could estimate it 1292 # not from test/train datasets explicitely, see 1293 # `ConfusionBasedError`, where confusion is simply 1294 # bound to classifiers confusion matrix 1295 states = self.states 1296 if states.isEnabled('confusion'): 1297 confusion = clf._summaryClass( 1298 #labels=self.labels, 1299 targets=testdataset.labels, 1300 predictions=predictions, 1301 values=clf.states.get('values', None)) 1302 try: 1303 confusion.labels_map = testdataset.labels_map 1304 except: 1305 pass 1306 states.confusion = confusion 1307 1308 if states.isEnabled('samples_error'): 1309 samples_error = [] 1310 for i, p in enumerate(predictions): 1311 samples_error.append(self.__errorfx(p, testdataset.labels[i])) 1312 1313 states.samples_error = dict(zip(testdataset.origids, samples_error)) 1314 1315 # compute error from desired and predicted values 1316 error = self.__errorfx(predictions, testdataset.labels) 1317 1318 return error
1319 1320
1321 - def _postcall(self, vdata, wdata=None, error=None):
1322 """ 1323 """ 1324 # estimate the NULL distribution when functor and training data is 1325 # given 1326 if not self.__null_dist is None and not wdata is None: 1327 # we need a matching transfer error instances (e.g. same error 1328 # function), but we have to disable the estimation of the null 1329 # distribution in that child to prevent infinite looping. 1330 null_terr = copy.copy(self) 1331 null_terr.__null_dist = None 1332 self.__null_dist.fit(null_terr, wdata, vdata) 1333 1334 1335 # get probability of error under NULL hypothesis if available 1336 if not error is None and not self.__null_dist is None: 1337 self.null_prob = self.__null_dist.p(error)
1338 1339 1340 @property
1341 - def errorfx(self): return self.__errorfx
1342 1343 @property
1344 - def null_dist(self): return self.__null_dist
1345
1346 1347 1348 -class ConfusionBasedError(ClassifierError):
1349 """For a given classifier report an error based on internally 1350 computed error measure (given by some `ConfusionMatrix` stored in 1351 some state variable of `Classifier`). 1352 1353 This way we can perform feature selection taking as the error 1354 criterion either learning error, or transfer to splits error in 1355 the case of SplitClassifier 1356 """ 1357
1358 - def __init__(self, clf, labels=None, confusion_state="training_confusion", 1359 **kwargs):
1360 """Initialization. 1361 1362 :Parameters: 1363 clf : Classifier 1364 Either trained or untrained classifier 1365 confusion_state 1366 Id of the state variable which stores `ConfusionMatrix` 1367 labels : list 1368 if provided, should be a set of labels to add on top of the 1369 ones present in testdata 1370 """ 1371 ClassifierError.__init__(self, clf, labels, **kwargs) 1372 1373 self.__confusion_state = confusion_state 1374 """What state to extract from""" 1375 1376 if not clf.states.isKnown(confusion_state): 1377 raise ValueError, \ 1378 "State variable %s is not defined for classifier %s" % \ 1379 (confusion_state, `clf`) 1380 if not clf.states.isEnabled(confusion_state): 1381 if __debug__: 1382 debug('CERR', "Forcing state %s to be enabled for %s" % 1383 (confusion_state, `clf`)) 1384 clf.states.enable(confusion_state)
1385 1386 1387 __doc__ = enhancedDocString('ConfusionBasedError', locals(), 1388 ClassifierError) 1389 1390
1391 - def _call(self, testdata, trainingdata=None):
1392 """Extract transfer error. Nor testdata, neither trainingdata is used 1393 """ 1394 confusion = self.clf.states.getvalue(self.__confusion_state) 1395 self.confusion = confusion 1396 return confusion.error
1397