%PARZENML Optimum smoothing parameter in Parzen density estimation.% %   H = PARZENML(A,FID)% % INPUT	%   A    input dataset%   FID  File ID to write progress to (default [], see PRPROGRESS)%% OUTPUT%   H    scalar smoothing parameter%% DESCRIPTION% Maximum likelihood estimation for the smoothing parameter H in the % Parzen denstity estimation of the data in A. A leave-one out % maximum likelihood estimation is used. %% This routine does not use class information and computes a single% smoothing parameter. It may be profitable to scale the data before% calling it. eg. WS = SCALEM(A,'variance'); A = A*WS; If desired,% remove unlabeled objects first, e.g. by SELDAT.% % SEE ALSO% DATASETS, MAPPINGS, SCALEM, SELDAT, PARZENM, PARZENDC, PRPROGRESS% 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: parzenml.m,v 1.8 2006/03/07 16:15:54 duin Exp $function h = parzenml(A,fid)	prtrace(mfilename);		if nargin < 2, fid = []; end		[m,k] = size(A);	DD= distm(+A) + diag(1e70*ones(1,m));	E = min(DD);		h1 = sqrt(max(E));    % initial estimate of h	F1 = derl(DD,E,h1,k); % derivative	prprogress(fid,'\nparzenml: ML Smoothing Parameter Optimization\n')	prprogress(fid,'  h = %5.3f   F = %8.3e \n',h1,F1);	if abs(F1) < 1e-70 		h = h1;		prwarning(4,'jump out\n');		return;	end		a1 = (F1+m*k)*h1*h1;	h2 = sqrt(a1/(m*k));  % second guess	F2 = derl(DD,E,h2,k); % derivative	prprogress(fid,'  h = %5.3f   F = %8.3e \n',h2,F2);	if (abs(F2) < 1e-70) | (abs(1e0-h1/h2) < 1e-6) 		h = h2;		prwarning(4,'jump out\n');		return	end		% find zero-point of derivative to optimize h^2	% stop if improvement is small, or h does not change significantly		alf = 1;	while abs(1e0-F2/F1) > 1e-4 & abs(1e0-h2/h1) > 1e-3 & abs(F2) > 1e-70		h3 = (h1*h1*h2*h2)*(F2-F1)/(F2*h2*h2-F1*h1*h1);		if h3 < 0 % this should not happen			h3 = sqrt((F2+m*k)*h2*h2/(m*k));		else			h3 = sqrt(h3);		end		h3 = h2 +alf*(h3-h2);		F3 = derl(DD,E,h3,k);		prprogress(fid,'  h = %5.3f   F = %8.3e \n',h3,F3);		F1 = F2; F2 = F3;		h1 = h2; h2 = h3;		alf = alf*0.99; % decrease step size	end	h = h2;	prprogress(fid,'parzenml finished')returnfunction F = derl(DD,E,h,k)	% computation of the likelihood derivative for Parzen density	% given distances D and their object minima E (for increased accuracy)	m = size(DD,1);	Y = (DD-repmat(E,m,1))/(2*h*h); % correct for minimum distance to save accuracy	IY = find(Y<20);                % take small distance only, others don't contribute	P = zeros(m,m);	P(IY) = exp(-Y(IY));	PP = sum(P,2)';	FU = repmat(realmax,1,m);	J = find(PP~=0); 	FU(J) = 1./PP(J);	FF = sum(DD.*P,2);	F = (FU*FF)./(h*h) - m*k;return