python EM算法1

原文:https://blog.****.net/hx14301009/article/details/81228692
https://blog.****.net/zouxy09/article/details/8537620

假设我们需要调查我们学校的男生和女生的身高分布。在校园里随便地活捉了100个男生和100个女生，他们共200个人（也就是200个身高的样本数据）。

一 高斯模型

你开始喊：“男的左边，女的右边，其他的站中间！”。然后你就先统计抽样得到的100个男生的身高。假设他们的身高是服从高斯分布$N(\mu,\sigma)$的。但是这个分布的均值$\mu$和方差$\sigma^2$我们不知道，这两个参数就是我们要估计的。记作$\theta=[\mu,\sigma]^T$。

用数学的语言来说就是：在学校那么多男生（身高）中，我们独立地按照概率密度$p(x|\theta)$$（高斯分布$N(\mu,\sigma)$N(μ,σ)的形式）抽取100了个（身高），组成样本集X，我们想通过样本集X来估计出未知参数θ。抽到的样本集是$X=\{x_1,x_2,\cdots,x_N\}$X={x1​,x2​,⋯,xN​}$，其中`$，其中‘x_i$`表示抽到的第i个人的身高，这里抽到的样本个数N就是100。

因为这些男生（的身高）是服从同一个高斯分布$p(x|\theta)$的。那么我抽到男生A（的身高）的概率是$p(x_A|\theta)$$，抽到男生B的概率是$p(x_B|\theta)$p(xB​∣θ)$，那因为他们是独立的，所以很明显，我同时抽到男生A和男生B的概率是`$，那因为他们是独立的，所以很明显，我同时抽到男生A和男生B的概率是‘p(x_A|\theta)*p(x_B|\theta)$$，同理，我同时抽到这100个男生的概率就是他们各自概率的乘积了。即从分布是`$，同理，我同时抽到这100个男生的概率就是他们各自概率的乘积了。即从分布是‘p(x|\theta)$`的总体样本中抽取到这100个样本的概率，也就是样本集X中各个样本的联合概率，用下式表示：

L(\theta)=L(x_1,\cdots,x_n;\theta)=\prod\limits_{i=1}^n p(x_i;\theta),\theta\in \Theta

有时，可以看到$L(\theta)$是连乘的，所以为了便于分析，还可以定义对数似然函数，将其变成连加的：

\ln L(\theta)=\ln \prod\limits_{i=1}^n p(x_i;\theta)=\sum\limits_{i=1}^n \ln p(x_i;\theta),\theta\in \Theta

通过抽取得到的那100个男生的身高和已知其身高服从高斯分布，我们通过最大化其似然函数，就可以得到了对应高斯分布的参数$\theta=[\mu,\sigma]^T$了。对于女生的身高分布也可以用同样的方法得到。

二 高斯混合模型

再回到例子本身，如果没有“男的左边，女的右边，其他的站中间！”这个步骤，或者说我抽到这200个人都戴了面具。那现在这200个人已经混到一起了，这时候，你从这200个人（的身高）里面随便给我指一个人（的身高），我都无法确定这个人（的身高）是男生（的身高）还是女生（的身高）。也就是说你不知道抽取的那200个人里面的每一个人到底是从男生的那个身高分布里面抽取的，还是女生的那个身高分布抽取的。用数学的语言就是，抽取得到的每个样本都不知道是从哪个分布抽取的。

这个时候，对于每一个样本或者你抽取到的人，就有两个东西需要猜测或者估计的了，一是这个人是男的还是女的？二是男生和女生对应的身高的高斯分布的参数是多少？

为了解决这个你依赖我，我依赖你的循环依赖问题，总得有一方要先打破僵局，先随便整一个值出来，然后我再根据你的变化调整我的变化，然后如此迭代着不断互相推导，最终就会收敛到一个解。

2.1 EM解高斯混合模型

其实这个思想无处在不。

例如，小时候，老妈给一大袋糖果给你，叫你和你姐姐等分，然后你懒得去点糖果的个数，所以你也就不知道每个人到底该分多少个。咱们一般怎么做呢？先把一袋糖果目测的分为两袋，然后把两袋糖果拿在左右手，看哪个重，如果右手重，那很明显右手这代糖果多了，然后你再在右手这袋糖果中抓一把放到左手这袋，然后再感受下哪个重，然后再从重的那袋抓一小把放进轻的那一袋，继续下去，直到你感觉两袋糖果差不多相等了为止。呵呵，然后为了体现公平，你还让你姐姐先选了。

EM算法就是这样，假设我们想估计知道A和B两个参数，在开始状态下二者都是未知的，但如果知道了A的信息就可以得到B的信息，反过来知道了B也就得到了A。可以考虑首先赋予A某种初值，以此得到B的估计值，然后从B的当前值出发，重新估计A的取值，这个过程一直持续到收敛为止。

EM的意思是“Expectation Maximization”，在我们上面这个问题里面，我们是先随便猜一下男生（身高）的正态分布的参数：如均值和方差是多少。例如男生的均值是1米7，方差是0.1米（当然了，刚开始肯定没那么准），然后计算出每个人更可能属于第一个还是第二个正态分布中的（例如，这个人的身高是1米8，那很明显，他最大可能属于男生的那个分布），这个是属于Expectation一步。

有了每个人的归属，或者说我们已经大概地按上面的方法将这200个人分为男生和女生两部分，我们就可以根据之前说的最大似然那样，通过这些被大概分为男生的n个人来重新估计第一个分布的参数，女生的那个分布同样方法重新估计。这个是Maximization。

然后，当我们更新了这两个分布的时候，每一个属于这两个分布的概率又变了，那么我们就再需要调整E步……如此往复，直到参数基本不再发生变化为止。

这里把每个人（样本）的完整描述看做是三元组$y_i=\{x_i,z_{i1},z_{i2}\}$$，其中，$x_i$xi​是第i个样本的观测值，也就是对应的这个人的身高，是可以观测到的值。$z_{i1}$zi1​和$z_{i2}$zi2​表示男生和女生这两个高斯分布中哪个被用来产生值$x_i$xi​，就是说这两个值标记这个人到底是男生还是女生（的身高分布产生的）。这两个值我们是不知道的，是隐含变量。确切的说，$z_{ij}$zij​在$x_i$xi​由第j个高斯分布产生时值为1，否则为0。例如一个样本的观测值为1.8，然后他来自男生的那个高斯分布，那么我们可以将这个样本表示为$\{1.8, 1, 0\}${1.8,1,0}$。如果zi1和zi2的值已知，也就是说每个人我已经标记为男生或者女生了，那么我们就可以利用上面说的最大似然算法来估计他们各自高斯分布的参数。但是它们未知，因此我们只能用EM算法。

假设已经知道这个隐含变量了，那么直接按上面说的最大似然估计求解那个分布的参数就好了。我们可以先给这个分布弄一个初始值，然后求这个隐含变量的期望，当成是这个隐含变量的已知值，那么现在就可以用最大似然求解那个分布的参数了吧，那假设这个参数比之前的那个随机的参数要好，它更能表达真实的分布，那么我们再通过这个参数确定的分布去求这个隐含变量的期望，然后再最大化，得到另一个更优的参数，……迭代，就能得到一个皆大欢喜的结果了。

三、EM算法

3.1 Jensen不等式

凸函数：优化理论中，设f是定义域为实数的函数，如果对于所有的实数x，$f''(x)\geq 0$，那么f是凸函数。当x是向量时，如果其hessian矩阵H是半正定的（$H\geq 0$），那么f是凸函数。如果$f''(x)\leq 0$或者$H< 0$，那么称f是严格凸函数。 当f是（严格）凹函数当且仅当f是（严格）凸函数。比如$\log(x)$是凹函数。

Note: 发现好多地方凸函数定义可能是相反的，国内外定义好像也不一样。这里只要记得，正常的碗就是凸函数！

Jensen不等式：如果f是凸函数，X是随机变量，那么$E[f(X)]\geq f(EX)$$。特别地，如果f是严格凸函数，那么$E[f(X)]= f(EX)$E[f(X)]=f(EX)当且仅当$p(x=E[X])=1$p(x=E[X])=1，也就是说X是常量。Jensen不等式应用于凹函数时，不等号方向反向，也就是$E[f(X)]\leq f(EX)$E[f(X)]≤f(EX)$。这里我们将`$。这里我们将‘f(E[X])$`简写为`$‘简写为‘f(EX)$$。

如下所示$E[f(X)]\geq f(EX)$成立:

3.2 EM算法推导

假设我们有一个样本集`$\{x^{(1)}\cdots,x^{(m)}\}${x(1)⋯,x(m)}$，包含m个独立的样本。但每个样本i对应的类别z(i)是未知的（相当于聚类），也即隐含变量。故我们需要估计概率模型p(x,z)的参数θ，但是由于里面包含隐含变量z，所以很难用最大似然求解，但如果z知道了，那我们就很容易求解了.

对于参数估计，我们本质上还是想获得一个使似然函数最大化的那个参数θ，现在与最大似然不同的只是似然函数式中多了一个未知的变量z，见下式（1）。也就是说我们的目标是找到适合的θ和z让L(θ)最大。那我们也许会想，你就是多了一个未知的变量而已啊，我也可以分别对未知的θ和z分别求偏导，再令其等于0，求解出来不也一样吗？

\begin{array}{lllc}
\sum\limits_i \log p(x^{(i)},\theta)&=&\sum\limits_i \log\sum\limits_{z^{(i)}}p(x^{(i)},z^{(i)};\theta)&(1)\\
&=&\sum\limits_i\log\sum\limits_{z^{(i)}}Q_i(z^{(i)})\frac{p(x^{(i)},z^{(i)};\theta)}{Q_i(z^{(i)}}&(2)\\
&\geq&\sum\limits_i\sum\limits_{z^{(i)}}Q_i(z^{(i)})\log\frac{p(x^{(i)},z^{(i)};\theta)}{Q_i(z^{(i)})}&(3)
\end{array}

本质上我们是需要最大化（1）式（对（1）式，我们回忆下联合概率密度下某个变量的边缘概率密度函数的求解，注意这里z也是随机变量。对每一个样本i的所有可能类别z求等式右边的联合概率密度函数和，也就得到等式左边为随机变量x的边缘概率密度），也就是似然函数，但是可以看到里面有“和的对数”，求导后形式会非常复杂（自己可以想象下log(f1(x)+ f2(x)+ f3(x)+…)复合函数的求导），所以很难求解得到未知参数z和θ。那OK，我们可否对（1）式做一些改变呢？我们看（2）式，（2）式只是分子分母同乘以一个相等的函数，还是有“和的对数”啊，还是求解不了，那为什么要这么做呢？咱们先不管，看（3）式，发现（3）式变成了“对数的和”，那这样求导就容易了。我们注意点，还发现等号变成了不等号，为什么能这么变呢？这就是Jensen不等式的大显神威的地方。

2）到（3）利用了Jensen不等式，考虑到$\log(x)$是凹函数（二阶导数小于0），而且

$$\sum\limits_{z^{(i)}}Q_i(z^{(i)})\left[\frac{p(x^{(i)},z^{(i)};\theta)}{Q_i(z^{(i)})}\right]$z(i)∑​Qi​(z(i))[Qi​(z(i))p(x(i),z(i);θ)​]就是$\frac{p(x^{(i)},z^{(i)};\theta)}{Q_i(z^{(i)})}$Qi​(z(i))p(x(i),z(i);θ)​的期望,（考虑到E(X)=∑x*p(x)，f(X)是X的函数，则E(f(X))=∑f(x)*p(x)）,又$\sum\limits_{z^{(i)}}Q_i(z^{(i)})=1$z(i)∑​Qi​(z(i))=1`,所以就可以得到公式（3）的不等式了:

`$f\left(E_{z^{(i)}}\sim Q_i\left[ \frac{p(x^{(i)},z^{(i)};\theta)}{Q_i(z^{(i)})}\right]\right)\geq E_{z^{(i)}\sim Q_i}\left[f\left(\frac{p(x^{(i)},z^{(i)};\theta)}{Q_i(z^{(i)})}\right)\right]$f(Ez(i)​∼Qi​[Qi​(z(i))p(x(i),z(i);θ)​])≥Ez(i)∼Qi​​[f(Qi​(z(i))p(x(i),z(i);θ)​)]$.

OK，到这里，现在式（3）就容易地求导了，但是式（2）和式（3）是不等号啊，式（2）的最大值不是式（3）的最大值啊，而我们想得到式（2）的最大值，那怎么办呢？

现在我们就需要一点想象力了，上面的式（2）和式（3）不等式可以写成：似然函数L(θ)>=J(z,Q)，那么我们可以通过不断的最大化这个下界J，来使得L(θ)不断提高，最终达到它的最大值。

见上图，我们固定θ，调整Q(z)使下界J(z,Q)上升至与L(θ)在此点θ处相等（绿色曲线到蓝色曲线），然后固定Q(z)，调整θ使下界J(z,Q)达到最大值（KaTeX parse error: Can't use function '$' in math mode at position 9: \theta^t$̲`到`$\theta^{t+1…），然后再固定θ，调整Q(z)……直到收敛到似然函数L(θ)的最大值处的θ*。这里有两个问题：什么时候下界J(z,Q)与L(θ)在此点θ处相等？为什么一定会收敛？

首先第一个问题，在Jensen不等式中说到，当自变量X是常数的时候，等式成立。而在这里，即:

\frac{p(x^{(i)},z^{(i)};\theta)}{Q_i(z^{(i)})}=c

再推导下，由于$\sum\limits_{z^{(i)}}Q_i(z^{(i)})=1$（因为Q是随机变量$z^{(i)}$的概率密度函数），则可以得到：分子的和等于c（分子分母都对所有$z^{(i)}$求和：多个等式分子分母相加不变，这里认为每个样例的两个概率比值都是c），则：

\begin{array}{lll}
Q_i(z^{(i)})&=&\frac{p(x^{(i)},z^{(i)};\theta)}{\sum\limits_{z^{(i)}}p(x^{(i)},z;\theta)}\\
&=&\frac{p(x^{(i)},z^{(i)};\theta)}{p(x^{(i)};\theta)}\\
&=&p(z^{(i)}\vert x^{(i)};\theta)
\end{array}

至此，我们推出了在固定参数θ后，使下界拉升的Q(z)的计算公式就是后验概率，解决了Q(z)如何选择的问题。这一步就是E步，建立L(θ)的下界。接下来的M步，就是在给定Q(z)后，调整θ，去极大化L(θ)的下界J（在固定Q(z)后，下界还可以调整的更大）。那么一般的EM算法的步骤如下：

EM算法的步骤

Note: 通过男女混合高斯模型来记！

初始化分布参数θ；
循环重复直到收敛 {
KaTeX parse error: Can't use function '$' in math mode at position 7: \qquad$̲`（E步:建立`$\ell(\…\qquad\qquad\quad$#E步骤：根据参数初始值或上一次迭代的模型参数来计算出隐性变量的后验概率，其实就是隐性变量的期望。作为隐藏变量的现估计值： $$\qquad\qquad\qquad Q_i(z^{(i)}):=p(z^{(i)}\vert x^{(i)};\theta)$$ $$\qquad$（M步:优化下界）计算
KaTeX parse error: Can't use function '$' in math mode at position 18: …quad\qquad\quad$̲`#将似然函数最大化以获得新的…\qquad\qquad\qquad \theta:=\underset{\theta}{\arg\max}\sum\limits_i\sum\limits_{z^{(i)}}Q_i(z{(i)})\log\frac{p(x^{(i)},z{(i)};\theta)}{Q_i(z^{(i)})}$$

这个不断的迭代，就可以得到使似然函数L(θ)最大化的参数θ了。那就得回答刚才的第二个问题了，它会收敛吗？

3.3 那么究竟怎么确保EM收敛？

假定$\theta^{(t)}$和$\theta^{(t+1)}$是EM第t次和t+1次迭代后的结果。

如果我们证明`$\theta^{(t)}\leq\theta^{(t+1)}$θ(t)≤θ(t+1)$，也就是说极大似然估计单调增加，那么最终我们会到达最大似然估计的最大值。

下面来证明，选定$\theta^{(t)}$后，我们得到E步

`$Q_{i}^{(t)}(z^{(i)}):=p(z^{(i)}\vert x^{(i)};\theta ^{(t)})$Qi(t)​(z(i)):=p(z(i)∣x(i);θ(t))$

这一步保证了在给定$\theta^{(t)}$时，Jensen不等式中的等式成立，也就是

$\ell (\theta^{(t)})=\sum\limits_i\sum\limits_{z^{(i)}}Q_i(z^{(i)})\log\frac{p(x^{(i)},z^{(i)};\theta^{(t)})}{Q_i(z^{(i)})}$ℓ(θ(t))=i∑​z(i)∑​Qi​(z(i))logQi​(z(i))p(x(i),z(i);θ(t))​

然后进行M步，固定$Q_{i}^{(t)}(z^{(i)})$$，并将$\theta^{(t)}$θ(t)视作变量，对上面的$\ell (\theta^{(t)})$ℓ(θ(t))求导后，得到$\theta^{(t+1)}$θ(t+1)$，这样经过一些推导会有以下式子成立：

\begin{array}{lllc}
\ell(\theta^{(t+1)})&\geq&\sum\limits_i\sum\limits_{z^{(i)}}Q_i(z^{(i)})\log\frac{p(x^{(i)},z^{(i)};\theta^{(t+1)})}{Q_i(z^{(i)})}&(4)\\
&\geq&\sum\limits_i\sum\limits_{z^{(i)}}Q_i(z^{(i)})\log\frac{p(x^{(i)},z^{(i)};\theta^{(t)})}{Q_i(z^{(i)})}&(5)\\
&=&\ell(\theta^{(t)})&(6)
\end{array}

解释第（4）步，得到$\theta^{(t+1)}$时，只是最大化$\theta^{(t+1)}$$，而没有使等式成立，等式成立只有是在固定$\theta$θ，并按E步得到$Q_i$Qi​`时才能成立。

第（5）步利用了M步的定义，M步就是将$\theta^{(t)}$调整到`$\theta^{(t+1)}$θ(t+1)$，使得下界最大化。

（6）是之前的等式结果.

如果我们定义:

`$J(Q,\theta)=\sum\limits_i\sum\limits_{z^{(i)}}Q_i(z^{(i)})\log\frac{p(x^{(i)},z^{(i)};\theta)}{Q_i(z^{(i)})}$J(Q,θ)=i∑​z(i)∑​Qi​(z(i))logQi​(z(i))p(x(i),z(i);θ)​$,

从前面的推导中我们知道$\ell(\theta)\geq J(Q,\theta)$$，EM可以看作是J的坐标上升法，E步固定$\theta$θ，优化Q，M步固定Q优化$\theta$θ`。

3.4 EM算法的具体实现

E步很简单，按照一般EM公式得到：

$\omega_j^{(i)}=Q_i(z^{(i)}=j)=P(z^{(i)}=j\vert x^{(i)};\phi,\mu,\Sigma)$ωj(i)​=Qi​(z(i)=j)=P(z(i)=j∣x(i);ϕ,μ,Σ).

简单解释就是每个样例i的隐含类别$z^{(i)}$为j的概率可以通过后验概率计算得到。

在M步中，我们需要$\ell(\theta),\theta=[\phi,\mu,\Sigma]$对$\theta,\mu,\Sigma$分别求偏导.

首先,

$\ell(\theta)=\sum\limits_{i=1}^{m}\sum\limits_{z^{(i)}}Q_i(z^{(i)})\log\frac{P(x^{(i)},z^{(i)};\phi,\mu,\Sigma)}{Q_i(z^{(i)})}$ℓ(θ)=i=1∑m​z(i)∑​Qi​(z(i))logQi​(z(i))P(x(i),z(i);ϕ,μ,Σ)​

$\qquad=\sum\limits_{i=1}{m}\sum\limits_{j=1}^k Q_i(z^{(i)}=j)\log\frac{P(x^{(i)}\vert z^{(i)}=j;\mu,\Sigma)p(z^{(i)}=j;\phi)}{Q_i(z^{(i)}=j)}$=i=1∑​mj=1∑k​Qi​(z(i)=j)logQi​(z(i)=j)P(x(i)∣z(i)=j;μ,Σ)p(z(i)=j;ϕ)​

$\qquad=\sum\limits_{i=1}^m\sum\limits_{j=1}^k \omega_j^{(i)}\log\frac{\frac{1}{(2\pi)^{n/2}|\Omega_j|^{1/2}}\exp(-\frac{1}{2}(x^{(i)}-\mu_j)^T\Omega_j^{-1}(x^{(i)}-\mu_j))\dot \phi_j}{\omega_j^{(i)}}$=i=1∑m​j=1∑k​ωj(i)​logωj(i)​(2π)n/2∣Ωj​∣1/21​exp(−21​(x(i)−μj​)TΩj−1​(x(i)−μj​))ϕ˙​j​​,

这里将$z^{(i)}$的k种情况展开,其中$\phi_j,\mu_j,\Omega_j$是未知参数.

其次,分别对$\phi_j,\mu_j,\Omega_j$求偏导.

(1) 固定$\phi_j$和$\Sigma_j$,对$\mu_j$求导得

$\nabla_{\mu_\ell}\sum\limits_{i=1}^m\sum\limits_{j=1}^k \omega_j^{(i)}\log\frac{\frac{1}{(2\pi)^{n/2}|\Omega_j|^{1/2}}\exp(-\frac{1}{2}(x^{(i)}-\mu_j)^T\Omega_j^{-1}(x^{(i)}-\mu_j))\dot \phi_j}{\omega_j^{(i)}}$∇μℓ​​i=1∑m​j=1∑k​ωj(i)​logωj(i)​(2π)n/2∣Ωj​∣1/21​exp(−21​(x(i)−μj​)TΩj−1​(x(i)−μj​))ϕ˙​j​​

$\qquad=-\nabla_{\mu_\ell}\sum\limits_{i=1}^m\sum\limits_{j=1}^k \omega_j^{(i)}\frac{1}{2}(x^{(i)}-\mu_j)^T \Omega_{j}^{-1}(x^{(i)}-\mu_j)$=−∇μℓ​​i=1∑m​j=1∑k​ωj(i)​21​(x(i)−μj​)TΩj−1​(x(i)−μj​)

$$\qquad=\frac{1}{2}\sum\limits_{i=1}^m \omega_{\ell}^{(i)}\nabla_{\mu_{\ell}}2\mu_{\ell}^T\Sigma_{\ell}^{-1}x^{(i)}-\mu_{\ell}^T\Sigma_{\ell}^{-1}\mu_{\ell}=0$=21​i=1∑m​ωℓ(i)​∇μℓ​​2μℓT​Σℓ−1​x(i)−μℓT​Σℓ−1​μℓ​=0`,

得到
$\qquad\mu_{\ell}:=\frac{\sum_{i=1}^m\omega_{\ell}^{(i)}x^{(i)}}{\sum_{i=1}^m\omega_{\ell}^{(i)}}$μℓ​:=∑i=1m​ωℓ(i)​∑i=1m​ωℓ(i)​x(i)​.

(2) 固定$\mu_j$和$\Sigma_j$,对$\phi_j$求导得

公式$\ell(\theta)=\sum\limits_{i=1}^{m}\sum\limits_{z^{(i)}}Q_i(z^{(i)})\log\frac{P(x^{(i)},z^{(i)};\phi,\mu,\Sigma)}{Q_i(z^{(i)})}$中,在$\mu_j$和$\Omega_j$确定后,分子上面的一串都是常数了，实际上需要优化的公式是：

$$\qquad\sum\limits_{i=1}^m\sum\limits_{j=1}^k\omega_j^{(i)}\log\phi_j$i=1∑m​j=1∑k​ωj(i)​logϕj​`.

需要知道的是，$$\phi_j$ϕj​还需要满足一定的约束条件就是$\sum\limits_{j=1}^k\phi_j=1$j=1∑k​ϕj​=1`。

这个优化问题我们很熟悉了，直接构造拉格朗日乘子如下:

$L(\phi)=\sum\limits_{i=1}^m\sum\limits_{j=1}^k\omega_j^{(i)}\log\phi_j+\beta(\sum\limits_{j=1}^k\phi_j-1)$$,这里$\phi_j\geq 0$ϕj​≥0`.

求导得:

KaTeX parse error: Can't use function '$' in math mode at position 105: …\phi_j}+\beta=0$̲`,得到`$\phi_j=\f…,再次使用$\sum\limits_{j=1}^k\phi_j=1$,得到$-\beta=\sum\limits_{i=1}^m\sum\limits_{j=1}^k\omega_{j=1}^{(i)}=\sum\limits_{i=1}^m 1=m$.

从而,得到M步中$\phi_j$的更新公式:

$\qquad\phi_j:=\frac{1}{m}\sum\limits_{i=1}^m \omega_j^{(i)}$ϕj​:=m1​i=1∑m​ωj(i)​.

(3) 固定$\mu_j$和$\phi_j$,对$\Sigma_j$求导得

$\nabla_{\Sigma_j}\sum\limits_{i=1}^m\sum\limits_{j=1}^k \omega_j^{(i)}\log\frac{\frac{1}{(2\pi)^{n/2}|\Sigma_j|^{1/2}}\exp(-\frac{1}{2}(x^{(i)}-\mu_j)^T\Sigma_j^{-1}(x^{(i)}-\mu_j))\dot \phi_j}{\omega_j^{(i)}}$∇Σj​​i=1∑m​j=1∑k​ωj(i)​logωj(i)​(2π)n/2∣Σj​∣1/21​exp(−21​(x(i)−μj​)TΣj−1​(x(i)−μj​))ϕ˙​j​​

$\qquad=\nabla_{\Sigma_j}\sum\limits_{i=1}^m\sum\limits_{j=1}^k \omega_j^{(i)}\left[\log\frac{1}{(2\pi)^{n/2}|\Sigma_j|^{1/2}}\exp(-\frac{1}{2}(x^{(i)}-\mu_j)^T\Sigma_j^{-1}(x^{(i)}-\mu_j))\right]$=∇Σj​​i=1∑m​j=1∑k​ωj(i)​[log(2π)n/2∣Σj​∣1/21​exp(−21​(x(i)−μj​)TΣj−1​(x(i)−μj​))]

$\qquad=\nabla_{\Sigma_j}\sum\limits_{i=1}^m\sum\limits_{j=1}^k \omega_j^{(i)}\left[-\log((2\pi)^{n/2}|\Sigma_j|^{1/2})+\log(\exp(-\frac{1}{2}(x^{(i)}-\mu_j)^T\Sigma_j^{-1}(x^{(i)}-\mu_j)))\right]$=∇Σj​​i=1∑m​j=1∑k​ωj(i)​[−log((2π)n/2∣Σj​∣1/2)+log(exp(−21​(x(i)−μj​)TΣj−1​(x(i)−μj​)))]

$\qquad=\nabla_{\Sigma_j}\sum\limits_{i=1}^m\sum\limits_{j=1}^k \omega_j^{(i)}\left[-\frac{1}{2}\log|\Sigma_j|-\frac{1}{2}(x^{(i)}-\mu_j)^T\Sigma_j^{-1}(x^{(i)}-\mu_j)\right]$=∇Σj​​i=1∑m​j=1∑k​ωj(i)​[−21​log∣Σj​∣−21​(x(i)−μj​)TΣj−1​(x(i)−μj​)]

KaTeX parse error: Can't use function '$' in math mode at position 7: \qquad$̲`下面用到了矩阵公式:\nabla_A|A|=|A|(A^{-1})T,\ \frac{\partial a^TX{-1}b}{\partial X}=-X^{-T}abTX^{-T}$$.

$\qquad=\sum\limits_{i=1}^m\sum\limits_{j=1}^k \omega_j^{(i)}\left[-\frac{1}{2}\frac{1}{|\Sigma_j|}|\Sigma_j|(\Sigma_j^{-1})^T+\frac{1}{2}\Sigma_j^{-T}(x^{(i)}-\mu_j)(x^{(i)}-\mu_j)^T\Sigma_j^{-T}\right]$=i=1∑m​j=1∑k​ωj(i)​[−21​∣Σj​∣1​∣Σj​∣(Σj−1​)T+21​Σj−T​(x(i)−μj​)(x(i)−μj​)TΣj−T​]

$$\qquad=\sum\limits_{i=1}^m\sum\limits_{j=1}^k \omega_j^{(i)}\left[-\frac{1}{2}+\frac{1}{2}(x^{(i)}-\mu_j)(x^{(i)}-\mu_j)^T\Sigma_j^{-T}\right]\Sigma_j^{-T}=0$=i=1∑m​j=1∑k​ωj(i)​[−21​+21​(x(i)−μj​)(x(i)−μj​)TΣj−T​]Σj−T​=0`

可得:

$\qquad\Sigma_j=\frac{\sum\limits_{i=1}^m\omega_j^{(i)}(x^{(i)}-\mu_j)(x^{(i)}-\mu_j)^T}{\sum\limits_{i=1}^m\omega_j^{(i)}}$Σj​=i=1∑m​ωj(i)​i=1∑m​ωj(i)​(x(i)−μj​)(x(i)−μj​)T​.

综上所述:

E步:

$\qquad\qquad\quad\, \omega_j^{(i)}=Q_i(z^{(i)}=j)=P(z^{(i)}=j\vert x^{(i)};\phi,\mu,\Sigma)$ωj(i)​=Qi​(z(i)=j)=P(z(i)=j∣x(i);ϕ,μ,Σ),

$\qquad\qquad\qquad\quad=\frac{\phi_j N(x^{(i)}|\mu_j,\Sigma_j)}{\sum\limits_{j=1}^kN(x^{(i)}|\mu_j,\Sigma_j)}$=j=1∑k​N(x(i)∣μj​,Σj​)ϕj​N(x(i)∣μj​,Σj​)​.

M步:

$\qquad\qquad\quad\,\mu_{j}:=\frac{\sum_{i=1}^m\omega_{\ell}^{(i)}x^{(i)}}{\sum_{i=1}^m\omega_{\ell}^{(i)}}$μj​:=∑i=1m​ωℓ(i)​∑i=1m​ωℓ(i)​x(i)​

$\qquad\qquad\quad\,\Sigma_j=\frac{\sum\limits_{i=1}^m\omega_j^{(i)}(x^{(i)}-\mu_j)(x^{(i)}-\mu_j)^T}{\sum\limits_{i=1}^m\omega_j^{(i)}}$Σj​=i=1∑m​ωj(i)​i=1∑m​ωj(i)​(x(i)−μj​)(x(i)−μj​)T​

$\qquad\qquad\quad\,\phi_j:=\frac{1}{m}\sum\limits_{i=1}^m \omega_j^{(i)}$ϕj​:=m1​i=1∑m​ωj(i)​

高斯混合聚类算法步骤的完整描述：

---

输入:样本集`$D=\{x^{(1)},x^{(2)},\cdots,x^{(m)}\};$D={x(1),x(2),⋯,x(m)};$

$$\quad\ $`高斯混合成分个数k.

过程:
1: 初始化高斯混合分布的模型参数$\{(\phi_j,\mu_j,\Sigma_j)\vert 1\leq j \leq k\}$$ 2: repeat 3: $$\qquad$for $j=1,2,\cdots,m$ do
4:KaTeX parse error: Can't use function '$' in math mode at position 13: \qquad\qquad$̲`计算由`$x^{(i)}$`…\qquad\qquad\quad, \omega_j^{(i)}=Q_i(z{(i)}=j)=P(z^{(i)}=j\vert x^{(i)};\phi,\mu,\Sigma)$$,

$\qquad\qquad\qquad\quad=\frac{\phi_j N(x^{(i)}|\mu_j,\Sigma_j)}{\sum\limits_{j=1}^kN(x^{(i)}|\mu_j,\Sigma_j)}$=j=1∑k​N(x(i)∣μj​,Σj​)ϕj​N(x(i)∣μj​,Σj​)​.
5:KaTeX parse error: Can't use function '$' in math mode at position 9: \qquad\ $̲`end for 6: \qquad$`for `$‘for‘j=1,2,\cdots,k$` do 7:$‘do7:$\qquad\qquad$计算新均值向量:$$\mu_{j}:=\frac{\sum_{i=1}^m\omega_{\ell}^{(i)}x^{(i)}}{\sum_{i=1}^m\omega_{\ell}^{(i)}}$$ 8:$$\qquad\qquad$计算新协方差向量:$\Sigma_j=\frac{\sum\limits_{i=1}^m\omega_j^{(i)}(x^{(i)}-\mu_j)(x^{(i)}-\mu_j)^T}{\sum\limits_{i=1}^m\omega_j^{(i)}}$Σj​=i=1∑m​ωj(i)​i=1∑m​ωj(i)​(x(i)−μj​)(x(i)−μj​)T​
9:KaTeX parse error: Can't use function '$' in math mode at position 13: \qquad\qquad$̲`计算新混合向量:\phi_j:=\frac{1}{m}\sum\limits_{i=1}^m \omega_j^{(i)}$10:$10:\qquad$`end for 11:$‘endfor11:$\qquad$将模型参数更新为$\{(\phi_j,\mu_j,\Sigma_j)\vert 1\leq j\leq k\}${(ϕj​,μj​,Σj​)∣1≤j≤k}$
12:until 满足停止条件
13:$C_j=\emptyset(1\leq j\leq k)$$ 14:for$i=1,2,\cdots,m$i=1,2,⋯,mdo 15: $$\qquad$将$x_i$的簇标记为$\lambda_i$;
16: KaTeX parse error: Can't use function '$' in math mode at position 7: \qquad$̲`将`$x_i$`划入相应的簇…
17:end for
输出:簇划分`$C=\{C_1,C_2,\cdots,C_k\}$C={C1​,C2​,⋯,Ck​}$

---