pca_dd.m

数据挖掘的工具箱,最新版的,希望对做这方面研究的人有用

%PCA_DD Principal Component data description
%
%       W = PCA_DD(A,FRACREJ,N)
%
% Traininig of a PCA, with N features (or explaining a fraction N of
% the variance).
%
% Default: N=0.9

% Copyright: D.M.J. Tax, R.P.W. Duin, [e-mail]
% Faculty of Applied Physics, Delft University of Technology
% P.O. Box 5046, 2600 GA Delft, The Netherlands

function W = pca_dd(a,fracrej,n)

if (nargin<3)
  n = 0.9;
end
if (nargin<2)
  fracrej = 0.05;
end
if (nargin<1)|isempty(a)
	W = mapping(mfilename,{fracrej,n});
	W = setname(W,'Principal Component Analysis data description');
	return
end

if ~ismapping(fracrej)           %training

	% remove the labels:
	a = target_class(a);     % only use the target class
	[m,k] = size(a);
	% Be careful with the mean:
	meana = repmat(mean(a),m,1);
	a = (a - meana);

	% Train it and compute the reconstruction error:
	if (n<0) % we requested the low- instead of high-variance directions!
		w = pca(a,k);
		w = w(:,end+n+1:end);
	else
		w = pca(a,n);
	end
	W = +w; % extract the projection matrix from the eigenvectors:
	W = W.rot;
	if size(W,2)==k
		warning('PCA_DD: Output dimensionality is equal to input dimensionality!');
	end
	Proj = W*inv(W'*W)*W';
	% project and find the distribution of the distance:
	dif = a - a*Proj;
	d = sum(dif.*dif,2);

	% obtain the threshold:
	thr = dd_threshold(d,1-fracrej);

	%and save all useful data:
	% (I know I just have to store W instead of Proj, but I do not like
	% to compute the inverse of W'*W over and over again, this uses just
	% some disk/memory space):
	W = [];  % W was already used, forget that one...
	W.P = Proj;
	W.m = meana(1,:);
	W.threshold = thr;
	W.scale = mean(d);
	W = mapping(mfilename,'trained',W,str2mat('target','outlier'),k,2);

else                               %testing

	W = getdata(fracrej);  % unpack
	m = size(a,1);

	%compute reconstruction error:
	dif = +a - repmat(W.m,m,1);
	dif = dif - dif*W.P;
	out = sum(dif.*dif,2);
	newout = -[out, repmat(W.threshold,m,1)];

	% map to probabilities:
	%newout = dist2dens(newout,W.scale);

	W = setdat(a,newout,fracrej);

end

return