baycopt.m

细胞生长结构可视化工具箱-MATLAB Toolbox1999.zip

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%     Calculate the performance as a classifier
function Gcs = BayCOpt(Gcs,test)

increase = 1;

tp=0;	%No True positives 
tn=0; 	%true negatives
fp=0;	% false positives
fn=0;	% false negatives

Act=[];

for i=1:Gcs.NoClasses
	eval(['act_' int2str(i) '=[];']);
	eval(['Gcs.pip_' int2str(i) '=[];']);
end

Gcs.class=[];
sens=[];
spec=[];
oldsens=[];
oldspec=[];

% Calculate the classification probabilities for each test sample
%disp('Calculating test sample classification probabilities');
for i=1:size(test,1)
	Ip = test(rem(i-1,size(test,1))+1,1:Gcs.n);	
	Gcs.class(i) =  test(rem(i-1,size(test,1))+1,Gcs.n+1);
	for k = 1:size(Gcs.wvis,1)
		Act(k) = exp(-((norm((Ip-Gcs.w(k,:)),Gcs.metric)^2)./(Gcs.sigmav(k)).^2));
	end
	act_total(i) = 0;
	for cindex=1:Gcs.NoClasses			% NEED TO ADJUST FOR PRIORS !!
		eval(['act_' int2str(cindex) '(i) = Act*Gcs.z' int2str(cindex) '*Gcs.Prior' int2str(cindex) ';']);
		eval(['act_total(i) = act_total(i) + act_' int2str(cindex) '(i);']);
	end
	for cindex=1:Gcs.NoClasses
		eval(['Gcs.pip_' int2str(cindex) '(i) = act_' int2str(cindex) '(i)/act_total(i);']);
	end
end

% Now find optimal threshold where sens = spec
cont = 1;
thres=0;
if Gcs.NoClasses == 2
	disp('Trying different thresholds')
	while cont == 1;
		%disp(['Threshold = ',num2str(thres)]);
		tp =0; fp=0; tn=0; fn=0;
		for i=1:size(test,1)
			if (Gcs.pip_1(i)>thres)
				if (Gcs.class(i) == 1)
					tp = tp+1;	% z2 is positive count
				else
					fp = fp+1;
				end
			else
				if (Gcs.class(i) == 2)
					tn = tn+1;
				else
					fn = fn+1;
				end
			end
		end	
		oldsens = sens;
		oldspec = spec;
		sens = tp/(tp+fn);
		spec = tn/(tn+fp);
		disp(['Sensitivity = ',num2str(sens)]);
		disp(['Specificity = ',num2str(spec)]);

		if (abs(increase)<0.0005)|( abs(sens-spec)<0.01)
			cont = 0;	% terminate search
		else
			if increase>0
				if (spec>=sens)	% Keep going until have passed equality point 
					increase = -increase/2;	% reverse direction of search and halve step size
				end
			else
				if (spec<=sens)
					increase = -increase/2;
				end
			end
			thres = thres+increase;
		end	
	end
	acc = (tp+tn)/(tp+tn+fp+fn);
	acc1 = tp/(tp+fn);
	acc2 = tn/(tn+fp);
else

	% Implement MAP classifier
	NoTrue=0;		% Total numbers across all classes
	NoFalse=0;
	for cindex=1:Gcs.NoClasses
		eval(['NoTrue' int2str(cindex) ' = 0;']);
		eval(['NoFalse' int2str(cindex) ' = 0;']);
	end
	for i=1:size(test,1)
		maxp = -1;
		maxclass=-1;
		for cindex=1:Gcs.NoClasses
			eval(['result = maxp < Gcs.pip_' int2str(cindex) '(i);']);
			if result == 1
				eval(['maxp = Gcs.pip_' int2str(cindex) '(i);']);
				maxclass = cindex;
			end
		end
		%disp(['class ' int2str(class(i)), 'prediction ' int2str(maxclass)])
		if maxclass == Gcs.class(i)
			NoTrue = NoTrue+1;
			eval(['NoTrue' int2str(maxclass) ' = NoTrue' int2str(maxclass) '+1;' ]);
		else
			NoFalse = NoFalse+1;
			eval(['NoFalse' int2str(Gcs.class(i)) ' = NoFalse' int2str(Gcs.class(i)) '+1;' ]);
		end
		
	end

	acc = NoTrue/(NoTrue+NoFalse);
	for cindex=1:Gcs.NoClasses
		eval(['acc' int2str(cindex) ' = NoTrue' int2str(cindex) '/(NoTrue' int2str(cindex) '+NoFalse' int2str(cindex) ');' ]);
	end
end
disp('Performance of Bayes Classifier');
disp(['Accuracy        ',num2str(acc)]);
%disp(['Threshold       ',num2str(thres)]);

Gcs.accv = [Gcs.accv,acc];
for cindex=1:Gcs.NoClasses
	eval(['Gcs.accv' int2str(cindex) ' = [Gcs.accv' int2str(cindex) ',acc' int2str(cindex) '];']);
end