Code covered by the BSD License

# EM algorithm for Gaussian mixture model with background noise

by

### Andrew (view profile)

16 May 2012 (Updated )

Standard EM algorithm to fit a GMM with the (optional) consideration of background noise.

GM_EM(X,k,bn_noise,reps,max_iters,tol)
```function [idx,mu] = GM_EM(X,k,bn_noise,reps,max_iters,tol)

% GM_EM - fit a Gaussian mixture model to N points located in n-dimensional
% space.
%
% Note: This function requires the Statistical Toolbox and, if you wish to
% plot (for k = 2), the function error_ellipse
%
% Elementary usage:
% GM_EM(X,k) -  fit a GMM to X, where X is N x n and k is the number of
%               clusters. Algorithm follows steps outlined in Bishop
%               (2009) 'Pattern Recognition and Machine Learning', Chapter 9.
%
%   bn_noise - allow for uniform background noise term ('T' or 'F',
%   default 'T'). If 'T', relevant classification uses the
%   (k+1)th cluster
%   reps - number of repetitions with different initial conditions
%   (default = 10). Note: only the best fit (in a likelihood sense) is
%   returned.
%   max_iters - maximum iteration number for EM algorithm (default = 100)
%   tol - tolerance value (default = 0.01)
%
% Outputs
%   idx - classification/labelling of data in X
%   mu - GM centres

% Preamble
if nargin < 6
tol = 0.01;
end
if nargin < 5
max_iters = 100;
end
if nargin < 4
reps = 10;
end
if nargin < 3
bn_noise = 'T';
end
if nargin < 2
k = 2;
end
if nargin  == 0
X = [mvnrnd([1;1],eye(2),20);mvnrnd([10;10],eye(2),20);10*rand(20,2)];
end

% Parameters
n = size(X,2);                      % dimensionality
N = size(X,1);                      % sample size
dom_area=  prod(max(X)-min(X));     % domain area
mu_save = zeros(n,k,reps);          % all mu
Sigma_save = zeros(n,n,k,reps);     % all Sigma
log_lk_save = zeros(reps,1);        % all log-likelihoods
if bn_noise == 'T'
pi_save = zeros(reps,k+1);      % all mixture weights
else
pi_save = zeros(reps,k);
end

for rep_num = 1:reps

% Initialise
if bn_noise == 'T'
pi_weight = repmat(0.99/k,1,k);
pi_weight(k+1) = 1-sum(pi_weight(1:k));
else
pi_weight = repmat(1/k,1,k);
end

Sigma = repmat(abs(max(max(X)))*eye(n),[1,1,k]);    % random initialisations
mu = max(max(X))*rand(n,k);
log_lk = zeros(max_iters,1);

% Until max_iters of convergence is reached
for m = 1:max_iters
mvnpdfs = zeros(N,1);
for j = 1:k
[~,err] = cholcov(Sigma(:,:,j),0);          % Check Cholesky
if err ~= 0
break
end
mvnpdfs(:,j) = mvnpdf(X,mu(:,j)',Sigma(:,:,j));
end
if err ~= 0                                     % If Cholesky fails, abandon this attempt
break
end
if bn_noise == 'T'
mvnpdfs(:,j+1) = 1/dom_area;                % Final "mixture" is a uniform distribution
end
log_lk(m) = sum(log(sum(repmat(pi_weight,N,1).*mvnpdfs,2)));  % Compute log-lk
if (m > 1) && (abs(log_lk(m) - log_lk(m-1)) < tol)
disp(['Rep ',num2str(rep_num),': Convergence reached in ',num2str(m),' iterations'])
break
end

den = sum(repmat(pi_weight,N,1).*mvnpdfs,2);        % Follow steps in Bishop, Ch. 9

% Find mu_k
gamma_znk = repmat(pi_weight(1:k),N,1).*mvnpdfs(:,1:k)./repmat(den,1,k);

for j = 1:k
Nk = sum(gamma_znk(:,j));
mu(:,j) =  sum(X.*repmat(gamma_znk(:,j),1,n))/Nk;  % Update mu

tot_sum = (X'-repmat(mu(:,j),1,N))*diag(gamma_znk(:,j))*(X'-repmat(mu(:,j),1,N))';
Sigma(:,:,j) = tot_sum/Nk;                         % Update Sigma
pi_weight(j) = Nk/N;                               % Update weights
end

if bn_noise == 'T'
pi_weight(k+1) = 1- sum(pi_weight(1:k));
end

end
if m == max_iters
disp(['Warning: Rep ',num2str(rep_num),': Maximum number of iterations reached'])
end
if err ~= 0
disp('Warning: chol failed, algorithm abandoned');
end
log_lk_save(rep_num) = log_lk(m);
mu_save(:,:,rep_num) = mu;
pi_save(rep_num,:) = pi_weight;
Sigma_save(:,:,:,rep_num) = Sigma;      % Save last values of M-step
end

%Choose best fit
best_rep = find(log_lk_save == max(log_lk_save));
mu = mu_save(:,:,best_rep);
Sigma = Sigma_save(:,:,:,best_rep);
pi_weight = pi_save(best_rep,:);

% Now that we have our underlying prob. distribution we can classify
if bn_noise == 'T'
response = zeros(N,k+1);
else
response = zeros(N,k);
end
for i = 1:k
response(:,i) = pi_weight(i)*mvnpdf(X,mu(:,i)',Sigma(:,:,i));
end
if bn_noise == 'T'
response(:,i+1) = pi_weight(i)/dom_area;
end

[~,idx] =max(response,[],2);

% Plot if in 1 or 2-dimensions
str_color = ['r','g','b','c','m','y'];
if (n == 1) && (k < 6);
figure
stem(X(:,1),ones(N,1),'k'); hold on;
for i = 1:k
stem(X(idx==i,1),ones(length(X(idx==i,1)),1),str_color(i)); hold on;
stem(mu(i),2,str_color(i),'Linewidth',3)
end
set(gca,'Ylim',[0 3])
end

if (n == 2) && (k < 6);
figure
str_color = ['r','g','b','c','m','y'];
scatter(X(:,1),X(:,2),'kx'); hold on;
for i = 1:k
scatter(X(idx==i,1),X(idx==i,2),[str_color(i),'x']);
plot(mu(1,i),mu(2,i),[str_color(i),'o'],'Linewidth',2)
try
error_ellipse(Sigma(:,:,i),mu(:,i),'style',str_color(i),'conf',0.95)
catch exception
if strcmp(exception.identifier,'MATLAB:UndefinedFunction')
end
end
end
end

% Return transpose of mu to be compatible with other algorithms (e.g.
% kmeans)
mu = mu';

```