pca_b.m

在MATLAB上所使用的PCA程序

function [eigvector, eigvalue] = pca(x)
%PCA	Principal component analysis
%	[EIGVECTOR, EIGVALUE] = PCA(X)
%	X: rows of vectors of data points
%	EIGVECTOR: each column is a eigenvector of X'*X
%	EIGVALUE: eigenvalues of X'*X 
%
%	Type "pca" for a self-demo.

%	Roger Jang, April-6-97

if nargin == 0,
	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];
	[v, d] = pca(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);
	return
end

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

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