EM算法在高斯混合模型学习中的应用

电玩女神 2022-12-10 07:15 288阅读 0赞

> 本篇文章是之前[期望极大算法(EM算法)][EM]文章的后续，有需要可以先看看那篇文章关于EM算法的推导。

### 高斯混合模型 ###

高斯混合模型是研究算法的人避不开的一个东西，其在非深度学习的远古时代经常被用到，比如图像处理任务的前背景提取，点云处理任务的点云聚类等等等等。

具体的，高斯混合模型是指具有如下形式的概率分布模型：  
 P ( y ∣ θ ) = ∑ k = 1 K α k ϕ ( y ∣ θ k ) P(y \\mid \\theta)=\\sum\_\{k=1\}^\{K\} \\alpha\_\{k\} \\phi\\left(y \\mid \\theta\_\{k\}\\right) P(y∣θ)=k=1∑Kαkϕ(y∣θk)  
其中， α k \\alpha\_\{k\} αk是系数， α k ⩾ 0 \\alpha\_\{k\} \\geqslant 0 αk⩾0, ∑ k = 1 K α k = 1 \\quad \\sum\_\{k=1\}^\{K\} \\alpha\_\{k\}=1 ∑k=1Kαk=1； ϕ ( y ∣ θ k ) \\phi\\left(y \\mid \\theta\_\{k\}\\right) ϕ(y∣θk)是高斯分布密度， θ k = ( μ k , σ k 2 ) \\theta\_\{k\}=\\left(\\mu\_\{k\}, \\sigma\_\{k\}^\{2\}\\right) θk=(μk,σk2)，  
 ϕ ( y ∣ θ k ) = 1 2 π σ k exp ⁡ ( − ( y − μ k ) 2 2 σ k 2 ) \\phi\\left(y \\mid \\theta\_\{k\}\\right)=\\frac\{1\}\{\\sqrt\{2 \\pi\} \\sigma\_\{k\}\} \\exp \\left(-\\frac\{\\left(y-\\mu\_\{k\}\\right)^\{2\}\}\{2 \\sigma\_\{k\}^\{2\}\}\\right) ϕ(y∣θk)=2πσk1exp(−2σk2(y−μk)2)  
称为第 k k k个分模型。

###  Q Q Q函数的一般表达 ###

在算法处理的过程中，将问题建模成高斯混合模型后，往往需要去解模型中的参数，这个时候就需要用到EM算法。

在[期望极大算法(EM算法)][EM]文章的分析中我们已经知道，要进行EM算法，得先得到 Q Q Q函数：  
 Q ( θ , θ ( i ) ) = ∑ Z log ⁡ P ( Y , Z ∣ θ ) P ( Z ∣ Y , θ ( i ) ) Q\\left(\\theta, \\theta^\{(i)\}\\right)=\\sum\_\{Z\} \\log P(Y, Z \\mid \\theta)P\\left(Z \\mid Y, \\theta^\{(i)\}\\right) Q(θ,θ(i))=Z∑logP(Y,Z∣θ)P(Z∣Y,θ(i))  
在随后的抛硬币问题中我们也介绍了一种求解这个 Q Q Q函数的方法.但如果看过李老师的书的人会发现，本人没有用书中给出的看起来更具总结性的 Q Q Q函数形式来求解问题，原因在于当时本人也不明白那个式子是怎么来的。。。

在李航老师的书上， Q Q Q函数是这样定义的：**完全数据的对数似然函数 log ⁡ P ( Y , Z ∣ θ ) \\log P(Y, Z \\mid \\theta) logP(Y,Z∣θ)关于在给定观测数据 Y Y Y和当前参数 θ ( i ) \\theta^\{(i)\} θ(i)下对未观测数据 Z Z Z的条件概率分布 P ( Z ∣ Y , θ ( i ) ) P\\left(Z \\mid Y, \\theta^\{(i)\}\\right) P(Z∣Y,θ(i))的期望称为 Q Q Q函数，即：**  
 Q ( θ , θ ( i ) ) = E Z \[ log ⁡ P ( Y , Z ∣ θ ) ∣ Y , θ ( i ) \] Q\\left(\\theta, \\theta^\{(i)\}\\right)=E\_\{Z\}\\left\[\\log P(Y, Z \\mid \\theta) \\mid Y, \\theta^\{(i)\}\\right\] Q(θ,θ(i))=EZ\[logP(Y,Z∣θ)∣Y,θ(i)\]  
第一次看到这个定义和下面的公式感觉整个人都不好了！实在不知道他们之间为什么是个相等的关系。在经过一两个小时的发呆、无助、掉头发后突然看懂了！

假设 log ⁡ P ( Y , Z ∣ θ ) \\log P(Y, Z \\mid \\theta) logP(Y,Z∣θ)是只与变量 Z Z Z相关的函数，则可以把其写成 f ( Z ) f(Z) f(Z)，当 Z Z Z取得一个定值的时候，其就是一个固定的数值。如果这个时候对它取期望，就有：  
 E Z ( f ( Z ) ) = ∑ Z \[ f ( Z ) P ( Z ) \] E\_Z(f(Z))=\\sum\_\{Z\}\\left\[f(Z)P\\left(Z\\right)\\right\] EZ(f(Z))=Z∑\[f(Z)P(Z)\]  
而如果Z的取值本身受到别的参数 x x x， y y y影响，而这些参数都已经给出，则原式可以写成：  
 E ( f ( Z ) ∣ x , y ) = ∑ Z \[ f ( Z ) P ( Z ∣ x , y ) \] E(f(Z)\\mid x,y)=\\sum\_\{Z\}\\left\[f(Z)P\\left(Z\\mid x,y\\right )\\right\] E(f(Z)∣x,y)=Z∑\[f(Z)P(Z∣x,y)\]  
代入我们已知的量，等式得证。

由证明的过程也可以知道，这里的给定观测数据 Y Y Y和当前参数两个已知量影响的只有 Z Z Z这一变量。由于他们与 log ⁡ P ( Y , Z ∣ θ ) \\log P(Y, Z \\mid \\theta) logP(Y,Z∣θ)中的两个参数看起来似乎有关系，因此才大大增加了理解的难度。

有了 Q Q Q函数的一般表达，从式子的形式我们知道关键的一步是把高斯混合模型的 log ⁡ P ( Y , Z ∣ θ ) \\log P(Y, Z \\mid \\theta) logP(Y,Z∣θ)给列出来，也就是把其完全数据的似然函数的对数列出来。

### 高斯混合模型参数估计的EM算法 ###

假设数据 y 1 , y 2 , ⋯   , y N y\_\{1\},y\_\{2\}, \\cdots, y\_\{N\} y1,y2,⋯,yN由高斯混合模型生成，  
 P ( y ∣ θ ) = ∑ k = 1 K α k ϕ ( y ∣ θ k ) P(y \\mid \\theta)=\\sum\_\{k=1\}^\{K\} \\alpha\_\{k\} \\phi\\left(y \\mid \\theta\_\{k\}\\right) P(y∣θ)=k=1∑Kαkϕ(y∣θk)  
其中， θ = ( α 1 , α 2 , ⋯   , α K ; θ 1 , θ 2 , ⋯   , θ K ) \\theta=\\left(\\alpha\_\{1\}, \\alpha\_\{2\}, \\cdots, \\alpha\_\{K\} ; \\theta\_\{1\}, \\theta\_\{2\}, \\cdots, \\theta\_\{K\}\\right) θ=(α1,α2,⋯,αK;θ1,θ2,⋯,θK)，求高斯混合模型的参数，就是用EM算法估计高斯混合模型的参数 θ \\theta θ。

回顾观测数据 y j , j = 1 , 2 , ⋯   , N , y\_\{j\},j=1,2, \\cdots, N, yj,j=1,2,⋯,N, 的产生过程：

1.  依概率 α k \\alpha\_\{k\} αk选择第 k k k个高斯分布分模型 ϕ ( y ∣ θ k ) \\phi\\left(y \\mid \\theta\_\{k\}\\right) ϕ(y∣θk)；
2.  依第 k k k个分模型的概率分布 ϕ ( y ∣ θ k ) \\phi\\left(y \\mid \\theta\_\{k\}\\right) ϕ(y∣θk)生成观测数据 y j y\_\{j\} yj。

在高斯混合模型建模中，结果通常是已知的—观测数据 y j , j = 1 , 2 , ⋯   , N , y\_\{j\},j=1,2, \\cdots, N, yj,j=1,2,⋯,N, 是已知的；

而结果产生自哪个分模型未知—反映观测数据 y j y\_\{j\} yj来自第 k k k个分模型的数据未知， k = 1 , 2 , ⋯   , K ， k=1,2, \\cdots, K ， k=1,2,⋯,K，以隐变量 γ j k \\gamma\_\{j k\} γjk表示，可定义如下:  
 γ j k = \{ 1 , 第 j 个观测来自第 k 个分模型 0 , 否则 \\gamma\_\{j k\}=\\left\\\{\\begin\{array\}\{ll\} 1, & \\text \{ 第 \} j \\text \{ 个观测来自第 \} k \\text \{ 个分模型 \} \\\\ 0, & \\text \{ 否则 \} \\end\{array\}\\right. γjk=\{ 1,0,第j个观测来自第k个分模型否则

j = 1 , 2 , ⋯   , N ; k = 1 , 2 , ⋯   , K j=1,2, \\cdots, N ; \\quad k=1,2, \\cdots, K j=1,2,⋯,N;k=1,2,⋯,K

有了观测数据 y j y\_j yj及未观测数据 γ j k \\gamma\_\{j k\} γjk，那么完全数据是  
 ( y j , γ j 1 , γ j 2 , ⋯   , γ j K ) , j = 1 , 2 , ⋯   , N \\left(y\_\{j\}, \\gamma\_\{j 1\}, \\gamma\_\{j 2\}, \\cdots, \\gamma\_\{j K\}\\right), \\quad j=1,2, \\cdots, N (yj,γj1,γj2,⋯,γjK),j=1,2,⋯,N  
于是，可以写出完全数据的似然函数：  
 P ( y , γ ∣ θ ) = ∏ j = 1 N P ( y j , γ j 1 , γ j 2 , ⋯   , γ j K ∣ θ ) = ∏ k = 1 K ∏ j = 1 N \[ α k ϕ ( y j ∣ θ k ) \] γ j k = ∏ k = 1 K α k n ∏ j = 1 N \[ ϕ ( y j ∣ θ k ) \] γ j k = ∏ k = 1 K α k n k ∏ j = 1 N \[ 1 2 π σ k exp ⁡ ( − ( y j − μ k ) 2 2 σ k 2 ) \] γ j k \\begin\{array\}\{l\} P(y, \\gamma \\mid \\theta)=\\prod\_\{j=1\}^\{N\} P\\left(y\_\{j\}, \\gamma\_\{j 1\}, \\gamma\_\{j 2\}, \\cdots, \\gamma\_\{j K\} \\mid \\theta\\right)\\\\ =\\prod\_\{k=1\}^\{K\} \\prod\_\{j=1\}^\{N\}\\left\[\\alpha\_\{k\} \\phi\\left(y\_\{j\} \\mid \\theta\_\{k\}\\right)\\right\]^\{\\gamma\_\{jk\}\} \\\\ =\\prod\_\{k=1\}^\{K\} \\alpha\_\{k\}^\{n\} \\prod\_\{j=1\}^\{N\}\\left\[\\phi\\left(y\_\{j\} \\mid \\theta\_\{k\}\\right)\\right\]^\{\\gamma\_\{jk\}\} \\\\ =\\prod\_\{k=1\}^\{K\} \\alpha\_\{k\}^\{n\_\{k\}\} \\prod\_\{j=1\}^\{N\}\\left\[\\frac\{1\}\{\\sqrt\{2 \\pi\} \\sigma\_\{k\}\} \\exp \\left(-\\frac\{\\left(y\_\{j\}-\\mu\_\{k\}\\right)^\{2\}\}\{2 \\sigma\_\{k\}^\{2\}\}\\right)\\right\]^\{\\gamma\_\{jk\}\} \\end\{array\} P(y,γ∣θ)=∏j=1NP(yj,γj1,γj2,⋯,γjK∣θ)=∏k=1K∏j=1N\[αkϕ(yj∣θk)\]γjk=∏k=1Kαkn∏j=1N\[ϕ(yj∣θk)\]γjk=∏k=1Kαknk∏j=1N\[2πσk1exp(−2σk2(yj−μk)2)\]γjk  
式中， n k = ∑ j = 1 N γ j k n\_\{k\}=\\sum\_\{j=1\}^\{N\} \\gamma\_\{j k\} nk=∑j=1Nγjk， ∑ k = 1 K n k = N \\sum\_\{k=1\}^\{K\} n\_\{k\}=N ∑k=1Knk=N。

那么，完全数据的对数似然函数为：  
 log ⁡ P ( y , γ ∣ θ ) = ∑ k = 1 K n k log ⁡ α k + ∑ j = 1 N γ j k \[ log ⁡ ( 1 2 π σ k ) − ( y j − μ k ) 2 2 σ k 2 \] = ∑ k = 1 K n k log ⁡ α k + ∑ j = 1 N γ j k \[ log ⁡ ( 1 2 π ) − log ⁡ σ k − ( y j − μ k ) 2 2 σ k 2 \] \\log P(y, \\gamma \\mid \\theta)=\\sum\_\{k=1\}^\{K\} n\_\{k\} \\log \\alpha\_\{k\}+\\sum\_\{j=1\}^\{N\} \\gamma\_\{j k\}\\left\[\\log \\left(\\frac\{1\}\{\\sqrt\{2 \\pi\}\\sigma\_\{k\}\}\\right)-\\frac\{\\left(y\_\{j\}-\\mu\_\{k\}\\right)^\{2\}\}\{2 \\sigma\_\{k\}^\{2\}\}\\right\]\\\\ =\\sum\_\{k=1\}^\{K\} n\_\{k\} \\log \\alpha\_\{k\}+\\sum\_\{j=1\}^\{N\} \\gamma\_\{j k\}\\left\[\\log \\left(\\frac\{1\}\{\\sqrt\{2 \\pi\}\}\\right)-\\log\{\\sigma\_\{k\}-\}\\frac\{\\left(y\_\{j\}-\\mu\_\{k\}\\right)^\{2\}\}\{2 \\sigma\_\{k\}^\{2\}\}\\right\]\\\\ logP(y,γ∣θ)=k=1∑Knklogαk\+j=1∑Nγjk\[log(2πσk1)−2σk2(yj−μk)2\]=k=1∑Knklogαk\+j=1∑Nγjk\[log(2π1)−logσk−2σk2(yj−μk)2\]  
整个 Q Q Q函数为：  
 Q ( θ , θ ( i ) ) = E \[ log ⁡ P ( y , γ ∣ θ ) ∣ y , θ ( i ) \] = E \{ ∑ k = 1 K n k log ⁡ α k + ∑ j = 1 N γ j k \[ log ⁡ ( 1 2 π ) − log ⁡ σ k − ( y j − μ k ) 2 2 σ k 2 \] \} = ∑ k = 1 K \{ ∑ j = 1 N ( E γ j k ) log ⁡ α k + ∑ j = 1 N ( E γ j k ) \[ log ⁡ ( 1 2 π ) − log ⁡ σ k − 1 2 σ k 2 ( y j − μ k ) 2 \] \} \\begin\{aligned\} Q\\left(\\theta, \\theta^\{(i)\}\\right) &=E\\left\[\\log P(y, \\gamma \\mid \\theta) \\mid y, \\theta^\{(i)\}\\right\] \\\\ &=E\\left\\\{\\sum\_\{k=1\}^\{K\} n\_\{k\} \\log \\alpha\_\{k\}+\\sum\_\{j=1\}^\{N\} \\gamma\_\{j k\}\\left\[\\log \\left(\\frac\{1\}\{\\sqrt\{2 \\pi\}\}\\right)-\\log \\sigma\_\{k\}-\\frac\{\\left(y\_\{j\}-\\mu\_\{k\}\\right)^\{2\}\}\{2 \\sigma\_\{k\}^\{2\}\}\\right\]\\right\\\} \\\\ &=\\sum\_\{k=1\}^\{K\}\\left\\\{\\sum\_\{j=1\}^\{N\}\\left(E \\gamma\_\{j k\}\\right) \\log \\alpha\_\{k\}+\\sum\_\{j=1\}^\{N\}\\left(E \\gamma\_\{j k\}\\right)\\left\[\\log \\left(\\frac\{1\}\{\\sqrt\{2 \\pi\}\}\\right)-\\log \\sigma\_\{k\}-\\frac\{1\}\{2 \\sigma\_\{k\}^\{2\}\}\\left(y\_\{j\}-\\mu\_\{k\}\\right)^\{2\}\\right\]\\right\\\} \\end\{aligned\} Q(θ,θ(i))=E\[logP(y,γ∣θ)∣y,θ(i)\]=E\{ k=1∑Knklogαk\+j=1∑Nγjk\[log(2π1)−logσk−2σk2(yj−μk)2\]\}=k=1∑K\{ j=1∑N(Eγjk)logαk\+j=1∑N(Eγjk)\[log(2π1)−logσk−2σk21(yj−μk)2\]\}  
这里出现了 E γ j k E \\gamma\_\{j k\} Eγjk，依据我们前面的推导它是已知观测数据和当前参数的情况下，隐函数的似然。可以写成： E ( γ j k ∣ y , θ ) E\\left(\\gamma\_\{j k\} \\mid y, \\theta\\right) E(γjk∣y,θ)，记为 γ ^ j k \\hat\{\\gamma\}\_\{j k\} γ^jk。有  
 γ ^ j k = E ( γ j k ∣ y , θ ) = P ( γ j k = 1 ∣ y , θ ) = P ( γ j k = 1 , y j ∣ θ ) ∑ k = 1 K P ( γ j k = 1 , y j ∣ θ ) = P ( y j ∣ γ j k = 1 , θ ) P ( γ j k = 1 ∣ θ ) ∑ k = 1 K P ( y j ∣ γ j k = 1 , θ ) P ( γ j k = 1 ∣ θ ) = α k ϕ ( y j ∣ θ k ) ∑ k = 1 K α k ϕ ( y j ∣ θ k ) , j = 1 , 2 , ⋯   , N ; k = 1 , 2 , ⋯   , K \\begin\{aligned\} \\hat\{\\gamma\}\_\{j k\} &=E\\left(\\gamma\_\{j k\} \\mid y, \\theta\\right)=P\\left(\\gamma\_\{j k\}=1 \\mid y, \\theta\\right) \\\\ &=\\frac\{P\\left(\\gamma\_\{j k\}=1, y\_\{j\} \\mid \\theta\\right)\}\{\\sum\_\{k=1\}^\{K\} P\\left(\\gamma\_\{j k\}=1, y\_\{j\} \\mid \\theta\\right)\} \\\\ &=\\frac\{P\\left(y\_\{j\} \\mid \\gamma\_\{j k\}=1, \\theta\\right) P\\left(\\gamma\_\{j k\}=1 \\mid \\theta\\right)\}\{\\sum\_\{k=1\}^\{K\} P\\left(y\_\{j\} \\mid \\gamma\_\{j k\}=1, \\theta\\right) P\\left(\\gamma\_\{j k\}=1 \\mid \\theta\\right)\} \\\\ &=\\frac\{\\alpha\_\{k\} \\phi\\left(y\_\{j\} \\mid \\theta\_\{k\}\\right)\}\{\\sum\_\{k=1\}^\{K\} \\alpha\_\{k\} \\phi\\left(y\_\{j\} \\mid \\theta\_\{k\}\\right)\}, \\quad j=1,2, \\cdots, N ; \\quad k=1,2, \\cdots, K \\end\{aligned\} γ^jk=E(γjk∣y,θ)=P(γjk=1∣y,θ)=∑k=1KP(γjk=1,yj∣θ)P(γjk=1,yj∣θ)=∑k=1KP(yj∣γjk=1,θ)P(γjk=1∣θ)P(yj∣γjk=1,θ)P(γjk=1∣θ)=∑k=1Kαkϕ(yj∣θk)αkϕ(yj∣θk),j=1,2,⋯,N;k=1,2,⋯,K  
推到这可以知道 γ ^ j k \\hat\{\\gamma\}\_\{j k\} γ^jk等于当前模型参数下第 j j j个观测数据来自第 k k k个分模型的概率，称为分模型 k k k对观测数据 y j y\_j yj的响应度。代入 Q Q Q函数可以得到:  
 Q ( θ , θ ( i ) ) = ∑ k = 1 K \{ n k log ⁡ α k + ∑ j = 1 N γ ^ j k \[ log ⁡ ( 1 2 π ) − log ⁡ σ k − 1 2 σ k 2 ( y j − μ k ) 2 \] \} Q\\left(\\theta, \\theta^\{(i)\}\\right) =\\sum\_\{k=1\}^\{K\}\\left\\\{n\_\{k\} \\log \\alpha\_\{k\}+\\sum\_\{j=1\}^\{N\}\\hat\{\\gamma\}\_\{j k\}\\left\[\\log \\left(\\frac\{1\}\{\\sqrt\{2 \\pi\}\}\\right)-\\log \\sigma\_\{k\}-\\frac\{1\}\{2 \\sigma\_\{k\}^\{2\}\}\\left(y\_\{j\}-\\mu\_\{k\}\\right)^\{2\}\\right\]\\right\\\} Q(θ,θ(i))=k=1∑K\{ nklogαk\+j=1∑Nγ^jk\[log(2π1)−logσk−2σk21(yj−μk)2\]\}  
到此就得到了只含有模型参数的 Q Q Q函数，真的不容易啊！要是没有李航老师的书做参考，估计推到吐血都推不出来！

有了 Q Q Q函数相当于 E E E步就有了， M M M步很简单，就是 Q Q Q函数取相应模型参数的偏导然后求其极值点也就是等于0的点即可。求得结果如下：  
 μ ^ k = ∑ j = 1 N γ ^ j k y j ∑ j = 1 N γ ^ j k , k = 1 , 2 , ⋯   , K \\hat\{\\mu\}\_\{k\}=\\frac\{\\sum\_\{j=1\}^\{N\} \\hat\{\\gamma\}\_\{j k\} y\_\{j\}\}\{\\sum\_\{j=1\}^\{N\} \\hat\{\\gamma\}\_\{j k\}\}, \\quad k=1,2, \\cdots, K μ^k=∑j=1Nγ^jk∑j=1Nγ^jkyj,k=1,2,⋯,K

σ ^ k 2 = ∑ j = 1 N γ ^ j k ( y j − μ k ) 2 ∑ j = 1 N γ ^ j k , k = 1 , 2 , ⋯   , K \\hat\{\\sigma\}\_\{k\}^\{2\}=\\frac\{\\sum\_\{j=1\}^\{N\} \\hat\{\\gamma\}\_\{j k\}\\left(y\_\{j\}-\\mu\_\{k\}\\right)^\{2\}\}\{\\sum\_\{j=1\}^\{N\} \\hat\{\\gamma\}\_\{j k\}\}, \\quad k=1,2, \\cdots, K σ^k2=∑j=1Nγ^jk∑j=1Nγ^jk(yj−μk)2,k=1,2,⋯,K

α ^ k = n k N = ∑ j = 1 N γ ^ j k N , k = 1 , 2 , ⋯   , K \\hat\{\\alpha\}\_\{k\}=\\frac\{n\_\{k\}\}\{N\}=\\frac\{\\sum\_\{j=1\}^\{N\} \\hat\{\\gamma\}\_\{j k\}\}\{N\}, \\quad k=1,2, \\cdots, K α^k=Nnk=N∑j=1Nγ^jk,k=1,2,⋯,K

可以看到其只包含有 γ ^ j k \\hat\{\\gamma\}\_\{j k\} γ^jk这一与当前参数相关的变量，因此在实际计算过程中，我们只需要设定初值，然后在 E E E步计算出 γ ^ j k \\hat\{\\gamma\}\_\{j k\} γ^jk，在M步将计算得到的 γ ^ j k \\hat\{\\gamma\}\_\{j k\} γ^jk代入求得新的参数，一直重复直到收敛即可。

真的是推导要老命，编程3分钟啊！

### 高斯混合模型使用EM算法解决实际问题 ###

> 已知观测数据 -67，-48，6，8，14，16，23，24，28，29，41，49，56，60，75 试估计两个分量的高斯混合模型的5个参数。

由上面的推导已经知道了所有需要的信息，因此只要设定初值然后直接代公式即可，代码如下：

#include <iostream>
    #include <cmath>
    
    #define N 15
    #define pi 3.1415926535898
    
    class theta
    {
        
    public:
    	double alpha;
    	double mu;
    	double sigma;
    	void print()
    	{
        
    		std::cout << "--------------" << std::endl;
    		std::cout << "alpha:" << alpha << std::endl;
    		std::cout << "mu:" << mu << std::endl;
    		std::cout << "sigma:" << sigma << std::endl;
    		std::cout << "sigma平方" << sigma* sigma << std::endl;
    	}
    };
    
    theta mtheta[2];//两个高斯分模型参数
    double gamma[2][N];//E步结果gamma
    double y[N] = {
         -67,-48,6,8,14,16,23,24,28,29,41,49,56,60,75 };//观测结果
    
    double phi(theta& mtheta, double yj)
    {
        
    	return 1 / (sqrt(2 * pi) * mtheta.sigma) * exp(-pow((yj - mtheta.mu), 2) / (2 * pow(mtheta.sigma, 2)));
    }
    
    void EStep()
    {
        
    	for (size_t j = 0; j < N; j++)
    	{
        
    		gamma[0][j] = mtheta[0].alpha * phi(mtheta[0], y[j]);
    		gamma[1][j] = mtheta[1].alpha * phi(mtheta[1], y[j]);
    		double sum = gamma[0][j] + gamma[1][j];
    		gamma[0][j] = gamma[0][j] / sum;
    		gamma[1][j] = gamma[1][j] / sum;
    		//std::cout << "gamma0:" << gamma[0][j] << std::endl;
    		//std::cout << "gamma1:" << gamma[1][j] << std::endl;
    	}
    }
    
    void MStep()
    {
        
    	for (size_t k = 0; k < 2; k++)
    	{
        
    		double mu = 0;
    		double sigma = 0;
    		double sumgamma = 0;
    		for (size_t j = 0; j < N; j++)
    		{
        
    			mu += gamma[k][j] * y[j];
    			sigma += gamma[k][j] * pow((y[j] - mtheta[k].mu), 2);
    			sumgamma += gamma[k][j];
    		}
    		mtheta[k].mu = mu / sumgamma;
    		mtheta[k].sigma = sqrt(sigma / sumgamma);
    		mtheta[k].alpha = sumgamma / N;
    	}
    }
    
    int main()
    {
        
    	//初始化高斯分模型参数变量
    	mtheta[0].alpha = 0.5;
    	mtheta[0].mu = 30;
    	mtheta[0].sigma = sqrt(500);
    	mtheta[1].alpha = 0.5;
    	mtheta[1].mu = -30;
    	mtheta[1].sigma = sqrt(500);
    
    	for (size_t k = 0; k < 2; k++)
    	{
        
    		mtheta[k].print();
    	}
    	//迭代10次
    	for (size_t i = 0; i < 10; i++)
    	{
        
    		EStep();
    		MStep();
    		for (size_t k = 0; k < 2; k++)
    		{
        
    			mtheta[k].print();
    		}
    	}
    }

如果求出了模型参数后想知道观察结果在这一参数下的分类，再求一次 γ ^ j k \\hat\{\\gamma\}\_\{j k\} γ^jk即可！

### 参考文章 ###

https://zhuanlan.zhihu.com/p/32508410

[EM]: https://blog.csdn.net/weixinhum/article/details/108395402