function [label, model, llh] = emgm(X, init)% EM algorithm for Gaussian mixture model%% initializationfprintf('EM for Gaussian mixture: running ... ');R = initialization(X,init);tol = 1e-6;maxiter = 500;llh = -inf(1,maxiter);converged = false;t = 1;while ~converged && t < maxiter t = t+1; model = maximization(X,R); [R, llh(t)] = expectation(X,model); converged = llh(t)-llh(t-1) < tol*abs(llh(t));end[~,label(1,:)] = max(R,[],2);llh = llh(2:t);if converged fprintf('converged in %d steps.\n',t);else fprintf('not converged in %d steps.\n',maxiter);endfunction R = initialization(X, init)[d,n] = size(X);if isstruct(init) % initialize with model R = expectation(X,init);elseif length(init) == 1 % random initialization k = init; idx = randsample(n,k); m = X(:,idx); [~,label] = max(bsxfun(@minus,m'*X,sum(m.^2,1)'/2)); while k ~= unique(label) idx = randsample(n,k); m = X(:,idx); [~,label] = max(bsxfun(@minus,m'*X,sum(m.^2,1)'/2)); end R = full(sparse(1:n,label,1,n,k,n));elseif size(init,1) == 1 && size(init,2) == n % initialize with labels label = init; k = max(label); R = full(sparse(1:n,label,1,n,k,n));elseif size(init,1) == d && size(init,2) > 1 %initialize with only centers k = size(init,2); m = init; [~,label] = max(bsxfun(@minus,m'*X,sum(m.^2,1)'/2)); R = full(sparse(1:n,label,1,n,k,n));else error('ERROR: init is not valid.');endfunction [R, llh] = expectation(X, model)mu = model.mu;Sigma = model.Sigma;w = model.weight;n = size(X,2);k = size(mu,2);R = zeros(n,k);for i = 1:k R(:,i) = loggausspdf(X,mu(:,i),Sigma(:,:,i));endR = bsxfun(@plus,R,log(w));T = logsumexp(R,2);llh = sum(T)/n; % loglikelihoodR = bsxfun(@minus,R,T);R = exp(R);function model = maximization(X, R)[d,n] = size(X);k = size(R,2);sigma0 = eye(d)*(1e-6); % regularization factor for covariances = sum(R,1);w = s/n;mu = bsxfun(@rdivide, X*R, s);Sigma = zeros(d,d,k);for i = 1:k Xo = bsxfun(@minus,X,mu(:,i)); Xo = bsxfun(@times,Xo,sqrt(R(:,i)')); Sigma(:,:,i) = (Xo*Xo'+sigma0)/s(i);endmodel.mu = mu;model.Sigma = Sigma;model.weight = w;
