www.pudn.com > Classification_toolbox.part01.rar > Bayesian_Model_Comparison.m

function test_targets = Bayesian_Model_Comparison(train_patterns, train_targets, test_patterns, Ngaussians) 
 
% Classify using the Bayesian model comparison algorithm. This function accepts as inputs 
% the maximum number of Gaussians per class and returns a decision surface based on the  
% most likely number of Gaussians in each class 
% 
% Inputs: 
% 	training_patterns   - Train patterns 
%	training_targets	- Train targets 
%   test_patterns       - Test  patterns 
%	Ngaussians          - Number of redraws 
% 
% Outputs 
%	test_targets        - Predicted targets 
% 
% Strongly built for only two classes! 
 
[Ndim, M]	= size(train_patterns); 
P_D_given_h	= zeros(2, max(Ngaussians)+2); 
 
%Use holdout for diffrentiating between training data for finding the gaussian parameters 
%and the likelihood 
holdout			    = 0.5; 
i					= randperm(M); 
train_patterns      = train_patterns(:,i); 
train_targets	    = train_targets(i); 
EMindices		    = 1:floor(M*holdout); 
MLindices		    = floor(M*holdout)+1:M; 
i0					= MLindices(find(train_targets(MLindices) == 0)); 
i1					= MLindices(find(train_targets(MLindices) == 1)); 
Ni0				    = length(i0); 
Ni1				    = length(i1); 
 
for i = 1:max(Ngaussians), 
    
   %Calculate the error for each possible model. Use only EMindices data 
   G			= min([ones(1,length(Ngaussians))*i; Ngaussians]); 
   [D, param]	= EM(train_patterns(:,EMindices), train_targets(:,EMindices), train_patterns(:,EMindices), G); 
    
   %Calculate likelihood of the data given these Gaussians 
   %Use only the MLindices data 
   if (P_D_given_h(1, G(1)) == 0),	%Do it only if it wasn't already computed 
		P_D_given_h(1, G(1)) = computeML(train_patterns(:,i0), param(1).mu, param(1).sigma, param(1).w);       
   end 
   if (P_D_given_h(2, G(2)) == 0),	%Do it only if it wasn't already computed 
		P_D_given_h(2, G(2)) = computeML(train_patterns(:,i1), param(2).mu, param(2).sigma, param(2).w);       
   end 
    
end 
P_D_given_h(find(isnan(P_D_given_h)))	=	0; 
P_D_given_h = P_D_given_h./(eps+sum(P_D_given_h')'*ones(1,size(P_D_given_h,2))); %Normalize 
 
%Compute the Hessian for each class 
H 						= diff([zeros(2,2), P_D_given_h], 2, 2); 
 
P_D_given_h			= P_D_given_h .* (abs(H).^(-0.5)) .* (2*pi); 
 
likelihood = P_D_given_h(1,[1:Ngaussians(1)])' * P_D_given_h(2,[1:Ngaussians(2)]); 
 
%Choose the ML model as the one with the lowest error 
[i1, i2] = find(likelihood == max(max(likelihood))); 
if isempty(i1), 
    error('Could not find a likely pair.') 
end 
 
i1       = i1(1); i2 = i2(1); %To give preference for simpler models... 
 
test_targets = EM(train_patterns, train_targets, test_targets, [i1, i2]); 
 
disp(['FINAL SELECTION: Using ' num2str(i1) ' Gaussians for class 1 and ' num2str(i2) ' Gaussians for class 2']) 
 
 
function P = computeML(patterns, mu, sigma, w) 
 
M	= size(patterns,2); 
Ng	= size(mu,1); 
p	= zeros(Ng, M); 
 
warning off 
 
if Ng == 1, 
	for j = 1:M, 
   	x		= patterns(:,j); 
      p(j)	= w(1)/(2*pi*sqrt(det(sigma)))*exp(-0.5*(x-mu')'*inv(sigma)*(x-mu')); 
   end 
   P			= prod(p); 
else 
   for j = 1:M, 
   	x		= patterns(:,j); 
      for k = 1:length(w), 
         p(k, j) = w(k)/(2*pi*sqrt(det(squeeze(sigma(k,:,:)))))*... 
            exp(-0.5*(x-mu(k,:)')'*inv(squeeze(sigma(k,:,:)))*(x-mu(k,:)')); 
      end 
   end 
   P			= prod(sum(p));    
end