我们现在考虑一类具有隐变量的统计推断问题，用数学的语言就是说:
1) 我们的总体样本空间$mathcal{S}inmathbb{R}^{M} imesmathbb{R}^{L}$, 其分布的概率密度函数是$P(x,ymid heta)$，其中$xinmathbb{R}^{M},yinmathbb{R}^{L}$, $ heta$是未知的参数。
2) 我们在抽样的时候，由于条件限制，对每一个样本$(x,y)$，$xinmathbb{R}^{M},yinmathbb{R}^{L}$，后L个变量无法观测到，也就是向量$y$无法获悉，只能观测到$x$。我们称前M个变量为可观测变量，x为可观测向量，后L个变量为不可观测变量，y为不可观测向量。
问题1：
我们有数据集合$D=lbrace x_{1},...,x_{N} brace$, 其中$x_{i}inmathbb{R}^{M}$为第i个样本可观测向量，样本之间相互独立，试推断可能的参数$ heta$。
EM原算法:
为了方便起见，我们记大写字母$X=(x_{1},...,x_{N})$表示所有的样本的可观测向量的总体，而$Y=(y_{1},...,y_{N})$表示所有的样本的不可观测向量的总体, 我们自然地将概率密度表示为紧凑的形式：
$$P(X,Ymid heta) riangleq P((x_{1},y_{1}),...,(x_{N},y_{N})mid heta)=prod_{i=1}^{N}P((x_{1},y_{1}),...,(x_{N},y_{N})mid heta)$$
这个时候自然的，我们的似然函数就是：
egin{equation}P(Xmid heta)=int P(X,Ymid heta)dYend{equation}
当$Y$各个分量连续，或者是：
egin{equation}P(Xmid heta)=sum_{Y}P(X,Ymid heta)end{equation}
当$Y$各个分量离散，或者是还有混合的情形，Y部分分量离散部分连续。为了记号方便起见我们假定Y各个分量均离散，其实各种情形下计算过程基本无差异。
现在我们的目标是使得似然函数$P(Xmid heta)$尽量大。注意到一般来说由于表达式(2)中求和号的存在，直接求似然函数的极值比较困难，现在我们想办法设计一种迭代算法使得极大似然函数尽量大。假设我们第$i$步迭代已经得到参数$ heta^{(i)}$, 现在希望迭代更新使得似然函数尽量大。注意到：
egin{split}&frac{P(Xmid heta)}{P(Xmid heta^{(i)})} ewline =&sum_{Y}frac{P(X,Ymid heta)}{P(Xmid heta^{(i)})} ewline =&sum_{Y}frac{P(X,Ymid heta)P(Ymid X, heta^{(i)})}{P(X,Ymid heta^{(i)})}end{split}
所以两边取$log$并且运用Jensen不等式我们得到：
egin{split}log(frac{P(Xmid heta)}{P(Xmid heta^{(i)})})=&log(sum_{Y}frac{P(X,Ymid heta)P(Ymid X, heta^{(i)})}{P(X,Ymid heta^{(i)})}) ewline geq &sum_{Y}P(Ymid X, heta^{(i)})logfrac{P(X,Ymid heta)}{P(X,Ymid heta^{(i)})} ewline =&sum_{Y}[log P(X,Ymid heta)]P(Ymid X, heta^{(i)})-sum_{Y}[log P(X,Ymid heta^{(i)})]P(Ymid X, heta^{(i)}) ewline =&E_{Y}(log P(X,Ymid heta)ig | X, heta^{(i)})-E_{Y}(log P(X,Ymid heta^{(i)})ig | X, heta^{(i)})end{split}
我们令辅助函数:
egin{equation}Q( heta, heta^{(i)}) riangleq E_{Y}(log P(X,Ymid heta)ig | X, heta^{(i)})end{equation}
所以有：
egin{equation}log P(Xmid heta)-log P(Xmid heta^{(i)})geq Q( heta, heta^{(i)})-Q( heta^{(i)}, heta^{(i)})end{equation}
为了使得 $log p(Xmid heta)$显著增加，我们只需要选取一个:
egin{equation} heta^{(i+1)} riangleqmathop{ ext{argmax}}limits_{ heta} Q( heta, heta^{(i)})end{equation}
从上述推导我们看到，我们首先得显式的算出函数$Q$，它由一个求期望(Expectation)的计算给出，然后下一步是求Q的极大值使得Q极大化(Maximization)，所以我们的算法也会分为两步，一步是求期望(E步)，再求极大值（M步），这种方法于是被称为EM算法。
现在总结EM算法如下：
输入：观测变量$X$, 隐含变量$Y$, 联合分布$p(X,Ymid heta)$, $epsilon$
输出：模型参数$ heta$
Step1. 选择模型的初始参数$ heta^{(0)}$
Step2. E步：记$ heta^{(i)}$是第$i$步迭代得到的参数，则在第$i+1$步迭代中，计算函数：
egin{equation}Q( heta, heta^{(i)}) riangleq E_{Y}(log p(X,Ymid heta)ig | X, heta^{(i)})end{equation}
Step3 M步：更新 $ heta^{(i+1)} riangleqmathop{ ext{argmax}}limits_{ heta} Q( heta, heta^{(i)})$, 如果 $Q( heta^{(i+1)}, heta^{(i)})-Q( heta^{(i)}, heta^{(i)})geq epsilon$ 则回到Step2继续，否则停止，输出$ heta^{(i)}$.
高斯混合模型
所谓的高斯混合模型是对如下的隐变量模型应用EM算法统计推断：
问题2：
我们研究一类问题1的特例，即样本总体$mathcal{S}=mathbb{R}^{M} imeslbrace 1,2,...,K brace$, 分布密度函数$p(x,kmid heta)=pi_{k}P(xmidSigma_{k},mu_{k})$, 其中参数$ heta=(Sigma_{1},...,Sigma_{K},mu_{1},...,mu_{K},pi_{1},...,pi_{i})$ 满足对任意的k=1,...,K：
1)$Sigma_{k}$为对称正定$M imes M$矩阵；
2)$mu_{k}inmathbb{R}^{M}$；
3)$pi_{k}geq 0$,$sum_{k=1}^{K}pi_{i}=1$,
且我们用$P(xmidSigma,mu)$表示$M$维正态分布$mathcal{N}(mu_{k},Sigma_{k})$的概率密度函数。现在如果我们有数据集合$D=lbrace x_{1},...,x_{N} brace$, 其中$x_{i}inmathbb{R}^{M}$为第i个样本可观测向量，样本之间相互独立，试推断可能的参数$ heta$。
E步：
我们对上面的问题用EM算法求解，先开始E步，求函数Q。为了方便起见，对于任意$i=1,...,N$,我们定义函数： $l_{i}: SGL(M) imes mathbb{R}^{M}longrightarrow mathbb{R}:$
$$l_{i}(Sigma,mu)=log[P(x_{i}mid Sigma, mu)],$$
对任意的$Sigmain SGL(M), muinmathbb{R}^{M}$, 这里$SGL(M)$表示全体对称可逆的$M imes M$矩阵的集合。
这时：
egin{split}Q( heta, heta^{(s)})=& E_{Y}(log P(x_{1},...,x_{N},y_{1},...,y_{N}mid heta)ig | x_{1},...,x_{N}, heta^{(s)}) ewline =&sum_{i=1}^{N}E_{Y}(log P(x_{i},y_{i}mid heta)ig | X, heta^{(s)}) ewline =&sum_{i=1}^{N}E_{y_{i}}(log P(x_{i},y_{i}mid heta)ig | x_{i}, heta^{(s)}) ewline =&sum_{i=1}^{N}sum_{k=1}^{K}[log P(x_{i},kmid heta)]P(kmid x_{i}, heta^{(s)}) ewline =&sum_{i=1}^{N}sum_{k=1}^{K}P_{ki}^{(s)}log [P(x_{i}mid Sigma_{k},mu_{k})cdotpi_{k}] ewline =&sum_{i=1}^{N}sum_{k=1}^{K}P_{ki}^{(s)}l_{i}(Sigma_{k},mu_{k})+sum_{k=1}^{K}(sum_{i=1}^{N}P_{ki}^{(s)})log(pi_{k})end{split}
M步：
现在开始求Q的极大值，先得计算一下偏导数。下面我们先重点计算一下函数$l_{i}$的微分. 这时候 $P(x_{i}mid Sigma, mu)=(2pi)^{-frac{M}{2}}det(Sigma)^{-frac{1}{2}}explbrace -frac{1}{2}(x_{i}-mu)^{T}cdotSigma^{-1}cdot(x_{i}-mu) brace$, $$l_{i}(Sigma,mu)=-frac{1}{2}log(det(Sigma))-frac{1}{2}(x_{i}-mu)^{T}Sigma^{-1}(x_{i}-mu)$$.
我们注意到行列式求导公式：$d(log(detSigma))= ext{tr}(Sigma^{-1}dSigma)$所以容易得到：
egin{equation}dl_{i}=-frac{1}{2}trlbrace[Sigma^{-1}-Sigma^{-1}(x_{i}-mu)cdot(x_{i}-mu)^{T}Sigma^{-1}]dSigma brace+(x_{i}-mu)^{T}Sigma^{-1}dmuend{equation}
注意到我们要求极值问题：
egin{split} ext{max }& Q( heta, heta^{(s)}) ewline ext{s.t. }& pi_{k}in [0,1],& sum_{k=1}^{K}pi_{k}=1end{split}
只需要求某个模型参数$ heta$以及常数$lambdainmathbb{R}$ 满足：
egin{equation} dQ=lambda sum_{k=1}^{K}dpi_{k}end{equation}
结合(5)我们很容易得到：
$$dQ=-frac{1}{2} ext{tr}(A_{k}cdot dSigma_{k})+[sum_{i=1}^{N}P^{(s)}_{ki}x_{i}-(sum_{i=1}^{N}P_{ki})mu_{k}]d mu_{k}+sum_{k=1}^{K}(sum_{i=1}^{N}P_{ki}^{(s)})frac{d pi_{k}}{pi_{k}}$$其中矩阵$A_{k}=Sigma_{k}^{-1}(sum_{i=1}^{N}P_{ki}^{(s)})-Sigma_{k}^{-1}[sum_{i=1}^{N}P_{ki}^{(s)}(x_{i}-mu_{k})cdot (x_{i}-mu_{k})^{T}]Sigma_{k}^{-1}$
再结合（6）我们有：$A_{k}=0$, $sum_{i=1}^{N}P^{(s)}_{ki}x_{i}-(sum_{i=1}^{N}P_{ki})mu_{k}=0$, $sum_{i=1}^{N}P_{ki}^{(s)}=lambdapi_{k}$,
所以立即得到对任意$k=1,...,K$, 当$ heta$各个参数取如下值的时候$Q( heta, heta_{k})$将取极大值：
egin{equation}Sigma_{k}=frac{sum_{i=1}^{N}P_{ki}^{(s)}(x_{i}-mu_{k})cdot(x_{i}-mu_{k})^{T}}{sum_{i=1}^{N}P_{ki}^{(s)}},end{equation}
egin{equation}mu_{k}=frac{sum_{i=1}^{N}P_{ki}^{(s)}x_{i}}{sum_{i=1}^{N}P_{ki}^{(s)}},end{equation}
egin{equation}pi_{k}=frac{sum_{i=1}^{N}P^{(s)}_{ki}}{N},end{equation}
此时$M$步计算完毕。
此时我们可以总结一下Gauss混合模型的算法如下：
输入：数据集$D= lbrace x_{1},...,x_{N} brace $, 混合类别个数$K$，某停止条件
输出：高斯混合模型的参数$ heta=(Sigma_{1},...,Sigma_{K},mu_{1},...,mu_{K},pi_{1},...,pi_{K})$
Step1. 初始化某参数$ heta^{(0)}$
Step2. 对于$s=1,2,...$执行：
1. E-Step: 计算： egin{equation}P^{(s)}_{ki}=frac{P(x_{i}mid Sigma_{k}^{(s)},mu_{k}^{(s)})pi_{k}^{(s)}}{sum_{j=1}^{K} P(x_{i}mid Sigma_{j}^{(s)},mu_{j}^{(s)})pi_{j}^{(s)}}end{equation}
2. M-Step: 对任意的$k=1,...,K$计算：egin{equation}Sigma_{k}^{(s+1)}=frac{sum_{i=1}^{N}P_{ki}^{(s)}(x_{i}-mu_{k})cdot(x_{i}-mu_{k})^{T}}{sum_{i=1}^{N}P_{ki}^{(s)}},end{equation}
egin{equation}mu_{k}^{(s+1)}=frac{sum_{i=1}^{N}P_{ki}^{(s)}x_{i}}{sum_{i=1}^{N}P_{ki}^{(s)}},end{equation}
egin{equation}pi_{k}^{(s+1)}=frac{sum_{i=1}^{N}P^{(s)}_{ki}}{N},end{equation}
令$ heta^{(s+1)}=(Sigma_{1}^{(s+1)},...,Sigma_{K}^{(s+1)},mu_{1}^{(s+1)},...,mu_{K}^{(s+1)},pi_{1}^{(s+1)},...,pi_{K}^{(s+1)})$,
如果满足停止条件，则输出$ heta^{(s)}$
参考文献：
【1】李航：《统计学习方法》，北京，清华大学出版社，2012
【2】周志华：《机器学习》，北京，清华大学出版社，2016