# demand problem with different credal set ############################################################################### pmfs = c( 0.6, 0.4, 0.65, 0.35) rvars = c( 440, 260, 420, 300, 370, 370) getexpectations = getexpectationsfunc(2, pmfs) getexpectations(rvars) getlowerprevisions = getlowerprevisionsfunc(getexpectations) getupperprevisions = getupperprevisionsfunc(getexpectations) isgammamaximin = isgammamaxisomethingfunc(getlowerprevisions) isgammamaximax = isgammamaxisomethingfunc(getupperprevisions) isintervalmaximal = ismaximalfunc(getexpectations, intervalcompare) isrbayesmaximal = ismaximalfunc(getexpectations, rbayescompare) isrbayesadmissible = isrbayesadmissiblefunc(getexpectations) isgammamaximin(rvars) # overtime isgammamaximax(rvars) # overtime isintervalmaximal(rvars) # machinery or overtime isrbayesmaximal(rvars) # overtime isrbayesadmissible(rvars) # overtime # investment problem ############################################################################### pmfs = c( 0.0, 0.6, 0.4, 0.3, 0.3, 0.4) rvars = c( 100, 50, -25, 75, 50, 0, 60, 55, 10, 35, 35, 35) getexpectations = getexpectationsfunc(3, pmfs) getexpectations(rvars) getlowerprevisions = getlowerprevisionsfunc(getexpectations) getupperprevisions = getupperprevisionsfunc(getexpectations) isgammamaximin = isgammamaxisomethingfunc(getlowerprevisions) isgammamaximax = isgammamaxisomethingfunc(getupperprevisions) isintervalmaximal = ismaximalfunc(getexpectations, intervalcompare) isrbayesmaximal = ismaximalfunc(getexpectations, rbayescompare) isrbayesadmissible = isrbayesadmissiblefunc(getexpectations) isgammamaximin(rvars) # 3 isgammamaximax(rvars) # 3 isintervalmaximal(rvars) # 2 or 3 isrbayesmaximal(rvars) # 3 isrbayesadmissible(rvars) # 3 # environmental problem ############################################################################### # assuming precise expectation: option 2 has highest expected suitability # beta=1: as above # beta=0: maximin = 1, maximax = 1, interval dominance=1,2,3, maximality=1,2,3 # beta=0.5: # option | lower prev | upper prev # 1 | 5.7 | 7.8 # 2 | 7.9 | 10.3 # 3 | 6.8 | 8.6 # so for interval dominance = any option # maximality is a bit more tricky lowprev = function(beta) function(a0, a1, a2, a3) a0 + beta * (a1 * 6 + a2 * 10 + a3 * 8) + (1-beta) * (min(a1 * 5, a1 * 12) + min(a2 * 3, a2 * 11) + min(a3 * 4, a3 * 10)) # confirm lower and upper previsions above lowprev(0.7)(0, 1, 0, 0) lowprev(0.7)(0, 0, 1, 0) lowprev(0.7)(0, 0, 0, 1) -lowprev(0.7)(0, -1, 0, 0) -lowprev(0.7)(0, 0, -1, 0) -lowprev(0.7)(0, 0, 0, -1) # maximality lowprev(0.7)(0, -1, 1, 0) # only positive one! lowprev(0.7)(0, -1, 0, 1) lowprev(0.7)(0, 1, -1, 0) lowprev(0.7)(0, 0, -1, 1) lowprev(0.7)(0, 1, 0, -1) lowprev(0.7)(0, 0, 1, -1) # so 2 dominates 1, hence 1 is not maximal, but 2 and 3 are # more cleverly: because we know 1 is interval dominated, it cannot be # maximal, so we really only need to compare 2 and 3 # credal naive Bayes classifier based on maximality ############################################################################### source("classification.r") classifier.credal2 = function(s=2) { list( trainer=function(train, attribs, class) { list( attribs=attribs, class=class, c.count=table(train[class]), a.c.count=lapply( attribs, function(attrib) { table(train[c(class, attrib)]) }) ) }, tester=function(model, testrow) { dominates = function(c1, c2) { helperfunc = function(x) { k = length(model$attribs) rc = (model$c.count[c2] + x) / (model$c.count[c1] + s - x) ras = sapply( 1:k, function(i) { colid = model$attribs[i] a.level = testrow[1,colid] model$a.c.count[[i]][c1,a.level] / (model$a.c.count[[i]][c2,a.level] + x) }) (rc ** (k - 1)) * prod(ras) } val = optim(s/2, helperfunc, lower=0, upper=s, method="L-BFGS-B") val$value > 1 } c.levels = names(model$c.count) result = Filter( function(c1) !any(sapply(c.levels, function(c2) { dominates(c2, c1) })), c.levels) accuracy = (testrow[1, model$class] %in% result) list( acc=accuracy, deter=(length(result) == 1), singleacc=if (length(result) == 1) accuracy else NA, setsize=if (length(result) != 1) length(result) else NA, setacc=if (length(result) != 1) accuracy else NA) } ) } test.credal2 = function() { # set up the classifier myclassifier = classifier.composed( list( classifier.credal(2), # interval dominance classifier.credal2(2))) # maximality # test with just first 30 observations to get some interesting effects mammo = getdata()[1:30,] # predict severity (column 6) from all attributes print(kfcv.classifier(mammo, 1:5, 6, myclassifier)) # predict severity from all attributes except BIRADS (column 1) print(kfcv.classifier(mammo, 2:5, 6, myclassifier)) # predict severity from BIRADS only print(kfcv.classifier(mammo, 1, 6, myclassifier)) # test with full data mammo = getdata() # predict BIRADS from other attributes print(kfcv.classifier(mammo, 2:5, 1, myclassifier)) } test.credal2()