lans_pcurve1.m

模式识别工具包

%	lans_pcurve1	- Principal curve single iteration
%
%	[psurf,prsurf,prdata,mse,se] = lans_pcurve1(psurf,prsurf,prdata[,para])
%
%	_____OUTPUT_____________________________________________________________
%	psurf	new principal surface				(structure)
%		see lans_psurf.m
%	prsurf	projection of points onto new psurf		(structure)
%		see lans_psurf.m
%	prdata	ordered data w.r.t. prsurf.x			(col vectors)
%	mse	mse of data to prsurf.f				(scalar)
%		(NOT mean distance, but mean 'squared distance')
%	se	squared error(dist) of each data to prsurf.f	(row vector)
%
%	_____INPUT______________________________________________________________
%	psurf	current principal curve				(structure)
%		see lans_psurf.m
%	prsurf	projection curve
%		# knots == # data points MUST!!
%	prdata	data vectors					(col vectors)
%	para	see lans_paraget.m			(string)
%	
%	_____SEE ALSO___________________________________________________________
%	lans_psurf lans_expect lans_project lans_projdist	
%
%	_____NOTES______________________________________________________________
%	for demo, call function without parameters
%
%	- psurf will NOT be ordered if input psurf is NOT ordered
%	- Requires current projection of data onto the curve, i.e. prsurf
%	lans_project.m must be executed (at least once) before
%	calling this function
%	- #data = #knots
%	- only psurf.x and psurf.f is modified, preserving other components
%
%	(C) 2003.07.29 Kui-yu Chang
%	http://www.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______________________________________________________________
%	- Tibshirani algorithm algo=3 NOT working well
%		- knots tend to cluster together (unoptimized)
%		- the lone knot drifts away
function [psurf,prsurf,prdata,mse,se] = lans_pcurve1(psurf,prsurf,prdata,para)

if nargin>0
%__________ REGULAR ____________________________________________________________
if nargin<4
	para=[];
end

algo		= lans_paraget('-algo',para);
orient	= lans_paraget('-orient',para);

%-----	compute expectation
[psurf]			= lans_expect(prdata,prsurf,psurf,para);
if orient
	psurf	= lans_orient(psurf,para);
end

%-----	project data onto new curve (useless for EM algorithm 3)
[prsurf,prdata,mse,se]	= lans_project(prdata,psurf,para);

%__________ REGULAR ends _______________________________________________________
else
%__________ DEMO _______________________________________________________________
rseed	= 1;
clf;clc;
disp('running lans_pcurve1.m in demo mode');

%para	= '-algo 1 -clos 0 -recompute unitspeed -smoother 3 -smoopara 1 -span 5e-2 -showt 1';

para	= '-algo 1 -clos 1 -recompute unitspeed -smoother 1 -smoopara 2 -span 1 -showt 1';

%para	= '-algo 3 -clos 0 -covartype spherical -orient 0 -covarclip 1 -covarmin 1e-2 -covarmax 5e-2 -recompute unitspeed -pinittype 1 -showt 1 -smoother 3 -smoopara 1 -span 1e-2';

%bad
%para	= '-algo 3 -clos 0 -covartype spherical -orient 0 -covarclip 1 -covarmin 1e-2 -covarmax 5e-2 -recompute unitspeed -pinittype 1 -showt 1 -smoother 3 -smoopara 1 -span 1e-2';

% slow, '-clos 0' also works
%para = '-algo 1 -clos 1 -recompute unitspeed -smoother 1 -smoopara 2 -span 1 -showt 1';

% fast, '-clos 0' also works
para	= '-algo 2 -clos 1 -recompute unitspeed -smoother 3 -smoopara 1 -span 5e-2 -showt 1';




N	= 50;		% # of samples 
K	= 50;		% # of iterations
cres	= 20;		% resolution of covariances circle plot
cint	= 10;		% interval between plotting covars

x						= lans_gencircle(N,rseed);
%----------	Initial Principal Curve
k	= 0;
[psurf,prsurf,prdata,mse]	= lans_psurf(x,para);

algo	= lans_paraget('-algo',para);
ctype	= lans_paraget('-covartype', para, 0);
%----------	K iterations
for k=1:K
	[psurf,prsurf,prdata,mse]	= lans_pcurve1(psurf,prsurf,prdata,para);
	if ~mod(k,cint)
		clf;
		lans_plotmd(x,'mx',psurf.f,'k-',para);
	end
   if ~mod(k,cint)
      if (ctype~=0) 
         lans_plotcovars(psurf.covars,psurf.f,'m-',cres);
      end
		title(sprintf('[%2d] MSE = %f',k,mse));
		drawnow;
	end;
end

if algo>2
	lans_plotcovars(psurf.covars,psurf.f,'m-',cres);
end
%__________ DEMO ends __________________________________________________________
end