pca2.m

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

function [new_data, eigvector, eigvalue] = pca2(data, eigvector_n)
%PCA	Principal component analysis
%	Usage:
%	[NEW_DATA, EIGVECTOR, EIGVALUE] = PCA(DATA, EIG_N)
%
%	DATA: Rows of vectors of (zero-mean) data points
%	EIG_N: No. of selected eigenvectors
%	NEW_DATA: The data after projection
%	EIGVECTOR: Each column of this matrix is a eigenvector of
%		1. DATA'*DATA if size(DATA, 1) >= size(DATA, 2)
%		2. DATA*DATA' if size(DATA, 1) <  size(DATA, 2)
%	EIGVALUE: Eigenvalues of
%		1. DATA'*DATA if size(DATA, 1) >= size(DATA, 2)
%		2. DATA*DATA' if size(DATA, 1) <  size(DATA, 2)
%
%	Note that DATA must be zero-mean'ed before calling this function. 
%
%	Type "pca2" for a self-demo.

%	Roger Jang, 970406, 990612, 991215

R = data'*data;
[v, d] = eig(R);
eigvector = v;
eigvalue = diag(d);

% Sort based on descending order
[junk, index] = sort(-eigvalue);
eigvalue = eigvalue(index);
eigvector = eigvector(:, index);
new_data = data*eigvector;

if eigvector_n < size(data, 2),
	eigvalue = eigvalue(1:eigvector_n);
	eigvector = eigvector(:, 1:eigvector_n);
	new_data = new_data(:, 1:eigvector_n);
end

function selfdemo

	% ====== Demo for 2D data
	data_n = 1000;
	x = 5+randn(data_n, 1);
	y = 10+randn(data_n, 1)/2;
	x_mean = mean(x);
	y_mean = mean(y);
	theta = pi/6;
	z = ((x-x_mean) + i*(y-y_mean))*exp(i*theta)+x_mean+i*y_mean;
	x = real(z);
	y = imag(z);
	plot(x, y, '.');
	axis image;

	t = [x-x_mean y-y_mean];
	[new_data, v, d] = feval(mfilename, t);
	v1 = -v(:, 1);
	v2 = -v(:, 2);

	arrow_x = [-1 0 nan -0.1 0 -0.1]+1;
	arrow_y = [0 0 nan 0.1 0 -0.1];
	arrow = arrow_x + j*arrow_y;
	arrow1 = 2*arrow*(v1(1)+i*v1(2))*d(1)/data_n+x_mean+i*y_mean;
	arrow2 = 2*arrow*(v2(1)+i*v2(2))*d(2)/data_n+x_mean+i*y_mean;
	line(real(arrow1), imag(arrow1), 'color', 'c', 'linewidth', 2);
	line(real(arrow2), imag(arrow2), 'color', 'g', 'linewidth', 2);
	title('Axes for PCA');

	% ====== Demo for Iris data
	load iris.dat
	data = iris(:, 1:4);		% Take only the input part
	data_n = size(data, 1);
	data = data - ones(data_n, 1)*mean(data);	% Make data zero-mean
	[new_data, eigvector, eigvalue] = feval(mfilename, data);
	index1 = find(iris(:,5)==1);
	index2 = find(iris(:,5)==2);
	index3 = find(iris(:,5)==3);
	figure;
	plot(	new_data(index1, 1), new_data(index1, 2), '*', ...
		new_data(index2, 1), new_data(index2, 2), 'o', ...
		new_data(index3, 1), new_data(index3, 2), 'x');
	legend('Class 1', 'Class 2', 'Class 3');
	title('PCA projection of IRIS data onto the first 2 eigenvectors');
	axis equal; axis tight;

	figure;
	plot(	new_data(index1, 3), new_data(index1, 4), '*', ...
		new_data(index2, 3), new_data(index2, 4), 'o', ...
		new_data(index3, 3), new_data(index3, 4), 'x');
	legend('Class 1', 'Class 2', 'Class 3');
	title('PCA projection of IRIS data onto the last 2 eigenvectors');
	axis equal; axis tight;