MATLAB实现高斯混合分布的EM算法及二维时概率密度曲面、置信椭圆绘制

具体EM算法并不打算介绍太多，详细公式及证明可以参考文末链接的参考内容

EM算法

1 定义分量数目 K K K,对每个分量 k k k设置 π k , μ k , Σ k \pi_k,\boldsymbol{\mu}_{k},\boldsymbol{\Sigma}_{k} πk​,μk​,Σk​的初始值。 2 E step 根据当前的 π k , μ k , Σ k \pi_k,\boldsymbol{\mu}_{k},\boldsymbol{\Sigma}_{k} πk​,μk​,Σk​计算后验概率 γ ( z n k ) \gamma(z_{nk}) γ(znk​)。 γ ( z n k ) = π k N ( x n ∣ μ n , Σ n ) ∑ j = 1 K π j N ( x n ∣ μ j , Σ j ) \gamma\left(z_{n k}\right)=\frac{\pi_{k} \mathcal{N}\left(\boldsymbol{x}_{n} \mid \boldsymbol{\mu}_{n}, \boldsymbol{\Sigma}_{n}\right)}{\sum_{j=1}^{K} \pi_{j} \mathcal{N}\left(\boldsymbol{x}_{n} \mid \boldsymbol{\mu}_{j}, \boldsymbol{\Sigma}_{j}\right)} γ(znk​)=∑j=1K​πj​N(xn​∣μj​,Σj​)πk​N(xn​∣μn​,Σn​)​ 3 M step 根据 E step 中计算的 γ ( z n k ) \gamma(z_{nk}) γ(znk​)再计算新的 π k , μ k , Σ k \pi_k,\boldsymbol{\mu}_{k},\boldsymbol{\Sigma}_{k} πk​,μk​,Σk​： μ k n e w = 1 N k ∑ n = 1 N γ ( z n k ) x n Σ k n e w = 1 N k ∑ n = 1 N γ ( z n k ) ( x n − μ k n e w ) ( x n − μ k n e w ) T π k n e w = N k N \begin{aligned} \boldsymbol{\mu}_{k}^{n e w} &=\frac{1}{N_{k}} \sum_{n=1}^{N} \gamma\left(z_{n k}\right) \boldsymbol{x}_{n} \\ \boldsymbol{\Sigma}_{k}^{n e w} &=\frac{1}{N_{k}} \sum_{n=1}^{N} \gamma\left(z_{n k}\right)\left(\boldsymbol{x}_{n}-\boldsymbol{\mu}_{k}^{n e w}\right)\left(\boldsymbol{x}_{n}-\boldsymbol{\mu}_{k}^{n e w}\right)^{T} \\ \pi_{k}^{n e w} &=\frac{N_{k}}{N} \end{aligned} μknew​Σknew​πknew​​=Nk​1​n=1∑N​γ(znk​)xn​=Nk​1​n=1∑N​γ(znk​)(xn​−μknew​)(xn​−μknew​)T=NNk​​​ 其中： N k = ∑ n = 1 N γ ( z n k ) N_{k}=\sum_{n=1}^{N} \gamma\left(z_{n k}\right) Nk​=n=1∑N​γ(znk​) 4 检查参数是否收敛或对数似然函数是否收敛，若不收敛，则返回第2步。 对数似然函数： ln ⁡ p ( x ∣ π , μ , Σ ) = ∑ n = 1 N ln ⁡ { ∑ k = 1 K π k N ( x k ∣ μ k , Σ k ) } \ln p(\boldsymbol{x} \mid \boldsymbol{\pi}, \boldsymbol{\mu}, \boldsymbol{\Sigma})=\sum_{n=1}^{N} \ln \left\{\sum_{k=1}^{K} \pi_{k} \mathcal{N}\left(\boldsymbol{x}_{k} \mid \boldsymbol{\mu}_{k}, \boldsymbol{\Sigma}_{k}\right)\right\} lnp(x∣π,μ,Σ)=n=1∑N​ln{k=1∑K​πk​N(xk​∣μk​,Σk​)}

EM算法MATLAB实现（附带详细注释）

此EM算法代码利用大量矩阵运算，和反复转置，减小了中间变量的大小，显著提高效率。

function [Mu,Sigma,Pi,Class]=gaussKMeans(pntSet,K,initM)
% @author:slandarer
% ===============================================================
% pntSet  | NxD数组   | 点坐标集                                |
% K       | 数值      | 划分堆数量                              |
% --------+-----------+-----------------------------------------+
% Mu      | KxD数组   | 每一行为一类的坐标中心                  |
% Sigma   | DxDxK数组 | 每一层为一类的协方差矩阵                |
% Pi      | Kx1列向量 | 每一个数值为一类的权重(占比)            |
% Class   | Nx1列向量 | 每一个数值为每一个元素的标签(属于哪一类)|
% --------+-----------+-----------------------------------------+


[N,D]=size(pntSet); % N:元素个数 | D:维数

% 初始化数据===============================================================
if nargin