pca.m

一个关于数据聚类和模式识别的程序,在生物化学,化学中因该都可以用到.希望对大家有用,谢谢支持

function [DS2, eigVec, eigValue] = pca(DS, eigVecNum)
%pca: Principal component analysis
%	Usage: [DS2, eigVec, eigValue] = pca(DS, eigVecNum)
%		DS: DS.input is the data matrix, where each column is a data vector
%			Please try "DS=dcData(1)" to get an example of DS.
%		eigVecNum: No. of selected eigenvectors
%		DS2: output data set, where DS2.input is the data after projection
%		eigVec: Each column of this matrix is a eigenvector of DS.input*DS.input' sorted by its decending order of eigen values
%		eigValue: Eigenvalues of (DS.input*DS.input') corresponding to eigVec
%
%	Note that DS.input must be zero-mean'ed before calling this function. 
%
%	Type "pca" for a self-demo.

%	Roger Jang, 970406, 990612, 991215, 20060506

if nargin<1, selfdemo; return; end
if ~isstruct(DS)
	fprintf('Please try "DS=prData(1)" to get an example of DS.\n');
	error('The input DS should be a structure variable!');
end
if nargin<2, eigVecNum = min(size(DS.input)); end

m = size(DS.input,1);	% Dimension of data point
n = size(DS.input,2);	% No. of data point
A = DS.input;

if n>=m
	[eigVec, eigValue] = eig(A*A');
	eigValue = diag(eigValue);
	% ====== Sort based on descending order
	[junk, index] = sort(-eigValue);
	eigValue = eigValue(index);
	eigVec = eigVec(:, index);
	if eigVecNum<m
		eigValue = eigValue(1:eigVecNum);
		eigVec = eigVec(:, 1:eigVecNum);
	end
else	% This is an efficient method which computes the eigvectors of A*A' when size(A,1)>size(A,2)
	% A*A'*x=lambda*x ===> A'*A*A'*x=lambda*A'*x ===> eigVec of A'*A is A'*x, multiply these eigVec by A, we have A*A'*x=lambda*x ===> Got it!
	[eigVec, eigValue] = eig(A'*A);
	eigValue = diag(eigValue);
	% ====== Sort based on descending order
	[junk, index] = sort(-eigValue);
	eigValue = eigValue(index);
	eigVec = eigVec(:, index);		% Eigenvectors of A'*A
	eigVec = A*eigVec;			% Eigenvectors of A*A'
	eigVec = eigVec*diag(1./(sum(eigVec.^2).^0.5)); % Normalization
	if eigVecNum<n
		eigValue = eigValue(1:eigVecNum);
		eigVec = eigVec(:, 1:eigVecNum);
	end
end

DS2=DS;
DS2.input=eigVec'*A;


% ====== Self demo
function selfdemo
	% ====== Demo for 2D data
	dataNum = 1000;
	data = randn(1,dataNum)+j*randn(1,dataNum)/3;
	data = data*exp(j*pi/6);	% 臂锣30