CODEWITHSUNDEEP
algorithmmatlabgaussianmixture-model
evolved
evolved

Reputation: 2210

Matlab: EM for Gaussian Mixture Models without gmdistribution

I have to train a Gaussian Mixture model using four components on a given dataset. The set is three dimensional and contains 300 samples.

The problem is that I cannot check for convergence using the log-likelihood because it is -Inf. This results from rounded zero values while evaluating the gaussian in the formula of the responsibilities (see E-step).

Can you tell me if my implementation of the EM algorithm is correct so far? And how to account for the problem with the rounded zero values?

Here is my implementation of the EM algorithm (one iteration):

First I initialized the means and the covariance of the components using kmeans:

load('data1.mat');

X = Data'; % 300x3 data set
D = size(X,2); % dimension
N = size(X,1); % number of samples
K = 4; % number of Gaussian Mixture components

% Initialization
p = [0.2, 0.3, 0.2, 0.3]; % arbitrary pi
[idx,mu] = kmeans(X,K); % initial means of the components

% compute the covariance of the components
sigma = zeros(D,D,K);
for k = 1:K
    sigma(:,:,k) = cov(X(idx==k,:));
end

For the E-step I am using the following formula to calculate the responsibilities responsibility

Here is the corresponding code:

gm = zeros(K,N); % gaussian component in the nominator - 
                 % some values evaluate to zero
sumGM = zeros(N,1); % denominator of responsibilities
% E-step:  Evaluate the responsibilities using the current parameters
% compute the nominator and denominator of the responsibilities
for k = 1:K
    for i = 1:N
        % HERE values evalute to zero e.g. exp(-746.6228) = -Inf
        gm(k,i) = p(k)/sqrt(det(sigma(:,:,k))*(2*pi)^D)*exp(-0.5*(X(i,:)-mu(k,:))*inv(sigma(:,:,k))*(X(i,:)-mu(k,:))');
        sumGM(i) = sumGM(i) + gm(k,i);
    end
end
res = zeros(K,N); % responsibilities
Nk = zeros(4,1);
for k = 1:K
    for i = 1:N
         res(k,i) = gm(k,i)/sumGM(i);
    end
    Nk(k) = sum(res(k,:));
end

Nk(k) is computed using the formula given in the M-step.

M-step

reestimate parameters using current responsibilities

% M-step: Re-estimate the parameters using the current responsibilities
mu = zeros(K,3);
for k = 1:K
    for i = 1:N
        mu(k,:) = mu(k,:) + res(k,i).*X(k,:);
        sigma(:,:,k) = sigma(:,:,k) + res(k,i).*(X(k,:)-mu(k,:))*(X(k,:)-mu(k,:))';
    end
    mu(k,:) = mu(k,:)./Nk(k);
    sigma(:,:,k) = sigma(:,:,k)./Nk(k);
    p(k) = Nk(k)/N;
end

Now in order to check for convergence the log-likelihood is computed using this formula: log-likelihood

% Evaluate the log-likelihood and check for convergence of either 
% the parameters or the log-likelihood. If not converged, go to E-step.
loglikelihood = 0;
for i = 1:N
    for k = 1:K
         loglikelihood = loglikelihood + log(gm(k,i));
    end
end

loglikelihood is -Inf because some gm(k,i) values in the E-step are zero. Therefore the log is obviously negative infinity.

How can I solve this problem?

Could it be solved by increasing Matlab's precision?

Or is there something wrong with my implementation?

Upvotes: 1

Views: 1056

Answers (1)

Eskapp
Eskapp

Reputation: 3635

According to the formula, you should compute the logarithm of the sum of the gm quantities. (so, log(sum(gm(i,:))) ). Within the k components, at least one will have a likelihood greater than 0. This will solve your problem hopefully.

Another very general remark, when numbers are too big for functions as the exponential one, and when you are sure you are using the right formula, you can always try to work with the log of the quantities. But you should not need to do this here as 0 is a good approximation of exp(-746) ;)

Upvotes: 2

Related Questions