lans_posterior.m

模式识别工具包

%	lans_posterior	- Compute posterior responsibility matrix/log likelihood
%
%	[R,Lcomp,Lavg]	= lans_posterior(y,bas,Rold)
%
%	_____OUTPUTS____________________________________________________________
%	R	posterior responsibility matrix (M x N)		(matrix)
%		P(m|n)	= p(y_n | m)*P(m)/sum_h[p(y_n | h)*P(h)]
%			M = # components
%			N = # samples
%	Lcomp	Complete log likelihood				(scalar)
%		sum_n log{sum_m [p(y_n | m)*P(m)]}
%	Lavg	Average log likelihood				(scalar)
%		sum_m sum_n {R_mn^OLD * ln [p(y_n | m)]}
%		
%	_____INPUTS_____________________________________________________________
%	y	data points					(col vectors)
%	bas	component basis	up to M components		(structure)
%		bas.name		bas.para
%		variance		. x M
%		prior			1 x M
%	Rold	Old posterior for computing Lavg		(matrix)
%
%	_____EXAMPLE____________________________________________________________
%
%	_____NOTES______________________________________________________________
%	- NO regularization 
%	- Rold must be specified to correctly compute Lavg, otherwise an
%	  estimate using 'new' R will be used
%	- summing R across each row gives 1 i.e.
%	  sum(R)	= ones(1,N), i.e. contribution from all components to
%	  point n should sum to 1
%	- small +ve numbers < PREC=realmin replaced by 1
%
%	_____SEE ALSO___________________________________________________________
%	lans_basis	lans_pdf
%
%	(C) 1999.12.14 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______________________________________________________________

function	[R,Lcomp,Lavg] = lans_posterior(y,bas,Rold)

PREC	= realmin;		% SMALLEST Positive number on computer

N	= size(y,2);
M	= size(bas.cen,2);

wprior	= lans_findcell('prior',bas.name);
prior	= bas.para{wprior};				% M x 1 col vector

%__________	p(y_n | x_m) for all m,n		% M x N
p_ygivenx= lans_pdf(y,bas);

%__________	p(y_n,x_m) for all m,n			% M x N
if sum(abs(prior-prior(1)))<eps
	% unequal priors
	p_yx	= p_ygivenx.*(prior'*ones(1,N));
else
	% equal priors
	p_yx	= p_ygivenx/M;
end

%__________	p(y_n) for all n			% 1 x N
p_y	= sum(p_yx);
						
%__________	p(x_m | y_n) Posterior/Responsibility	% M x N 
den	= ones(M,1)*p_y;			% M x N	(identical rows)
den	= den + (den==0);
%den	= lans_replace(PREC,1,den);		% replace small values by 1
R	= p_yx./den;

% ok to replace small values with 1 because Lcomp, Lavg is just an output of
% the simulation, which will not affect the workings of the simulation
if nargout>1
	Lcomp	= sum(log(p_y+(p_y<PREC)));
	if nargout>2
		if nargin<3
			% estimate
			Lavg	= sum(sum(R.*log(p_yx+(p_yx<PREC))));
		else
			Lavg	= sum(sum(Rold.*log(p_yx+(p_yx<PREC))));
		end
	end
end