【轉】EM算法MATLAB代碼及詳細注解
版權聲明：本文為博主原創文章，遵循 CC 4.0 BY-SA 版權協議，轉載請附上原文出處鏈接和本聲明。
本文鏈接：https://blog.csdn.net/qq_23968185/article/details/70940197
作者：鼴鼠的胡須
覺得有用的話,歡迎一起討論相互學習~
下面代碼為PRML所附的基於混合高斯（MoG）的代碼，個人認為編碼可讀性和風格都值得借鑒。
function [label, model, llh] = mixGaussEm(X, init)
% Perform EM algorithm for fitting the Gaussian mixture model.
% Input:
% X: d x n data matrix
% init: k (1 x 1) number of components or label (1 x n, 1<=label(i)<=k) or model structure
% Output:
% label: 1 x n cluster label
% model: trained model structure
% llh: loglikelihood
% Written by Mo Chen (sth4nth@gmail.com).
%% init
fprintf('EM for Gaussian mixture: running ... \n');
tol = 1e-6;
maxiter = 500;
llh = -inf(1,maxiter);
R = initialization(X,init);
for iter = 2:maxiter
[~,label(1,:)] = max(R,[],2); %label表示1*n的類別向量
R = R(:,unique(label)); % remove empty clusters，unique輸出label向量中不重復的元素，表示非空類別向量
model = maximization(X,R); %EM算法的M step，X表示數據矩陣，R表示類別矩陣，model是結構體，表示模型，其屬性是模型的參數
[R, llh(iter)] = expectation(X,model); %EM算法的E step，X表示數據矩陣，model表示模型結構體，R表示返回的隸屬度矩陣，llh表示似然函數的目標值
if abs(llh(iter)-llh(iter-1)) < tol*abs(llh(iter)); break; end;
end
llh = llh(2:iter);
function R = initialization(X, init) %X是數據矩陣，init用於初始化MoG的成分，R返回的是一個n行k列的矩陣，第ij個元素表示第i個樣本由第j個成分生成的概率
n = size(X,2); %n是樣本個數
if isstruct(init) % init with a model %isstruct判斷輸入是否是一個matlab結構體
R = expectation(X,init); %如果init是一個結構體，直接用該模型進行E step
elseif numel(init) == 1 %如果init是一個整數
k = init; %用init表示混合成分的個數，即類別個數
label = ceil(k*rand(1,n)); %ceil用於向數軸的正方向取整，初始化樣本的label
R = full(sparse(1:n,label,1,n,k,n)); %sparse通過記錄稀疏矩陣非負元素的索引和值來節省內存，full是一相反作用；R是n行k列矩陣，n表示樣本個數，k表示類別數，每一行
%是一個one-hot向量，表示該樣本屬於哪一類
elseif all(size(init)==[1,n]) % init with labels %若init是一個一行n列的向量，則為樣本類別的向量
label = init;
k = max(label);
R = full(sparse(1:n,label,1,n,k,n));
else
error('ERROR: init is not valid.');
end
%EM算法的E step，X表示數據矩陣，model表示模型結構體，R表示返回的隸屬度矩陣，llh表示似然函數的目標值
function [R, llh] = expectation(X, model)
mu = model.mu;
Sigma = model.Sigma;
w = model.w; %w為MoG的混合系數向量
n = size(X,2); %n為樣本個數
k = size(mu,2); %k為MoG混合成分的個數，即類別個數
R = zeros(n,k); %R隸屬度矩陣，行數為樣本個數，列數為類別個數，第ij個元素表示第i個樣本由第j個成分生成的概率
for i = 1:k %計算樣本的每個gauss概率的對數
R(:,i) = loggausspdf(X,mu(:,i),Sigma(:,:,i));
end
R = bsxfun(@plus,R,log(w)); %計算隸屬度（未歸一化）矩陣的對數
T = logsumexp(R,2); %對R取指數加和再取對數
llh = sum(T)/n; % loglikelihood %似然函數的均值
R = exp(bsxfun(@minus,R,T)); %計算隸屬度矩陣
%EM算法的M step，X表示數據矩陣，R表示隸屬度矩陣，第ij個元素表示第i個樣本由第j個成分生成的概率，model是結構體，表示模型，其屬性是模型的參數
function model = maximization(X, R)
[d,n] = size(X); %d表示樣本維數，n表示樣本個數
k = size(R,2); %k表示MoG成分的個數
nk = sum(R,1); %nk表示求隸屬度矩陣R的列和
w = nk/n; %w表示混合成分系數
mu = bsxfun(@times, X*R, 1./nk); %mu是一個m行k列的矩陣，表示k個高斯成分的期望，每個都是m元隨機變量
Sigma = zeros(d,d,k); %Sigma是一個三維張量，表示第k個高斯成分的協方差矩陣是d*d的
r = sqrt(R);
for i = 1:k %循環計算每個成分的協方差
Xo = bsxfun(@minus,X,mu(:,i));
Xo = bsxfun(@times,Xo,r(:,i)');
Sigma(:,:,i) = Xo*Xo'/nk(i)+eye(d)*(1e-6);
end
model.mu = mu;
model.Sigma = Sigma;
model.w = w;
function y = loggausspdf(X, mu, Sigma) %計算Gauss概率分布函數的對數的函數，輸入變量分別為數據X，期望mu，期望協方差Sigma
d = size(X,1); %d表示樣本維數
X = bsxfun(@minus,X,mu); %樣本與均值作差
[U,p]= chol(Sigma); %chol表示將協方差矩陣Sigma進行一個上三角矩陣分解，U表示上三角因子矩陣，Sigma=U'的逆與U作積（將協方差矩陣分解求逆加快計算效率）
if p ~= 0 %如果p不為0則Sigma不是正定矩陣，報錯
error('ERROR: Sigma is not PD.');
end
Q = U'\X; %Q=U'的逆與X的乘積
q = dot(Q,Q,1); % quadratic term (M distance) %dot表示點乘之后求列和
c = d*log(2*pi)+2*sum(log(diag(U))); % normalization constant
y = -(c+q)/2;
function s = logsumexp(X, dim)
% Compute log(sum(exp(X),dim)) while avoiding numerical underflow.
% By default dim = 1 (columns).
% Written by Mo Chen (sth4nth@gmail.com).
if nargin == 1,
% Determine which dimension sum will use
dim = find(size(X)~=1,1);
if isempty(dim), dim = 1; end
end
% subtract the largest in each dim
y = max(X,[],dim);
s = y+log(sum(exp(bsxfun(@minus,X,y)),dim)); % TODO: use log1p
i = isinf(y);
if any(i(:))
s(i) = y(i);
end