lans_psurf.m

模式识别工具包

%	lans_psurf	- Principal curve generation/initialization
%
%	[psurf,prsurf,prdata,mse]	= lans_psurf(data) =
%	lans_psurf(data,options)
%
%	_____OUTPUTS____________________________________________________________
%	psurf	principal curve structure			(structure)
%		.x 	knots in lower/latent dimension		(col vectors)
%		.f	psurf values f(x) in original dimension	(col vectors)
%		.covars	covariance of data about .f		(2/3-D array)
%		.covartype covariance type			(string)
%			'spherical'	.covars = row	vector	(1xM)
%			'diagonal'	.covars = col vectors	(DxM)
%			'full' 		.covars = 3-D array	(DxDxM)
%		.priors		prior probability of each knot	(row vector)
%		.M	# of knots/latent vectors
%		.D	dimensionality of original space
%		.Q	dimensionality of manifold
%		prsurf	projection of points onto current psurf	(structure)
%			see lans_psurf.m
%		prdata	data sorted with respect to prsurf.x	(col vectors)
%		mse	Mean squared error of prdata to prsurf	(scalar)
%
%	_____INPUTS_____________________________________________________________
%	data	unormalized data in original dimension		(col vectors)
%	options	see lanspara.m					(string)
%		-clos						(binary)
%		-covartype	spherical diagonal full
%		-pinittype					(integer)
%		1	1st principle component
%		2	random values about mean with .1 stdev
%		-M	# knots %				(integer)
%			default = # data points
%
%	_____NOTES______________________________________________________________
%	- uses randn.m
%	- psurf	is arc-length parameterized but NOT NECESSARY
%	ordered except when initialized to first principal component
%	- code added to ensure that first dimension of eigenvector is >0 (+ve)
%	this is
%	to take care of the eigs.m discrepencies between matlab and matcom
%
%	_____SEE ALSO___________________________________________________________
%	lans_psurfinit	lans_project	lans_1speed
%
%	(C) 2001.12.30 Kui-yu Chang
%	http://lans.ece.utexas.edu/~kuiyu

%	This program is free software; you can redistribute it and/or modify
%	it under the terms of the GNU General Public License as published by
%	the Free Software Foundation; either version 2 of the License, or
%	(at your option) any later version.
%
%	This program is distributed in the hope that it will be useful,
%	but WITHOUT ANY WARRANTY; without even the implied warranty of
%	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%	GNU General Public License for more details.
%
%	You should have received a copy of the GNU General Public License
%	along with this program; if not, write to the Free Software
%	Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307, USA
%	or check
%			http://www.gnu.org/

%	_____TO DO_________________________________________________________ 
%	- M (variable # of knots)
%	- closed PCA initialization

function [psurf,prsurf,prdata,mse]= lans_psurf(data,options)
if nargin<2
	options	= [];
	if nargin<1
		clc;
		help lans_psurf;
		return;
	end
end

[D,N]		= size(data);

algo		= lans_paraget('-algo',options);
beta		= lans_paraget('-beta',options);
	beta	= beta+(beta==0);		% set to 1 if unspecified
clos		= lans_paraget('-clos',options);
covartype	= lans_paraget('-covartype',options);
M		= lans_paraget('-M',options);
orient	= lans_paraget('-orient',options);
pinittype	= lans_paraget('-pinittype',options,1);

%----- same # of knots as # data points if M unspecified
if M==0|algo<4
	M	= N;
end
	
%----------	set up mixture parameters (needed later for sorting)
psurf.M		= M;
psurf.D		= D;
psurf.Q		= 1;
psurf.beta	= beta;
if algo>2
psurf.covartype	= covartype;
psurf.priors	= ones(1,M)/M;
switch covartype
	case 'spherical'
		psurf.covars	= ones(1, M)/beta;
	case 'diagonal'
    		psurf.covars	= ones(D, M)/beta;
    	case 'full'
    		psurf.covars	= repmat(eye(D), [1 1 M])/beta;  
end
end
%----------	standardize data
[ndata,xmean,xstd]	= lans_stand(data);
switch(pinittype)
	case 1				% 1st principal component direction
	%-----	recover variance
	ndata			= ndata.*(xstd*ones(1,N));

	%-----	compute principal axis
	[pc_data,pvar,paxis]	= lans_pca(ndata);
	paxis1			= paxis(:,1);

	%-----	project(linear) data on principal axis
	psurf.x	= (ndata'*paxis1)';		% gives component value
	psurf.f	= paxis1*psurf.x;		% gives value of curve vector f

	%-----	sort points w.r.t. projections
	psurf.x		= psurf.x - min(psurf.x);	% sort w.r.t. x
	[dumb,sidx]	= sort(psurf.x);
	psurf		= lans_porder(psurf,sidx,options);

	%-----	recover mean
	psurf.f		= psurf.f + xmean*ones(1,N);

	%-----	add one point to make a triangle for closed curve
	if clos~=0

	end;

	case 2,				%small random values about mean (UNORDERED) 
	randompoints	= xmean*ones(1,N)+.1*xstd*ones(1,N).*randn(size(data));
	psurf.x		= lans_1speed(randompoints);
	psurf.f		= randompoints;

	case 3,				%small random values perturbed about data
	varvec		= .1*sqrt(diag(cov(data')))*ones(1,N);
	psurf.f		= data+varvec.*randn(size(data));
	psurf.x		= lans_1speed(psurf.f);
end;

%__________	orient covariances here if specified
if orient
	psurf	= lans_orient(psurf,options);
end
psurf
%----------	compute projected surfaces and corresponding sorted data
[prsurf,prdata,mse]	= lans_project(data,psurf);