bayesian_model_comparison.m

数据挖掘matlab源码

function D = Bayesian_Model_Comparison(train_features, train_targets, Ngaussians, region)

% 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:
% 	features	- Train features
%	targets	    - Train targets
%	Ngaussians  - Number of redraws
%	region	    - Decision region vector: [-x x -y y number_of_points]
%
% Outputs
%	D			- Decision sufrace
%
% Strongly built for only two classes!

[Ndim, M]	= size(train_features);
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_features = train_features(:,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 decision region and error for each possible model. Use only EMindices data
   G						= min([ones(1,length(Ngaussians))*i; Ngaussians]);
   [D, param]	= EM(train_features(:,EMindices), train_targets(:,EMindices), G, region);
   
   %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_features(:,i0), param.m0, param.s0, param.w0);      
   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_features(:,i1), param.m1, param.s1, param.w1);      
   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...

D        = EM(train_features, train_targets, [i1, i2], region);

disp(['FINAL SELECTION: Using ' num2str(i1) ' Gaussians for class 1 and ' num2str(i2) ' Gaussians for class 2'])


function P = computeML(features, mu, sigma, w)

M	= size(features,2);
Ng	= size(mu,1);
p	= zeros(Ng, M);

warning off

if Ng == 1,
	for j = 1:M,
   	x		= features(:,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		= features(:,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