function [mask P M S] = gmm_classifier (x,K)

%PUBLISHED UNDER THE CREATIVE COMMONS LICENSE (SEE MY WEBPAGE FOR MORE DETAILS)

%INPUT:
%x is the input vector where clustering is performed, x in R
%K is the number of clusters

%OUTPUT:
%mask is the classification of points of x => mask(i) is the category of point x(i), categories in {1,2,3,...,K}
%P: priors of mixture model
%M: centers of clusters
%S: covariances of clusters

%to provide clusters' centers' initialization set the global variable MU_INIT
%MU_INIT is an Kx1 matrix
global MU_INIT;

EPSILON = 1E-06;
ITER    = 20; %EM ITTERATIONS UPPER BOUND

N = size (x,1);

p1 = zeros (N,K); p_post = zeros (N,K);

P0 = zeros (K,1);
M0 = zeros (K,1);
S0 = zeros (K,1);

if isempty(MU_INIT)
    %INITIALIZE WITH K-MEANS
    mask = kmeans (x, K, 'emptyaction','singleton','display','off',...
        'start','cluster','replicates',20);
    for i=1:K
        ind = find (mask == i);
        if isempty (ind)
            P0(i) = eps;
            M0(i) = 0;
            S0(i) = 1000;
        else
            P0(i) = length(ind) / N;
            M0(i) = mean (x(ind));
            S0(i) = std (x(ind));
        end
    end
else
    M0 = MU_INIT;
    S0 = ones(K,1);
    P0 = (1/K) * ones (K,1);
end

%EM Loop
Qnew = inf;
for iter=1:ITER
    %Expectation
    for j=1:K
        p1(:,j) = normpdf (x, M0(j), S0(j))*P0(j);
    end
    px = sum (p1,2);
    px (find(px == 0)) = eps;

    for j=1:K
        p_post (:,j) = p1(:,j) ./ px;
    end
    
    P = P0;
    M = M0;
    S = S0;
    
    %convergence condition
    Qold = Qnew; Qnew = - sum (log(px));
    if Qold - Qnew < EPSILON
        break;
    end
    
    %Maximization
    P0 = (sum(p_post,1))';
    P0(find(P0 == 0)) = eps;
    for j=1:K
        M0(j) = sum(p_post(:,j).*x,1) / P0(j);
        S0(j) = sqrt (sum (p_post(:,j).*((x - M0(j)) .^ 2)) /P0(j));
        if S0(j) < 1.0
            P0(j) = eps;
            S0(j) = 1000;
        end
    end
    P0 = P0 / N;
end


% normalization of priors
P = P / sum(P);

for j=1:K
    p1(:,j) = normpdf (x, M(j), S(j)) * P(j);
end
px = sum (p1,2);
px (find(px == 0)) = eps;

for j=1:K
    p_post (:,j) = p1(:,j) ./ px;
end
[q mask] = max (p_post,[],2);