www.gusucode.com > 超全的模式识别matlab源码程序 > code/Gibbs.m
function test_targets = Gibbs(train_patterns, train_targets, test_patterns, Ndiv) % Classify using the Gibbs algorithm % Inputs: % train_patterns - Train patterns % train_targets - Train targets % test_patterns - Test patterns % Ndiv - Resolution of the distribution of Theta given the data % % Outputs % test_targets - Predicted targets [Ni, M] = size(train_patterns); Uclasses = unique(train_targets); Nc = length(Uclasses); mu = zeros(Nc, Ni, Ndiv); sigma = zeros(Nc, Ni, Ni, Ndiv); %We will assume diagonal matrices for simplicity for i = 1:Nc, in = find(train_targets == Uclasses(i)); %Find the range of the data min_range = min(train_patterns(:,in)')'; max_range = max(train_patterns(:,in)')'; %Make a "guess" on theta based on this range of data for j = 1:Ni, mu(i, j, :) = linspace(min_range(j), max_range(j),Ndiv); sigma(i,j,j,:) = linspace(min_range(j), max_range(j),Ndiv); end end %Now test the probability of the test set given these parameters (P(D|theta)) P_D_given_Theta = zeros(Nc, Ndiv, Ndiv); X = train_patterns; for i = 1:Nc, in = find(train_targets == Uclasses(i)); x = X(:,in); pD_given_theta = zeros(1, length(in)); m = length(in); for j = 1:Ndiv, for k = 1:Ndiv, temp_mu = squeeze(mu(i, :, j))'; temp_si = squeeze(sigma(i, :, :, k)); pD_given_theta = 1/((2*pi*det(temp_si))^(Ni/2)*exp(-0.5*diag((x-temp_mu*ones(1,m))'*inv(temp_si)*(x-temp_mu*ones(1,m))))); P_D_given_Theta(i,j,k) = sum(pD_given_theta); %P_D_given_Theta(i,j,k) = prod(pD_given_theta(find(pD_given_theta~=0))); end end end %Find P(Theta|D), and draw one using it P_Theta_given_D = zeros(size(P_D_given_Theta)); Theta_mu = zeros(Nc, Ni, 1); Theta_sigma = zeros(Nc, Ni, Ni); for i = 1:Nc, P_Theta_given_D(i,:,:) = squeeze(P_D_given_Theta(i,:,:))/sum(sum(squeeze(P_D_given_Theta(i,:,:)))); %Make a draw n = rand(1); seq = cumsum(reshape(squeeze(P_Theta_given_D(i,:,:)),Ndiv^2,1)); loc = max(find(n > seq))+1; if isempty(loc) indices(i) = 1; else indices(i) = loc; end seq = zeros(size(seq)); seq(loc)=1; [j,k] = find(reshape(seq, Ndiv, Ndiv)); Theta_mu(i, :, 1) = squeeze(mu(i, :, j))'; Theta_sigma(i, :, :) = squeeze(sigma(i, :, :, k)); end %Classify test patterns p = hist(train_targets, Uclasses) / length(train_targets); for i = 1:Nc, param_struct(i).mu = squeeze(Theta_mu(i, :, :)); param_struct(i).sigma = squeeze(Theta_sigma(i, :, :)); param_struct(i).p = p(i); param_struct(i).w = 1/Nc; param_struct(i).type = 'Gaussian'; end test_targets = classify_paramteric(param_struct, test_patterns);