www.gusucode.com > 超全的模式识别matlab源码程序 > code/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