% results = bm([contingency matrix]) % % contingency matrix -> contingency table of network prediction and training targets % ie. classification x target' % % results <- structure of results with the following fields % contingencyMatrix: As provided as argument 1 to function % recall: Ratio of times class correctly predicted % to total times it was actually that class % precision: Ratio of times class correctly predicted % to total times it was predicted to be that class % weightedAverage: Ratio of correct predicted classes % to total number of cases % F: Harmonic mean % Fall: F for all classes % G: Geometric mean % Gall: G for all classes % bookmakerMatrix: Matrix of bookmaker results per class % bookmaker: Bookmaker result for all classes % % matlab: sf 16/03/03 % octave: dp 11/04/03 - also extended to weight k>2 combinations of bm, F and G correctly % octave: dp 30/12/07 - simplified and extended to markedness, correlation and detmination function [results,vectors] = bm(cm) if (size(cm,1) ~= size(cm,2)) error('Contingency matrix must be square'); else k = size(cm,1); end N = sum(sum(cm)); rprob = sum(cm) ./ N; pprob = sum(cm') ./ N; recall = diag(cm)' ./ sum(cm); precision = diag(cm')' ./ sum(cm'); wav = sum(diag(cm)) ./ N; F = (2.*precision.*recall) ./ (precision + recall); G = sqrt(precision.*recall); % Fall - use weighted harmonic mean % Fall = k / sum(ones(1,k)./F); %(assume equiprobable) Fall = 1 / sum(rprob./F); %(assume real distribution) % Gall - use weighted geometric mean % Gall = (prod(G)^(1/k)); %(assume equiprobable) Gall = prod(G.^rprob); %(assume real distribution) % mask = diag(ones(1,k)); % eye(k) for compatability % maskc = reshape(mask,k*k,1); % ind = find(maskc==0); % maskc(ind) = -1; % mask = reshape(maskc,k,k); prob = rprob(ones(k,1),:)' - (ones(k,k) - eye(k)); bmcm = cm ./ prob; bms = sum(bmcm') / N; bm = bms * pprob'; clear bmcm, prob; prob = pprob(ones(k,1),:)' - (ones(k,k) - eye(k)); mkcm = cm' ./ prob; mks = sum(mkcm') / N; mk = mks * rprob'; clear mkcn, prob; %ipprob = ones(1,k) - pprob; %irprob = ones(1,k) - rprob; %ppprob = prod(pprob); %prprob = prod(rprob); dof = (k-1)*(k-1); pow = 1-1/((k-1)*k); eav = rprob*pprob'; excm = pprob'*rprob; ll = sum(sum(cm.*log((cm+ones(k,k)*0.00001)/N./excm))); g2 = 2*ll; pg = 1-chi2cdf(g2,dof); % approx probabilty of this contingency or better %chi2 = sum(sum((abs(cm-excm*N)).^2./(excm*N))); % uncorrected %chi2 = sum(sum((abs(cm-excm*N)-0.5).^2./(excm*N))); % Yates correction chi2 = sum(sum((abs(cm-excm*N)-pow).^2./(excm*N))); % Powers correction px = 1-chi2cdf(chi2,dof); % approx probabilty of this contingency or better % Fisher exact test - point probability [prod(realmar!)prod(predmar!)]/[N!prod(cell!)] % gammaln(N+1) = ln(N!) logmarfact = sum(gammaln(N.*[pprob,rprob]+ones(1,2*k))); logcellNfact = sum(sum(gammaln(cm+ones(k,k))))+gammaln(N+1); pf1 = exp(logmarfact-logcellNfact); % Fisher hypergeometric probability of this contingency % Powers correction to Fisher exact test approximates the cdf from slope of pdf pwr = log(N)/N; % or gammaln(N)/N^2 logpwrfact = sum(gammaln(N.*[pprob,rprob]+(1+pwr).*ones(1,2*k))); pfp = exp(logpwrfact-logcellNfact); % approx probabilty of this contingency or better if (pfp>0.5) pfp=pf1; end %results.contingencyMatrix = cm; results.K = k; %1,48 results.N = N; results.degreesoffreedom = dof; results.evenar = sum(rprob)/k; %4,51 results.evengr = prod(rprob)^(1/k); results.evenhr = k/sum(1./rprob); results.evenap = sum(pprob)/k; results.evengp = prod(pprob)^(1/k); results.evenhp = k/sum(1./pprob); results.rratio = results.evenhr*results.evenar/results.evengr^2; %10,57 results.pratio = results.evenhp*results.evenap/results.evengp^2; results.gratio = sqrt(results.pratio*results.rratio); results.williamsr = 1+((k/results.evenhr-1))/(6*N*dof); results.williamsp = 1+((k/results.evenhp-1))/(6*N*dof); results.williamsg = 1+((k/results.evenhp-1)*(k/results.evenhr-1))^0.5/(6*N*dof); results.detroot = abs(det(cm))^(1/k)/N; %16,63 results.randAverage = wav; results.Fall = Fall; results.Gall = Gall; results.kappa = (wav - eav)/(1 - eav); %20,67 results.bookmaker = bm; %21,- can be negative results.markedness = mk; %22,- can be negative results.correlation = sqrt(abs(mk*bm)); %23,68 bm and mk can have opposite sign results.dicorrelation = (results.detroot/results.evengr)^2; %24,69 results.dmcorrelation = (results.detroot/results.evengp)^2; results.cramervgcorr = sqrt(g2/N/(k-1)); results.cramervxcorr = sqrt(chi2/N/(k-1)); results.determination = (mk*bm); %28,73 results.kdetermination = results.kappa^2; results.bdetermination = bm^2; results.mdetermination = mk^2; results.cdeterminationa = (results.detroot^2/results.evenap/results.evenar)^2; results.cdeterminationg = (results.detroot^2/results.evengp/results.evengr)^2; results.cdeterminationh = (results.detroot^2/results.evenhp/results.evenhr)^2; results.cramervgsqdet = results.cramervgcorr^2; results.cramervxsqdet = results.cramervxcorr^2; if(results.determination==0) results.detratioa = 0; %37,82 results.detratiog = 0; results.detratioh = 0; else results.detratioa = results.cdeterminationa/results.determination; %37,82 results.detratiog = results.cdeterminationg/results.determination; results.detratioh = results.cdeterminationh/results.determination; end results.loglikelihood = ll; %40,85 results.gsq = g2; results.mutualinformation = ll/N; results.pearsonxsq = chi2; results.pg = pg; %44,89 results.px = px; results.fisherpf1 = pf1; results.fisherpfp = pfp; %47,92 vectors.recall = recall; vectors.precision = precision; vectors.F = F; vectors.G = G; vectors.bookmakerSum = bms; vectors.prevalence = rprob; vectors.mardednessSum = mks; vectors.bias = pprob;