lans_plotcovars.m

来自「模式识别工具包」· M 代码 · 共 238 行

%	lans_plotcovars	- Plot 2-D covariance matrices outlines
%
%	[]	= lans_plotcovars(covars,[mu,[,ptype[,divs]]]);
%
%	_____INPUTS_____________________________________________________________
%	covars		Co-variances matrix
%		'full'		dxdxM	array			(3-D array)
%		'diagonal'	dxM	array			(col vectors)
%		'spherical'	1xM	array			(row vector)
%	mu		centres					(col vectors)
%			= 0 if all centered @ origin		(scalar)
%	ptype		plot type (default)			(string)
%			'contour' selects contour plot with divs
%			specifiying # of lines (defaults 20)
%	divs		# of segments to plot			(scalar)
%			defaults	10
%
%	_____NOTES______________________________________________________________
%	- for demo, call function without parameters
%	- to plot contours, set ptype to 'contour'
%	- Only handles 2-D data
%	- the outline denotes 1 standard deviation, set internally by ndev
%
%	_____SEE ALSO___________________________________________________________
%	lans_plotmd gmm
%
%	(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/

function	[]	= lans_plotcovars(covars,mu,ptype,divs);

if nargin>0
%--------------------	regular
plotcontour	= 0;
ndevcontour	= 2.5;
ndev		= 1;

if nargin<4
	divs	= 20;
	if nargin<3
		ptype	= 'b-';
		if nargin<2
			mu	= 0;
		end
	end
end

if strcmp('contour',ptype)
	plotcontour	= 1;
else
	t		= [0:(2*pi)/(divs-1):2*pi];
	cir1		= [sin(t);cos(t)];
	options		= '-hold 1';
end

[D,M]	= size(mu);
smatrix	= size(covars);

if (length(mu)==1)
	if length(smatrix)>2		% 3-D covars
		D	= smatrix(2);
		M	= smatrix(3);
	elseif length(smatrix)==2	% 1-D or 2-D covars, might be full or diag
		if smatrix(1)==smatrix(2)		% assume full
			M = 1;	
			D = smatrix(1);
		else					% diagonal/isotropic
			D	= smatrix(1);
			M	= smatrix(2);
		end
	end
	mu	= mu*zeros(D,M);
end

if M>1		% more than 1 centers
	if length(smatrix)>2
		ctype = 'full';
	else
		if smatrix(1)==1
			ctype	= 'spherical';
		else
			ctype	= 'diagonal';
		end
	end
else		% only 1 center
	switch prod(smatrix)
		case 1,
			ctype	= 'spherical';
		case 2,
			ctype	= 'diagonal';
		otherwise
			ctype	= 'full';
	end
end

switch ctype
	case 'full',
		for m=1:M
			if M>1
				C	= covars(:,:,m);
			else
				C	= covars;
			end
			%find main diagonal
%			[E,D]	= eigtuned(C);
			[E,D]	= eig(C);
			e1	= E(:,1);
			clkang	= acos(e1(2));
			if (e1(1,:)>0)
				clkang	= -clkang;
			end
	
			%-----	scale principal axes length			
			sd		= sqrt(diag(D));
			d1		= ndevcontour*sd(1,1);
			%-----	make mesh
			if plotcontour
				ndiv	= divs+1;
				int	= d1*2/divs;				
				mesh1	= -d1:int:d1;
				xmesh1	= mesh1+mu(1,m);
				ymesh1	= mesh1+mu(2,m);
				X	= ones(ndiv,1)*xmesh1;
				Y	= ymesh1'*ones(1,ndiv,1);
				xvec	= reshape(X,ndiv*ndiv,1);
				yvec	= reshape(Y,ndiv*ndiv,1);
				F	= gauss(mu(:,m)',C,[xvec yvec]);
		
				hold on;
				[c,h]	= contour3(X,Y,reshape(F,ndiv,ndiv));
%				clabel(c,h);
			else
				cir(1,:)	= cir1(1,:)*ndev*sd(2);
				cir(2,:)	= cir1(2,:)*ndev*sd(1);

				%-----	rotate contours
				c	= cos(clkang);
				s	= sin(clkang);
				R	= [c -s;s c];
				cir	= R*cir;
				cir(1,:)= cir(1,:)+mu(1,m);
				cir(2,:)= cir(2,:)+mu(2,m);
				
				lans_plotmd(cir,ptype,options);			
			end
		end
	case 'diagonal',
		for m=1:M
			if plotcontour
				ndiv	= divs+1;
				d1	= ndevcontour*max(sqrt(covars(:,m)));
				int	= d1*2/divs;				
				mesh1	= -d1:int:d1;
				xmesh1	= mesh1+mu(1,m);
				ymesh1	= mesh1+mu(2,m);
				X	= ones(ndiv,1)*xmesh1;
				Y	= ymesh1'*ones(1,ndiv,1);
			gX	= ones(ndiv,1)*gauss(mu(1,m)',covars(1,m),xmesh1')';
			gY	= gauss(mu(2,m)',covars(2,m),ymesh1')*ones(1,ndiv,1);
				F	= gX.*gY;
				hold on;
				[c,h]	= contour3(X,Y,F);
			else
			cir(1,:)	= cir1(1,:)*ndev*sqrt(covars(2,m))+mu(1,m);
			cir(2,:)	= cir1(2,:)*ndev*sqrt(covars(1,m))+mu(2,m);
				lans_plotmd(cir,ptype,options);
			end
		end
	case 'spherical',
		for m=1:M
			if plotcontour
				ndiv	= divs+1;
				d1	= ndevcontour*sqrt(covars(m));
				int	= d1*2/divs;				
				mesh1	= -d1:int:d1;
				xmesh1	= mesh1+mu(1,m);
				ymesh1	= mesh1+mu(2,m);
				X	= ones(ndiv,1)*xmesh1;
				Y	= ymesh1'*ones(1,ndiv,1);
			gX	= ones(ndiv,1)*gauss(mu(1,m)',covars(m),xmesh1')';
			gY	= gauss(mu(2,m)',covars(m),ymesh1')*ones(1,ndiv,1);
				F	= gX.*gY;
				hold on;
				[c,h]	= contour3(X,Y,F);
			else
				cir(1,:)	= cir1(1,:)*ndev*sqrt(covars(m))+mu(1,m);
				cir(2,:)	= cir1(2,:)*ndev*sqrt(covars(m))+mu(2,m);
				lans_plotmd(cir,ptype,options);
			end
		end
end

axis equal;

%verifies that contour lies where prob = const*exp(-.5) 
%gauss(mu',covars,cir')
%gauss(mu',covars,mu')*exp(-.5);

%--------------------	regular ends
else
%--------------------	demo
clf;clc;
disp('running lans_plotcovars.m in demo mode');
%randn('seed',0);
x	= randn(2,50);
%x(2,:)	= x(2,:)*.5;

mu	= mean(x')';
C	= cov(x');	% full
%C	= [1;4]*.1;	% diagonal
%C	= 1;		% spherical

[pc_data,pvar,paxis]	= lans_pca(x);
scaledaxis		= paxis.*(ones(2,1)*sqrt(pvar)');
lans_plotmd(x,'go');
lans_plotaxis(scaledaxis,mu,'r-','-hold 1');
lans_plotcovars(C,mu,'contour',50);
%axis equal;

%--------------------	demo ends
end