function [test_targets, P] = Marginalization(train_patterns, train_targets, test_patterns, params)
% Classify data with missing features using the marginal distribution
%
% Inputs:
% train_patterns - Training patterns
% train_targets - Training targets
% test_patterns - Test patterns
% params - [The number of the missing feature, Number of bins for the data]
%
% Outputs
% targets - Output targets
% P - P(w|x) (Probability of a class given the good features)
%warning off;
[missing, N] = process_params(params);
[d, L] = size(train_patterns);
if (missing > d),
error(['The number of the missing pattern must be between 1 and ' num2str(d)])
end
not_missing = find(~ismember(1:d, missing));
%Calculate the marginal distribution using histograms
Nbins = max(3,floor(size(train_patterns,2).^(1/3)));
[p, b, region] = high_histogram(train_patterns,Nbins);
classes = unique(train_targets);
Ptrain = zeros(length(classes), Nbins^(d-1));
for i = 1:length(classes),
indices = find(train_targets == classes(i));
g_i = high_histogram(train_patterns(:,indices),Nbins,region);
Pc = squeeze(sum(g_i.*p,missing)./sum(p,missing));
bad = find(~isfinite(Pc));
Pc(bad) = 0;
Ptrain(i,:) = Pc(:)';
end
if (d > 2)
Ptrain = reshape(Ptrain, [length(classes) ones(1,d-1)*d]);
end
% Classify the marginal probabilities for all the patterns
[temp, binned_patterns] = high_histogram(test_patterns, Nbins, region);
P = zeros(length(classes), size(test_patterns,2));
for i = 1:size(P,2),
s = ':,';
for j = 1:d-1,
s = [s, num2str(binned_patterns(not_missing(j),i))];
if (j ~= d-1)
s = [s, ','];
end
end
s = [s, ''];
P(:,i) = eval(['Ptrain(', s, ')']);
end
[m, I] = max(P);
test_targets = classes(I);