roc.m

一个关于数据聚类和模式识别的程序,在生物化学,化学中因该都可以用到.希望对大家有用,谢谢支持

function [threshold, fp, fn, mu1, var1, mu2, var2, a, b]=roc(data1, data2, plotOpt);
%roc
%	Usage: [threshold, fp, fn, mu1, var1, mu2, var2, a, b]=roc(data1, data2, plotOpt);
%		data1: vector for data of negative set
%		data2: vector for data of positive set
%		threshold: threhold (using Baysian rule where priors are multiplied)
%		fp: false positive error rate (using Baysian rule where priors are multiplied)
%		fn: false negative error rate (using Baysian rule where priors are multiplied)
%		mu1: mu of class 1
%		var1: variance of class 1
%		mu2: mu of class 2
%		var2: variance of class 2
%		a, b: the fitting parameters of y=1/(1+exp(-a*(x-b))) for class-2 conditional probability (priors are multiplied)
%
%	For example:
%		data1=randn(1,100)/2;
%		data2=randn(1,800)+3;
%		roc(data1, data2, 1);

%	Roger Jang, 20050717

if nargin<1; selfdemo; return; end
if nargin<3, plotOpt=0; end

data1=data1(:)';
data2=data2(:)';
gaussianParam1=gaussianMle(data1);
gaussianParam2=gaussianMle(data2);
n1=length(data1);
n2=length(data2);

x=linspace(min([data1, data2]), max([data1, data2]));
g1=gaussian(x, gaussianParam1);
g2=gaussian(x, gaussianParam2);

absDiff1=abs(n1*g1-n2*g2);

%temp=absDiff1;
%[junk, id1]=max(n1*g1);
%[junk, id2]=max(n2*g2);
%temp=inf*absDiff1;
%temp(id1:id2)=absDiff1(id1:id2);
%[junk, minIndex]=min(temp);	% The threshold must fall between x(id1) and x(id2)
%threshold=x(minIndex);
%if mu1<threshold && threshold < mu2
%	fp=sum(data1>=threshold)/(n1+n2);
%	fn=sum(data2<threshold)/(n1+n2);
%elseif mu2<threshold && threshold < mu1
%	fp=sum(data1<=threshold)/(n1+n2);
%	fn=sum(data2>threshold)/(n1+n2);
%else
%	fp=inf;
%	fn=inf;
%	fprintf('Something wrong in roc.m!\n');
%end

thresholds=x;
%allData=[data1(:); data2(:)]';
%sorted=sort(allData);
%thresholds=([sorted(1)-1, sorted]+[sorted, sorted(end)+1])/2;

fps=0*thresholds;
fns=0*thresholds;
for i=1:length(thresholds);
	th=thresholds(i);
	fps(i)=sum(data1>=th)/(n1+n2);
	fns(i)=sum(data2< th)/(n1+n2);
end
[minValue, minIndex]=min(fps+fns);
index=find((fps+fns)==minValue);	 % Check if there is more than one minValue

if length(index)>1	% If there is more than one minValue (usually a plateau of perfect separation), taking weighted average
	start=thresholds(index(1));
	stop=thresholds(index(end));
	w1=sum(abs(start-data1));
	w2=sum(abs(data2-stop));
	minIndex=round((w1*index(end)+w2*index(1))/(w1+w2));
end

fp=fps(minIndex);
fn=fns(minIndex);
threshold=thresholds(minIndex);

absDiff2=abs(n1*g1./(n1*g1+n2*g2)-n2*g2./(n1*g1+n2*g2));

% 璸衡だ计Ρ絬