功能描述：
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
联系：highspeedlogic
QQ ：1224848052
微信：HuangL1121
邮箱：1224848052@qq.com
网站：http://www.mat7lab.com/
网站：http://www.hslogic.com/
微信扫一扫：
|
