function model = rsde(X,options)
% RSDE Reduced Set Density Estimator.
%
% Synopsis:
% model = rsde(X,options)
%
% Description:
% This function implements the Reduced Set Density Estimator
% [Girol03] which provides kernel density estimate optimal
% in the L2 sense. The density is modeled as the weighted sum
% of Gaussians (RBF kernel) centered in seleceted subset of
% training data.
%
% Input:
% X [dim x num_data] Input data sample.
% options [struct] Control parameters:
% .arg [1x1] Standard deviation of the Gaussian kernel.
% .solver [string] Quadratic Programming solver (default 'quadprog');
%
% Output:
% model [struct] Output density model:
% .Alpha [nsv x 1] Weights of the kernel functions.
% .sv.X [dim x nsv] Selected centers of kernel functions.
% .nsv [1x1] Number of selected centers.
% .options.arg = options.arg.
% .options.ker = 'rbf'
%
% Example:
% gnd = struct('Mean',[-2 3],'Cov',[1 0.5],'Prior',[0.4 0.6]);
% sample = gmmsamp( gnd, 200 );
% figure; hold on; ppatterns(sample.X);
% plot([-4:0.1:8], pdfgmm([-4:0.1:8],gnd),'r');
% model = rsde(sample.X,struct('arg',0.7));
% x = linspace(-4,8,100);
% plot(x,kernelproj(x,model),'g');
% ppatterns(model.sv.X,'ob',13);
%
% See also
% KERNELPROJ.
%
% About: Statistical Pattern Recognition Toolbox
% (C) 1999-2003, Written by Vojtech Franc and Vaclav Hlavac
% <a href="http://www.cvut.cz">Czech Technical University Prague</a>
% <a href="http://www.feld.cvut.cz">Faculty of Electrical Engineering</a>
% <a href="http://cmp.felk.cvut.cz">Center for Machine Perception</a>
% Modifications:
% 17-sep-2004, VF, revised
[dim,num_data] = size(X);
if nargin < 2, options = []; else options = c2s(options); end
if ~isfield(options,'arg'), options.arg = 1; end
if ~isfield(options,'solver'), options.solver = 'quadprog'; end
options.ker = 'rbf';
switch options.solver,
case 'quadprog'
G2h = 1/((2*pi)^(dim/2)*(sqrt(2)*options.arg)^dim)*...
kernel(X,options.ker,options.arg*sqrt(2));
Ph = -1/((2*pi)^(dim/2)*options.arg^dim)*...
sum(kernel(X,options.ker,options.arg),2)/num_data;
Aeq = ones(1,num_data);
beq = 1;
LB = zeros(num_data,1);
UB = ones(num_data,1);
[Alpha,fval]= quadprog(G2h,Ph,[],[],Aeq,beq,LB,UB);
otherwise
error('Unknown quadratic programming solver.');
end
inx = find(Alpha > 0);
model.Alpha = Alpha(inx)/((2*pi)^(dim/2)*options.arg^dim);
model.b = 0;
model.sv.X = X(:,inx);
model.nsv = length(inx);
model.options = options;
model.fun = 'kernelproj';
return;