lans_gausspdf.m

模式识别工具包

%	lans_gausspdf	- Evaluates Gaussian pdf @ multiple vector points 
%
%	[pdf]	= lans_gausspdf(x,mu,sig)
%		= lans_gausspdf(x,mu,sig,nonorm)
%		= lans_gausspdf(x,mu,sig,alpha,EET)
%		= lans_gausspdf(x,mu,sig,nonorm,alpha,EET)
%
%	_____OUTPUTS____________________________________________________________
%	pdf	pdf @ x						(row vector)
%
%	_____INPUTS_____________________________________________________________
%	x	locations					(col vectors)
%	mu	mean						(vector)
%	sig	covariance				
%			spherical				(scalar)
%			diagonal	1xD or Dx1 vector	(vector)
%			full		DxD			(matrix)
%
%	nonorm	does not compute normalization factor if defined (dummy)
%	alpha	specifies amount of clamping			(scalar)
%	EET	space span by tangential orthonormal matrix	(DxQ)
%		EET	= E*E'
%
%	_____EXAMPLE____________________________________________________________
%
%
%	_____NOTES______________________________________________________________
%	- optimized for each covariance format
%	- in the case of a specially clamped covariance matrix in the R^Q
%	  submanifold spanned by the column space of E, the covariance can be
%	  specified as
%
%	  	sigfull	= B*eye(D)+(S-B)EET;
%
%	  with the following simplified inverse and determinant
%
%	  	inv(sigfull)	= eye(D)/B - ((S-B)/(B*S))*EET;
%	  	det(sigfull)	= (S^Q)*(B^(D-Q))
%	where
%		B	= (D-alpha*Q)*sig/(D-Q)
%		S	= alpha*sig;
%
%	in this case, the sig provided to this function is simply a scalar 
%	- if D=Q, then attenuate eye(Q) by alpha
%	- if nonorm specified, then the Gaussian activation (not a pdf) is
%	  computed instead, speeding up Radial Basis Function (RBF) computation.
%
%	_____SEE ALSO___________________________________________________________
%	lans_gausspdf1
%
%	(C) 1999.10.13 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______________________________________________________________
%	auto-detection of diagonal and spherical by examining full matrix

function	[pdf]	= lans_gausspdf(x,mu,sig,nonorm,alpha,EET)

if nargin>0
%__________ REGULAR ____________________________________________________________
[D N]	= size(x);

xm	= x-mu*ones(1,N);			% mean centered
xm1	= xm';					% transposed
if nargin<5
	if length(sig)==1				% 1 spherical
		if nargin<4
			fac	= 1/sqrt((2*pi*sig)^D);
		end
		raised	= vdist2(xm)/sig;	
	elseif min(size(sig))==1			% 2 diagonal
		if nargin<4
			fac	= 1/sqrt((2*pi)^D*prod(sig));
		end
		isig	= diag(diag(1./sqrt(sig)));
		raised	= vdist2(xm.*(isig*ones(1,N)));	
	else						% 3 full
		if nargin<4
			fac	= 1/sqrt((2*pi)^D*det(sig));
		end
		isig	= inv(sig);
		raised	= sum((xm1*isig).*xm1,2)';
	end
else							% 4 clamped covariance
	if nargin==5
		EET	= alpha;
		alpha	= nonorm;
	end
	
	Q	= rank(EET);

	if D>Q
		B		= (D-alpha*Q)*sig/(D-Q);
		S		= alpha*sig;
		sigfull		= B*eye(D) + (S-B)*EET;
		sigfulldet	= (S^Q)*(B^(D-Q));
		isig		= eye(D)/B - (S-B)*EET/(B*S);
		if nargin==5
			fac		= 1/sqrt((2*pi)^D*sigfulldet);
		end
		raised		= sum((xm1*isig).*xm1,2)';
	else	% reverts to unit spherical variance
		if nargin==5
			fac	= 1/sqrt((2*pi)^D);
		end
		raised	= vdist2(xm);	
	end
end
if (nargin==3)|(nargin==6)
	pdf	= fac*exp(-.5*raised);
else
	pdf	= exp(-.5*raised);
end

%__________ REGULAR ends _______________________________________________________
else
%__________ DEMO _______________________________________________________________
clc
disp('running lans_gausspdf.m in demo mode');
N	= 10;
D	= 4;
mu	= [1:D]';
sig	= 2*eye(D);
x	= rand(D,N);

flops(0)
lans_gausspdf(x,mu,sig)
dstr	= sprintf('Full used %d flops',flops);
disp(dstr);

flops(0)
lans_gausspdf(x,mu,diag(sig))
dstr	= sprintf('Diagonal used %d flops',flops);
disp(dstr);

flops(0)
lans_gausspdf(x,mu,sig(1))
dstr	= sprintf('Spherical used %d flops',flops);
disp(dstr);

%__________ DEMO ends __________________________________________________________
end