Code covered by the BSD License  

Highlights from
EM algorithm for Gaussian mixture model with background noise

image thumbnail

EM algorithm for Gaussian mixture model with background noise

by

 

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.
%
% Additional inputs:
%   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')
                disp('Warning: ellipse not plotted (please download error_ellipse.m)')
            end
        end
    end
end

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

Contact us