%LDC Linear Bayes Normal Classifier (BayesNormal_1)
%
% [W.R,S,M] = LDC(A,R,S,M)
% W = A*LDC([],R,S,M);
%
% INPUT
% A Dataset
% R,S Regularization parameters, 0 <= R,S <= 1
% (optional; default: no regularization, i.e. R,S = 0)
% M Dimension of subspace structure in covariance matrix (default: K,
% all dimensions)
%
% OUTPUT
% W Linear Bayes Normal Classifier mapping
% R Value of regularization parameter R as used
% S Value of regularization parameter S as used
% M Value of regularization parameter M as used
%
% DESCRIPTION
% Computation of the linear classifier between the classes of the dataset A
% by assuming normal densities with equal covariance matrices. The joint
% covariance matrix is the weighted (by a priori probabilities) average of
% the class covariance matrices. R and S (0 <= R,S <= 1) are regularization
% parameters used for finding the covariance matrix G by:
%
% G = (1-R-S)*G + R*diag(diag(G)) + S*mean(diag(G))*eye(size(G,1))
%
% This covariance matrix is then decomposed as G = W*W' + sigma^2 * eye(K),
% where W is a K x M matrix containing the M leading principal components
% and sigma^2 is the mean of the K-M smallest eigenvalues.
%
% The use of soft labels is supported. The classification A*W is computed
% by NORMAL_MAP.
%
% If R, S or M is NaN the regularisation parameter is optimised by REGOPTC.
% The best result are usually obtained by R = 0, S = NaN, M = [], or by
% R = 0, S = 0, M = NaN (which is for problems of moderate or low dimensionality
% faster). If no regularisation is supplied a pseudo-inverse of the
% covariance matrix is used in case it is close to singular.
%
% Note that A*(KLMS([],N)*NMC) performs a similar operation by first
% pre-whitening the data in an N-dimensional space, followed by the
% nearest mean classifier. The regularization controlled by N is different
% from the above in LDC as it entirely removes small variance directions.
%
% To some extend LDC is also similar to FISHERC.
%
% EXAMPLES
% See PREX_PLOTC.
%
% REFERENCES
% 1. R.O. Duda, P.E. Hart, and D.G. Stork, Pattern classification, 2nd edition,
% John Wiley and Sons, New York, 2001.
% 2. A. Webb, Statistical Pattern Recognition, John Wiley & Sons, New York, 2002.
% 3. C. Liu and H. Wechsler, Robust Coding Schemes for Indexing and Retrieval
% from Large Face Databases, IEEE Transactions on Image Processing, vol. 9,
% no. 1, 2000, 132-136.
%
% SEE ALSO
% MAPPINGS, DATASETS, REGOPTC, NMC, NMSC, LDC, UDC, QUADRC, NORMAL_MAP, FISHERC
% Copyright: R.P.W. Duin, r.p.w.duin@prtools.org
% Faculty EWI, Delft University of Technology
% P.O. Box 5031, 2600 GA Delft, The Netherlands
% $Id: ldc.m,v 1.9 2008/01/25 10:16:23 duin Exp $
function [W,r,s,dim] = weighted_ldc(a, r, s, dim, bUseWeight_Vector)
prtrace(mfilename);
if (nargin < 5) || isempty(bUseWeight_Vector)
prwarning(4,'Weighting parameter W not given, assuming equal weights');
bUseWeight_Vector = false;
end
if (nargin < 4) || isempty(dim)
prwarning(4,'subspace dimensionality M not given, assuming K');
dim = [];
end
if (nargin < 3) || isempty(s)
prwarning(4,'Regularisation parameter S not given, assuming 0.');
s = 0;
end
if (nargin < 2) || isempty(r)
prwarning(4,'Regularisation parameter R not given, assuming 0.');
r = 0;
end
if (nargin < 1) | (isempty(a)) % No input arguments:
W = prmapping(mfilename,{r,s,dim,bUseWeight_Vector}); % return an untrained mapping.
elseif any(isnan([r,s,dim])) % optimize regularisation parameters
defs = {0,0,[]};
parmin_max = [1e-8,9.9999e-1;1e-8,9.9999e-1;1,size(a,2)];
[W,r,s,dim] = regoptc(a,mfilename,{r,s,dim,bUseWeight_Vector},defs,[3 2 1],parmin_max,testc([],'soft'),[1 1 0]);
else % training
islabtype(a,'crisp','soft');
isvaldfile(a,2,2); % at least 2 object per class, 2 classes
[m,k,c] = getsize(a);
X = (+a)';
% not_nan_idx = find(all(~isnan(X)));
% m = length(not_nan_idx);
if (nargin < 2) | isempty(bUseWeight_Vector)
prwarning(4,'Weighting parameter W not given, assuming equal weights');
weight_matrix = ones(c,m);
bUseWeight_Vector = false;
else
if( bUseWeight_Vector )
feature_vector_weight = getident(a, 'feature_vector_weight');
end
weight_matrix = zeros(c,m);
for ii = 1:c
class_indexes = findnlab(a,ii);
if( bUseWeight_Vector )
weight_matrix(ii,class_indexes) = feature_vector_weight(class_indexes);
else
weight_matrix(ii,class_indexes) = 1;
end
end
if(any(isnan(weight_matrix)))
error('Something VERY wrong with W');
end
end
% Averiguo los features direccionales.
cFeaturesDomain = getfeatdom(a);
DirectionalFeatures = [];
for ii = 1:length(cFeaturesDomain)
if(~isempty(cFeaturesDomain{ii}))
DirectionalFeatures = [DirectionalFeatures ii];
end
end
Allfeatures_idx = 1:k;
if( ~isempty(DirectionalFeatures) )
OtherFeatures = setdiff(Allfeatures_idx,DirectionalFeatures);
else
OtherFeatures = Allfeatures_idx;
end
U = nan(k,c);
% if( ~isempty(OtherFeatures) )
% U(OtherFeatures,:)=X(OtherFeatures,not_nan_idx)*weight_matrix(:,not_nan_idx)'/K;
for ii = rowvec(OtherFeatures)
not_nan_idx = find(~isnan(X(ii,:)));
K=diag(sum(weight_matrix(:,not_nan_idx),2));
U(ii,:)=X(ii,not_nan_idx)*weight_matrix(:,not_nan_idx)'/K;
end
% end
% if( ~isempty(DirectionalFeatures) )
% U(DirectionalFeatures,:)=angle(exp(1i*X(DirectionalFeatures,not_nan_idx))*weight_matrix(:,not_nan_idx)'/K);
for ii = rowvec(DirectionalFeatures)
not_nan_idx = find(~isnan(X(ii,:)));
K=diag(sum(weight_matrix(:,not_nan_idx),2));
U(ii,:)=angle(exp(1i*X(ii,not_nan_idx))*weight_matrix(:,not_nan_idx)'/K);
end
% end
% U = nan(k,c);
% for ii = 1:c
% class_indexes = findnlab(a,ii);
% U(:,ii)=median(X(:,class_indexes),2);
% end
iAux = nan(k, m);
not_nan_idx = 1:m;
weight_matrix_sqrt = sqrt(weight_matrix(:,not_nan_idx));
K=diag(sum(weight_matrix,2));
if( ~isempty(OtherFeatures) )
iAux(OtherFeatures,:) = ( ( X(OtherFeatures,not_nan_idx).*(ones(length(OtherFeatures),c)*weight_matrix_sqrt) ) - U(OtherFeatures,:)*weight_matrix_sqrt );
end
if( ~isempty(DirectionalFeatures) )
for ii = DirectionalFeatures
iAux1 = ( X(ii,:).*(ones(1,c)*weight_matrix_sqrt) );
iAux2 = U(ii,:)*weight_matrix_sqrt;
iAux3 = [ iAux1 - iAux2; 2*pi + iAux1 - iAux2; iAux1 - 2*pi - iAux2 ];
[dummy, iMinIndex] = min(abs(iAux3), [], 1);
iMinIndex = sub2ind(size(iAux3),iMinIndex,1:m);
iAux(ii,:) = iAux3(iMinIndex);
end
end
% G = iAux * iAux' / sum(sum(K));
G = zeros(k,k);
for ii = 1:k
for jj = 1:k
aux_val = iAux(ii,:) .* iAux(jj,:);
bAux = ~isnan(aux_val);
if( ~any(bAux) ); continue; end;
K = sum(sum(weight_matrix(:,bAux)));
G(ii,jj) = sum(aux_val(bAux)) / K;
end
end
%Cuando no separaba features direccionales de no direccionales.
% G=X*( ( X.*(ones(k,c)*weight_matrix) ) - U*weight_matrix )' / sum(sum(K));
% Regularize
if (s > 0) || (r > 0)
G = (1-r-s)*G + r * diag(diag(G)) + s*mean(diag(G))*eye(size(G,1));
end
if (dim < k)
dim = min(rank(G)-1,dim);
[eigvec,eigval] = eig(G); eigval = diag(eigval);
[dummy,ind] = sort(-eigval);
% Estimate sigma^2 as avg. eigenvalue outside subspace.
sigma2 = mean(eigval(ind(dim+1:end)));
% Subspace basis: first DIM eigenvectors * sqrt(eigenvalues).
G = eigvec(:,ind(1:dim)) * diag(eigval(ind(1:dim))) * eigvec(:,ind(1:dim))' + ...
sigma2 * eye(k);
end
w.mean = U';
w.prior = getprior(a);
w.cov = G;
W = normal_map_new(w,getlablist(a),k,c);
W = setcost(W,a);
W = setuser(W,cFeaturesDomain);
end
W = setname(W,'Bayes-Normal-1');
return