Changeset 7396:f4da44a215c5 in orange
 Timestamp:
 02/04/11 10:19:06 (3 years ago)
 Branch:
 default
 Convert:
 ce33b68b15f2cea16c0f122ab2edc6490267cb2e
 Location:
 orange
 Files:

 1 added
 2 edited
Legend:
 Unmodified
 Added
 Removed

orange/Orange/evaluation/scoring.py
r7388 r7396 1 """ 2 3 .. index: scoring 4 5 ================================= 6 Orange Statistics for Predictions 7 ================================= 8 9 This module contains various measures of quality for classification and 10 regression. Most functions require an argument named res, an instance of 11 :class:Orange.evaluation.testing.ExperimentResults as computed by functions from orngTest and which contains 12 predictions obtained through crossvalidation, leave oneout, testing on 13 training data or test set examples. 14 """ 15 16 import statc 17 from operator import add 18 import orngMisc, orngTest 19 import numpy 20 21 #### Private stuff 22 23 def log2(x): 24 """Calculate logarithm in base 2.""" 25 return math.log(x)/math.log(2) 26 27 def checkNonZero(x): 28 if x==0.0: 29 raise ValueError, "Cannot compute the score: no examples or sum of weights is 0.0." 30 31 def gettotweight(res): 32 totweight = reduce(lambda x, y: x+y.weight, res.results, 0) 33 if totweight==0.0: 34 raise ValueError, "Cannot compute the score: sum of weights is 0.0." 35 return totweight 36 37 def gettotsize(res): 38 if len(res.results): 39 return len(res.results) 40 else: 41 raise ValueError, "Cannot compute the score: no examples." 42 43 44 def splitByIterations(res): 45 if res.numberOfIterations < 2: 46 return [res] 47 48 ress = [orngTest.ExperimentResults(1, res.classifierNames, res.classValues, res.weights, classifiers=res.classifiers, loaded=res.loaded) 49 for i in range(res.numberOfIterations)] 50 for te in res.results: 51 ress[te.iterationNumber].results.append(te) 52 return ress 53 54 55 def classProbabilitiesFromRes(res, **argkw): 56 probs = [0.0] * len(res.classValues) 57 if argkw.get("unweighted", 0) or not res.weights: 58 for tex in res.results: 59 probs[int(tex.actualClass)] += 1.0 60 totweight = gettotsize(res) 61 else: 62 totweight = 0.0 63 for tex in res.results: 64 probs[tex.actualClass] += tex.weight 65 totweight += tex.weight 66 checkNonZero(totweight) 67 return [prob/totweight for prob in probs] 68 69 70 def statisticsByFolds(stats, foldN, reportSE, iterationIsOuter): 71 # remove empty folds, turn the matrix so that learner is outer 72 if iterationIsOuter: 73 if not stats: 74 raise ValueError, "Cannot compute the score: no examples or sum of weights is 0.0." 75 numberOfLearners = len(stats[0]) 76 stats = filter(lambda (x, fN): fN>0.0, zip(stats,foldN)) 77 stats = [ [x[lrn]/fN for x, fN in stats] for lrn in range(numberOfLearners)] 78 else: 79 stats = [ [x/Fn for x, Fn in filter(lambda (x, Fn): Fn > 0.0, zip(lrnD, foldN))] for lrnD in stats] 80 81 if not stats: 82 raise ValueError, "Cannot compute the score: no classifiers" 83 if not stats[0]: 84 raise ValueError, "Cannot compute the score: no examples or sum of weights is 0.0." 85 86 if reportSE: 87 return [(statc.mean(x), statc.sterr(x)) for x in stats] 88 else: 89 return [statc.mean(x) for x in stats] 90 91 def ME(res, **argkw): 92 MEs = [0.0]*res.numberOfLearners 93 94 if argkw.get("unweighted", 0) or not res.weights: 95 for tex in res.results: 96 MEs = map(lambda res, cls, ac = float(tex.actualClass): 97 res + abs(float(cls)  ac), MEs, tex.classes) 98 totweight = gettotsize(res) 99 else: 100 for tex in res.results: 101 MEs = map(lambda res, cls, ac = float(tex.actualClass), tw = tex.weight: 102 res + tw*abs(float(cls)  ac), MEs, tex.classes) 103 totweight = gettotweight(res) 104 105 return [x/totweight for x in MEs] 106 107 MAE = ME 108 109 ######################################################################### 110 # PERFORMANCE MEASURES: 111 # Scores for evaluation of numeric predictions 112 113 def checkArgkw(dct, lst): 114 """checkArgkw(dct, lst) > returns true if any items have nonzero value in dct""" 115 return reduce(lambda x,y: x or y, [dct.get(k, 0) for k in lst]) 116 117 def regressionError(res, **argkw): 118 """regressionError(res) > regression error (default: MSE)""" 119 if argkw.get("SE", 0) and res.numberOfIterations > 1: 120 # computes the scores for each iteration, then averages 121 scores = [[0.0] * res.numberOfIterations for i in range(res.numberOfLearners)] 122 if argkw.get("normabs", 0) or argkw.get("normsqr", 0): 123 norm = [0.0] * res.numberOfIterations 124 125 nIter = [0]*res.numberOfIterations # counts examples in each iteration 126 a = [0]*res.numberOfIterations # average class in each iteration 127 for tex in res.results: 128 nIter[tex.iterationNumber] += 1 129 a[tex.iterationNumber] += float(tex.actualClass) 130 a = [a[i]/nIter[i] for i in range(res.numberOfIterations)] 131 132 if argkw.get("unweighted", 0) or not res.weights: 133 # iterate accross test cases 134 for tex in res.results: 135 ai = float(tex.actualClass) 136 nIter[tex.iterationNumber] += 1 137 138 # compute normalization, if required 139 if argkw.get("normabs", 0): 140 norm[tex.iterationNumber] += abs(ai  a[tex.iterationNumber]) 141 elif argkw.get("normsqr", 0): 142 norm[tex.iterationNumber] += (ai  a[tex.iterationNumber])**2 143 144 # iterate accross results of different regressors 145 for i, cls in enumerate(tex.classes): 146 if argkw.get("abs", 0): 147 scores[i][tex.iterationNumber] += abs(float(cls)  ai) 148 else: 149 scores[i][tex.iterationNumber] += (float(cls)  ai)**2 150 else: # unweighted<>0 151 raise NotImplementedError, "weighted error scores with SE not implemented yet" 152 153 if argkw.get("normabs") or argkw.get("normsqr"): 154 scores = [[x/n for x, n in zip(y, norm)] for y in scores] 155 else: 156 scores = [[x/ni for x, ni in zip(y, nIter)] for y in scores] 157 158 if argkw.get("R2"): 159 scores = [[1.0  x for x in y] for y in scores] 160 161 if argkw.get("sqrt", 0): 162 scores = [[math.sqrt(x) for x in y] for y in scores] 163 164 return [(statc.mean(x), statc.std(x)) for x in scores] 165 166 else: # single iteration (testing on a single test set) 167 scores = [0.0] * res.numberOfLearners 168 norm = 0.0 169 170 if argkw.get("unweighted", 0) or not res.weights: 171 a = sum([tex.actualClass for tex in res.results]) \ 172 / len(res.results) 173 for tex in res.results: 174 if argkw.get("abs", 0): 175 scores = map(lambda res, cls, ac = float(tex.actualClass): 176 res + abs(float(cls)  ac), scores, tex.classes) 177 else: 178 scores = map(lambda res, cls, ac = float(tex.actualClass): 179 res + (float(cls)  ac)**2, scores, tex.classes) 180 181 if argkw.get("normabs", 0): 182 norm += abs(tex.actualClass  a) 183 elif argkw.get("normsqr", 0): 184 norm += (tex.actualClass  a)**2 185 totweight = gettotsize(res) 186 else: 187 # UNFINISHED 188 for tex in res.results: 189 MSEs = map(lambda res, cls, ac = float(tex.actualClass), 190 tw = tex.weight: 191 res + tw * (float(cls)  ac)**2, MSEs, tex.classes) 192 totweight = gettotweight(res) 193 194 if argkw.get("normabs", 0) or argkw.get("normsqr", 0): 195 scores = [s/norm for s in scores] 196 else: # normalize by number of instances (or sum of weights) 197 scores = [s/totweight for s in scores] 198 199 if argkw.get("R2"): 200 scores = [1.0  s for s in scores] 201 202 if argkw.get("sqrt", 0): 203 scores = [math.sqrt(x) for x in scores] 204 205 return scores 206 207 def MSE(res, **argkw): 208 """MSE(res) > meansquared error""" 209 return regressionError(res, **argkw) 210 211 def RMSE(res, **argkw): 212 """RMSE(res) > root meansquared error""" 213 argkw.setdefault("sqrt", True) 214 return regressionError(res, **argkw) 215 216 def MAE(res, **argkw): 217 """MAE(res) > mean absolute error""" 218 argkw.setdefault("abs", True) 219 return regressionError(res, **argkw) 220 221 def RSE(res, **argkw): 222 """RSE(res) > relative squared error""" 223 argkw.setdefault("normsqr", True) 224 return regressionError(res, **argkw) 225 226 def RRSE(res, **argkw): 227 """RRSE(res) > root relative squared error""" 228 argkw.setdefault("normsqr", True) 229 argkw.setdefault("sqrt", True) 230 return regressionError(res, **argkw) 231 232 def RAE(res, **argkw): 233 """RAE(res) > relative absolute error""" 234 argkw.setdefault("abs", True) 235 argkw.setdefault("normabs", True) 236 return regressionError(res, **argkw) 237 238 def R2(res, **argkw): 239 """R2(res) > Rsquared""" 240 argkw.setdefault("normsqr", True) 241 argkw.setdefault("R2", True) 242 return regressionError(res, **argkw) 243 244 def MSE_old(res, **argkw): 245 """MSE(res) > meansquared error""" 246 if argkw.get("SE", 0) and res.numberOfIterations > 1: 247 MSEs = [[0.0] * res.numberOfIterations for i in range(res.numberOfLearners)] 248 nIter = [0]*res.numberOfIterations 249 if argkw.get("unweighted", 0) or not res.weights: 250 for tex in res.results: 251 ac = float(tex.actualClass) 252 nIter[tex.iterationNumber] += 1 253 for i, cls in enumerate(tex.classes): 254 MSEs[i][tex.iterationNumber] += (float(cls)  ac)**2 255 else: 256 raise ValueError, "weighted RMSE with SE not implemented yet" 257 MSEs = [[x/ni for x, ni in zip(y, nIter)] for y in MSEs] 258 if argkw.get("sqrt", 0): 259 MSEs = [[math.sqrt(x) for x in y] for y in MSEs] 260 return [(statc.mean(x), statc.std(x)) for x in MSEs] 261 262 else: 263 MSEs = [0.0]*res.numberOfLearners 264 if argkw.get("unweighted", 0) or not res.weights: 265 for tex in res.results: 266 MSEs = map(lambda res, cls, ac = float(tex.actualClass): 267 res + (float(cls)  ac)**2, MSEs, tex.classes) 268 totweight = gettotsize(res) 269 else: 270 for tex in res.results: 271 MSEs = map(lambda res, cls, ac = float(tex.actualClass), tw = tex.weight: 272 res + tw * (float(cls)  ac)**2, MSEs, tex.classes) 273 totweight = gettotweight(res) 274 275 if argkw.get("sqrt", 0): 276 MSEs = [math.sqrt(x) for x in MSEs] 277 return [x/totweight for x in MSEs] 278 279 def RMSE_old(res, **argkw): 280 """RMSE(res) > root meansquared error""" 281 argkw.setdefault("sqrt", 1) 282 return MSE_old(res, **argkw) 283 284 285 ######################################################################### 286 # PERFORMANCE MEASURES: 287 # Scores for evaluation of numeric predictions 288 289 def checkArgkw(dct, lst): 290 """checkArgkw(dct, lst) > returns true if any items have nonzero value in dct""" 291 return reduce(lambda x,y: x or y, [dct.get(k, 0) for k in lst]) 292 293 def regressionError(res, **argkw): 294 """regressionError(res) > regression error (default: MSE)""" 295 if argkw.get("SE", 0) and res.numberOfIterations > 1: 296 # computes the scores for each iteration, then averages 297 scores = [[0.0] * res.numberOfIterations for i in range(res.numberOfLearners)] 298 if argkw.get("normabs", 0) or argkw.get("normsqr", 0): 299 norm = [0.0] * res.numberOfIterations 300 301 nIter = [0]*res.numberOfIterations # counts examples in each iteration 302 a = [0]*res.numberOfIterations # average class in each iteration 303 for tex in res.results: 304 nIter[tex.iterationNumber] += 1 305 a[tex.iterationNumber] += float(tex.actualClass) 306 a = [a[i]/nIter[i] for i in range(res.numberOfIterations)] 307 308 if argkw.get("unweighted", 0) or not res.weights: 309 # iterate accross test cases 310 for tex in res.results: 311 ai = float(tex.actualClass) 312 nIter[tex.iterationNumber] += 1 313 314 # compute normalization, if required 315 if argkw.get("normabs", 0): 316 norm[tex.iterationNumber] += abs(ai  a[tex.iterationNumber]) 317 elif argkw.get("normsqr", 0): 318 norm[tex.iterationNumber] += (ai  a[tex.iterationNumber])**2 319 320 # iterate accross results of different regressors 321 for i, cls in enumerate(tex.classes): 322 if argkw.get("abs", 0): 323 scores[i][tex.iterationNumber] += abs(float(cls)  ai) 324 else: 325 scores[i][tex.iterationNumber] += (float(cls)  ai)**2 326 else: # unweighted<>0 327 raise NotImplementedError, "weighted error scores with SE not implemented yet" 328 329 if argkw.get("normabs") or argkw.get("normsqr"): 330 scores = [[x/n for x, n in zip(y, norm)] for y in scores] 331 else: 332 scores = [[x/ni for x, ni in zip(y, nIter)] for y in scores] 333 334 if argkw.get("R2"): 335 scores = [[1.0  x for x in y] for y in scores] 336 337 if argkw.get("sqrt", 0): 338 scores = [[math.sqrt(x) for x in y] for y in scores] 339 340 return [(statc.mean(x), statc.std(x)) for x in scores] 341 342 else: # single iteration (testing on a single test set) 343 scores = [0.0] * res.numberOfLearners 344 norm = 0.0 345 346 if argkw.get("unweighted", 0) or not res.weights: 347 a = sum([tex.actualClass for tex in res.results]) \ 348 / len(res.results) 349 for tex in res.results: 350 if argkw.get("abs", 0): 351 scores = map(lambda res, cls, ac = float(tex.actualClass): 352 res + abs(float(cls)  ac), scores, tex.classes) 353 else: 354 scores = map(lambda res, cls, ac = float(tex.actualClass): 355 res + (float(cls)  ac)**2, scores, tex.classes) 356 357 if argkw.get("normabs", 0): 358 norm += abs(tex.actualClass  a) 359 elif argkw.get("normsqr", 0): 360 norm += (tex.actualClass  a)**2 361 totweight = gettotsize(res) 362 else: 363 # UNFINISHED 364 for tex in res.results: 365 MSEs = map(lambda res, cls, ac = float(tex.actualClass), 366 tw = tex.weight: 367 res + tw * (float(cls)  ac)**2, MSEs, tex.classes) 368 totweight = gettotweight(res) 369 370 if argkw.get("normabs", 0) or argkw.get("normsqr", 0): 371 scores = [s/norm for s in scores] 372 else: # normalize by number of instances (or sum of weights) 373 scores = [s/totweight for s in scores] 374 375 if argkw.get("R2"): 376 scores = [1.0  s for s in scores] 377 378 if argkw.get("sqrt", 0): 379 scores = [math.sqrt(x) for x in scores] 380 381 return scores 382 383 def MSE(res, **argkw): 384 """MSE(res) > meansquared error""" 385 return regressionError(res, **argkw) 386 387 def RMSE(res, **argkw): 388 """RMSE(res) > root meansquared error""" 389 argkw.setdefault("sqrt", True) 390 return regressionError(res, **argkw) 391 392 def MAE(res, **argkw): 393 """MAE(res) > mean absolute error""" 394 argkw.setdefault("abs", True) 395 return regressionError(res, **argkw) 396 397 def RSE(res, **argkw): 398 """RSE(res) > relative squared error""" 399 argkw.setdefault("normsqr", True) 400 return regressionError(res, **argkw) 401 402 def RRSE(res, **argkw): 403 """RRSE(res) > root relative squared error""" 404 argkw.setdefault("normsqr", True) 405 argkw.setdefault("sqrt", True) 406 return regressionError(res, **argkw) 407 408 def RAE(res, **argkw): 409 """RAE(res) > relative absolute error""" 410 argkw.setdefault("abs", True) 411 argkw.setdefault("normabs", True) 412 return regressionError(res, **argkw) 413 414 def R2(res, **argkw): 415 """R2(res) > Rsquared""" 416 argkw.setdefault("normsqr", True) 417 argkw.setdefault("R2", True) 418 return regressionError(res, **argkw) 419 420 def MSE_old(res, **argkw): 421 """MSE(res) > meansquared error""" 422 if argkw.get("SE", 0) and res.numberOfIterations > 1: 423 MSEs = [[0.0] * res.numberOfIterations for i in range(res.numberOfLearners)] 424 nIter = [0]*res.numberOfIterations 425 if argkw.get("unweighted", 0) or not res.weights: 426 for tex in res.results: 427 ac = float(tex.actualClass) 428 nIter[tex.iterationNumber] += 1 429 for i, cls in enumerate(tex.classes): 430 MSEs[i][tex.iterationNumber] += (float(cls)  ac)**2 431 else: 432 raise ValueError, "weighted RMSE with SE not implemented yet" 433 MSEs = [[x/ni for x, ni in zip(y, nIter)] for y in MSEs] 434 if argkw.get("sqrt", 0): 435 MSEs = [[math.sqrt(x) for x in y] for y in MSEs] 436 return [(statc.mean(x), statc.std(x)) for x in MSEs] 437 438 else: 439 MSEs = [0.0]*res.numberOfLearners 440 if argkw.get("unweighted", 0) or not res.weights: 441 for tex in res.results: 442 MSEs = map(lambda res, cls, ac = float(tex.actualClass): 443 res + (float(cls)  ac)**2, MSEs, tex.classes) 444 totweight = gettotsize(res) 445 else: 446 for tex in res.results: 447 MSEs = map(lambda res, cls, ac = float(tex.actualClass), tw = tex.weight: 448 res + tw * (float(cls)  ac)**2, MSEs, tex.classes) 449 totweight = gettotweight(res) 450 451 if argkw.get("sqrt", 0): 452 MSEs = [math.sqrt(x) for x in MSEs] 453 return [x/totweight for x in MSEs] 454 455 def RMSE_old(res, **argkw): 456 """RMSE(res) > root meansquared error""" 457 argkw.setdefault("sqrt", 1) 458 return MSE_old(res, **argkw) 459 460 461 ######################################################################### 462 # PERFORMANCE MEASURES: 463 # Scores for evaluation of classifiers 464 465 def CA(res, reportSE = False, **argkw): 466 if res.numberOfIterations==1: 467 if type(res)==ConfusionMatrix: 468 div = nm.TP+nm.FN+nm.FP+nm.TN 469 checkNonZero(div) 470 ca = [(nm.TP+nm.TN)/div] 471 else: 472 CAs = [0.0]*res.numberOfLearners 473 if argkw.get("unweighted", 0) or not res.weights: 474 totweight = gettotsize(res) 475 for tex in res.results: 476 CAs = map(lambda res, cls: res+(cls==tex.actualClass), CAs, tex.classes) 477 else: 478 totweight = 0. 479 for tex in res.results: 480 CAs = map(lambda res, cls: res+(cls==tex.actualClass and tex.weight), CAs, tex.classes) 481 totweight += tex.weight 482 checkNonZero(totweight) 483 ca = [x/totweight for x in CAs] 484 485 if reportSE: 486 return [(x, x*(1x)/math.sqrt(totweight)) for x in ca] 487 else: 488 return ca 489 490 else: 491 CAsByFold = [[0.0]*res.numberOfIterations for i in range(res.numberOfLearners)] 492 foldN = [0.0]*res.numberOfIterations 493 494 if argkw.get("unweighted", 0) or not res.weights: 495 for tex in res.results: 496 for lrn in range(res.numberOfLearners): 497 CAsByFold[lrn][tex.iterationNumber] += (tex.classes[lrn]==tex.actualClass) 498 foldN[tex.iterationNumber] += 1 499 else: 500 for tex in res.results: 501 for lrn in range(res.numberOfLearners): 502 CAsByFold[lrn][tex.iterationNumber] += (tex.classes[lrn]==tex.actualClass) and tex.weight 503 foldN[tex.iterationNumber] += tex.weight 504 505 return statisticsByFolds(CAsByFold, foldN, reportSE, False) 506 507 508 # Obsolete, but kept for compatibility 509 def CA_se(res, **argkw): 510 return CA(res, True, **argkw) 511 512 513 def AP(res, reportSE = False, **argkw): 514 if res.numberOfIterations == 1: 515 APs=[0.0]*res.numberOfLearners 516 if argkw.get("unweighted", 0) or not res.weights: 517 for tex in res.results: 518 APs = map(lambda res, probs: res + probs[tex.actualClass], APs, tex.probabilities) 519 totweight = gettotsize(res) 520 else: 521 totweight = 0. 522 for tex in res.results: 523 APs = map(lambda res, probs: res + probs[tex.actualClass]*tex.weight, APs, tex.probabilities) 524 totweight += tex.weight 525 checkNonZero(totweight) 526 return [AP/totweight for AP in APs] 527 528 APsByFold = [[0.0]*res.numberOfLearners for i in range(res.numberOfIterations)] 529 foldN = [0.0] * res.numberOfIterations 530 if argkw.get("unweighted", 0) or not res.weights: 531 for tex in res.results: 532 APsByFold[tex.iterationNumber] = map(lambda res, probs: res + probs[tex.actualClass], APsByFold[tex.iterationNumber], tex.probabilities) 533 foldN[tex.iterationNumber] += 1 534 else: 535 for tex in res.results: 536 APsByFold[tex.iterationNumber] = map(lambda res, probs: res + probs[tex.actualClass] * tex.weight, APsByFold[tex.iterationNumber], tex.probabilities) 537 foldN[tex.iterationNumber] += tex.weight 538 539 return statisticsByFolds(APsByFold, foldN, reportSE, True) 540 541 542 def BrierScore(res, reportSE = False, **argkw): 543 """Computes Brier score""" 544 # Computes an average (over examples) of sum_x(t(x)  p(x))^2, where 545 # x is class, 546 # t(x) is 0 for 'wrong' and 1 for 'correct' class 547 # p(x) is predicted probabilty. 548 # There's a trick: since t(x) is zero for all classes but the 549 # correct one (c), we compute the sum as sum_x(p(x)^2)  2*p(c) + 1 550 # Since +1 is there for each example, it adds 1 to the average 551 # We skip the +1 inside the sum and add it just at the end of the function 552 # We take max(result, 0) to avoid 0.0000x due to rounding errors 553 554 if res.numberOfIterations == 1: 555 MSEs=[0.0]*res.numberOfLearners 556 if argkw.get("unweighted", 0) or not res.weights: 557 totweight = 0.0 558 for tex in res.results: 559 MSEs = map(lambda res, probs: 560 res + reduce(lambda s, pi: s+pi**2, probs, 0)  2*probs[tex.actualClass], MSEs, tex.probabilities) 561 totweight += tex.weight 562 else: 563 for tex in res.results: 564 MSEs = map(lambda res, probs: 565 res + tex.weight*reduce(lambda s, pi: s+pi**2, probs, 0)  2*probs[tex.actualClass], MSEs, tex.probabilities) 566 totweight = gettotweight(res) 567 checkNonZero(totweight) 568 if reportSE: 569 return [(max(x/totweight+1.0, 0), 0) for x in MSEs] ## change this, not zero!!! 570 else: 571 return [max(x/totweight+1.0, 0) for x in MSEs] 572 573 BSs = [[0.0]*res.numberOfLearners for i in range(res.numberOfIterations)] 574 foldN = [0.] * res.numberOfIterations 575 576 if argkw.get("unweighted", 0) or not res.weights: 577 for tex in res.results: 578 BSs[tex.iterationNumber] = map(lambda rr, probs: 579 rr + reduce(lambda s, pi: s+pi**2, probs, 0)  2*probs[tex.actualClass], BSs[tex.iterationNumber], tex.probabilities) 580 foldN[tex.iterationNumber] += 1 581 else: 582 for tex in res.results: 583 BSs[tex.iterationNumber] = map(lambda res, probs: 584 res + tex.weight*reduce(lambda s, pi: s+pi**2, probs, 0)  2*probs[tex.actualClass], BSs[tex.iterationNumber], tex.probabilities) 585 foldN[tex.iterationNumber] += tex.weight 586 587 stats = statisticsByFolds(BSs, foldN, reportSE, True) 588 if reportSE: 589 return [(x+1.0, y) for x, y in stats] 590 else: 591 return [x+1.0 for x in stats] 592 593 def BSS(res, **argkw): 594 return [1x/2 for x in apply(BrierScore, (res, ), argkw)] 595 596 ##def _KL_div(actualClass, predicted): 597 ## 598 ##def KL(res, **argkw): 599 ## KLs = [0.0]*res.numberOfLearners 600 ## 601 ## if argkw.get("unweighted", 0) or not res.weights: 602 ## for tex in res.results: 603 ## KLs = map(lambda res, predicted: res+KL(tex.actualClass, predicted), KLs, tex.probabilities) 604 ## totweight = gettotsize(res) 605 ## else: 606 ## for tex in res.results: 607 ## ## TEGA SE NISI! 608 ## CAs = map(lambda res, cls: res+(cls==tex.actualClass and tex.weight), CAs, tex.classes) 609 ## totweight = gettotweight(res) 610 ## 611 ## return [x/totweight for x in CAs] 612 613 614 ##def KL_se(res, **argkw): 615 ## # KullbackLeibler divergence 616 ## if res.numberOfIterations==1: 617 ## if argkw.get("unweighted", 0) or not res.weights: 618 ## totweight = gettotsize(res) 619 ## else: 620 ## totweight = gettotweight(res) 621 ## return [(x, x*(1x)/math.sqrt(totweight)) for x in apply(CA, (res,), argkw)] 622 ## else: 623 ## KLsByFold = [[0.0]*res.numberOfIterations for i in range(res.numberOfLearners)] 624 ## foldN = [0.0]*res.numberOfIterations 625 ## 626 ## if argkw.get("unweighted", 0) or not res.weights: 627 ## for tex in res.results: 628 ## for lrn in range(res.numberOfLearners): 629 ## CAsByFold[lrn][tex.iterationNumber] += 630 ## foldN[tex.iterationNumber] += 1 631 ## else: 632 ## for tex in res.results: 633 ## for lrn in range(res.numberOfLearners): 634 ## CAsByFold[lrn][tex.iterationNumber] += 635 ## foldN[tex.iterationNumber] += tex.weight 636 ## 637 ## newFolds = [] 638 ## for lrn in range(res.numberOfLearners): 639 ## newF = [] 640 ## for fold in range(res.numberOfIterations): 641 ## if foldN[fold]>0.0: 642 ## newF.append(CAsByFold[lrn][fold]/foldN[fold]) 643 ## newFolds.append(newF) 644 ## 645 ## checkNonZero(len(newFolds)) 646 ## return [(statc.mean(cas), statc.sterr(cas)) for cas in newFolds] 647 ## 648 649 def IS_ex(Pc, P): 650 "Pc aposterior probability, P aprior" 651 if (Pc>=P): 652 return log2(P)+log2(Pc) 653 else: 654 return (log2(1P)+log2(1Pc)) 655 656 def IS(res, apriori=None, reportSE = False, **argkw): 657 if not apriori: 658 apriori = classProbabilitiesFromRes(res) 659 660 if res.numberOfIterations==1: 661 ISs = [0.0]*res.numberOfLearners 662 if argkw.get("unweighted", 0) or not res.weights: 663 for tex in res.results: 664 for i in range(len(tex.probabilities)): 665 cls = tex.actualClass 666 ISs[i] += IS_ex(tex.probabilities[i][cls], apriori[cls]) 667 totweight = gettotsize(res) 668 else: 669 for tex in res.results: 670 for i in range(len(tex.probabilities)): 671 cls = tex.actualClass 672 ISs[i] += IS_ex(tex.probabilities[i][cls], apriori[cls]) * tex.weight 673 totweight = gettotweight(res) 674 if reportSE: 675 return [(IS/totweight,0) for IS in ISs] 676 else: 677 return [IS/totweight for IS in ISs] 678 679 680 ISs = [[0.0]*res.numberOfIterations for i in range(res.numberOfLearners)] 681 foldN = [0.] * res.numberOfIterations 682 683 # compute info scores for each fold 684 if argkw.get("unweighted", 0) or not res.weights: 685 for tex in res.results: 686 for i in range(len(tex.probabilities)): 687 cls = tex.actualClass 688 ISs[i][tex.iterationNumber] += IS_ex(tex.probabilities[i][cls], apriori[cls]) 689 foldN[tex.iterationNumber] += 1 690 else: 691 for tex in res.results: 692 for i in range(len(tex.probabilities)): 693 cls = tex.actualClass 694 ISs[i][tex.iterationNumber] += IS_ex(tex.probabilities[i][cls], apriori[cls]) * tex.weight 695 foldN[tex.iterationNumber] += tex.weight 696 697 return statisticsByFolds(ISs, foldN, reportSE, False) 698 699 700 def Friedman(res, statistics, **argkw): 701 sums = None 702 for ri in splitByIterations(res): 703 ranks = statc.rankdata(apply(statistics, (ri,), argkw)) 704 if sums: 705 sums = sums and [ranks[i]+sums[i] for i in range(k)] 706 else: 707 sums = ranks 708 k = len(sums) 709 N = res.numberOfIterations 710 k = len(sums) 711 T = sum([x*x for x in sums]) 712 F = 12.0 / (N*k*(k+1)) * T  3 * N * (k+1) 713 return F, statc.chisqprob(F, k1) 714 715 716 def Wilcoxon(res, statistics, **argkw): 717 res1, res2 = [], [] 718 for ri in splitByIterations(res): 719 stats = apply(statistics, (ri,), argkw) 720 if (len(stats) != 2): 721 raise TypeError, "Wilcoxon compares two classifiers, no more, no less" 722 res1.append(stats[0]) 723 res2.append(stats[1]) 724 return statc.wilcoxont(res1, res2) 725 726 def rankDifference(res, statistics, **argkw): 727 if not res.results: 728 raise TypeError, "no experiments" 729 730 k = len(res.results[0].classes) 731 if (k<2): 732 raise TypeError, "nothing to compare (less than two classifiers given)" 733 if k==2: 734 return apply(Wilcoxon, (res, statistics), argkw) 735 else: 736 return apply(Friedman, (res, statistics), argkw) 737 738 class ConfusionMatrix: 739 def __init__(self): 740 self.TP = self.FN = self.FP = self.TN = 0.0 741 742 def addTFPosNeg(self, predictedPositive, isPositive, weight = 1.0): 743 if predictedPositive: 744 if isPositive: 745 self.TP += weight 746 else: 747 self.FP += weight 748 else: 749 if isPositive: 750 self.FN += weight 751 else: 752 self.TN += weight 753 754 755 def confusionMatrices(res, classIndex=1, **argkw): 756 tfpns = [ConfusionMatrix() for i in range(res.numberOfLearners)] 757 758 if classIndex<0: 759 numberOfClasses = len(res.classValues) 760 if classIndex < 1 or numberOfClasses > 2: 761 cm = [[[0.0] * numberOfClasses for i in range(numberOfClasses)] for l in range(res.numberOfLearners)] 762 if argkw.get("unweighted", 0) or not res.weights: 763 for tex in res.results: 764 trueClass = int(tex.actualClass) 765 for li, pred in enumerate(tex.classes): 766 predClass = int(pred) 767 if predClass < numberOfClasses: 768 cm[li][trueClass][predClass] += 1 769 else: 770 for tex in enumerate(res.results): 771 trueClass = int(tex.actualClass) 772 for li, pred in tex.classes: 773 predClass = int(pred) 774 if predClass < numberOfClasses: 775 cm[li][trueClass][predClass] += tex.weight 776 return cm 777 778 elif res.baseClass>=0: 779 classIndex = res.baseClass 780 else: 781 classIndex = 1 782 783 cutoff = argkw.get("cutoff") 784 if cutoff: 785 if argkw.get("unweighted", 0) or not res.weights: 786 for lr in res.results: 787 isPositive=(lr.actualClass==classIndex) 788 for i in range(res.numberOfLearners): 789 tfpns[i].addTFPosNeg(lr.probabilities[i][classIndex]>cutoff, isPositive) 790 else: 791 for lr in res.results: 792 isPositive=(lr.actualClass==classIndex) 793 for i in range(res.numberOfLearners): 794 tfpns[i].addTFPosNeg(lr.probabilities[i][classIndex]>cutoff, isPositive, lr.weight) 795 else: 796 if argkw.get("unweighted", 0) or not res.weights: 797 for lr in res.results: 798 isPositive=(lr.actualClass==classIndex) 799 for i in range(res.numberOfLearners): 800 tfpns[i].addTFPosNeg(lr.classes[i]==classIndex, isPositive) 801 else: 802 for lr in res.results: 803 isPositive=(lr.actualClass==classIndex) 804 for i in range(res.numberOfLearners): 805 tfpns[i].addTFPosNeg(lr.classes[i]==classIndex, isPositive, lr.weight) 806 return tfpns 807 808 809 # obsolete (renamed) 810 computeConfusionMatrices = confusionMatrices 811 812 813 def confusionChiSquare(confusionMatrix): 814 dim = len(confusionMatrix) 815 rowPriors = [sum(r) for r in confusionMatrix] 816 colPriors = [sum([r[i] for r in confusionMatrix]) for i in range(dim)] 817 total = sum(rowPriors) 818 rowPriors = [r/total for r in rowPriors] 819 colPriors = [r/total for r in colPriors] 820 ss = 0 821 for ri, row in enumerate(confusionMatrix): 822 for ci, o in enumerate(row): 823 e = total * rowPriors[ri] * colPriors[ci] 824 if not e: 825 return 1, 1, 1 826 ss += (oe)**2 / e 827 df = (dim1)**2 828 return ss, df, statc.chisqprob(ss, df) 829 830 831 def sens(confm): 832 """Return sensitivity (recall rate) over the given confusion matrix.""" 833 if type(confm) == list: 834 return [sens(cm) for cm in confm] 835 else: 836 tot = confm.TP+confm.FN 837 if tot < 1e6: 838 import warnings 839 warnings.warn("Can't compute sensitivity: one or both classes have no instances") 840 return 1 841 842 return confm.TP/tot 843 844 def recall(confm): 845 """Return recall rate (sensitivity) over the given confusion matrix.""" 846 return sens(confm) 847 848 849 def spec(confm): 850 """Return specificity over the given confusion matrix.""" 851 if type(confm) == list: 852 return [spec(cm) for cm in confm] 853 else: 854 tot = confm.FP+confm.TN 855 if tot < 1e6: 856 import warnings 857 warnings.warn("Can't compute specificity: one or both classes have no instances") 858 return 1 859 return confm.TN/tot 860 861 862 def PPV(confm): 863 """Return positive predictive value (precision rate) over the given confusion matrix.""" 864 if type(confm) == list: 865 return [PPV(cm) for cm in confm] 866 else: 867 tot = confm.TP+confm.FP 868 if tot < 1e6: 869 import warnings 870 warnings.warn("Can't compute PPV: one or both classes have no instances") 871 return 1 872 return confm.TP/tot 873 874 875 def precision(confm): 876 """Return precision rate (positive predictive value) over the given confusion matrix.""" 877 return PPV(confm) 878 879 880 def NPV(confm): 881 """Return negative predictive value over the given confusion matrix.""" 882 if type(confm) == list: 883 return [NPV(cm) for cm in confm] 884 else: 885 tot = confm.FN+confm.TN 886 if tot < 1e6: 887 import warnings 888 warnings.warn("Can't compute NPV: one or both classes have no instances") 889 return 1 890 return confm.TN/tot 891 892 def F1(confm): 893 """Return F1 score (harmonic mean of precision and recall) over the given confusion matrix.""" 894 if type(confm) == list: 895 return [F1(cm) for cm in confm] 896 else: 897 p = precision(confm) 898 r = recall(confm) 899 if p + r > 0: 900 return 2. * p * r / (p + r) 901 else: 902 import warnings 903 warnings.warn("Can't compute F1: P + R is zero or not defined") 904 return 1 905 906 def Falpha(confm, alpha=1.0): 907 """Return the alphamean of precision and recall over the given confusion matrix.""" 908 if type(confm) == list: 909 return [Falpha(cm, alpha=alpha) for cm in confm] 910 else: 911 p = precision(confm) 912 r = recall(confm) 913 return (1. + alpha) * p * r / (alpha * p + r) 914 915 def MCC(confm): 916 ''' 917 Return Mattew correlation coefficient over the given confusion matrix. 918 919 MCC is calculated as follows: 920 MCC = (TP*TN  FP*FN) / sqrt( (TP+FP)*(TP+FN)*(TN+FP)*(TN+FN) ) 921 922 [1] Matthews, B.W., Comparison of the predicted and observed secondary 923 structure of T4 phage lysozyme. Biochim. Biophys. Acta 1975, 405, 442451 924 925 code by Boris Gorelik 926 ''' 927 if type(confm) == list: 928 return [MCC(cm) for cm in confm] 929 else: 930 truePositive = confm.TP 931 trueNegative = confm.TN 932 falsePositive = confm.FP 933 falseNegative = confm.FN 934 935 try: 936 r = (((truePositive * trueNegative)  (falsePositive * falseNegative))/ 937 math.sqrt( (truePositive + falsePositive) * 938 ( truePositive + falseNegative ) * 939 ( trueNegative + falsePositive ) * 940 ( trueNegative + falseNegative ) ) 941 ) 942 except ZeroDivisionError: 943 # Zero difision occurs when there is either no true positives 944 # or no true negatives i.e. the problem contains only one 945 # type of classes. 946 import warnings 947 warnings.warn("Can't compute MCC: TP or TN is zero or not defined") 948 r = None 949 950 return r 951 952 def scottsPi(confm, bIsListOfMatrices=True): 953 """Compute Scott's Pi for measuring interrater agreement for nominal data 954 955 http://en.wikipedia.org/wiki/Scott%27s_Pi 956 Scott's Pi is a statistic for measuring interrater reliability for nominal 957 raters. 958 959 @param confm: confusion matrix, or list of confusion matrices. To obtain 960 nonbinary confusion matrix, call 961 orngStat.computeConfusionMatrices and set the 962 classIndex parameter to 2. 963 @param bIsListOfMatrices: specifies whether confm is list of matrices. 964 This function needs to operate on nonbinary 965 confusion matrices, which are represented by python 966 lists, therefore one needs a way to distinguish 967 between a single matrix and list of matrices 968 """ 969 970 if bIsListOfMatrices: 971 try: 972 return [scottsPi(cm, bIsListOfMatrices=False) for cm in confm] 973 except TypeError: 974 # Nevermind the parameter, maybe this is a "conventional" binary 975 # confusion matrix and bIsListOfMatrices was specified by mistake 976 return scottsPiSingle(confm, bIsListOfMatrices=False) 977 else: 978 if isinstance(confm, ConfusionMatrix): 979 confm = numpy.array( [[confm.TP, confm.FN], 980 [confm.FP, confm.TN]], dtype=float) 981 else: 982 confm = numpy.array(confm, dtype=float) 983 984 marginalSumOfRows = numpy.sum(confm, axis=0) 985 marginalSumOfColumns = numpy.sum(confm, axis=1) 986 jointProportion = (marginalSumOfColumns + marginalSumOfRows)/ \ 987 (2.0 * numpy.sum(confm, axis=None)) 988 # In the eq. above, 2.0 is what the Wikipedia page calls 989 # the number of annotators. Here we have two annotators: 990 # the observed (true) labels (annotations) and the predicted by 991 # the learners. 992 993 prExpected = numpy.sum(jointProportion ** 2, axis=None) 994 prActual = numpy.sum(numpy.diag(confm), axis=None)/numpy.sum(confm, axis=None) 995 996 ret = (prActual  prExpected) / (1.0  prExpected) 997 return ret 998 999 def AUCWilcoxon(res, classIndex=1, **argkw): 1000 import corn 1001 useweights = res.weights and not argkw.get("unweighted", 0) 1002 problists, tots = corn.computeROCCumulative(res, classIndex, useweights) 1003 1004 results=[] 1005 1006 totPos, totNeg = tots[1], tots[0] 1007 N = totPos + totNeg 1008 for plist in problists: 1009 highPos, lowNeg = totPos, 0.0 1010 W, Q1, Q2 = 0.0, 0.0, 0.0 1011 for prob in plist: 1012 thisPos, thisNeg = prob[1][1], prob[1][0] 1013 highPos = thisPos 1014 W += thisNeg * (highPos + thisPos/2.) 1015 Q2 += thisPos * (lowNeg**2 + lowNeg*thisNeg + thisNeg**2 /3.) 1016 Q1 += thisNeg * (highPos**2 + highPos*thisPos + thisPos**2 /3.) 1017 1018 lowNeg += thisNeg 1019 1020 W /= (totPos*totNeg) 1021 Q1 /= (totNeg*totPos**2) 1022 Q2 /= (totPos*totNeg**2) 1023 1024 SE = math.sqrt( (W*(1W) + (totPos1)*(Q1W**2) + (totNeg1)*(Q2W**2)) / (totPos*totNeg) ) 1025 results.append((W, SE)) 1026 return results 1027 1028 AROC = AUCWilcoxon # for backward compatibility, AROC is obsolote 1029 1030 def compare2AUCs(res, lrn1, lrn2, classIndex=1, **argkw): 1031 import corn 1032 return corn.compare2ROCs(res, lrn1, lrn2, classIndex, res.weights and not argkw.get("unweighted")) 1033 1034 compare2AROCs = compare2AUCs # for backward compatibility, compare2AROCs is obsolote 1035 1036 1037 def computeROC(res, classIndex=1): 1038 import corn 1039 problists, tots = corn.computeROCCumulative(res, classIndex) 1040 1041 results = [] 1042 totPos, totNeg = tots[1], tots[0] 1043 1044 for plist in problists: 1045 curve=[(1., 1.)] 1046 TP, TN = totPos, 0.0 1047 FN, FP = 0., totNeg 1048 for prob in plist: 1049 thisPos, thisNeg = prob[1][1], prob[1][0] 1050 # thisPos go from TP to FN 1051 TP = thisPos 1052 FN += thisPos 1053 # thisNeg go from FP to TN 1054 TN += thisNeg 1055 FP = thisNeg 1056 1057 sens = TP/(TP+FN) 1058 spec = TN/(FP+TN) 1059 curve.append((1spec, sens)) 1060 results.append(curve) 1061 1062 return results 1063 1064 ## TC's implementation of algorithms, taken from: 1065 ## T Fawcett: ROC Graphs: Notes and Practical Considerations for Data Mining Researchers, submitted to KDD Journal. 1066 def ROCslope((P1x, P1y, P1fscore), (P2x, P2y, P2fscore)): 1067 if (P1x == P2x): 1068 return 1e300 1069 return (P1y  P2y) / (P1x  P2x) 1070 1071 def ROCaddPoint(P, R, keepConcavities=1): 1072 if keepConcavities: 1073 R.append(P) 1074 else: 1075 while (1): 1076 if len(R) < 2: 1077 R.append(P) 1078 return R 1079 else: 1080 T = R.pop() 1081 T2 = R[1] 1082 if ROCslope(T2, T) > ROCslope(T, P): 1083 R.append(T) 1084 R.append(P) 1085 return R 1086 return R 1087 1088 def TCcomputeROC(res, classIndex=1, keepConcavities=1): 1089 import corn 1090 problists, tots = corn.computeROCCumulative(res, classIndex) 1091 1092 results = [] 1093 P, N = tots[1], tots[0] 1094 1095 for plist in problists: 1096 ## corn gives an increasing by scores list, we need a decreasing by scores 1097 plist.reverse() 1098 TP = 0.0 1099 FP = 0.0 1100 curve=[] 1101 fPrev = 10e300 # "infinity" score at 0.0, 0.0 1102 for prob in plist: 1103 f = prob[0] 1104 if f <> fPrev: 1105 if P: 1106 tpr = TP/P 1107 else: 1108 tpr = 0.0 1109 if N: 1110 fpr = FP/N 1111 else: 1112 fpr = 0.0 1113 curve = ROCaddPoint((fpr, tpr, fPrev), curve, keepConcavities) 1114 fPrev = f 1115 thisPos, thisNeg = prob[1][1], prob[1][0] 1116 TP += thisPos 1117 FP += thisNeg 1118 if P: 1119 tpr = TP/P 1120 else: 1121 tpr = 0.0 1122 if N: 1123 fpr = FP/N 1124 else: 1125 fpr = 0.0 1126 curve = ROCaddPoint((fpr, tpr, f), curve, keepConcavities) ## ugly 1127 results.append(curve) 1128 1129 return results 1130 1131 ## returns a list of points at the intersection of the tangential isoperformance line and the given ROC curve 1132 ## for given values of FPcost, FNcost and pval 1133 def TCbestThresholdsOnROCcurve(FPcost, FNcost, pval, curve): 1134 m = (FPcost*(1.0  pval)) / (FNcost*pval) 1135 1136 ## put the isoperformance line in point (0.0, 1.0) 1137 x0, y0 = (0.0, 1.0) 1138 x1, y1 = (1.0, 1.0 + m) 1139 d01 = math.sqrt((x1  x0)*(x1  x0) + (y1  y0)*(y1  y0)) 1140 1141 ## calculate and find the closest point to the line 1142 firstp = 1 1143 mind = 0.0 1144 a = (x0*y1  x1*y0) 1145 closestPoints = [] 1146 for (x, y, fscore) in curve: 1147 d = ((y0  y1)*x + (x1  x0)*y + a) / d01 1148 d = abs(d) 1149 if firstp or d < mind: 1150 mind, firstp = d, 0 1151 closestPoints = [(x, y, fscore)] 1152 else: 1153 if abs(d  mind) <= 0.0001: ## close enough 1154 closestPoints.append( (x, y, fscore) ) 1155 return closestPoints 1156 1157 def frange(start, end=None, inc=None): 1158 "A range function, that does accept float increments..." 1159 1160 if end == None: 1161 end = start + 0.0 1162 start = 0.0 1163 1164 if inc == None or inc == 0: 1165 inc = 1.0 1166 1167 L = [start] 1168 while 1: 1169 next = start + len(L) * inc 1170 if inc > 0 and next >= end: 1171 L.append(end) 1172 break 1173 elif inc < 0 and next <= end: 1174 L.append(end) 1175 break 1176 L.append(next) 1177 1178 return L 1179 1180 ## input ROCcurves are of form [ROCcurves1, ROCcurves2, ... ROCcurvesN], 1181 ## where ROCcurvesX is a set of ROC curves, 1182 ## where a (one) ROC curve is a set of (FP, TP) points 1183 ## 1184 ## for each (sub)set of input ROC curves 1185 ## returns the average ROC curve and an array of (vertical) standard deviations 1186 def TCverticalAverageROC(ROCcurves, samples = 10): 1187 def INTERPOLATE((P1x, P1y, P1fscore), (P2x, P2y, P2fscore), X): 1188 if (P1x == P2x) or ((X > P1x) and (X > P2x)) or ((X < P1x) and (X < P2x)): 1189 raise ValueError, "assumptions for interpolation are not met: P1 = %f,%f P2 = %f,%f X = %f" % (P1x, P1y, P2x, P2y, X) 1190 dx = float(P2x)  float(P1x) 1191 dy = float(P2y)  float(P1y) 1192 m = dy/dx 1193 return P1y + m*(X  P1x) 1194 1195 def TP_FOR_FP(FPsample, ROC, npts): 1196 i = 0 1197 while i < npts  1: 1198 (fp, _, _) = ROC[i + 1] 1199 if (fp <= FPsample): 1200 i += 1 1201 else: 1202 break 1203 (fp, tp, _) = ROC[i] 1204 if fp == FPsample: 1205 return tp 1206 elif fp < FPsample and i + 1 < len(ROC): 1207 return INTERPOLATE(ROC[i], ROC[i+1], FPsample) 1208 elif fp < FPsample and i + 1 == len(ROC): # return the last 1209 return ROC[i][1] 1210 raise ValueError, "cannot compute: TP_FOR_FP in TCverticalAverageROC" 1211 #return 0.0 1212 1213 average = [] 1214 stdev = [] 1215 for ROCS in ROCcurves: 1216 npts = [] 1217 for c in ROCS: 1218 npts.append(len(c)) 1219 nrocs = len(ROCS) 1220 1221 TPavg = [] 1222 TPstd = [] 1223 for FPsample in frange(0.0, 1.0, 1.0/samples): 1224 TPsum = [] 1225 for i in range(nrocs): 1226 TPsum.append( TP_FOR_FP(FPsample, ROCS[i], npts[i]) ) ##TPsum = TPsum + TP_FOR_FP(FPsample, ROCS[i], npts[i]) 1227 TPavg.append( (FPsample, statc.mean(TPsum)) ) 1228 if len(TPsum) > 1: 1229 stdv = statc.std(TPsum) 1230 else: 1231 stdv = 0.0 1232 TPstd.append( stdv ) 1233 1234 average.append(TPavg) 1235 stdev.append(TPstd) 1236 1237 return (average, stdev) 1238 1239 ## input ROCcurves are of form [ROCcurves1, ROCcurves2, ... ROCcurvesN], 1240 ## where ROCcurvesX is a set of ROC curves, 1241 ## where a (one) ROC curve is a set of (FP, TP) points 1242 ## 1243 ## for each (sub)set of input ROC curves 1244 ## returns the average ROC curve, an array of vertical standard deviations and an array of horizontal standard deviations 1245 def TCthresholdlAverageROC(ROCcurves, samples = 10): 1246 def POINT_AT_THRESH(ROC, npts, thresh): 1247 i = 0 1248 while i < npts  1: 1249 (px, py, pfscore) = ROC[i] 1250 if (pfscore > thresh): 1251 i += 1 1252 else: 1253 break 1254 return ROC[i] 1255 1256 average = [] 1257 stdevV = [] 1258 stdevH = [] 1259 for ROCS in ROCcurves: 1260 npts = [] 1261 for c in ROCS: 1262 npts.append(len(c)) 1263 nrocs = len(ROCS) 1264 1265 T = [] 1266 for c in ROCS: 1267 for (px, py, pfscore) in c: 1268 ## try: 1269 ## T.index(pfscore) 1270 ## except: 1271 T.append(pfscore) 1272 T.sort() 1273 T.reverse() ## ugly 1274 1275 TPavg = [] 1276 TPstdV = [] 1277 TPstdH = [] 1278 for tidx in frange(0, (len(T)  1.0), float(len(T))/samples): 1279 FPsum = [] 1280 TPsum = [] 1281 for i in range(nrocs): 1282 (fp, tp, _) = POINT_AT_THRESH(ROCS[i], npts[i], T[int(tidx)]) 1283 FPsum.append(fp) 1284 TPsum.append(tp) 1285 TPavg.append( (statc.mean(FPsum), statc.mean(TPsum)) ) 1286 ## vertical standard deviation 1287 if len(TPsum) > 1: 1288 stdv = statc.std(TPsum) 1289 else: 1290 stdv = 0.0 1291 TPstdV.append( stdv ) 1292 ## horizontal standard deviation 1293 if len(FPsum) > 1: 1294 stdh = statc.std(FPsum) 1295 else: 1296 stdh = 0.0 1297 TPstdH.append( stdh ) 1298 1299 average.append(TPavg) 1300 stdevV.append(TPstdV) 1301 stdevH.append(TPstdH) 1302 1303 return (average, stdevV, stdevH) 1304 1305 ## Calibration Curve 1306 ## returns an array of (curve, yesClassPredictions, noClassPredictions) elements, where: 1307 ##  curve is an array of points (x, y) on the calibration curve 1308 ##  yesClassRugPoints is an array of (x, 1) points 1309 ##  noClassRugPoints is an array of (x, 0) points 1310 def computeCalibrationCurve(res, classIndex=1): 1311 import corn 1312 ## merge multiple iterations into one 1313 mres = orngTest.ExperimentResults(1, res.classifierNames, res.classValues, res.weights, classifiers=res.classifiers, loaded=res.loaded) 1314 for te in res.results: 1315 mres.results.append( te ) 1316 1317 problists, tots = corn.computeROCCumulative(mres, classIndex) 1318 1319 results = [] 1320 P, N = tots[1], tots[0] 1321 1322 bins = 10 ## divide interval between 0.0 and 1.0 into N bins 1323 1324 for plist in problists: 1325 yesClassRugPoints = [] 1326 noClassRugPoints = [] 1327 1328 yesBinsVals = [0] * bins 1329 noBinsVals = [0] * bins 1330 for (f, (thisNeg, thisPos)) in plist: 1331 yesClassRugPoints.append( (f, thisPos) ) #1.0 1332 noClassRugPoints.append( (f, thisNeg) ) #1.0 1333 1334 index = int(f * bins ) 1335 index = min(index, bins  1) ## just in case for value 1.0 1336 yesBinsVals[index] += thisPos 1337 noBinsVals[index] += thisNeg 1338 1339 curve = [] 1340 for cn in range(bins): 1341 f = float(cn * 1.0 / bins) + (1.0 / 2.0 / bins) 1342 yesVal = yesBinsVals[cn] 1343 noVal = noBinsVals[cn] 1344 allVal = yesVal + noVal 1345 if allVal == 0.0: continue 1346 y = float(yesVal)/float(allVal) 1347 curve.append( (f, y) ) 1348 1349 ## smooth the curve 1350 maxnPoints = 100 1351 if len(curve) >= 3: 1352 # loessCurve = statc.loess(curve, 3, 0.6) 1353 loessCurve = statc.loess(curve, maxnPoints, 0.5, 3) 1354 else: 1355 loessCurve = curve 1356 clen = len(loessCurve) 1357 if clen > maxnPoints: 1358 df = clen / maxnPoints 1359 if df < 1: df = 1 1360 curve = [loessCurve[i] for i in range(0, clen, df)] 1361 else: 1362 curve = loessCurve 1363 curve = [(c)[:2] for c in curve] ## remove the third value (variance of epsilon?) that suddenly appeared in the output of the statc.loess function 1364 results.append((curve, yesClassRugPoints, noClassRugPoints)) 1365 1366 return results 1367 1368 1369 ## Lift Curve 1370 ## returns an array of curve elements, where: 1371 ##  curve is an array of points ((TP+FP)/(P + N), TP/P, (th, FP/N)) on the Lift Curve 1372 def computeLiftCurve(res, classIndex=1): 1373 import corn 1374 ## merge multiple iterations into one 1375 mres = orngTest.ExperimentResults(1, res.classifierNames, res.classValues, res.weights, classifiers=res.classifiers, loaded=res.loaded) 1376 for te in res.results: 1377 mres.results.append( te ) 1378 1379 problists, tots = corn.computeROCCumulative(mres, classIndex) 1380 1381 results = [] 1382 P, N = tots[1], tots[0] 1383 for plist in problists: 1384 ## corn gives an increasing by scores list, we need a decreasing by scores 1385 plist.reverse() 1386 TP = 0.0 1387 FP = 0.0 1388 curve = [(0.0, 0.0, (10e300, 0.0))] 1389 for (f, (thisNeg, thisPos)) in plist: 1390 TP += thisPos 1391 FP += thisNeg 1392 if FP > N: 1393 import warnings 1394 warnings.warn("The sky is falling!!") 1395 curve.append( ((TP+FP)/(P + N), TP, (f, FP/(N or 1))) ) 1396 results.append(curve) 1397 1398 return P, N, results 1399 ### 1400 1401 class CDT: 1402 """ Stores number of concordant (C), discordant (D) and tied (T) pairs (used for AUC) """ 1403 def __init__(self, C=0.0, D=0.0, T=0.0): 1404 self.C, self.D, self.T = C, D, T 1405 1406 def isCDTEmpty(cdt): 1407 return cdt.C + cdt.D + cdt.T < 1e20 1408 1409 1410 def computeCDT(res, classIndex=1, **argkw): 1411 """Obsolete, don't use""" 1412 import corn 1413 if classIndex<0: 1414 if res.baseClass>=0: 1415 classIndex = res.baseClass 1416 else: 1417 classIndex = 1 1418 1419 useweights = res.weights and not argkw.get("unweighted", 0) 1420 weightByClasses = argkw.get("weightByClasses", True) 1421 1422 if (res.numberOfIterations>1): 1423 CDTs = [CDT() for i in range(res.numberOfLearners)] 1424 iterationExperiments = splitByIterations(res) 1425 for exp in iterationExperiments: 1426 expCDTs = corn.computeCDT(exp, classIndex, useweights) 1427 for i in range(len(CDTs)): 1428 CDTs[i].C += expCDTs[i].C 1429 CDTs[i].D += expCDTs[i].D 1430 CDTs[i].T += expCDTs[i].T 1431 for i in range(res.numberOfLearners): 1432 if isCDTEmpty(CDTs[0]): 1433 return corn.computeCDT(res, classIndex, useweights) 1434 1435 return CDTs 1436 else: 1437 return corn.computeCDT(res, classIndex, useweights) 1438 1439 ## THIS FUNCTION IS OBSOLETE AND ITS AVERAGING OVER FOLDS IS QUESTIONABLE 1440 ## DON'T USE IT 1441 def ROCsFromCDT(cdt, **argkw): 1442 """Obsolete, don't use""" 1443 if type(cdt) == list: 1444 return [ROCsFromCDT(c) for c in cdt] 1445 1446 C, D, T = cdt.C, cdt.D, cdt.T 1447 N = C+D+T 1448 if N < 1e6: 1449 import warnings 1450 warnings.warn("Can't compute AUC: one or both classes have no instances") 1451 return (1,)*8 1452 if N < 2: 1453 import warnings 1454 warnings.warn("Can't compute AUC: one or both classes have too few examples") 1455 1456 som = (CD)/N 1457 c = 0.5*(1+som) 1458 1459 if (C+D): 1460 res = (C/N*100, D/N*100, T/N*100, N, som, (CD)/(C+D), (CD)/(N*(N1)/2), 0.5*(1+som)) 1461 else: 1462 res = (C/N*100, D/N*100, T/N*100, N, som, 1.0, (CD)/(N*(N1)/2), 0.5*(1+som)) 1463 1464 if argkw.get("print"): 1465 print "Concordant = %5.1f Somers' D = %1.3f" % (res[0], res[4]) 1466 print "Discordant = %5.1f Gamma = %1.3f" % (res[1], res[5]>0 and res[5] or "N/A") 1467 print "Tied = %5.1f Taua = %1.3f" % (res[2], res[6]) 1468 print " %6d pairs c = %1.3f" % (res[3], res[7]) 1469 1470 return res 1471 1472 AROCFromCDT = ROCsFromCDT # for backward compatibility, AROCFromCDT is obsolote 1473 1474 1475 1476 # computes AUC using a specified 'cdtComputer' function 1477 # It tries to compute AUCs from 'ite' (examples from a single iteration) and, 1478 # if C+D+T=0, from 'all_ite' (entire test set). In the former case, the AUCs 1479 # are divided by 'divideByIfIte'. Additional flag is returned which is True in 1480 # the former case, or False in the latter. 1481 def AUC_x(cdtComputer, ite, all_ite, divideByIfIte, computerArgs): 1482 cdts = cdtComputer(*(ite, ) + computerArgs) 1483 if not isCDTEmpty(cdts[0]): 1484 return [(cdt.C+cdt.T/2)/(cdt.C+cdt.D+cdt.T)/divideByIfIte for cdt in cdts], True 1485 1486 if all_ite: 1487 cdts = cdtComputer(*(all_ite, ) + computerArgs) 1488 if not isCDTEmpty(cdts[0]): 1489 return [(cdt.C+cdt.T/2)/(cdt.C+cdt.D+cdt.T) for cdt in cdts], False 1490 1491 return False, False 1492 1493 1494 # computes AUC between classes i and j as if there we no other classes 1495 def AUC_ij(ite, classIndex1, classIndex2, useWeights = True, all_ite = None, divideByIfIte = 1.0): 1496 import corn 1497 return AUC_x(corn.computeCDTPair, ite, all_ite, divideByIfIte, (classIndex1, classIndex2, useWeights)) 1498 1499 1500 # computes AUC between class i and the other classes (treating them as the same class) 1501 def AUC_i(ite, classIndex, useWeights = True, all_ite = None, divideByIfIte = 1.0): 1502 import corn 1503 return AUC_x(corn.computeCDT, ite, all_ite, divideByIfIte, (classIndex, useWeights)) 1504 1505 1506 # computes the average AUC over folds using a "AUCcomputer" (AUC_i or AUC_ij) 1507 # it returns the sum of what is returned by the computer, unless at a certain 1508 # fold the computer has to resort to computing over all folds or even this failed; 1509 # in these cases the result is returned immediately 1510 def AUC_iterations(AUCcomputer, iterations, computerArgs): 1511 subsum_aucs = [0.] * iterations[0].numberOfLearners 1512 for ite in iterations: 1513 aucs, foldsUsed = AUCcomputer(*(ite, ) + computerArgs) 1514 if not aucs: 1515 return None 1516 if not foldsUsed: 1517 return aucs 1518 subsum_aucs = map(add, subsum_aucs, aucs) 1519 return subsum_aucs 1520 1521 1522 # AUC for binary classification problems 1523 def AUC_binary(res, useWeights = True): 1524 if res.numberOfIterations > 1: 1525 return AUC_iterations(AUC_i, splitByIterations(res), (1, useWeights, res, res.numberOfIterations)) 1526 else: 1527 return AUC_i(res, 1, useWeights)[0] 1528 1529 # AUC for multiclass problems 1530 def AUC_multi(res, useWeights = True, method = 0): 1531 numberOfClasses = len(res.classValues) 1532 1533 if res.numberOfIterations > 1: 1534 iterations = splitByIterations(res) 1535 all_ite = res 1536 else: 1537 iterations = [res] 1538 all_ite = None 1539 1540 # by pairs 1541 sum_aucs = [0.] * res.numberOfLearners 1542 usefulClassPairs = 0. 1543 1544 if method in [0, 2]: 1545 prob = classProbabilitiesFromRes(res) 1546 1547 if method <= 1: 1548 for classIndex1 in range(numberOfClasses): 1549 for classIndex2 in range(classIndex1): 1550 subsum_aucs = AUC_iterations(AUC_ij, iterations, (classIndex1, classIndex2, useWeights, all_ite, res.numberOfIterations)) 1551 if subsum_aucs: 1552 if method == 0: 1553 p_ij = prob[classIndex1] * prob[classIndex2] 1554 subsum_aucs = [x * p_ij for x in subsum_aucs] 1555 usefulClassPairs += p_ij 1556 else: 1557 usefulClassPairs += 1 1558 sum_aucs = map(add, sum_aucs, subsum_aucs) 1559 else: 1560 for classIndex in range(numberOfClasses): 1561 subsum_aucs = AUC_iterations(AUC_i, iterations, (classIndex, useWeights, all_ite, res.numberOfIterations)) 1562 if subsum_aucs: 1563 if method == 0: 1564 p_i = prob[classIndex] 1565 subsum_aucs = [x * p_i for x in subsum_aucs] 1566 usefulClassPairs += p_i 1567 else: 1568 usefulClassPairs += 1 1569 sum_aucs = map(add, sum_aucs, subsum_aucs) 1570 1571 if usefulClassPairs > 0: 1572 sum_aucs = [x/usefulClassPairs for x in sum_aucs] 1573 1574 return sum_aucs 1575 1576 1577 # Computes AUC, possibly for multiple classes (the averaging method can be specified) 1578 # Results over folds are averages; if some folds examples from one class only, the folds are merged 1579 def AUC(res, method = 0, useWeights = True): 1580 if len(res.classValues) < 2: 1581 raise ValueError("Cannot compute AUC on a singleclass problem") 1582 elif len(res.classValues) == 2: 1583 return AUC_binary(res, useWeights) 1584 else: 1585 return AUC_multi(res, useWeights, method) 1586 1587 AUC.ByWeightedPairs = 0 1588 AUC.ByPairs = 1 1589 AUC.WeightedOneAgainstAll = 2 1590 AUC.OneAgainstAll = 3 1591 1592 1593 # Computes AUC; in multivalued class problem, AUC is computed as one against all 1594 # Results over folds are averages; if some folds examples from one class only, the folds are merged 1595 def AUC_single(res, classIndex = 1, useWeights = True): 1596 if classIndex<0: 1597 if res.baseClass>=0: 1598 classIndex = res.baseClass 1599 else: 1600 classIndex = 1 1601 1602 if res.numberOfIterations > 1: 1603 return AUC_iterations(AUC_i, splitByIterations(res), (classIndex, useWeights, res, res.numberOfIterations)) 1604 else: 1605 return AUC_i( res, classIndex, useWeights)[0] 1606 1607 # Computes AUC for a pair of classes (as if there were no other classes) 1608 # Results over folds are averages; if some folds have examples from one class only, the folds are merged 1609 def AUC_pair(res, classIndex1, classIndex2, useWeights = True): 1610 if res.numberOfIterations > 1: 1611 return AUC_iterations(AUC_ij, splitByIterations(res), (classIndex1, classIndex2, useWeights, res, res.numberOfIterations)) 1612 else: 1613 return AUC_ij(res, classIndex1, classIndex2, useWeights) 1614 1615 1616 # AUC for multiclass problems 1617 def AUC_matrix(res, useWeights = True): 1618 numberOfClasses = len(res.classValues) 1619 numberOfLearners = res.numberOfLearners 1620 1621 if res.numberOfIterations > 1: 1622 iterations, all_ite = splitByIterations(res), res 1623 else: 1624 iterations, all_ite = [res], None 1625 1626 aucs = [[[] for i in range(numberOfClasses)] for i in range(numberOfLearners)] 1627 prob = classProbabilitiesFromRes(res) 1628 1629 for classIndex1 in range(numberOfClasses): 1630 for classIndex2 in range(classIndex1): 1631 pair_aucs = AUC_iterations(AUC_ij, iterations, (classIndex1, classIndex2, useWeights, all_ite, res.numberOfIterations)) 1632 if pair_aucs: 1633 for lrn in range(numberOfLearners): 1634 aucs[lrn][classIndex1].append(pair_aucs[lrn]) 1635 else: 1636 for lrn in range(numberOfLearners): 1637 aucs[lrn][classIndex1].append(1) 1638 return aucs 1639 1640 1641 def McNemar(res, **argkw): 1642 nLearners = res.numberOfLearners 1643 mcm = [] 1644 for i in range(nLearners): 1645 mcm.append([0.0]*res.numberOfLearners) 1646 1647 if not res.weights or argkw.get("unweighted"): 1648 for i in res.results: 1649 actual = i.actualClass 1650 classes = i.classes 1651 for l1 in range(nLearners): 1652 for l2 in range(l1, nLearners): 1653 if classes[l1]==actual: 1654 if classes[l2]!=actual: 1655 mcm[l1][l2] += 1 1656 elif classes[l2]==actual: 1657 mcm[l2][l1] += 1 1658 else: 1659 for i in res.results: 1660 actual = i.actualClass 1661 classes = i.classes 1662 for l1 in range(nLearners): 1663 for l2 in range(l1, nLearners): 1664 if classes[l1]==actual: 1665 if classes[l2]!=actual: 1666 mcm[l1][l2] += i.weight 1667 elif classes[l2]==actual: 1668 mcm[l2][l1] += i.weight 1669 1670 for l1 in range(nLearners): 1671 for l2 in range(l1, nLearners): 1672 su=mcm[l1][l2] + mcm[l2][l1] 1673 if su: 1674 mcm[l2][l1] = (abs(mcm[l1][l2]mcm[l2][l1])1)**2 / su 1675 else: 1676 mcm[l2][l1] = 0 1677 1678 for l1 in range(nLearners): 1679 mcm[l1]=mcm[l1][:l1] 1680 1681 return mcm 1682 1683 1684 def McNemarOfTwo(res, lrn1, lrn2): 1685 tf = ft = 0.0 1686 if not res.weights or argkw.get("unweighted"): 1687 for i in res.results: 1688 actual=i.actualClass 1689 if i.classes[lrn1]==actual: 1690 if i.classes[lrn2]!=actual: 1691 tf += i.weight 1692 elif i.classes[lrn2]==actual: 1693 ft += i.weight 1694 else: 1695 for i in res.results: 1696 actual=i.actualClass 1697 if i.classes[lrn1]==actual: 1698 if i.classes[lrn2]!=actual: 1699 tf += 1.0 1700 elif i.classes[lrn2]==actual: 1701 ft += 1.0 1702 1703 su = tf + ft 1704 if su: 1705 return (abs(tfft)1)**2 / su 1706 else: 1707 return 0 1708 1709 1710 def Friedman(res, stat=CA): 1711 """ Compares classifiers by Friedman test, treating folds as different examles. 1712 Returns F, p and average ranks 1713 """ 1714 res_split = splitByIterations(res) 1715 res = [stat(r) for r in res_split] 1716 1717 N = len(res) 1718 k = len(res[0]) 1719 sums = [0.0]*k 1720 for r in res: 1721 ranks = [kx+1 for x in statc.rankdata(r)] 1722 if stat==BrierScore: # reverse ranks for BrierScore (lower better) 1723 ranks = [k+1x for x in ranks] 1724 sums = [ranks[i]+sums[i] for i in range(k)] 1725 1726 T = sum([x*x for x in sums]) 1727 sums = [x/N for x in sums] 1728 1729 F = 12.0 / (N*k*(k+1)) * T  3 * N * (k+1) 1730 1731 return F, statc.chisqprob(F, k1), sums 1732 1733 1734 def WilcoxonPairs(res, avgranks, stat=CA): 1735 """ Returns a triangular matrix, where element[i][j] stores significance of difference 1736 between ith and jth classifier, as computed by Wilcoxon test. The element is positive 1737 if ith is better than jth, negative if it is worse, and 1 if they are equal. 1738 Arguments to function are ExperimentResults, average ranks (as returned by Friedman) 1739 and, optionally, a statistics; greater values should mean better results.append 1740 """ 1741 res_split = splitByIterations(res) 1742 res = [stat(r) for r in res_split] 1743 1744 k = len(res[0]) 1745 bt = [] 1746 for m1 in range(k): 1747 nl = [] 1748 for m2 in range(m1+1, k): 1749 t, p = statc.wilcoxont([r[m1] for r in res], [r[m2] for r in res]) 1750 if avgranks[m1]<avgranks[m2]: 1751 nl.append(p) 1752 elif avgranks[m2]<avgranks[m1]: 1753 nl.append(p) 1754 else: 1755 nl.append(1) 1756 bt.append(nl) 1757 return bt 1758 1759 1760 def plotLearningCurveLearners(file, allResults, proportions, learners, noConfidence=0): 1761 plotLearningCurve(file, allResults, proportions, [orngMisc.getobjectname(learners[i], "Learner %i" % i) for i in range(len(learners))], noConfidence) 1762 1763 def plotLearningCurve(file, allResults, proportions, legend, noConfidence=0): 1764 import types 1765 fopened=0 1766 if (type(file)==types.StringType): 1767 file=open(file, "wt") 1768 fopened=1 1769 1770 file.write("set yrange [0:1]\n") 1771 file.write("set xrange [%f:%f]\n" % (proportions[0], proportions[1])) 1772 file.write("set multiplot\n\n") 1773 CAs = [CA_dev(x) for x in allResults] 1774 1775 file.write("plot \\\n") 1776 for i in range(len(legend)1): 1777 if not noConfidence: 1778 file.write("'' title '' with yerrorbars pointtype %i,\\\n" % (i+1)) 1779 file.write("'' title '%s' with linespoints pointtype %i,\\\n" % (legend[i], i+1)) 1780 if not noConfidence: 1781 file.write("'' title '' with yerrorbars pointtype %i,\\\n" % (len(legend))) 1782 file.write("'' title '%s' with linespoints pointtype %i\n" % (legend[1], len(legend))) 1783 1784 for i in range(len(legend)): 1785 if not noConfidence: 1786 for p in range(len(proportions)): 1787 file.write("%f\t%f\t%f\n" % (proportions[p], CAs[p][i][0], 1.96*CAs[p][i][1])) 1788 file.write("e\n\n") 1789 1790 for p in range(len(proportions)): 1791 file.write("%f\t%f\n" % (proportions[p], CAs[p][i][0])) 1792 file.write("e\n\n") 1793 1794 if fopened: 1795 file.close() 1796 1797 1798 def printSingleROCCurveCoordinates(file, curve): 1799 import types 1800 fopened=0 1801 if (type(file)==types.StringType): 1802 file=open(file, "wt") 1803 fopened=1 1804 1805 for coord in curve: 1806 file.write("%5.3f\t%5.3f\n" % tuple(coord)) 1807 1808 if fopened: 1809 file.close() 1810 1811 1812 def plotROCLearners(file, curves, learners): 1813 plotROC(file, curves, [orngMisc.getobjectname(learners[i], "Learner %i" % i) for i in range(len(learners))]) 1814 1815 def plotROC(file, curves, legend): 1816 import types 1817 fopened=0 1818 if (type(file)==types.StringType): 1819 file=open(file, "wt") 1820 fopened=1 1821 1822 file.write("set yrange [0:1]\n") 1823 file.write("set xrange [0:1]\n") 1824 file.write("set multiplot\n\n") 1825 1826 file.write("plot \\\n") 1827 for leg in legend: 1828 file.write("'' title '%s' with lines,\\\n" % leg) 1829 file.write("'' title '' with lines\n") 1830 1831 for curve in curves: 1832 for coord in curve: 1833 file.write("%5.3f\t%5.3f\n" % tuple(coord)) 1834 file.write("e\n\n") 1835 1836 file.write("1.0\t1.0\n0.0\t0.0e\n\n") 1837 1838 if fopened: 1839 file.close() 1840 1841 1842 1843 def plotMcNemarCurveLearners(file, allResults, proportions, learners, reference=1): 1844 plotMcNemarCurve(file, allResults, proportions, [orngMisc.getobjectname(learners[i], "Learner %i" % i) for i in range(len(learners))], reference) 1845 1846 def plotMcNemarCurve(file, allResults, proportions, legend, reference=1): 1847 if reference<0: 1848 reference=len(legend)1 1849 1850 import types 1851 fopened=0 1852 if (type(file)==types.StringType): 1853 file=open(file, "wt") 1854 fopened=1 1855 1856 #file.write("set yrange [0:1]\n") 1857 #file.write("set xrange [%f:%f]\n" % (proportions[0], proportions[1])) 1858 file.write("set multiplot\n\n") 1859 file.write("plot \\\n") 1860 tmap=range(reference)+range(reference+1, len(legend)) 1861 for i in tmap[:1]: 1862 file.write("'' title '%s' with linespoints pointtype %i,\\\n" % (legend[i], i+1)) 1863 file.write("'' title '%s' with linespoints pointtype %i\n" % (legend[tmap[1]], tmap[1])) 1864 file.write("\n") 1865 1866 for i in tmap: 1867 for p in range(len(proportions)): 1868 file.write("%f\t%f\n" % (proportions[p], McNemarOfTwo(allResults[p], i, reference))) 1869 file.write("e\n\n") 1870 1871 if fopened: 1872 file.close() 1873 1874 defaultPointTypes=("{$\\circ$}", "{$\\diamond$}", "{$+$}", "{$\\times$}", "{$$}")+tuple([chr(x) for x in range(97, 122)]) 1875 defaultLineTypes=("\\setsolid", "\\setdashpattern <4pt, 2pt>", "\\setdashpattern <8pt, 2pt>", "\\setdashes", "\\setdots") 1876 1877 def learningCurveLearners2PiCTeX(file, allResults, proportions, **options): 1878 return apply(learningCurve2PiCTeX, (file, allResults, proportions), options) 1879 1880 def learningCurve2PiCTeX(file, allResults, proportions, **options): 1881 import types 1882 fopened=0 1883 if (type(file)==types.StringType): 1884 file=open(file, "wt") 1885 fopened=1 1886 1887 nexamples=len(allResults[0].results) 1888 CAs = [CA_dev(x) for x in allResults] 1889 1890 graphsize=float(options.get("graphsize", 10.0)) #cm 1891 difprop=proportions[1]proportions[0] 1892 ntestexamples=nexamples*proportions[1] 1893 xunit=graphsize/ntestexamples 1894 1895 yshift=float(options.get("yshift", ntestexamples/20.)) 1896 1897 pointtypes=options.get("pointtypes", defaultPointTypes) 1898 linetypes=options.get("linetypes", defaultLineTypes) 1899 1900 if options.has_key("numberedx"): 1901 numberedx=options["numberedx"] 1902 if type(numberedx)==types.IntType: 1903 if numberedx>0: 1904 numberedx=[nexamples*proportions[int(i/float(numberedx)*len(proportions))] for i in range(numberedx)]+[proportions[1]*nexamples] 1905 elif numberedx<0: 1906 numberedx = numberedx 1907 newn=[] 1908 for i in range(numberedx+1): 1909 wanted=proportions[0]+float(i)/numberedx*difprop 1910 best=(10, 0) 1911 for t in proportions: 1912 td=abs(wantedt) 1913 if td<best[0]: 1914 best=(td, t) 1915 if not best[1] in newn: 1916 newn.append(best[1]) 1917 newn.sort() 1918 numberedx=[nexamples*x for x in newn] 1919 elif type(numberedx[0])==types.FloatType: 1920 numberedx=[nexamples*x for x in numberedx] 1921 else: 1922 numberedx=[nexamples*x for x in proportions] 1923 1924 file.write("\\mbox{\n") 1925 file.write(" \\beginpicture\n") 1926 file.write(" \\setcoordinatesystem units <%10.8fcm, %5.3fcm>\n\n" % (xunit, graphsize)) 1927 file.write(" \\setplotarea x from %5.3f to %5.3f, y from 0 to 1\n" % (0, ntestexamples)) 1928 file.write(" \\axis bottom invisible\n")# label {#examples}\n") 1929 file.write(" ticks short at %s /\n" % reduce(lambda x,y:x+" "+y, ["%i"%(x*nexamples+0.5) for x in proportions])) 1930 if numberedx: 1931 file.write(" long numbered at %s /\n" % reduce(lambda x,y:x+y, ["%i " % int(x+0.5) for x in numberedx])) 1932 file.write(" /\n") 1933 file.write(" \\axis left invisible\n")# label {classification accuracy}\n") 1934 file.write(" shiftedto y=%5.3f\n" % yshift) 1935 file.write(" ticks short from 0.0 to 1.0 by 0.05\n") 1936 file.write(" long numbered from 0.0 to 1.0 by 0.25\n") 1937 file.write(" /\n") 1938 if options.has_key("default"): 1939 file.write(" \\setdashpattern<1pt, 1pt>\n") 1940 file.write(" \\plot %5.3f %5.3f %5.3f %5.3f /\n" % (0., options["default"], ntestexamples, options["default"])) 1941 1942 for i in range(len(CAs[0])): 1943 coordinates=reduce(lambda x,y:x+" "+y, ["%i %5.3f" % (proportions[p]*nexamples, CAs[p][i][0]) for p in range(len(proportions))]) 1944 if linetypes: 1945 file.write(" %s\n" % linetypes[i]) 1946 file.write(" \\plot %s /\n" % coordinates) 1947 if pointtypes: 1948 file.write(" \\multiput %s at %s /\n" % (pointtypes[i], coordinates)) 1949 1950 file.write(" \\endpicture\n") 1951 file.write("}\n") 1952 if fopened: 1953 file.close() 1954 file.close() 1955 del file 1956 1957 def legendLearners2PiCTeX(file, learners, **options): 1958 return apply(legend2PiCTeX, (file, [orngMisc.getobjectname(learners[i], "Learner %i" % i) for i in range(len(learners))]), options) 1959 1960 def legend2PiCTeX(file, legend, **options): 1961 import types 1962 fopened=0 1963 if (type(file)==types.StringType): 1964 file=open(file, "wt") 1965 fopened=1 1966 1967 pointtypes=options.get("pointtypes", defaultPointTypes) 1968 linetypes=options.get("linetypes", defaultLineTypes) 1969 1970 file.write("\\mbox{\n") 1971 file.write(" \\beginpicture\n") 1972 file.write(" \\setcoordinatesystem units <5cm, 1pt>\n\n") 1973 file.write(" \\setplotarea x from 0.000 to %5.3f, y from 0 to 12\n" % len(legend)) 1974 1975 for i in range(len(legend)): 1976 if linetypes: 1977 file.write(" %s\n" % linetypes[i]) 1978 file.write(" \\plot %5.3f 6 %5.3f 6 /\n" % (i, i+0.2)) 1979 if pointtypes: 1980 file.write(" \\put {%s} at %5.3f 6\n" % (pointtypes[i], i+0.1)) 1981 file.write(" \\put {%s} [lb] at %5.3f 0\n" % (legend[i], i+0.25)) 1982 1983 file.write(" \\endpicture\n") 1984 file.write("}\n") 1985 if fopened: 1986 file.close() 1987 file.close() 1988 del file 1989 1990 1991 def compute_friedman(avranks, N): 1992 """ 1993 Returns a tuple (friedman statistic, degrees of freedom) 1994 and (Iman statistic  Fdistribution, degrees of freedom) 1995 """ 1996 1997 k = len(avranks) 1998 1999 def friedman(N, k, ranks): 2000 return 12*N*(sum([rank**2.0 for rank in ranks])  (k*(k+1)*(k+1)/4.0) )/(k*(k+1)) 2001 2002 def iman(fried, N, k): 2003 return (N1)*fried/(N*(k1)  fried) 2004 2005 f = friedman(N, k, avranks) 2006 im = iman(f, N, k) 2007 fdistdof = (k1, (k1)*(N1)) 2008 2009 return (f, k1), (im, fdistdof) 2010 2011 def compute_CD(avranks, N, alpha="0.05", type="nemenyi"): 2012 """ 2013 if type == "nemenyi": 2014 critical difference for Nemenyi two tailed test. 2015 if type == "bonferronidunn": 2016 critical difference for BonferroniDunn test 2017 """ 2018 2019 k = len(avranks) 2020 2021 d = {} 2022 2023 d[("nemenyi", "0.05")] = [0, 0, 1.960, 2.343, 2.568, 2.728, 2.850, 2.949, 3.031, 3.102, 3.164 ] 2024 d[("nemenyi", "0.1")] = [0, 0, 1.645, 2.052, 2.291, 2.459, 2.589, 2.693, 2.780, 2.855, 2.920 ] 2025 d[("bonferronidunn", "0.05")] = [0, 0, 1.960, 2.241, 2.394, 2.498, 2.576, 2.638, 2.690, 2.724, 2.773 ] 2026 d[("bonferronidunn", "0.1")] = [0, 0, 1.645, 1.960, 2.128, 2.241, 2.326, 2.394, 2.450, 2.498, 2.539 ] 2027 2028 q = d[(type, alpha)] 2029 2030 cd = q[k]*(k*(k+1)/(6.0*N))**0.5 2031 2032 return cd 2033 2034 2035 def graph_ranks(filename, avranks, names, cd=None, cdmethod=None, lowv=None, highv=None, width=6, textspace=1, reverse=False, **kwargs): 2036 """ 2037 Draws a CD graph, which is used to display the differences in methods' 2038 performance. 2039 See Janez Demsar, Statistical Comparisons of Classifiers over 2040 Multiple Data Sets, 7(Jan):130, 2006. 2041 2042 Needs matplotlib to work. 2043 2044 Arguments: 2045 filename  Output file name (with extension). Formats supported 2046 by matplotlib can be used. 2047 avranks  List of average methods' ranks. 2048 names  List of methods' names. 2049 2050 Keyword arguments: 2051 cd  Critical difference. Used for marking methods that whose 2052 difference is not statistically significant. 2053 lowv  The lowest shown rank, if None, use 1. 2054 highv  The highest shown rank, if None, use len(avranks). 2055 width  Width of the drawn figure in inches, default 6 in. 2056 textspace  Space on figure sides left for the description 2057 of methods, default 1 in. 2058 reverse  If True, the lowest rank is on the right. Default: 2059 False. 2060 cdmethod  None by default. It can be an index of element in avranks or 2061 or names which specifies the method which should be marked 2062 with an interval. 2063 2064 Maintainer: Marko Toplak 2065 """ 2066 2067 width = float(width) 2068 textspace = float(textspace) 2069 2070 def nth(l,n): 2071 """ 2072 Returns only nth elemnt in a list. 2073 """ 2074 n = lloc(l,n) 2075 return [ a[n] for a in l ] 2076 2077 def lloc(l,n): 2078 """ 2079 List location in list of list structure. 2080 Enable the use of negative locations: 2081 1 is the last element, 2 second last... 2082 """ 2083 if n < 0: 2084 return len(l[0])+n 2085 else: 2086 return n 2087 2088 def mxrange(lr): 2089 """ 2090 Multiple xranges. Can be used to traverse matrices. 2091 This function is very slow due to unknown number of 2092 parameters. 2093 2094 >>> mxrange([3,5]) 2095 [(0, 0), (0, 1), (0, 2), (1, 0), (1, 1), (1, 2)] 2096 2097 >>> mxrange([[3,5,1],[9,0,3]]) 2098 [(3, 9), (3, 6), (3, 3), (4, 9), (4, 6), (4, 3)] 2099 2100 """ 2101 if len(lr) == 0: 2102 yield () 2103 else: 2104 #it can work with single numbers 2105 index = lr[0] 2106 if type(1) == type(index): 2107 index = [ index ] 2108 for a in range(*index): 2109 for b in mxrange(lr[1:]): 2110 yield tuple([a] + list(b)) 2111 2112 try: 2113 from matplotlib.figure import Figure 2114 from matplotlib.patches import Polygon 2115 from matplotlib.backends.backend_agg import FigureCanvasAgg 2116 except: 2117 import sys 2118 print >> sys.stderr, "Function requires matplotlib. Please install it." 2119 return 2120 2121 def printFigure(fig, *args, **kwargs): 2122 canvas = FigureCanvasAgg(fig) 2123 canvas.print_figure(*args, **kwargs) 2124 2125 sums = avranks 2126 2127 tempsort = sorted([ (a,i) for i,a in enumerate(sums) ], reverse=reverse) 2128 ssums = nth(tempsort, 0) 2129 sortidx = nth(tempsort, 1) 2130 nnames = [ names[x] for x in sortidx ] 2131 2132 if lowv == None: 2133 lowv = min(1, int(math.floor(min(ssums)))) 2134 if highv == None: 2135 highv = max(len(avranks), int(math.ceil(max(ssums)))) 2136 2137 cline = 0.4 2138 2139 k = len(sums) 2140 2141 lines = None 2142 sums = sorted(sums) 2143 2144 linesblank = 0 2145 scalewidth = width  2*textspace 2146 2147 def rankpos(rank): 2148 if not reverse: 2149 a = rank  lowv 2150 else: 2151 a = highv  rank 2152 return textspace+scalewidth/(highvlowv)*a 2153 2154 distanceh = 0.25 2155 2156 if cd and cdmethod == None: 2157 2158 #get pairs of non significant methods 2159 2160 def getLines(sums, hsd): 2161 2162 #get all pairs 2163 lsums = len(sums) 2164 allpairs = [ (i,j) for i,j in mxrange([[lsums], [lsums]]) if j > i ] 2165 #remove not significant 2166 notSig = [ (i,j) for i,j in allpairs if abs(sums[i]sums[j]) <= hsd ] 2167 #keep only longest 2168 2169 def noLonger((i,j), notSig): 2170 for i1,j1 in notSig: 2171 if (i1 <= i and j1 > j) or (i1 < i and j1 >= j): 2172 return False 2173 return True 2174 2175 longest = [ (i,j) for i,j in notSig if noLonger((i,j),notSig) ] 2176 2177 return longest 2178 2179 lines = getLines(ssums, cd) 2180 linesblank = 0.2 + 0.2 + (len(lines)1)*0.1 2181 2182 #add scale 2183 distanceh = 0.25 2184 cline += distanceh 2185 2186 #calculate height needed height of an image 2187 minnotsignificant = max(2*0.2, linesblank) 2188 height = cline + ((k+1)/2)*0.2 + minnotsignificant 2189 2190 fig = Figure(figsize=(width, height)) 2191 ax = fig.add_axes([0,0,1,1]) #reverse y axis 2192 ax.set_axis_off() 2193 2194 hf = 1./height # height factor 2195 wf = 1./width 2196 2197 def hfl(l): 2198 return [ a*hf for a in l ] 2199 2200 def wfl(l): 2201 return [ a*wf for a in l ] 2202 2203 """ 2204 Upper left corner is (0,0). 2205 """ 2206 2207 ax.plot([0,1], [0,1], c="w") 2208 ax.set_xlim(0, 1) 2209 ax.set_ylim(1, 0) 2210 2211 def line(l, color='k', **kwargs): 2212 """ 2213 Input is a list of pairs of points. 2214 """ 2215 ax.plot(wfl(nth(l,0)), hfl(nth(l,1)), color=color, **kwargs) 2216 2217 def text(x, y, s, *args, **kwargs): 2218 ax.text(wf*x, hf*y, s, *args, **kwargs) 2219 2220 line([(textspace, cline), (widthtextspace, cline)], linewidth=0.7) 2221 2222 bigtick = 0.1 2223 smalltick = 0.05 2224 2225 2226 import numpy 2227 2228 for a in list(numpy.arange(lowv, highv, 0.5)) + [highv]: 2229 tick = smalltick 2230 if a == int(a): tick = bigtick 2231 line([(rankpos(a), clinetick/2),(rankpos(a), cline)], linewidth=0.7) 2232 2233 for a in range(lowv, highv+1): 2234 text(rankpos(a), clinetick/20.05, str(a), ha="center", va="bottom") 2235 2236 k = len(ssums) 2237 2238 for i in range((k+1)/2): 2239 chei = cline+ minnotsignificant + (i)*0.2 2240 line([(rankpos(ssums[i]), cline), (rankpos(ssums[i]), chei), (textspace0.1, chei)], linewidth=0.7) 2241 text(textspace0.2, chei, nnames[i], ha="right", va="center") 2242 2243 for i in range((k+1)/2, k): 2244 chei = cline + minnotsignificant + (ki1)*0.2 2245 line([(rankpos(ssums[i]), cline), (rankpos(ssums[i]), chei), (textspace+scalewidth+0.1, chei)], linewidth=0.7) 2246 text(textspace+scalewidth+0.2, chei, nnames[i], ha="left", va="center") 2247 2248 if cd and cdmethod == None: 2249 2250 #upper scale 2251 if not reverse: 2252 begin, end = rankpos(lowv), rankpos(lowv+cd) 2253 else: 2254 begin, end = rankpos(highv), rankpos(highv  cd) 2255 2256 line([(begin, distanceh), (end, distanceh)], linewidth=0.7) 2257 line([(begin, distanceh + bigtick/2), (begin, distanceh  bigtick/2)], linewidth=0.7) 2258 line([(end, distanceh + bigtick/2), (end, distanceh  bigtick/2)], linewidth=0.7) 2259 text((begin+end)/2, distanceh  0.05, "CD", ha="center", va="bottom") 2260 2261 #non significance lines 2262 def drawLines(lines, side=0.05, height=0.1): 2263 start = cline + 0.2 2264 for l,r in lines: 2265 line([(rankpos(ssums[l])side, start), (rankpos(ssums[r])+side, start)], linewidth=2.5) 2266 start += height 2267 2268 drawLines(lines) 2269 2270 elif cd: 2271 begin = rankpos(avranks[cdmethod]cd) 2272 end = rankpos(avranks[cdmethod]+cd) 2273 line([(begin, cline), (end, cline)], linewidth=2.5) 2274 line([(begin, cline + bigtick/2), (begin, cline  bigtick/2)], linewidth=2.5) 2275 line([(end, cline + bigtick/2), (end, cline  bigtick/2)], linewidth=2.5) 2276 2277 printFigure(fig, filename, **kwargs) 2278 2279 if __name__ == "__main__": 2280 avranks = [3.143, 2.000, 2.893, 1.964] 2281 names = ["prva", "druga", "tretja", "cetrta" ] 2282 cd = compute_CD(avranks, 14) 2283 #cd = compute_CD(avranks, 10, type="bonferronidunn") 2284 print cd 2285 2286 print compute_friedman(avranks, 14) 2287 2288 #graph_ranks("test.eps", avranks, names, cd=cd, cdmethod=0, width=6, textspace=1.5) 
orange/doc/Orange/rst/index.rst
r7338 r7396 33 33 orange.projection.som 34 34 35 Orange.evaluation.scoring 36 35 37 Indices and tables 36 38 ==================
Note: See TracChangeset
for help on using the changeset viewer.