详解EM算法与混合高斯模型(Gaussian mixture model, GMM)

时间  2020-08-08 标签 详解 算法 混合 高斯 模型 gaussian mixture model gmm

 　　最近在看晓川老（shi）师(shu)的博士论文，接触了混合高斯模型（Gaussian mixture model, GMM）和EM(Expectation Maximization)算法，不由被论文中庞大的数学公式所吓退。本文经过查阅相关资料，在复杂巧妙的推理公式中融入了本身的理解，详细梳理了混合高斯模型和EM算法。html

1 单高斯模型(Gaussian single model, GSM）

 　　简单回顾一下几率论讲过的高斯模型。
　　高斯模型是一种经常使用的变量分布模型，在数理统计领域有着普遍的应用（……好吧读了这么多年书没白费，教科书般的话语已植入骨髓）。一维高斯分布的几率密度函数以下：
(1) f ( x ) = 1 2 π σ exp ⁡ ( − ( x − μ ) 2 2 σ 2 ) f(x) = \frac{1}{{\sqrt {2\pi } \sigma }}\exp ( - \frac{{{{(x - \mu )}^2}}}{{2{\sigma ^2}}})\tag{1} 　
　　 μ \mu 和 σ 2 {\sigma ^2} 分别是高斯分布的均值和方差。
　　譬如将男生身高视为变量X, 假设男生的身高服从高斯分布，则 X ∼ N ( μ , σ 2 ) X \sim N(\mu ,{\sigma ^2}) ，女生亦如此。只是男女生身高分布可能具备不一样的均值和方差。图1是从谷歌图库中搜索到的男女生身高分布图，来源不清，我的以为男生的均值身高虚高……四个记号分别表示0~3 σ \sigma 准则。
图1 男女生身高分布差别web

图1 男女生身高分布差别

 　　多维变量 X = ( x 1 , x 2 , . . . x n ) X = ({x_1},{x_2},...{x_n}) 的联合几率密度函数为：
(2) f ( X ) = 1 ( 2 π ) d / 2 ∣ Σ ∣ 1 / 2 exp ⁡ [ − 1 2 ( X − u ) T Σ − 1 ( X − u ) ] , X = ( x 1 , x 2 . . . x n ) f(X) = \frac{1}{{{{(2\pi )}^{d/2}}{{\left| \Sigma \right|}^{1/2}}}}\exp [ - \frac{1}{2}{(X - u)^T}{\Sigma ^{ - 1}}(X - u)],X = ({x_1},{x_2}...{x_n})\tag{2}
　　其中：
　　d：变量维度。对于二维高斯分布，有d=2;
　　 u = ( u 1 u 2 . . . u n ) u = \left( \begin{array}{l} {u_1}\\ {u_2}\\ ...\\ {u_n} \end{array} \right) ：各维变量的均值；
　　 Σ \Sigma ：协方差矩阵，描述各维变量之间的相关度。对于二维高斯分布，有：
(3) Σ = [ δ 11 δ 12 δ 21 δ 22 ] \Sigma= \left[ \begin{matrix} \delta _{11} & \delta _{12}\\ \delta _{21}& \delta _{22} \end{matrix} \right] \tag{3} 算法

这里写图片描述
图2 二维高斯数据分布

 　　图2是二维高斯分布产生的数据示例，参数设定为： u = ( 0 0 ) , Σ = [ 1 0.8 0.8 5 ] u = \left( \begin{array}{l} 0\\ 0 \end{array} \right),\Sigma= \left[ \begin{matrix} 1 & 0.8\\ 0.8&5 \end{matrix} \right] 。关于二维高斯分布的参数设定对为高斯曲面的影响，这篇文章二维高斯分布（Two-dimensional Gaussian distribution）的参数分析有说起，主要是为了下文理解混合高斯分布作铺垫。服从二维高斯分布的数据主要集中在一个椭圆内部，服从三维的数据集中在一个椭球内部。app

2 混合高斯模型（Gaussian mixture model, GMM）

2.1 为何要有混合高斯模型

先来看一组数据。　
　　　　　
　　ide

这里写图片描述
图3 混合高斯分布所产生数据

 　　若是咱们假设这组数据是由某个高斯分布产生的，利用极大似然估计（后文还会说起）对这个高斯分布作参数估计，获得一个最佳的高斯分布模型以下。svg

这里写图片描述
图4 用单高斯模型对样本做分析不合理示意

 　　有什么问题吗？通常来说越靠近椭圆的中心样本出现的几率越大，这是由几率密度函数决定的，可是这个高斯分布的椭圆中心的样本量却极少。显然样本服从单高斯分布的假设并不合理。单高斯模型没法产生这样的样本。
　　实际上，这是用两个不一样的高斯分布模型产生的数据。函数

这里写图片描述
图5 混合高斯模型对样本做分析示意

 　　正当单高斯模型抓耳挠腮的时候，混合高斯模型就大摇大摆地进场了。它经过求解两个高斯模型，并经过必定的权重将两个高斯模型融合成一个模型，即最终的混合高斯模型。这个混合高斯模型能够产生这样的样本。
　　更通常化的描述为：假设混合高斯模型由K个高斯模型组成（即数据包含K个类），则GMM的几率密度函数以下：
(3) p ( x ) = ∑ k = 1 K p ( k ) p ( x ∣ k ) = ∑ k = 1 K π k N ( x ∣ u k , Σ k ) p(x)=\sum\limits_{k = 1}^K {p(k)p(x|k) = } \sum\limits_{k = 1}^K {{\pi _k}N(x|{u_k},{\Sigma _k})} \tag{3}
　　其中， p ( x ∣ k ) = N ( x ∣ u k , Σ k ) p(x|k) = N(x|{u_k},{\Sigma _k}) 是第k个高斯模型的几率密度函数，能够当作选定第k个模型后，该模型产生x的几率； p ( k ) = π k p(k) = {\pi _k} 是第k个高斯模型的权重，称做选择第k个模型的先验几率，且知足 ∑ k = 1 K π k = 1 \sum\limits_{k = 1}^K {{\pi _k}} = 1 。
　　因此，混合高斯模型并非什么新奇的东西，它的本质就是融合几个单高斯模型，来使得模型更加复杂，从而产生更复杂的样本。理论上，若是某个混合高斯模型融合的高斯模型个数足够多，它们之间的权重设定得足够合理，这个混合模型能够拟合任意分布的样本。优化

2.2 直观上理解混合高斯模型

 　　下面经过几张图片来帮助理解混合高斯模型。
　　首先从简单的一维混合高斯模型提及。
　spa

这里写图片描述
图6 一维混合高斯模型

 　　在图6中，y1,y2和y3分别表示三个一维高斯模型，他们的参数设定如图所示。y4表示将三个模型的几率密度函数直接相加，注意的是这并非一个混合高斯模型，由于不知足 ∑ k = 1 K π k = 1 \sum\limits_{k = 1}^K {{\pi _k}} = 1 的条件。而y5和y6分别是由三个相同的高斯模型融合生成的不一样混合模型。因而可知，调整权重将极大影响混合模型的几率密度函数曲线。另外一方面也能够直观地理解混合高斯模型能够更好地拟合样本的缘由：它有更复杂更多变的几率密度函数曲线。理论上，混合高斯模型的几率密度函数曲线能够是任意形状的非线性函数。
　　下面再给出一个二维空间3个高斯模型混合的例子。
　.net

这里写图片描述
(a) 3个类别高斯分布截面轮廓线

这里写图片描述
(b) 混合高斯分布截面轮廓线

这里写图片描述
© 二维混合高斯分布几率密度函数图

图7 二维混合高斯模型

 　　(a) 图表示的是3个高斯模型的截面轮廓图，3个模型的权重系数已在图中注明，由截面轮廓图可知3个模型之间存在不少重叠区域。其实这也正是混合高斯模型所但愿的。由于若是它们之间的重叠区域较少，那么生成的混合高斯模型通常较为简单，难以生成较为复杂的样本。
　　设定好了3个高斯模型和它们之间的权重系数以后，就能够肯定二维混合高斯分布的几率密度函数曲面，如图©所示。图(b)是对于图©几率密度曲面的截面轮廓线。从图7也能够看出，整个混合高斯分布曲面相对比于单高斯分布曲面已经异常复杂。实际上，经过调整混合高斯分布的系数 ( π , μ , Σ ) (\pi ,\mu ,\Sigma ) ，可使得图©的几率密度曲面去拟合任意的三维曲面，从而采样生成所须要的数据样本。

3 极大似然估计(Maximum Likehood Estimate, MLE)（最大化对数似然函数）

● 最大化对数似然函数（log-likelihood function）的意义

 　　首先直观化地解释一下最大化对数似然函数要解决的是什么问题。
　　假设咱们采样获得一组样本 y t {y_t} ,并且咱们知道变量Y服从高斯分布（本文只说起高斯分布，其余变量分布模型相似），数学形式表示为 Y ∼ N ( μ , Σ ) Y \sim N(\mu ,\Sigma ) 。采样的样本如图8所示，咱们的目的就是找到一个合适的高斯分布（也就是肯定高斯分布的参数 μ , Σ \mu ,\Sigma ），使得这个高斯分布能产生这组样本的可能性尽量大。
　

这里写图片描述
图8 最大化似然函数的意义

 　　那怎么找到这个合适的高斯分布呢（在图8中的表示就是1~4哪一个分布较为合适）？这时候似然函数就闪亮登场了。
　　似然函数数学化：设有样本集 Y = y 1 , y 2 . . . y N Y = {y_1},{y_2}...{y_N} 。 p ( y n ∣ μ , Σ ) p({y_n}|\mu ,\Sigma ) 是高斯分布的几率分布函数，表示变量 Y = y n Y = {y_n} 的几率。假设样本的抽样是独立的，那么咱们同时抽到这N个样本的几率是抽到每一个样本几率的乘积，也就是样本集Y的联合几率。此联合几率即为似然函数：
(4) L ( μ , Σ ) = L ( y 1 , y 2 . . . y N ; μ , Σ ) = ∏ n = 1 N p ( y n ; μ , Σ ) L(\mu ,\Sigma ) = L({y_1},{y_2}...{y_N};\mu ,\Sigma ) = \prod\limits_{n = 1}^N {p({y_n};\mu ,\Sigma )}\tag{4}
　　对式子(4)进行求导并令导数为0（即最大化似然函数，通常还会先转化为对数似然函数再最大化），所求出的参数就是最佳的高斯分布对应的参数。
　　因此最大化似然函数的意义就是：经过使得样本集的联合几率最大来对参数进行估计，从而选择最佳的分布模型。
　　对于图8产生的样本用最大化似然函数的方法，最终能够获得序号1对应的高斯分布模型是最佳的模型。

4 EM算法（最大化Q函数）

4.1 为何要有EM算法（EM算法与极大似然估计分别适用于什么问题）

● 尝试用极大似然估计的方法来解GMM模型

 　　解GMM模型，实际上就是肯定GMM模型的参数 ( μ , Σ , π ) (\mu ,\Sigma ,\pi ) ，使得由这组参数肯定的GMM模型最有可能产生采样的样本。
　　先试试看用极大似然估计的方法来解GMM模型会出现什么样的问题。
　　如第3小节所述，要利用极大似然估计求解模型最重要的一步就是求出似然函数，即样本集出现的联合几率。而对于混合高斯模型，如何求解某个样本 y t {y_t} 的几率？显然咱们得先知道这个样原本源于哪一类高斯模型，而后求这个高斯模型生成这个样本的几率 p ( y t ) p({y_t}) 。
　　可是问题来了：咱们只有样本。不知道样本到底来源于哪一类的高斯模型。那么如何求解样本的生成几率 p ( y t ) p({y_t}) ？
　　先引入一个隐变量 γ \gamma 。它是一个K维二值随机变量，在它的K维取值中只有某个特定的元素 γ k {\gamma _k} 的取值为1，其它元素的取值为0。实际上，隐变量描述的就是：每一次采样，选择第k个高斯模型的几率，故有：
(5) p ( γ k = 1 ) = π k p({\gamma _k} = 1) = {\pi _k}\tag{5}
　　当给定了 γ \gamma 的一个特定的值以后（也就是知道了这个样本从哪个高斯模型进行采样），能够获得样本y的条件分布是一个高斯分布，知足：
(6) p ( y ∣ γ k = 1 ) = N ( y ∣ μ k , Σ k ) p(y|{\gamma _k} = 1) = N(y|{\mu _k},{\Sigma _k})\tag{6}
　　而实际上，每一个样本究竟是从这K个高斯模型中哪一个模型进行采样的，是都有可能的。故样本y的几率为：
(7) p ( y ) = ∑ γ p ( γ ) p ( y ∣ γ ) = ∑ k = 1 K π k N ( y ∣ μ k , Σ k ) p(y) = \sum\nolimits_\gamma {p(\gamma )} p(y|\gamma ){\rm{ = }}\sum\limits_{{\rm{k}} = 1}^K {{\pi _k}N(y|{\mu _k},{\Sigma _k})} \tag{7}
　　样本集Y(n个样本点)的联合几率为：
(8) L ( μ , Σ , π ) = L ( y 1 , y 2 . . . y N ; μ , Σ , π ) = ∏ n = 1 N p ( y n ; μ , Σ , π ) = ∏ n = 1 N ∑ k = 1 K π k N ( y n ∣ μ k , Σ k ) L(\mu ,\Sigma ,\pi ) = L({y_1},{y_2}...{y_N};\mu ,\Sigma ,\pi ) = \prod\limits_{n = 1}^N {p({y_n};\mu ,\Sigma ,\pi )} = \prod\limits_{n = 1}^N {\sum\limits_{{\rm{k}} = 1}^K {{\pi _k}N({y_n}|{\mu _k},{\Sigma _k})} } \tag{8}
　　对数似然函数表示为：
(9) ln ⁡ L ( μ , Σ , π ) = ∑ n = 1 N ln ⁡ ∑ k = 1 K π k N ( y n ∣ μ k , Σ k ) \ln L(\mu ,\Sigma ,\pi ) = \sum\limits_{n = 1}^N {\ln \sum\limits_{{\rm{k}} = 1}^K {{\pi _k}N({y_n}|{\mu _k},{\Sigma _k})} } \tag{9}
　　好了，而后求导，令导数为0，获得模型参数 ( μ , Σ , π ) (\mu ,\Sigma ,\pi ) 。
　　貌似问题已经解决了，喜大普奔。
　　然而仔细观察能够发现，对数似然函数里面，对数里面还有求和。实际上没有办法经过求导的方法来求这个对数似然函数的最大值。
　　MLE（极大似然估计）略显沮丧。这时候EM算法走过来，安慰着说：兄弟别哭，老哥帮你。

● 极大似然估计与EM算法适用问题分析

 　　下面先阐述一下极大似然估计与EM算法分别适用于解决什么样的问题。
　

这里写图片描述
图9 极大似然估计适用问题

这里写图片描述
图10 EM算法适用问题

 　　若是咱们已经清楚了某个变量服从的高斯分布，并且经过采样获得了这个变量的样本数据，想求高斯分布的参数，这时候极大似然估计能够胜任这个任务；而若是咱们要求解的是一个混合模型，只知道混合模型中各个类的分布模型（譬如都是高斯分布）和对应的采样数据，而不知道这些采样数据分别来源于哪一类（隐变量），那这时候就能够借鉴EM算法。EM算法能够用于解决数据缺失的参数估计问题（隐变量的存在实际上就是数据缺失问题，缺失了各个样原本源于哪一类的记录）。
　　下面将介绍EM算法的两个步骤：E-step(expectation-step，指望步)和M-step(Maximization-step,最大化步);

4.2 E-step

咱们现有样本集 Y = ( y 1 , y 2 . . . y T ) Y = ({y_1},{y_2}...{y_T}) ，经过隐变量 γ t , k {\gamma _{t,k}} （表示 y t {y_t} 这个样原本源于第k个模型）的引入，能够将数据展开成彻底数据：
( y t , γ t , 1 , γ t , 2 . . . γ t , K ) , t = 1 , 2... T ({y_t},{\gamma _t}_{,1},{\gamma _{t,2}}...{\gamma _{t,K}}),t = 1,2...T

所谓的彻底数据，就是不缺失的数据。只有样本集 Y = ( y 1 , y 2 . . . y T ) Y = ({y_1},{y_2}...{y_T}) 的数据是不完整的，存在信息缺失的。若 y t {y_t} 由第1类采样而来，则有 γ t , 1 = 1 , γ t , 2 = 0... γ t , K = 0 {\gamma _t}_{,1} = 1,{\gamma _{t,2}} = 0...{\gamma _{t,K}} = 0 ，表示为 ( y t , 1 , 0 , . . . 0 ) ({y_t},1,0,...0) 。
　　因此要求能采到这组数据的可能性，须要分两步走：①第t个样本由哪一类产生？②若是第t个样本由第k类产生，那么第k类产生第t个样本的几率为多少？
　　综合考虑上面两步，有了彻底数据的似然函数：
(10) p ( y , γ ∣ μ , Σ , π ) = ∏ t = 1 T p ( y t , γ t , 1 , γ t , 2 . . . γ t , K ∣ μ , Σ , π ) 　 　 　 　 　 　 　 = ∏ t = 1 T ∏ k = 1 K ( π k N ( y t ; μ k , Σ k ) ) γ t , k 　 　 　 　 　 　 　 = ∏ k = 1 K π k ∑ t = 1 T γ t , k ∏ t = 1 T ( N ( y t ; μ k , Σ k ) ) γ t , k \begin{array}{l} p(y,\gamma |\mu ,\Sigma ,\pi ) = \prod\limits_{t = 1}^T {p({y_t},{\gamma _t}_{,1},{\gamma _{t,2}}...{\gamma _{t,K}}|\mu ,\Sigma ,\pi )} \\ {\rm{ 　　　　　　　}} = \prod\limits_{t = 1}^T {\prod\limits_{k = 1}^K {{{({\pi _k}N({y_t};{\mu _k},{\Sigma _k}))}^{{\gamma _{t,k}}}}} } \\ {\rm{ 　　　　　　　}} = \prod\limits_{k = 1}^K {\pi _k^{\sum\nolimits_{t = 1}^T {^{{\gamma _{t,k}}}} }} \prod\limits_{t = 1}^T {{{(N({y_t};{\mu _k},{\Sigma _k}))}^{{\gamma _{t,k}}}}} \end{array} \tag{10}
　　第1个等号到第2个等号的理解：若 y t {y_t} 由第1类采样而来，则有 γ t , 1 = 1 , γ t , 2 = 0... γ t , K = 0 {\gamma _t}_{,1} = 1,{\gamma _{t,2}} = 0...{\gamma _{t,K}} = 0 ，
(11) p ( y t , γ t , 1 , γ t , 2 . . . γ t , K ∣ μ , Σ , π ) = ∏ k = 1 K ( π k N ( y t ; μ k , Σ k ) ) γ t , k 　 　 　 　 　 　 　 　 　 　 　 　 　 　 = ( π 1 N ( y t ; μ 1 , Σ 1 ) ) γ t , 1 ( π 2 N ( y t ; μ 2 , Σ 2 ) ) γ t , 2 . . . ( π K N ( y t ; μ K , Σ K ) ) γ t , K 　 　 　 　 　 　 　 　 　 　 　 　 　 　 = ( π 1 N ( y t ; μ 1 , Σ 1 ) ) 1 ( π 2 N ( y t ; μ 2 , Σ 2 ) ) 0 . . . ( π K N ( y t ; μ K , Σ K ) ) 0 　 　 　 　 　 　 　 　 　 　 　 　 　 　 = π 1 N ( y t ; μ 1 , Σ 1 ) \begin{array}{l} p({y_t},{\gamma _t}_{,1},{\gamma _{t,2}}...{\gamma _{t,K}}|\mu ,\Sigma ,\pi ){\rm{ = }}\prod\limits_{k = 1}^K {{{({\pi _k}N({y_t};{\mu _k},{\Sigma _k}))}^{{\gamma _{t,k}}}}} \\ {\rm{ 　　　　　　　　　　　　　　= }}{({\pi _1}N({y_t};{\mu _1},{\Sigma _1}))^{{\gamma _{t,1}}}}{({\pi _2}N({y_t};{\mu _2},{\Sigma _2}))^{{\gamma _{t,2}}}}...{({\pi _K}N({y_t};{\mu _K},{\Sigma _K}))^{{\gamma _{t,K}}}}\\ {\rm{ 　　　　　　　　　　　　　　 }} = {({\pi _1}N({y_t};{\mu _1},{\Sigma _1}))^1}{({\pi _2}N({y_t};{\mu _2},{\Sigma _2}))^0}...{({\pi _K}N({y_t};{\mu _K},{\Sigma _K}))^0}\\ {\rm{ 　　　　　　　　　　　　　　 }} = {\pi _1}N({y_t};{\mu _1},{\Sigma _1}) \end{array} \tag{11}
　　注意式子(11)与式子(7)的差异：若是求 p ( y t ) p({y_t}) 则须要考虑 y t {y_t} 有可能来源于k个类；而若是求的是 p ( y t , γ t , 1 , γ t , 2 . . . γ t , K ) p({y_t},{\gamma _t}_{,1},{\gamma _{t,2}}...{\gamma _{t,K}}) 则已经限定了 y t {y_t} 只会来源于某个类。
　　第2个等式到第3个等式的理解：先交换累乘符号。因为 π k {\pi _k} 与t无关，故而能够从内部的累乘符号中提取出来。
　　实际上彻底数据的似然函数描述的就是采集到这些样本的可能性。
　　彻底数据的对数似然函数为：
(12) ln ⁡ p ( y , γ ∣ μ , Σ , π ) = ∑ k = 1 K ( ∑ t = 1 T γ t , k ) ln ⁡ π k + ∑ t = 1 T γ t , k ( − ln ⁡ ( 2 π ) − 1 2 ln ⁡ ∣ Σ k ∣ − 1 2 ( y t − μ t ) T ( Σ k ) − 1 ( y t − μ t ) ) \ln p(y,\gamma |\mu ,\Sigma ,\pi ) = \sum\limits_{k = 1}^K {(\sum\limits_{t = 1}^T {{\gamma _{t,k}}} )\ln {\pi _k}} + \sum\limits_{t = 1}^T {{\gamma _{t,k}}} ( - \ln (2\pi ) - \frac{1}{2}\ln \left| {{\Sigma _k}} \right| - \frac{1}{2}{({y_t} - {\mu _t})^T}{({\Sigma _k})^{ - 1}}({y_t} - {\mu _t})) \tag{12}
　　这一步应该没啥问题吧。。。注意的是，此处考虑的是二维高斯分布的状况，对应于式子(2)中的d=2。
　　咱们的目标就是找出一组参数 ( μ ∗ , Σ ∗ , π ∗ ) (\mu *,\Sigma *,\pi *) 使得 ln ⁡ p ( y , γ ∣ μ , Σ , π ) \ln p(y,\gamma |\mu ,\Sigma ,\pi ) 最大。
　　那么问题来了： ln ⁡ p ( y , γ ∣ μ , Σ , π ) \ln p(y,\gamma |\mu ,\Sigma ,\pi ) 中含有隐变量 γ \gamma ， γ \gamma 的存在使得咱们无法最大化 ln ⁡ p ( y , γ ∣ μ , Σ , π ) \ln p(y,\gamma |\mu ,\Sigma ,\pi ) 。若是咱们知道了 γ \gamma ，那么最大化 ln ⁡ p ( y , γ ∣ μ , Σ , π ) \ln p(y,\gamma |\mu ,\Sigma ,\pi ) 就显得水到渠成。
　　可是坑爹的就是：咱们只有采样数据 y t {y_t} ， γ \gamma 未知。
　　那么怎么办呢？对 γ \gamma 来一个估计。
　　猜测咱们给了一组起始参数 ( μ 0 , Σ 0 , π 0 ) ({\mu ^0},{\Sigma ^0},{\pi ^0}) 或者优化过的第i次迭代的参数 ( μ i , Σ i , π i ) ({\mu ^i},{\Sigma ^i},{\pi ^i}) ，也就是说每个高斯分布的参数咱们都有了， γ \gamma 作的事不就是决定每一个样本由哪个高斯分布产生的嘛，有了每一个高斯分布的参数那咱们就能够猜测每一个样本最有可能来源于哪一个高斯分布没错吧！Done!
　　为此咱们不最大化 ln ⁡ p ( y , γ ∣ μ , Σ , π ) \ln p(y,\gamma |\mu ,\Sigma ,\pi ) （也没法最大化它），而是最大化Q函数。Q函数以下：
Q ( μ , Σ , π , μ i , Σ i , π i ) = E γ [ ln ⁡ p ( y , γ ∣ μ , Σ , π ) ∣ Y , μ i , Σ i , π i ] = E γ [ ∑ k = 1 K ( ∑ t = 1 T γ t , k ∣ y t , μ i , Σ i , π i ) ln ⁡ π k + ∑ t = 1 T ( γ t , k ∣ y t , μ i , Σ i , π i ) ( − ln ⁡ ( 2 π ) − 1 2 ln ⁡ ∣ Σ k ∣ − 1 2 ( y t − μ t ) T ( Σ k ) − 1 ( y t − μ t ) ) ] = ∑ k = 1 K ( ∑ t = 1 T E ( γ t , k ∣ y t , μ i , Σ i , π i ) ln ⁡ π k + ∑ t = 1 T E ( γ t , k ∣ y t , μ i , Σ i , π i ) ( − ln ⁡ ( 2 π ) − 1 2 ln ⁡ ∣ Σ k ∣ − 1 2 ( y t − μ t ) T ( Σ k ) − 1 ( y t − μ t ) ) ) \begin{array}{l} Q(\mu ,\Sigma ,\pi ,{\mu ^i},{\Sigma ^i},{\pi ^i}) = {E_\gamma }[\ln p(y,\gamma |\mu ,\Sigma ,\pi )|Y,{\mu ^i},{\Sigma ^i},{\pi ^i}]\\ = {E_\gamma }[\sum\limits_{k = 1}^K {(\sum\limits_{t = 1}^T {{\gamma _{t,k}}|{y_t},{\mu ^i},{\Sigma ^i},{\pi ^i}} )\ln {\pi _k}} + \sum\limits_{t = 1}^T {({\gamma _{t,k}}|{y_t},{\mu ^i},{\Sigma ^i},{\pi ^i})} ( - \ln (2\pi ) - \frac{1}{2}\ln \left| {{\Sigma _k}} \right| - \frac{1}{2}{({y_t} - {\mu _t})^T}{({\Sigma _k})^{ - 1}}({y_t} - {\mu _t}))]\\ = \sum\limits_{k = 1}^K {(\sum\limits_{t = 1}^T {E({\gamma _{t,k}}|{y_t},{\mu ^i},{\Sigma ^i},{\pi ^i})} \ln {\pi _k}} + \sum\limits_{t = 1}^T {E({\gamma _{t,k}}|{y_t},{\mu ^i},{\Sigma ^i},{\pi ^i})} ( - \ln (2\pi ) - \frac{1}{2}\ln \left| {{\Sigma _k}} \right| - \frac{1}{2}{({y_t} - {\mu _t})^T}{({\Sigma _k})^{ - 1}}({y_t} - {\mu _t}))) \end{array}
　　其中， E ( γ t , k ∣ y t , μ i , Σ i , π i ) E({\gamma _{t,k}}|{y_t},{\mu ^i},{\Sigma ^i},{\pi ^i}) 就是对 γ \gamma 的估计：
(14) E ( γ t , k ∣ y t , μ i , Σ i , π i ) = p ( γ t , k = 1 ∣ y t , μ i , Σ i , π i ) 　 　 　 　 　 　 　 　 　 　 = p ( γ t , k = 1 , y t ∣ μ i , Σ i , π i ) p ( y t ) 　 　 　 　 　 　 　 　 　 　 = p ( γ t , k = 1 , y t ∣ μ i , Σ i , π i ) ∑ k = 1 K p ( γ t , k = 1 , y t ∣ μ i , Σ i , π i ) 　 　 　 　 　 　 　 　 　 　 = p ( y t ∣ γ t , k = 1 , μ i , Σ i , π i ) p ( γ t , k = 1 ∣ μ i , Σ i , π i ) ∑ k = 1 K p ( y t ∣ γ t , k = 1 , μ i , Σ i , π i ) p ( γ t , k = 1 ∣ μ i , Σ i , π i ) 　 　 　 　 　 　 　 　 　 　 = π k i N ( y t ; μ k i , Σ k i ) ∑ k = 1 K π k i N ( y t ; μ k i , Σ k i ) \begin{array}{l} E({\gamma _{t,k}}|{y_t},{\mu ^i},{\Sigma ^i},{\pi ^i}){\rm{ = }}p({\gamma _{t,k}} = 1|{y_t},{\mu ^i},{\Sigma ^i},{\pi ^i})\\ 　　　　　　　　　　= \frac{{p({\gamma _{t,k}} = 1,{y_t}|{\mu ^i},{\Sigma ^i},{\pi ^i})}}{{p({y_t})}}\\ 　　　　　　　　　　= \frac{{p({\gamma _{t,k}} = 1,{y_t}|{\mu ^i},{\Sigma ^i},{\pi ^i})}}{{\sum\nolimits_{k = 1}^K {p({\gamma _{t,k}} = 1,{y_t}|{\mu ^i},{\Sigma ^i},{\pi ^i})} }}\\ 　　　　　　　　　　= \frac{{p({y_t}|{\gamma _{t,k}} = 1,{\mu ^i},{\Sigma ^i},{\pi ^i})p({\gamma _{t,k}} = 1|{\mu ^i},{\Sigma ^i},{\pi ^i})}}{{\sum\nolimits_{k = 1}^K {p({y_t}|{\gamma _{t,k}} = 1,{\mu ^i},{\Sigma ^i},{\pi ^i})p({\gamma _{t,k}} = 1|{\mu ^i},{\Sigma ^i},{\pi ^i})} }}\\ 　　　　　　　　　　= \frac{{\pi _k^iN({y_t};\mu _k^i,\Sigma _k^i)}}{{\sum\nolimits_{k = 1}^K {\pi _k^iN({y_t};\mu _k^i,\Sigma _k^i)} }} \end{array} \tag{14}
　　这公式是否是很可怕？？别急，带上几点声明，再去看公式就很好理解了！
　　① Q函数描述的其实就是在给定 ( μ i , Σ i , π i ) ({\mu ^i},{\Sigma ^i},{\pi ^i}) 参数下，先对样本Y作一个最有可能的划分（每一个样原本源于各个类的可能性，即对 γ \gamma 的估计 E ( γ t , k ∣ y t , μ i , Σ i , π i ) E({\gamma _{t,k}}|{y_t},{\mu ^i},{\Sigma ^i},{\pi ^i}) ），再描述可以产生这组样本的可能性（Q函数）；
　　② 有了对于 γ \gamma 的估计以后，Q函数只和样本有关（传统意义上的似然函数亦如此，彻底数据的似然函数还与 γ \gamma 有关），而再也不含有隐变量，从而使得最大化Q函数成为可能；
　　③ 最大化Q函数的过程实则就是使得可以产生这组样本的可能性最大，与最大化似然函数的思路一模一样。

4.3 M-step

有个Q函数，就能够对Q函数进行最大化，获得下一次迭代的模型参数了，即：
(15) μ i + 1 , Σ i + 1 , π i + 1 = arg ⁡ max ⁡ Q ( μ , Σ , π , μ i , Σ i , π i ) {\mu ^{i{\rm{ + }}1}},{\Sigma ^{i{\rm{ + }}1}},{\pi ^{i{\rm{ + }}1}}{\rm{ = }}\arg \max Q(\mu ,\Sigma ,\pi ,{\mu ^i},{\Sigma ^i},{\pi ^i}) \tag{15}
　　对Q函数进行求导，并另其导数为0，可得：
(16) μ k i + 1 = ∑ t = 1 T π k i N ( y t ; μ k i , Σ k i ) ∑ k = 1 K π k i N ( y t ; μ k i , Σ k i ) y t E ( γ t , k ∣ y t , μ i , Σ i , π i ) , k = 1 , 2... K \mu _k^{i + 1} = \frac{{\sum\nolimits_{t = 1}^T {\frac{{\pi _k^iN({y_t};\mu _k^i,\Sigma _k^i)}}{{\sum\nolimits_{k = 1}^K {\pi _k^iN({y_t};\mu _k^i,\Sigma _k^i)} }}} {y_t}}}{{E({\gamma _{t,k}}|{y_t},{\mu ^i},{\Sigma ^i},{\pi ^i})}},k = 1,2...K \tag{16}
(17) Σ k i + 1 = ∑ t = 1 T π k i N ( y t ; μ k i , Σ k i ) ∑ k = 1 K π k i N ( y t ; μ k i , Σ k i ) ( y t − μ k i ) 2 E ( γ t , k ∣ y t , μ i , Σ i , π i ) , k = 1 , 2... K \Sigma _k^{i + 1} = \frac{{\sum\nolimits_{t = 1}^T {\frac{{\pi _k^iN({y_t};\mu _k^i,\Sigma _k^i)}}{{\sum\nolimits_{k = 1}^K {\pi _k^iN({y_t};\mu _k^i,\Sigma _k^i)} }}} {{({y_t} - \mu _k^i)}^2}}}{{E({\gamma _{t,k}}|{y_t},{\mu ^i},{\Sigma ^i},{\pi ^i})}},k = 1,2...K\tag{17}
(18) π k i + 1 = E ( γ t , k ∣ y t , μ i , Σ i , π i ) T , k = 1 , 2... K \pi _k^{i + 1} = \frac{{E({\gamma _{t,k}}|{y_t},{\mu ^i},{\Sigma ^i},{\pi ^i})}}{T},k = 1,2...K\tag{18}
　　其中 μ k i + 1 , Σ k i + 1 , π k i + 1 \mu _k^{i + 1},\Sigma _k^{i + 1},\pi _k^{i + 1} 分别表示第（i+1）次迭代，第k个类的均值，协方差矩阵和所占的权重。

4.4 一个例子梳理EM算法的整个过程

 　　EM算法的核心思想是：经过迭代的过程来找到一组最优的参数 ( μ ∗ , Σ ∗ , π ∗ ) (\mu *,\Sigma *,\pi *) ，使得这组参数表示的模型最有可能产生现有的采样数据。每次迭代的过程就是参数矫正的过程。
　

这里写图片描述
图11 EM算法参数优化过程

 　　现假设初始化一组参数 ( μ 0 , Σ 0 , π 0 ) ({\mu ^0},{\Sigma ^0},{\pi ^0}) 。在这组参数下，2类二维高斯分布如图11绿色椭圆所示。而后利用现有的参数，E-step开始对样本数据进行划分（对 γ \gamma 进行估计）。蓝色的样本大多都被划分给第1类模型，橘黄色的样本大多都被划分给第2类模型。可是第1类模型还有优化空间：第1类模型还不能使得蓝色样本出现的联合几率达到最大。第2类模型也是如此。M-step便优化了2类模型的参数，获得新的参数 ( μ 1 , Σ 1 , π 1 ) ({\mu ^1},{\Sigma ^1},{\pi ^1}) ，使得优化后2类高斯分布如图11红色椭圆所示。其中，第1类模型主要优化的是模型均值（即椭圆的中心），第二类模型主要优化的是模型协方差矩阵（即椭圆的长轴、短轴和长短轴的方向）。而后重复进行E-step和M-step，直到参数 ( μ , Σ , π ) (\mu ,\Sigma ,\pi ) 收敛。
　　最后谈谈混合高斯模型的参数 π \pi 。
　　混合高斯模型的参数 μ , Σ \mu ,\Sigma 比较好理解，用于描述各个高斯分布的形状，对于它们的调整也比较直观：使得本高斯分布可以更好地接纳被划分到这类分布的样本。而为何要有参数 π \pi ？它描述的是各个高斯分布所占的比重，若是不加“歧视”的话（样原本源于各个高斯分布的可能性一致），则有 π k = 1 / K {\pi _k} = 1/K ；而若是对于某一类高斯分布（即为i）有侧重的话，则相应的 π i {\pi _i} 较大，体如今图11中就是被分配给各个类的样本数占样本总数的比例。若是一轮优化后，某一类高斯分布又接纳了更多样本，则其 π i {\pi _i} 变大，反之变小（因此图11从绿色椭圆调整为红色椭圆实际上两个类所对应的权重也被优化了）。
　　而从本质上来看参数 π \pi ，则是为了混合高斯模型能有更好的曲面拟合能力。当参数 π \pi 退化为某一类高斯分布的权重远远大于其余类高斯分布的时候，混合高斯模型就退化成了单高斯模型！

5 总结

 　　图12和图13梳理了高斯分布和混合高斯分布参数估计的逻辑流程。
　

这里写图片描述
图12 高斯分布参数估计逻辑流程

这里写图片描述
图13 混合高斯分布参数估计逻辑流程

 　　相对比于高斯分布的参数估计，混合高斯分布的参数估计更加复杂。主要缘由在于隐变量的存在。而为何混合高斯分布的参数估计须要屡次迭代循环进行？是由于EM算法中对于 γ \gamma 的估计利用的是初始化或者第i步迭代的参数 ( μ i , Σ i , π i ) ({\mu ^i},{\Sigma ^i},{\pi ^i}) ，这对于样本的分类划分是有偏差的。因此它只能经过屡次迭代优化寻找更佳的参数来抵消这一偏差。 　　终于把这篇文章梳理完了。世界杯要结束了，伪球迷也想见证一下冠军诞生。至此，本文结束。

联系我们 沪ICP备13005482号-10