%QUADRC Quadratic Discriminant Classifier% %		W = QUADRC(A,R,S)% % INPUT% 	A 		Dataset%		R,S		0 <= R,S <= 1, regularization parameters (default: R = 0, S = 0)% % OUTPUT% 	W     Quadratic Discriminant Classifier mapping%%	DESCRIPTION  %	Computation of the quadratic classifier between the classes of the dataset%	A assuming normal densities. R and S are regularization parameters used%	for finding the covariance matrix as%%		G = (1-R-S)*G + R*diag(diag(G)) + S*mean(diag(G))*eye(size(G,1))%%	NOTE% This routine differs from QDC; instead of using the densities, it is based%	on the class covariances. The multi-class problem is solved by multiple%	two-class quadratic discriminants. It is, thereby, a quadratic equivalent%	of FISHERC.% %	SEE ALSO%	MAPPINGS, DATASETS, FISHERC, NMC, NMSC, LDC, UDC, QDC% Copyright: R.P.W. Duin, duin@ph.tn.tudelft.nl% Faculty of Applied Sciences, Delft University of Technology% P.O. Box 5046, 2600 GA Delft, The Netherlands% $Id: quadrc.m,v 1.7 2004/04/24 22:16:48 duin Exp $function w = quadrc(a,arg2,s)	prtrace(mfilename);	if (nargin < 3)		prwarning(5,'regularisation parameter S not given, assuming 0');		s = 0; 	end	if (nargin < 2)		prwarning(5,'regularisation parameter R not given, assuming 0');		arg2 = 0; 	end	if (nargin < 1) | (isempty(a))		r = arg2;		w = mapping(mfilename,{r,s});		w = setname(w,'Quadr');		return	end	islabtype(a,'crisp');	isvaldset(a,2,2); % at least 2 objects per class, 2 classes	[m,k,c] = getsize(a);	if (~isa(arg2,'mapping'))					% Second argument is not a mapping: train.		r = arg2;		if (min(classsizes(a)) < 2)			error('Classes should contain more than one vector.')		end					if (c == 2)			% 2-class case: calculate quadratic discriminant parameters.			p = getprior(a); pa = p(1); pb = p(2);			JA = findnlab(a,1); JB = findnlab(a,2);			ma = mean(a(JA,:)); mb = mean(a(JB,:));			GA = covm(a(JA,:)); GB = covm(a(JB,:));			GA = (1-r-s) * GA + r * diag(diag(GA)) + ...															s * mean(diag(GA))*eye(size(GA,1));			GB = (1-r-s) * GB + r * diag(diag(GB)) + ...															s*mean(diag(GB))*eye(size(GB,1));			DGA = det(GA); DGB = det(GB);						GA = inv(GA);			GB = inv(GB);			par1 = 2*ma*GA-2*mb*GB; par2 = GB - GA; 			% If either covariance matrix is nearly singular, substitute FISHERC.			% Otherwise construct the mapping.			if (DGA <= 0) | (DGB <= 0)				warning('Covariance matrix nearly singular, regularization needed; using FISHERC instead')				w = fisherc(a);			else				par0 = (mb*GB*mb'-ma*GA*ma') + 2*log(pa/pb) + log(DGB) -log(DGA);				w = mapping(mfilename,'trained',{par0,par1',par2},getlablist(a),k,2);				w = cnormc(w,a);				w = setname(w,'Quadr');				w = setcost(w,a);			end		else			% For C > 2 classes, recursively call this function, using MCLASSC.			pars = feval(mfilename,[],r,s);			w = mclassc(a,pars);		end	else  														% Second argument is a trained mapping: test. Note that we can only		% have a 2-class case here. W's output will be [D, -D], as the distance		% of a sample to a class is the negative distance to the other class.		v = arg2; pars = getdata(v);		d = +sum((a*pars{3}).*a,2) + a*pars{2} + ones(m,1)*pars{1};		w = setdat(a,[d, -d],v);	endreturn	function p = logdet(cova,covb)%% Numerically feasible computation of log(det(cova)/det(covb)))% (seems not to help)k = size(cova,1);sa = mean(diag(cova));sb = mean(diag(covb));deta = det(cova./sa)detb = det(covb./sb)p = k*log(sa) - k*log(sb) + log(deta) - log(detb);