lans_expect.m

模式识别工具包

%	lans_expect	- Principal curve local averaging E(f|lambda)
%
%	[psurf]	= lans_expect(prdata,prsurf,psurf[,para])
%
%	_____OUTPUT_____________________________________________________________
%	psurf	new principal surface				(structure)
%		see lans_psurf.m
%
%	_____INPUT______________________________________________________________
%	prdata	data sorted w.r.t. prsurf.x			(col vectors)
%	prsurf	projection of points onto current psurf		(structure)
%		see lans_psurf.m
%	psurf	current principal surface			(structure)
%		see lans_psurf.m
%	para	see para.m paraget.m				(string)
%
%	_____NOTES______________________________________________________________
%	- requires NETLAB	: gmmpost.m
%	- duplicate knots aggregated via dupavg.m
%	- prior to calling this routine, lans_project needs to be run to obtain
%	  projections prsurf
%	- wlimit must be used to limit wi for EM algorithm
%	- psurf only used by EM algorithm
%	- non-EM algorithms DO NOT use psurf
%
%	_____SEE ALSO___________________________________________________________
%	lans_psurf lans_project dupavg
%
%
%	(C) 1999.05.01 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______________________________________________________________
%	- get rid of knots concentration problem
%		concentration seems unavoidable due to clustering nature
%		unless compute x by projection?
%	- minimize situations where wi goes to zero in smoother
%	- optimized knot selection
%	- set maximum allowed variance
%	- get rid of dependency on NETLAB

function	[psurf]	= lans_expect(prdata,prsurf,psurf,para)

if nargin<4
	para	= [];
	if nargin<1
		clc;
		help lans_expect;
		break;
	end
end

[d N]	= size(prdata);

%-----	delta value added to make singular covariances non-singular
delta	= 1e-10;

%----- get parameters
algo		= paraget('-algo',para);
covarclip	= paraget('-covarclip',para);
covarmax	= paraget('-covarmax',para);
covarmin	= paraget('-covarmin',para);
fixwhat		= paraget('-fix',para);
minprior	= paraget('-minprior',para);
orient		= paraget('-orient',para);
recompute	= paraget('-recompute',para);
smoother	= paraget('-smoother',para);
span		= paraget('-span',para);

%----------	compute span
sp	= round(span*N);	%nearest points excluding myself

%----------	FIRST get rid of duplicates (in k) for non EM methods
%		duplication for EM ok since # of centers is fixed
if algo<3
	[udf,x,w]	= dupavg([prdata;prsurf.f],prsurf.x);
	prdata		= udf(1:d,:);
	f		= udf(d+1:d+d,:);

%	The following works poorly
%	x	= lans_diversify(prsurf.x);
%	f	= prsurf.f;
%	w	= ones(1,length(x));
end

switch (algo)
	%	Hastie & Stuetzle algorithm
	case 1
		psurf.f		= smooth(x,prdata,w,para);
		psurf.x		= lans_1speed(psurf.f);

	%	Banfied and Raftery algorithm
	case 2
		nresidual	= smooth(x,prdata-f,w,para);
		psurf.f		= f + nresidual;
		psurf.x		= lans_1speed(psurf.f);

	%	Tibshirani EM algorithm (smoothed)
	%	orientation works only for spherical covariances
	case 3
		wlimit	= minprior/psurf.M;

		y	= prdata;
		%----------	(1) Compute Posterior	R_{M,N} = P(m|x)
		%		R is M x N
		%		M : # components
		%		N : # samples

		mix	= lans_psurf2mix(psurf);

		if orient
			% -psurf.beta contains the current spherical covariances
			% -psurf.covars contains the oriented full covariances
			% to be used in computing the posterior
			mix.covar_type	= 'full';
		end
		R	= gmmpost(mix,prdata')';
		w	= sum(R');

		%-----	Replace minuscule weights with .1 (IMPORTANT)
		w	= lans_replace(wlimit,wlimit,w);
		wi	= 1./w;
		wid	= ones(d,1)*wi;	

		%----------	(2) Recompute knots with others fixed
		% * assumes ordering exists in psurf.f!!!! false if psurf
		% NOT ordered initially e.g.
		% a) initialized using random points and/or
		% b) EM version does not use projection explicitly
		%
		% a can be taken care of if span is large enough, for b,
		% we have to reset psurf.f to the values of the prsurf
		% in the right order

		switch recompute
			case 'unitinterval',
				x=0:length(psurf.x)-1;
				x=x/length(psurf.x);
			case 'unitspeed',
				x	= lans_1speed(psurf.f);
				x	= lans_diversify(x);
			case 'projection',
%invalid because centers are determined via cubic spline smoothed EM algorithm
%i.e. the centers do not correspond to projections
				x	= prsurf.x-min(prsurf.x);
				x	= lans_diversify(x);
			case 'optimized',
				% x = *;
				% unimplemented;
		end

		%---------	(3) update psurf.f,psurf.covars,psurf.priors
		%---------	(3a)	psurf.f
		if ~strcmp(fixwhat,'centers')
			Diy		= (wid'.*(R*y'))';
			if smoother
				psurf.f	= smooth(x,Diy,w,para);
			else
				psurf.f	= Diy;			% no smoothing
			end
		end
		psurf.x	= lans_1speed(psurf.f);

		%---------	(3b)	psurf.covars
	if ~orient
		if strcmp(psurf.covartype,'spherical')
			for m=1:psurf.M
				d2	= vdist2(y,psurf.f(:,m)*ones(1,N))';
				C	= R(m,:)*d2*wi(m)/d;	
				%-----	only update if new covars is within range
				if covarclip
					C	= lans_clip(C,covarmin,covarmax);
					psurf.covars(m)	= C;
				else
					if (C>covarmin)&(C<covarmax)
						psurf.covars(m)	= C;
					end 	
				end
			end
		elseif strcmp(psurf.covartype,'diagonal')
			for m=1:psurf.M
				d1	= y-psurf.f(:,m)*ones(1,N);
				C	= (d1.*d1)*R(m,:)'*wi(m);
				%-----	only update if new covars is within range
				if covarclip
					C	= lans_clip(C,covarmin,covarmax);
					psurf.covars(:,m)	= C;
				else
					if (min(C)>covarmin)&(max(C)<covarmax)
						psurf.covars(:,m)	= C;
					end 	
				end

			end
		elseif strcmp(psurf.covartype,'full')
			for m=1:psurf.M
				d1	= y-psurf.f(:,m)*ones(1,N);
				d1R	= d1.*sqrt(ones(d,1)*R(m,:));
				C	= (d1R*d1R')*wi(m);
				[V,E]	= eigtuned(C);
				Es	= diag(E);
				%-----	only update if new covars is within range
				if covarclip
					Es	= lans_clip(Es,covarmin,covarmax);
					psurf.covars(:,:,m)	= V*diag(Es)*V';
				else
					if (Es(psurf.D)>covarmin)&(Es(1)<covarmax)
						psurf.covars(:,:,m)	= C;
					else
%tstr	= sprintf('[%d] added delta',m);
%disp(tstr);
						psurf.covars(:,:,m)	= C+delta*eye(psurf.D);
					end 	

				end


			end
		end
	else
	%assumes spherical/unoriented inversed covariances in psurf.beta
		for m=1:psurf.M
			d2	= vdist2(y,psurf.f(:,m)*ones(1,N))';
			C	= R(m,:)*d2*wi(m)/d;	
			%-----	only update if new covars is within range
			if covarclip
				C		= lans_clip(C,covarmin,covarmax);
				psurf.beta(m)	= 1/C;
			else
				if (C>covarmin)&(C<covarmax)
					psurf.beta(m)	= 1/C;
				end 	
			end
		end

	end
		%---------	(3c)	psurf.priors
		psurf.priors		= w/N;

	otherwise
		error('UNKNOWN Principal Curve algorithm number');
end