%PARZENC Optimisation of the Parzen classifier
%
% [W,H] = PARZENC(A,H)
% [W,H] = A*PARZENC([],H)
% [W,H] = A*PARZENC(H)
%
% INPUT
% A dataset
% H smoothing parameter (may be scalar, vector of per-class
% parameters, or matrix with parameters for each class (rows) and
% dimension (columns))
%
% OUTPUT
% W trained mapping
% H estimated smoothing (scalar value)
%
% DESCRIPTION
% Computation of the optimum smoothing parameter H for the Parzen
% classifier between the classes in the dataset A. The leave-one-out
% Lissack & Fu estimate is used for the classification error E. The
% final classifier is stored as a mapping in W. It may be converted
% into a classifier by W*CLASSC. PARZENC cannot be used for density
% estimation.
%
% In case smoothing H is specified, no learning is performed, just the
% discriminant W is produced for the given smoothing parameters H.
% Smoothing parameters may be scalar, vector of per-class parameters, or
% a matrix with individual smoothing for each class (rows) and feature
% directions (columns)
%
% EXAMPLES
% See PREX_DENSITY for densities and PREX_PARZEN for differences between
% PARZENC, PARZENDC and PARZENM.
%
% REFERENCES
% T. Lissack and K.S. Fu, Error estimation in pattern recognition via
% L-distance between posterior density functions, IEEE Trans. Inform.
% Theory, vol. 22, pp. 34-45, 1976.
%
% SEE ALSO (<a href="http://37steps.com/prtools">PRTools Guide</a>)
% DATASETS, MAPPINGS, PARZEN_MAP, PARZENML, PARZENDC, PARZENM, CLASSC,
% PREX_PARZEN
% Copyright: R.P.W. Duin, r.p.w.duin@37steps.com
% Faculty EWI, Delft University of Technology
% P.O. Box 5031, 2600 GA Delft, The Netherlands
% $Id: parzenc.m,v 1.6 2008/07/03 09:11:44 duin Exp $
function [W,h] = parzenc(varargin)
mapname = 'ParzenC';
argin = shiftargin(varargin,'scalar');
argin = setdefaults(argin,[],[]);
if mapping_task(argin,'definition')
W = define_mapping(argin,'untrained',mapname);
elseif mapping_task(argin,'training') % Train a mapping.
[a,h] = deal(argin{:});
islabtype(a,'crisp','soft');
isvaldfile(a,2,2); % at least 2 objects per class, 2 classes
a = testdatasize(a);
a = testdatasize(a,'objects');
[m,k,c] = getsize(a);
nlab = getnlab(a);
if ~isempty(h) % take user setting for smoothing parameter
if size(h,1) == 1, h = repmat(h,c,1); end
if size(h,2) == 1, h = repmat(h,1,k); end
if any(size(h) ~= [c,k])
error('Array with smoothing parameters has wrong size');
end
W = prmapping('parzen_map','trained',{a,h},getlablist(a),k,c);
W = setname(W,'Parzen Classifier');
return
end
% compute all object distances
% make diagonal inf to exclude objects own contribution
D = +distm(a) + diag(inf*ones(1,m));
% find object frequencies
if islabtype(a,'crisp')
csize = classsizes(a);
of = csize(nlab);
else
csize = sum(gettargets(a),1);
end
% find object weights q
p = getprior(a);
a = setprior(a,p);
q = p(nlab)./csize(nlab);
% initialise
h = max(std(a)); % for sure a too high value
L = -inf;
Ln = 0;
z = 0.1^(1/k); % initial step size
% iterate
iter = 0;
prwaitbar(100,'parzenc: Optimizing smoothing parameter',m > 100);
while abs(Ln-L) > 0.001 & z < 1
% In L we store the best performance estimate found so far.
% Ln is the actual performance (for the actual h)
% If Ln > L we improve the bound L, and so we rest it.
if Ln > L, L = Ln; end
iter = iter+1;
prwaitbar(100,100-100*exp(-iter/10));
r = -0.5/(h^2);
F = q(ones(1,m),:)'.*exp(D*r); % density contributions
FS = sum(F)*((m-1)/m); IFS = find(FS>0); % joint density distribution
if islabtype(a,'crisp');
G = sum(F .* (nlab(:,ones(1,m)) == nlab(:,ones(1,m))'));
G = G.*(of-1)./of; % true-class densities
else
% here we are for soft labels (stored in targets)
G = zeros(1,m);
for j=1:c
G = G + sum(F .* (a.targets(:,j) * a.targets(:,j)'));
% to be corrected for bias?
end
end
% performance estimate
en = max(p)*ones(1,m);
en(IFS) = (G(IFS))./FS(IFS);
Ln = exp(sum(log(en))/m);
if Ln < L % compute next estimate
z = sqrt(z); % adjust stepsize up (recall: 0 < z < 1)
h = h / z; % if we don't improve, increase h (approach h_opt from below)
else
h = h * z; % if we improve, decrease h (approach h_opt from above)
end
end
prwaitbar(0);
W = prmapping('parzen_map','trained',{a,repmat(h,c,k);},getlablist(a),k,c);
W = setname(W,mapname);
W = setcost(W,a);
end
return