【轉】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 ([email protected]).
%% 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 ([email protected]).
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