Kmeans和GMM参数学习的EM算法原理和Matlab实现

本文整理自JerryLead的博文“《K-means聚类算法》 ”，“《（EM算法）The EM Algorithm 》”，“《混合高斯模型（Mixtures of Gaussians）和EM算法 》”，以及自己编写的关于GMM的Matlab实现。
好的博文还有：
pluskid的《漫谈 Clustering (3): Gaussian Mixture Model》和《漫谈 Clustering (番外篇): Expectation Maximization》。

K-means和EM

K-means是聚类算法中最简单的一种，但是里面包含的思想却是不一般。聚类属于无监督学习，朴素贝叶斯、SVM等都是有类别标签$y$的，即已经给出了样本的分类。而聚类的样本中却没有给定$y$，只有特征$x$。聚类的目的是找到每个样本$x$潜在的类别$y$，并将同类别$y$的样本$x$放在一起。
在聚类问题中，给定训练样本${\{x_i\}}_{i=1}^{N}$，每个$x_i \in R^n$，没有类别标签$y$。
K-means算法是将样本聚类成$K$个簇（cluster），具体算法描述如下：

1、随机选取$K$个聚类质心点为${\{\mu_k\}}_{k=1}^{K},\mu_k\in R^n$ 。
2、重复下面过程直到收敛{
对于每一个样例$i$，计算其应该属于的类

$$z_i=arg\mathop{\min}\limits_j||x_i-\mu_j||^2$$
对于每一个类$j$，重新计算该类的质心
$$\mu_j=\frac{\sum_{i=1}^{N}{1{\{z_i=j\}x_i}}}{\sum_{i=1}^{N}1{\{z_i=j}\}}$$
}

其中，$K$是事先给定的聚类数，$z_i \in {1,2,..,K}$代表样例$i$与$K$个类中距离最近的那个类，$\mu_j$代表对属于同一个类的样本的质心。要将所有的样本聚成$K$类，首先随机选取$K$个点作为$K$个类的质心，第一步对于每一个样本计算其到$K$个质心的距离，选取距离最近的那个类作为$z_i$（经过第一步每个样本都有了所属的类）；第二步对于每一个类，重新计算它的质心$\mu_j$（对里面所有的样本坐标求平均）。重复迭代第一步和第二步直到质心不变或者变化很小。
K-means如何保证收敛?
K-means算法中强调结束条件就是收敛，可以证明的是K-means完全可以保证收敛性。下面定性的描述一下收敛性，定义畸变函数：

$$J(z,\mu)=\mathop {\sum}_{i=1}^{N}||x_i-\mu_{z_{i}}||^2$$
函数 $J$表示每个样本点到其质心的距离平方和。K-means就是要将$J$调整到最小。假设当前 $J$没有达到最小值，那么可以首先固定每个类的质心$\mu_j$，调整每个样例的所属类别 $z_i$来让 $J$函数减少，同样固定$z_i$，调整 $\mu_j$也可以使 $J$减小。这两个过程就是内循环中使$J$单调递减的过程。当 $J$递减到最小时，$\mu$和 $z$也同时收敛。由于畸变函数$J$是非凸函数，意味着不能保证取得的最小值是全局最小值，也就是说K-means对质心初始位置的选取比较敏感，但一般情况下K-means达到的局部最优已经满足需求。
K-means与EM的关系
K-means算法目的是将样本分成 $K$个类，即求每个样例$x$的隐含类别 $y$，然后利用隐含类别将$x$归类。由于事先不知道类别 $y$，那么可以先对每个样例假定一个$y$，但是怎么知道假定的对不对呢？怎么评价假定的好不好呢？在此，可以使用样本的极大似然估计来度量，即 $x$和$y$的联合分布 $p(x,y)$。如果找到的 $y$能够使 $p(x,y)$最大，那么找到的 $y$就是样例$x$的最佳类别， $x$同时就聚类了。但是第一次指定的$y$不一定会让 $p(x,y)$最大，而且 $p(x,y)$还依赖于其他未知参数。当然在给定 $y$的情况下，可以调整其他参数让$p(x,y)$最大；调整完参数后，发现有更好的 $y$可以指定，那么重新指定$y$。反复迭代直至没有更好的 $y$可以指定。
这个过程有几个难点：第一，怎么假定$y$？是每个样例硬指派一个 $y$还是不同的$y$有不同的概率，概率如何度量？第二，如何估计 $p(x,y)$， $p(x,y)$还可能依赖很多其他参数，如何调整里面的参数让 $p(x,y)$最大。
这可以采用EM的思想解决，E步就是估计隐含类别 $y$的期望值，M步调整其他参数使得在给定类别$y$的情况下，极大似然估计 $p(x,y)$能够达到极大值。然后在其他参数确定的情况下，重新估计 $y$，周而复始，直至收敛。
对应于K-means，最开始对于每个样例$x_i$随便指定一个 $z_i$，然后为了让 $p(x,y)$最大（让 $J$最小），求给定$z_i$情况下最小时的 $\mu_j$（前面提到的其他未知参数），然而此时发现，可以有更好的 $z_i$（质心与样例 $x_i$距离最小的类别）指定给样例 $x_i$，那么 $z_i$得到重新调整，上述过程就开始重复了，直到没有更好的 $z_i$指定。可见K-means里体现了EM思想，E步是确定隐含类别变量 $z$，M步更新其他参数$\mu$来使 $J$最小化。K-means隐含类别变量指定方法比较特殊，属于硬指定，即从$K$个类别中硬选出一个给样例，而不是对每个类别赋予不同的概率。总体思想还是一个迭代优化过程，有目标函数，也有参数变量，只是多了个隐含变量，确定其他参数估计隐含变量，再确定隐含变量估计其他参数，直至目标函数最优。

EM算法

下面介绍EM的整个推导过程：

Jensen不等式

设$f$是定义域为实数的函数，如果对于所有的实数$X$，$f''(x) \ge 0$，那么$f$是凸函数。当$x$是向量时，如果其Hessian矩阵$H$是半正定的，那么$f$是凸函数。如果$f''(x) > 0$，那么称$f$是严格凸函数。
Jensen不等式表述如下：如果$f$是凸函数，$X$是随机变量，那么

$$E[f(x)] \ge f(E[X])$$
特别地，如果$f$是严格凸函数，那么$E[f(x)] = f(E[X])$当且仅当$p(x=E(X))=1$，也就是说$X$是常量。
如果用图表示会很清晰：

图中实线$f$是凸函数，$X$是随机变量，有0.5的概率是a，有0.5的概率是b，$X$ 的期望值就是a和b的中值了，图中可以看到$E[f(x)] \ge f(E[X])$成立。
当$f$是（严格）凹函数当且仅当$-f$是（严格）凸函数。
Jensen不等式应用于凹函数时，不等号方向反向，即$E[f(x)] \le f(E[X])$ 。

EM算法

给定独立同分布训练样本${\{x_i\}}_{i=1}^{N}$，寻找到每个样例隐含的类别$z$，能使得$p(x,z)$最大。$p(x,z)$的最大似然估计如下：

$$L(\theta)=\mathop {\sum}_{i=1}^{N}logp(x_i;\theta) \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; =\mathop {\sum}_{i=1}^{N}log\sum_{z_i}p(x_i,z_i;\theta)$$
该式直接求 $\theta$一般比较困难，因为有隐藏变量 $z$存在，但是一般确定了$z$后，求解就容易了。
EM是一种解决存在隐含变量优化问题的有效方法。既然不能直接最大化 $L(\theta)$，那么就不断地建立 $L(\theta)$的下界（E步），然后优化下界（M步）。即对于每一个样例 $i$，让$Q(z_i)$表示该样例隐含变量 $z_i$的某种分布，且 $\sum_{z_i}Q(z_i)=1$， $Q(z_i) \ge 0$。由Jensen不等式得：
$$\begin{array}{l} L(\theta ) = \sum\limits_{i = 1}^N {\log \sum\limits_{{z_i}} {p({x_i},{z_i};\theta )} } \\ \;\;\;\;\;\;\;\; = \sum\limits_{i = 1}^N {\log \sum\limits_{{z_i}} {Q({z_i})\frac{{p({x_i},{z_i};\theta )}}{{Q({z_i})}}} } \\ \;\;\;\;\;\;\;\; \ge \sum\limits_{i = 1}^N {\sum\limits_{{z_i}} {Q({z_i})\log \frac{{p({x_i},{z_i};\theta )}}{{Q({z_i})}}} } \end{array}$$
其中， $log(x)$是凹函数且 $\sum\limits_{{z_i}} {Q({z_i})\frac{{p({x_i},{z_i};\theta )}}{{Q({z_i})}}}$
是 $[p({x_i},{z_i};\theta )/Q({z_i})]$的期望。
这个过程可以看作是对 $L(\theta)$求了下界。对于 $Q$的选择，有多种可能，那种更好的？假设$\theta$已经给定，那么 $l(\theta)$的值就决定于 $Q(z_i)$和 $p(x_i,z_i)$了。通过调整这两个概率使下界不断上升，以逼近 $L(\theta)$的真实值，那么什么时候算是调整好了呢？当不等式变成等式时，说明调整后的概率能够等价于 $L(\theta)$。按照这个思路，要找到等式成立的条件。根据Jensen不等式，要想让等式成立，需要让随机变量变成常数值，即：
$$\frac{{p({x_i},{z_i};\theta )}}{{Q({z_i})}} = c$$
其中， $c$为常数，不依赖于$z_i$。由于 $\sum\nolimits_z {Q(z)} = 1$，所以 $\sum\nolimits_{{z_i}} {p({x_i},{z_i};\theta )} = 1$，（多个等式分子分母相加不变，这个认为每个样例的两个概率比值都是 $c$），那么有下式：
$$\begin{array}{l} Q({z_i}) = \frac{{p({x_i},{z_i};\theta )}}{{\sum\nolimits_z {p({x_i},z;\theta )} }}\\ \;\;\;\;\;\;\;\;\; = \frac{{p({x_i},{z_i};\theta )}}{{p({x_i};\theta )}}\\ \;\;\;\;\;\;\;\;\; = p({z_i}|{x_i};\theta ) \end{array}$$
至此，导出了在固定其他参数 $\theta$后， $Q(z_i)$就是后验概率，解决了 $Q(z_i)$如何选择的问题。这一步就是E步，建立 $L(\theta)$的下界。
接下来的M步，就是在给定 $Q(z_i)$后，调整 $\theta$，去极大化 $L(\theta)$的下界。那么一般的EM算法的步骤如下：

循环重复直到收敛 {
（E步）对于每一个$i$，计算

$$Q({z_i}) = p({z_i}|{x_i};\theta )$$
（M步）计算

$$\theta = \arg \mathop {\max }\limits_\theta \sum\limits_{i = 1}^N {\sum\limits_{{z_i}} {Q({z_i})\log \frac{{p({x_i},{z_i};\theta )}}{{Q({z_i})}}} }$$

那么究竟怎么确保EM收敛？假定$\theta^{t}$和$\theta^{t+1}$是EM第$t$次和第$t+1$次迭代后的结果。如果$L({\theta ^t}) \le L({\theta ^{t + 1}})$，也就是说极大似然估计单调增加，且$L(\theta)$有界，那么迭代会到达最大似然估计的最大值。下面证明$L({\theta ^t}) \le L({\theta ^{t + 1}})$：
固定$\theta^{t}$后，E步

$${Q^t}({z_i}) = p({z_i}|{x_i};{\theta ^t})$$
该步保证了在给定 $\theta^{t}$时，Jensen不等式中的等式成立，即
$$L({\theta ^t}) = \sum\limits_{i = 1}^N {\sum\limits_{{z_i}} {{Q^t}({z_i})\log \frac{{p({x_i},{z_i};{\theta ^t})}}{{{Q^t}({z_i})}}} }$$
固定 ${Q^t}({z_i})$，M步：
将 $\theta^{t}$视作变量，对 $L(\theta^{t})$求导等于零可得到 $\theta^{t+1}$，则下式成立：
$$\begin{array}{l} L({\theta ^t}) = \sum\limits_{i = 1}^N {\sum\limits_{{z_i}} {{Q^t}({z_i})\log \frac{{p({x_i},{z_i};{\theta ^t})}}{{{Q^t}({z_i})}}} } \\ \;\;\;\;\;\;\;\; \le \sum\limits_{i = 1}^N {\sum\limits_{{z_i}} {{Q^t}({z_i})\log \frac{{p({x_i},{z_i};{\theta ^{t + 1}})}}{{{Q^t}({z_i})}}} } \\ \;\;\;\;\;\;\;\; \le \sum\limits_{i = 1}^N {\sum\limits_{{z_i}} {{Q^{t + 1}}({z_i})\log \frac{{p({x_i},{z_i};{\theta ^{t + 1}})}}{{{Q^{t + 1}}({z_i})}}} } \\ \;\;\;\;\;\;\;\; = L({\theta ^{t + 1}}) \end{array}$$
证毕。
如果定义
$$J(Q,\theta ) = \sum\limits_{i = 1}^N {\sum\limits_{{z_i}} {Q({z_i})\log \frac{{p({x_i},{z_i};\theta )}}{{Q({z_i})}}} } ,$$
则 $J(Q,\theta )$是 $L(\theta )$的下界。EM可以看作是 $J$的坐标上升法，E步固定$\theta$，优化 $Q$，M步固定$Q$优化 $\theta$。

EM和GMM

给定训练样本$\{ {x_i}\} _{i = 1}^N$，隐含类别标签用$z_i$表示。与K-means硬指定不同，GMM认为${z_i} \in \{ 1,2,...,K\}$满足多项式分布${z_i} \sim M(\phi )$，其中$p({z_i} = j) = {\phi _j}$，${\phi _j} \ge 0$，$\sum\nolimits_{j = 1}^K {{\phi _j}} = 1$。假定在给定$z_i$的条件下$x_i$满足多值高斯分布，即 $({x_i}|{z_i} = j) \sim N({\mu _j},{\sigma _j})$，则联合分布$p({x_i},{z_i}) = p({x_i}|{z_i})p({z_i})$。
GMM：

$$\sum\limits_{j = 1}^K {{\phi _j}} N({\mu _j},{\sigma _j})$$
对于每个样例 $x_i$，先从 $K$个类别中按多项式分布抽取一个 $z_i$，然后根据 $z_i$对应高斯分布生成样例 $x_i$，该过程称作混合高斯模型。注意： $z_i$是GMM的隐含随机变量， $\theta = (\phi ,\mu ,\sigma )$是GMM的参数，其中： ${\phi _j}$是样本类别中 ${z_i} = j$的比率， ${\mu _j}$是类别为 $j$的样本特征均值，${\sigma _j}$是类别为 $j$的样本特征方差。对数似然函数：
$$\begin{array}{l} L(\phi ,\mu ,\sigma ) = \log \prod\limits_{i = 1}^N {p({x_i};\phi ,\mu ,\sigma )} \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; = \sum\limits_{i = 1}^N {\log } p({x_i};\phi ,\mu ,\sigma )\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; = \sum\limits_{i = 1}^N {\log } \sum\limits_{{z_i} = 1}^K {p({x_i},{z_i};\phi ,\mu ,\sigma )} \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; = \sum\limits_{i = 1}^N {\log } \sum\limits_{{z_i} = 1}^K {p({x_i}|{z_i};\mu ,\sigma )p({z_i};\phi )} \end{array}$$

该式不存在闭合解。考虑EM的思想，第一步是猜测隐含类别变量 ，第二步是更新其他参数，以获得最大的最大似然估计。
E步：固定$\theta = ({\phi _j},{\mu _j},{\sigma _j})$，求

$$\begin{array}{l} {w_i}(j) = Q({z_i} = j;\theta)= p({z_i} = j|{x_i};\theta )\\ \;\;\;\;\;\;\;\;\; = \frac{{p({x_i},{z_i} = j;\theta )}}{{p({x_i};\theta )}}\\ \;\;\;\;\;\;\;\;\; = \frac{{p({x_i}|{z_i} = j;\mu ,\sigma )p({z_i} = j;\phi )}}{{\sum\nolimits_{l = 1}^K {p({x_i}|{z_i} = l;\mu ,\sigma )p({z_i} = l;\phi )} }}\\ \;\;\;\;\;\;\;\;\; = \frac{{\frac{1}{{\sqrt {2\pi } {\sigma _j}}}\exp \left( { - \frac{{{{({x_i} - {\mu _j})}^2}}}{{2\sigma _j^2}}} \right) \cdot {\phi _j}}}{{\sum\limits_{k = 1}^K {\frac{1}{{\sqrt {2\pi } {\sigma _k}}}\exp \left( { - \frac{{{{({x_i} - {\mu _k})}^2}}}{{2\sigma _k^2}}} \right) \cdot {\phi _k}} }} \end{array}$$

即每个样例$x_i$的隐含类别$z_i$为$j$的概率可以通过后验概率计算得到。对比K-means发现，每个样例分配的类别$z_i$是有一定的概率的，每个样例$i$都要计算属于每一个类别$j$的概率。
M步：固定${w_i}(j)$，最大化

$$\begin{array}{l} J(\theta ) = \sum\limits_{i = 1}^N {\sum\limits_{{z_i}} {Q({z_i})\log \frac{{p({x_i},{z_i};\theta )}}{{Q({z_i})}}} } \\ \;\;\;\;\;\;\;\; = \sum\limits_{i = 1}^N {\sum\limits_{j = 1}^K {Q({z_i} = j)\log \frac{{p({x_i},{z_i};\theta )}}{{Q({z_i} = j)}}} } \\ \;\;\;\;\;\;\;\; = \sum\limits_{i = 1}^N {\sum\limits_{j = 1}^K {{w_i}(j)\log \frac{{\frac{1}{{\sqrt {2\pi } {\sigma _j}}}\exp \left( { - \frac{{{{({x_i} - {\mu _j})}^2}}}{{2\sigma _j^2}}} \right) \cdot {\phi _j}}}{{{w_i}(j)}}} } \end{array}$$

求参数$\theta = ({\phi _j},{\mu _j},{\sigma _j})$。
对${\mu _j},{\sigma _j}$求导，得

$$\begin{array}{l} \frac{{\partial J(\theta )}}{{\partial {\mu _j}}} = \frac{\partial }{{\partial {\mu _j}}}\left( {\sum\limits_{i = 1}^N {\sum\limits_{j = 1}^K {{w_i}(j)\left[ {\log \frac{1}{{{w_i}(j)}} + \log \frac{1}{{\sqrt {2\pi } {\sigma _j}}} - \frac{{{{({x_i} - {\mu _j})}^2}}}{{2\sigma _j^2}} + \log {\phi _j}} \right]} } } \right)\\ \;\;\;\;\;\;\;\;\;\;\; = \sum\limits_{i = 1}^N {{w_i}(j)\left[ {\frac{1}{{\sigma _j^2}}({\mu _j} - {x_i}){\mu _j}} \right]} \end{array}$$
$$\begin{array}{l} \frac{{\partial J(\theta )}}{{\partial {\sigma _j}}} = \frac{\partial }{{\partial {\sigma _j}}}\left( {\sum\limits_{i = 1}^N {\sum\limits_{j = 1}^K {{w_i}(j)\left[ {\log \frac{1}{{{w_i}(j)}} + \log \frac{1}{{\sqrt {2\pi } }} + \log \frac{1}{{{\sigma _j}}} - \frac{{{{({x_i} - {\mu _j})}^2}}}{{2\sigma _j^2}} + \log {\phi _j}} \right]} } } \right)\\ \;\;\;\;\;\;\;\;\;\;\; = \sum\limits_{i = 1}^N {{w_i}(j)\left[ { - \frac{1}{{{\sigma _j}}} + \frac{1}{{\sigma _j^3}}{{({x_i} - {\mu _j})}^2}} \right]} \end{array}$$
令其等于0，解得：
$$\begin{array}{l} {\mu _j} = \frac{{\sum\nolimits_{i = 1}^N {{w_i}(j){x_i}} }}{{\sum\nolimits_{i = 1}^N {{w_i}(j)} }}\\ \sigma _j^2 = \frac{{\sum\nolimits_{i = 1}^N {{w_i}(j){{({x_i} - {\mu _j})}^2}} }}{{\sum\nolimits_{i = 1}^N {{w_i}(j)} }} \end{array}$$
同时，由于 ${\phi _j} \ge 0$，且 $\sum\nolimits_{j = 1}^K {{\phi _j}} = 1$，故建立拉格朗日函数
$$L({\phi _j}) = \sum\limits_{i = 1}^N {{w_i}(j)\log } {\phi _j} + \beta \left( {\sum\nolimits_{j = 1}^K {{\phi _j}} - 1} \right)$$
其中， ${\phi _j}$相关的常数被省略了。对 ${\phi _j}$求导等于0，结合 $\sum\nolimits_{j = 1}^K {{\phi _j}} = 1$可解得：
$$\beta = - \frac{{\sum\limits_{i = 1}^N {{w_i}(j)} }}{{{\phi _j}}} = - \frac{{\sum\limits_{j = 1}^K {\sum\limits_{i = 1}^N {{w_i}(j)} } }}{{\sum\limits_{j = 1}^K {{\phi _j}} }} = - \sum\limits_{i = 1}^N {\sum\limits_{j = 1}^K {{w_i}(j)} } = - \sum\limits_{i = 1}^N 1 = - N$$
即：
$${\phi _j} = \frac{1}{N}\sum\limits_{i = 1}^N {{w_i}(j)}$$
算法过程如下：

循环下面步骤，直到收敛：{
（E步）对于每一个$i$和$j$计算

$${w_i}(j) = \frac{{\frac{1}{{\sqrt {2\pi } {\sigma _j}}}\exp \left( { - \frac{{{{({x_i} - {\mu _j})}^2}}}{{2\sigma _j^2}}} \right) \cdot {\phi _j}}}{{\sum\limits_{k = 1}^K {\frac{1}{{\sqrt {2\pi } {\sigma _k}}}\exp \left( { - \frac{{{{({x_i} - {\mu _k})}^2}}}{{2\sigma _k^2}}} \right) \cdot {\phi _k}} }}$$
（M步），更新参数：
$$\begin{array}{l} {\phi _j} = \frac{1}{N}\sum\limits_{i = 1}^N {{w_i}(j)} \\ {\mu _j} = \frac{{\sum\nolimits_{i = 1}^N {{w_i}(j){x_i}} }}{{\sum\nolimits_{i = 1}^N {{w_i}(j)} }}\\ \sigma _j^2 = \frac{{\sum\nolimits_{i = 1}^N {{w_i}(j){{({x_i} - {\mu _j})}^2}} }}{{\sum\nolimits_{i = 1}^N {{w_i}(j)} }} \end{array}$$
}

E步中，将其他参数$\phi ,\mu ,\sigma$看作常量，计算$z_i$的后验概率，也就是估计隐含类别变量。估计好后，利用上面的公式重新计算其他参数，计算好后发现最大化最大似然估计时，${w_i}(j)$值又不对了，需要重新计算，周而复始，直至收敛。

GMM示例

给定由三个高斯的组成的混合模型（$K = 3$），其中：第一个高斯的均值$\mu_1=5$、方差$\sigma_1=3$ ，第二个高斯的均值$\mu_2=20$ 、方差$\sigma_2=5$ ，第三个高斯的均值$\mu_3=50$、方差$\sigma_3=10$ ，且三个高斯出现的概率为$\phi_1=0.2$、$\phi_1=0.4$ 、$\phi_1=0.4$。由该GMM生成10000个样本（ $N=10000$），并根据样本和给定的初始参数$\phi=[0.33,0.33,0.34]$ 、$\mu=[0,5,10]$、$\sigma=[5,5,5]$，迭代学习50次结果后收敛到$\phi=[0.1969,0.4082,0.3950]$ 、$\mu=[4.9686,20.0338,50.0925]$、$\sigma=[3.0639,5.0977,10.1096]$ 。下图给出了10000个样本的分布和50次迭代学习参数变化情况。

Matlab代码：

% 生成过程
phi1 = 0.2; mu1 = 5; sigma1 = 3;
phi2 = 0.4; mu2 = 20; sigma2 = 5;
phi3 = 0.4; mu3 = 50; sigma3 = 10;

N = 10000;
x = zeros(N,1);
for i = 1 : N
    rate = rand;
    if rate <= phi1
        x(i) = normrnd(mu1,sigma1);
    elseif rate <= phi1+phi2
        x(i) = normrnd(mu2,sigma2);
    else
        x(i) = normrnd(mu3,sigma3);
    end
end

figure(1); subplot(2,2,1); hist(x,1000);

% 学习过程
mu = [0, 5, 10];
sigma = [5, 5, 5];
phi = [0.33, 0.33, 0.34];
w = zeros(N,3);

T = 50;
mu_ = zeros(T+1,3);
sigma_ = zeros(T+1,3);
phi_ = zeros(T+1,3);
mu_(1,:) = mu;
sigma_(1,:) = sigma;
phi_(1,:) = phi;
for t = 1 : T
    % Expectation
    for k = 1 : 3
        w(:,k) = phi(k)*normpdf(x,mu(k),sigma(k));
    end
    w = w./repmat(sum(w,2),[1 3]);

    % Maximization
    for k = 1 : 3
        mu(k) = w(:,k)'*x / sum(w(:,k));
        sigma(k) = sqrt(w(:,k)'*((x-mu(k)).*(x-mu(k))) / sum(w(:,k)));
        phi(k) = sum(w(:,k)) / N;
    end
    mu_(t+1,:) = mu;
    sigma_(t+1,:) = sigma;
    phi_(t+1,:) = phi;
end
figure(1); subplot(2,2,2); plot(phi_); title('\phi');
figure(1); subplot(2,2,3); plot(mu_); title('\mu');
figure(1); subplot(2,2,4); plot(sigma_); title('\sigma');