rsde.m~

用matlab实现的统计模式识别工具箱

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, 100 );
%  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.75));
%  x = linspace(-4,8,100);
%  plot(x,kernelproj(x,model),'g'); 
%  ppatterns(model.sv.X,'ob',13);
%
% See also 
%  KERNELPROJ.
%

% Input arguments
%-------------------------------------------------------
[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';

% Solve associted QP task
%-------------------------------------------------------

switch options.solver,
  case 'quadprog'
   G2h = kernel(X,options.ker,options.arg*2);
   Ph = -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);
model.b = 0;
model.sv.X = X(:,inx);
model.nsv = length(inx);
model.options = options;
model.fun = 'kernelproj';

return;
% EOF