evalifgt.m

Non-parametric density estimation

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% evalIFGT  Evaluate the density estimate using the "improved" Fast Gauss Transform
%
%   [e,b] = evalIFGT(X,Y,N [,Nc,rC]) -- eval likelihood ("e") of the points Y under
%                       the density estimate X using N coefficients of the
%                       "improved" Fast Gauss Transform; the value "b" is the bound
%                       on the (absolute) error which could arise.
%
%  Optional arguments:
%    Nc  -- # of clusters to use for "X", default is sqrt(Npoints) 
%    rC  -- Cutoff radius (in std deviations) to exclude contributions, default 3
%
% See: Yang, Duraiswami, Gumerov; "Improved Fast Gauss Transform", submitted to 
%         the Siam Journal of Scientific Computing, 2004
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [estimate,errbound] = evalIFGT(pp,q,Ncoeff,Nclusters,rCutoff)
  p = kde(pp);	% copy constructor to dodge later rescaling...
  if (p.type ~= 0)  
    error('Sorry -- FGT = fast Gauss transform; it needs Gaussian kernels');
  end;
  if (size(p.bandwidth,2)>2*p.N)
    error('Sorry -- IFGT currently supports only uniform bandwidths');
  end;
  if (nargin<4) Nclusters = round(sqrt(getNpts(p))); end;
  if (nargin<5) rCutoff = 3; end;  
  if (isa(q,'kde')) qpts = getPoints(q); else qpts = q; end;
  
  BW = getBW(p,1); BWorig = BW;
  if (any( BW - BW(1) ))   % CONVERT TO SINGLE, SCALAR BW:
    p = rescale(p, 1./BW); %  if differ in dimensions, need to rescale
    qpts = qpts .* repmat(1./BW,[1,size(qpts,2)]);
    BW = 1;
  else BW = BW(1);         % already scalar; can just drop other dim's
  end;

  [c,cPts,cWts,cWt,cRad] = fpClusterK(p,Nclusters);
  %[c,cPts,cWts,cWt,cRad] = fpClusterR(p,sqrt(2)*BW);
  coeff = findCoeff(c,cPts,cWts,cRad,BW,Ncoeff);
  [estimate,errbound] = evalCoeff( qpts, c,coeff,Ncoeff,cWt,BW,cRad,rCutoff);

  % Change norm. constant (due to rescaling operation)
  scale = p.D*log(BW) - sum(log(BWorig));
  estimate = estimate * exp(scale); errbound = errbound * exp(scale);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% fpCluster -- fast, "farthest point" clustering method
%   cluster points of "p" into K clusters, described by "centers",
%    "clusters" (cell array of pts to each cluster),
%    cWeight (weight per cluster), and maximum radius of any cluster.
%
function [centers, clusters, weights, cWeight, radius] = fpClusterK(p, K)
  points = getPoints(p); wts = getWeights(p);
  [D,N] = size(points);
  centers = zeros(D,K); clusters = cell(1,K); weights = cell(1,K);
  assign = ones(1,N); dmin = zeros(1,N)+inf;
  next = fix(rand(1)*N)+1;  % choose 1st center at random
  for i=1:K
    centers(:,i) = points(:, next);
    d = points - repmat(centers(:,i),[1,N]);
    d = sqrt(sum(d.^2,1));
    F=find(d<dmin); dmin(F)=d(F); assign(F) = i;
    [radius, next] = max(dmin); % next center is a farthest point
  end;
  cWeight = zeros(1,K);
  for i=1:K
    clusters{i}=points(:, find(assign == i) );
    weights{i}=wts(:, find(assign == i) );
    cWeight(i) = sum(weights{i}); %size(clusters{i},2) / N;
  end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Same thing but cluster until radius < rMax
function [centers, clusters, weights, cWeight, radius] = fpClusterR(p, rMax)
  points = getPoints(p); wts = getWeights(p);
  [D,N] = size(points); K = N; centers = zeros(D,K);
  assign = ones(1,N); dmin = zeros(1,N)+inf;
  next = fix(rand(1)*N)+1;  % choose 1st center at random
  i=0; radius = inf;
  while (radius > rMax),
    i=i+1; centers(:,i) = points(:, next);
    d = points - repmat(centers(:,i),[1,N]);
    d = sqrt(sum(d.^2,1));
    F=find(d<dmin); dmin(F)=d(F); assign(F) = i;
    [radius, next] = max(dmin); % next center is a farthest point
  end;
  K=i; centers = centers(:,1:K);
  cWeight = zeros(1,K); clusters = cell(1,K); weights = cell(1,K);
  for i=1:K
    clusters{i}=points(:, find(assign == i) );
    weights{i}=wts(:, find(assign == i) );
    cWeight(i) = sum(weights{i}); %size(clusters{i},2) / N;
  end;



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% findCoeff -- find the Taylor series coefficients of the Gaussian sum
%   described by each cluster.
%
function coeff = findCoeff(centers, points, weights, radius, h, Nterms)
  h = sqrt(2)*h;        % stupid transform...

  Npart = length(points); coeff = cell(1,Npart); D = size(centers,1);
  NptsTotal = 0;
  Npts = zeros(1,Npart); for i=1:Npart, Npts(i) = size(points{i},2); end;
  NptsTotal = sum(Npts);

  for i=1:Npart
    NptsI = Npts(i);
    vals = ( points{i}-repmat(centers(:,i),[1,NptsI]) )/h;

    Ncoeff = round(exp( sum(log(Nterms:Nterms+D-1))-sum(log(1:Nterms)) ));
    coeffI = zeros(NptsI, Ncoeff);

    start = 0; startNew = 1;
    coeffI(:, start+1) = exp( -sum(vals.^2,1) )';
    pos = ones(1,D); alpha = zeros(D,1);

    for j=2:Nterms
      Nprev = startNew - start;
      Nadd = sum(Nprev-pos+1); alphaNew = zeros(D,Nadd);
      m = 1; posNew(1) = m;
      for k=1:D
        for l=pos(k):Nprev
          alphaNew(:,m) = alpha(:,l); alphaNew(k,m) = alphaNew(k,m)+1;
          constFactor = 2  ./ prod(max(alphaNew(:,m),1));
          coeffI(:,startNew+m) = vals(k,:)' .* coeffI(:,start+l) * constFactor;
          m = m+1;
        end;
        if (k ~= D) posNew(k+1) = m; end;
      end;
      pos = posNew; alpha = alphaNew; start = startNew; startNew = start+Nadd;
    end;
    %coeffI = sum(coeffI,1)/NptsTotal;
    coeffI = weights{i}*coeffI;
    coeff{i} = coeffI;
  end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% evalCoeff -- evaluate the Taylor series at a number of new locations
%   also returns an upper bound on the incurred error
%
function [est, err] = evalCoeff(locations,centers,coeff,Nterms,cWt,h,cRad,rCutoff)

  h = sqrt(2)*h;	% stupid transformation

  Npts = size(locations,2); Npart = size(centers,2); D = size(locations,1);
  est = zeros(1,Npts);
  err = zeros(1,Npts);
  for i=1:Npart  
      coeffI = coeff{i};
      vals = ( locations - repmat(centers(:,i),[1,Npts]) ) / h;
      distance2 = sum(vals.^2,1);
      PTS = find( distance2 < rCutoff^2);
      PTSN = find( distance2 >= rCutoff^2);

      start = 0; startNew = 1;
      terms = zeros(length(PTS),size(coeffI,2));
      terms(:,start+1) = exp( - distance2(PTS) )';
      pos = ones(1,D); 
      for j=2:Nterms
        Nprev = startNew - start;
	Nadd = sum(Nprev-pos+1); 
	m=1; posNew(1)=m;
	for k=1:D
	  for l=pos(k):Nprev
	    terms(:,startNew+m) = vals(k,PTS)' .* terms(:,start+l);
	    m = m+1;
	  end;
	  if (k~=D) posNew(k+1) = m; end;
	end;
	pos = posNew; start = startNew; startNew = start+Nadd;
      end;
      est(PTS) = est(PTS) + (coeffI * terms');

      % error bound addition for included points...  + Qin * 2^p/p! rhox^p rhoy^p
      err(PTS)=err(PTS) + cWt(i)*exp( Nterms*log(2*rCutoff)- sum(log(1:Nterms)) + Nterms*log(cRad/h) );

      % error bound addition for excluded points...  + Qin * exp(-rhoy^2+rhox^2)
      err(PTSN)=err(PTSN) + cWt(i)*exp( - rCutoff^2 + (cRad/h)^2 );

 end; 
 h = h / sqrt(2);
 
 est = est ./ (2*pi*h^2)^(D/2);
 err = err ./ (2*pi*h^2)^(D/2);