Code covered by the BSD License  

Highlights from
Nipals algorithm for Principle Component Analysis

image thumbnail

Nipals algorithm for Principle Component Analysis

by

 

10 Jun 2013 (Updated )

This function is written largely based on nipals function from R chemometrics package.

nipals(X,a,it,tol )
function [ T,P,pcvar ] = nipals(X,a,it,tol )
%Nipals algorithm for Principle Component Analysis
%
% Input:
%   X: a matrix with rows obseravations and columns variables. 
%   a: number of components desired in the output. 
%   it: maximum number of iterations allowed to find one component. default 1000.
%   tol: tolerance value to assese convergence. default 1e-4
%
% Output:
%   T: scores or transformed coordinates.
%   P: loadings or coefficients.
%   pcvar: variance captured in each output component.
%
% If X is m*n matrix, then T is m*a, P is n*a
%
% This function is written largely based on nipals function from R chemometrics
% package. It was frustrating that Matlab did not has a robust nipals
% function so I wrote one for my project and share here.
%
% Author: Qiaonan Duan, 6/7/2013, MSSM.
%
% Updates:
% Line 48 (6/10/2013): prec = thnew-th; --> prec = (thnew-th)'*(thnew-th);
%
if nargin == 2
		it = 1000;
		tol = 1e-4;
	elseif nargin == 3
		tol = 1e-4;
	end

	[obsCount,varCount] = size(X);
	Xh = X - repmat( mean(X,1), obsCount, 1 );
	T = zeros(obsCount,a);
	P = zeros(varCount,a);
	pcvar = zeros(varCount,1);
	varTotal = sum(var(Xh));
	currVar = varTotal;
	nr = 0;

	for h = 1:a
		th = Xh(:,1);
		ende = false;

		while(~ende)
			nr = nr+1;
			ph = Xh'*th/(th'*th);
			ph = ph/norm(ph);
			thnew = Xh*ph/(ph'*ph);
			prec = (thnew-th)'*(thnew-th);
			th = thnew;

			if prec <= tol^2
				ende = true;
			elseif it <= nr
				ende = true;
				disp('Iteration stops without convergence')
			end
		end

		Xh = Xh-th*ph';
		T(:,h) = th;
		P(:,h) = ph;
		oldVar = currVar;
		currVar = sum(var(Xh));
		pcvar(h) = ( oldVar - currVar )/varTotal;
		nr = 0;

end

Contact us