【机器学习第7章——贝叶斯分类器】

7.贝叶斯分类器

7.1贝叶斯决策论

• 一个医疗判断问题

  • 有两个可选的假设：病人有癌症；病人无癌症

  • 可用数据来自化验结果：正+和负-

  • 有先验知识

    • 在所有人口中，患病率是0.008
    • 对确实有病的患者的化验准确率是0.98
    • 对确实无病的患者的化验准确率是0.97、

  • 总结如下
    $$\begin{gathered} P(cancer)=0.008,p(\neg cancer)0.992 \\ P(+|cancer)=0.98,p(-|cancer)0.02 \\ P(+|\neg cancer)=0.08,p(-|cancer)0.97 \end{gathered}$$

  • 问题：假定有一个新病人，化验结果为正，是否应将病人断定为有癌症？求后验概率
    $$P(cancer|+),p(\neg cancer|+)$$

• 贝叶斯定理的意义在于，我们在生活中经常遇到这种情况：

  • 我们可以很容易直接得出P(A|B)，P(B|A)则很难直接得出，但我们更关心P(B|A)。

• 贝叶斯定理
  $$P(B|A)=\frac{P(A|B)P(B)}{P(A)}$$

• 例子：有一家企业M提供产品和服务，企业K考虑是否进入该市场。同时企业M会阻止K进入该市场，而企业K能否进入该市场完全取决于M为组织其进入所花费的成本大小。

  • 假设K并不知道M是属于高阻挠成本类型还是低阻挠成本类型
    • 如果M属于高阻挠成本类型，K进入市场时，M进行阻挠的概率是20%
    • 如果M属于低阻挠成本类型，K进入市场时，M进行阻挠的概率是100%
  • 现设K认为M属于高阻挠成本企业的概率为70%，K进入市场后，M进行了阻挠，判断企业M为高阻挠成本类型的概率。

• 综上，M为高阻挠成本类型的概率为0.32，也就是说M为低阻挠成本类型的概率为0.68，即M更有可能是低阻挠成本类型

7.2 朴素贝叶斯分类器

给定类标号y，朴素贝叶斯分类器在估计类条件概率时假设属性之间条件独立。条件独立假设可以形式化的表达如下:
$$\begin{gathered} P(X|Y=y)=\prod _{i=1}^{n}P(x_i|Y=y) \\ 其中每个训练样本可用一个属性向量X=(x_1,x_2,x_3,...,x_n)表示，各个属性之间条件独立 \end{gathered}$$

• 假设给定的训练样本数据，计算
  $$x=\left(\begin{array}{c}outlook=Sunny\\ temperature=Cool\\ humdity=High\\ wind=Strong\end{array}\right),y=[yes,no]$$

Day	Outlook	Temperate	Humidity	Wind	Play Tennis
D1	Sunny	Hot	High	Weak	No
D2	Sunny	Hot	High	Strong	No
D3	Overcast	Hot	High	Weak	Yes
D4	Rain	Mild	High	Weak	Yes
D5	Rain	Cool	Normal	Weak	Yes
D6	Rain	Cool	Normal	Strong	No
D7	Overcast	Cool	Normal	Strong	Yes
D8	Sunny	Mild	High	Weak	No
D9	Sunny	Cool	Normal	Weak	Yes
D10	Rain	Mild	Normal	Weak	Yes
D11	Sunny	Mild	Normal	Strong	Yes
D12	Overcast	Mild	High	Strong	Yes
D13	Overcast	Hot	Normal	Weak	Yes
D14	Rain	Mild	High	Strong	No

• 对于如果存在连续属性来说，假定概率服从高斯分布，则
  $${\mu }_{c,i}和{\sigma }_{c,i}^2表示第c类样本在第i个属性上取值的均值和方差$$

  $$p(x_i|c)=\frac{1}{\sqrt{2\pi }{\sigma }_{c,i}}{e}^{-\frac{(x_i-{\mu }_{c,i}{)}^2}{2{\sigma }_{c,i}^2}}$$

  • 假如有一个西瓜的密度为0.697，样本中好瓜密度的均值为0.574，好瓜标准差为0.129
    $$p(密度=0.697|好瓜=是)=\frac{1}{\sqrt{2\pi }\times 0.129}{e}^{-\frac{(0.697-0.574{)}^2}{2\times 0.129^2}}\approx 1.959$$
    剩下的属性计算步骤与上例相同

• 但是这样存在一个问题，假设来了一个新样本
  $$x=\left(\begin{array}{c}outlook=Cloudy\\ temperature=Cool\\ humdity=High\\ wind=Strong\end{array}\right)$$
  进行计算会发现
  $$\begin{gathered} P(outlook=Cloudy|Y=yes)=\frac{0}{9}=0 \\ P(outlook=Cloudy|Y=no)=\frac{0}{5}=0 \end{gathered}$$
  这样到最后会使整个式子等于0，没有记录到并不等于出现的概率是0，所以提出了条件概率的m估计

条件概率的m估计

当训练样本不能覆盖那么多的属性值时，都会出现上述的窘境。简单的使用样本比例来估计类条件概率的方法太脆弱了，尤其是当训练样本少而属性数目又很大时。解决方法是使用m估计方法来估计条件概率
$$P(X_i|Y)=\frac{n_c+mp}{n+m}$$

$$\begin{gathered} n是Y中的样本总数 \\ {n}_c是Y中取值x_i的样本数 \\ m是称为等价样本大小的参数 \\ p是用户指定的参数 \end{gathered}$$

也就是说，在分子和分母上同时加上一点儿数，使得分数不要为0

• 多项式模型

  • 基本原理：在多项式模型，设某文档
    $$\begin{gathered} d=(t_1,t_2,...,t_k) \\ {t}_k是出现过的单词，允许重复 \end{gathered}$$

  • 先验概率
    $$p(c)=\frac{类c下单词总数}{整个训练本的单词总数}$$

  • 条件概率
    $$\begin{gathered} P(t_k|c)=\frac{类c下单词t_k在各个文档中出现的次数+1}{类c下单词总数+|v|} \\ v是训练样本的单词表（即抽取单词，单词出现多次，只算一个） \\ 此时m=|v|,p=\frac{1}{|v|} \end{gathered}$$

  • 给定一个新样本
    $$doc:ChineseChineseChineseTokyoJapan$$
    用属性向量表示为
    $$d=(Chinese,Chinese,Chinese,Tokyo,Japan)$$

id	doc	类别 in China？
1	Chinese Beijing Chinese	yes
2	Chinese Chinese Shanghai	yes
3	Chinese Macao	yes
4	Tokyo Japan Chinese	no

• 伯努利模型

  • 先验概率
    $$p(c)=\frac{类c下文件数}{整个训练本的文件总数}$$

  • 条件概率
    $$\begin{gathered} P(t_k|c)=\frac{类c下单词t_k的文件数+1}{类c下文件数+2} \\ 此时m=2,p=\frac{1}{2} \end{gathered}$$

    $$p(t_c|c=yes)=\prod p(t_c|c=yes)(1-p(t_i|c=yes))$$

  • 给定一个新样本
    $$doc:ChineseChineseChineseTokyoJapan$$
    用属性向量表示为
    $$d=(Chinese,Chinese,Chinese,Tokyo,Japan)$$

id	doc	类别 in China？
1	Chinese Beijing Chinese	yes
2	Chinese Chinese Shanghai	yes
3	Chinese Macao	yes
4	Tokyo Japan Chinese	no

• 模型比较
  • 二者计算粒度不同，多项式以单词为粒度，而伯努利以文件为粒度
  • 在后验概率计算中，多项式只会计算文档d出现过的单词的概率，而伯努利则会把所有单词的概率都算进去（有出现在文档d中的概率为p，没出现的为1-p）
• 运用贝叶斯定理所需概率
  • P(B)（先验）
  • P(A|B)（类条件密度）
• 但是大部分情况我们并不能获得到如此完善的信息
• 如果可以将类条件密度参数化，则可以显著降低难度
  • 例如：P(A|B)的正态性，就将概率密度估计问题转化为参数估计问题
• 所以常见的估计有两种
  • 极大似然估计（MLE）
    • 将参数看作是确定的量，只是其值是未知，通过极大化所观察的样本概率得到最优的参数，接着通过分析得到参数的值
  • 贝叶斯估计
    • 把参数当成服从某种先验概率分布的随机变量，对样本进行观测的过程，就是把先验概率密度转化成为后验概率密度，使得对于每个新样本，后验概率密度函数在待估参数的真实值附近形成最大尖峰。

7.3 极大似然估计

优点

• 当样本数目增加时，收敛性质会更好
• 比其他可选择的技术更加简单

基本原理

$$D_c$$

表示训练集D中第c类样本组合的集合，样本之间是独立同分布的，则参数
$${\theta }_c$$
对于数据集的似然是
$$P(D_c|{\theta }_c)=\prod _{x\in D_c}P(x|{\theta }_c)$$
但这样会造成数据的下溢，需要对其取对数
$$LL({\theta }_c)=\log P(D_c|{\theta }_c)$$
此时我们对对数似然函数的
$${\theta }_c$$
求导，并令导数为0，可求解出
$${\hat{\theta }}_c$$
即为使对数似然最大的值，记为
$${\hat{\theta }}_c=arg\underset{{\theta }_c}{\max }LL({\theta }_c)$$

• 以高斯情况为例，x为标量的似然函数计算如下

• 机器学习中特征向量是多维的，这时候密度函数需要发生变化

$$\begin{gathered} X=[x_1,x_2,..,x_n{]}^T,\mu =[{\mu }_1,{\mu }_2,...,{\mu }_n{]}^T \\ 多维情况下将方差替换成协方差矩阵(\sum ) \end{gathered}$$

$$f(X)=\frac{1}{(\sqrt{2\pi }{)}^n|\sum {|}^{\frac{1}{2}}}{e}^{-\frac{(X-\mu {)}^T(X-\mu )}{2\sum }}$$

$$\begin{gathered} \ln f(X)=-\frac{1}{2}\ln [(2\pi {)}^n|\sum |]-\frac{1}{2}(X-\mu {)}^T\sum ^{-1}(X-\mu ) \\ \sum ^{-1}:协方差矩阵的逆 \end{gathered}$$

7.4 贝叶斯网络

• 贝叶斯网络是一个有向无环图，包含两个部分

  • 节点：表示随机变量
  • 节点间的有向边：表示节点间的互相关系

• 贝叶斯网络B（D，P），D表示一个有向无环图，P表示条件概率分布的集合

  • 联合概率分布定义
    $$\begin{gathered} P(X)=\prod _{i=1}^{n}P(X_i|{\pi }_i) \\ {X}_i=(X_1,X_2,...,X_n) \\ {\pi }_i:X_i的父节点 \end{gathered}$$

7.5 半朴素贝叶斯分类器

朴素贝叶斯分类器采用的是属性条件独立性假设，但是在现实任务中这个假设往往很难成立，所以就产生了半朴素贝叶斯分类器

基本思想：适当考虑一部分属性间的相互依赖信息，从而既不需要进行完全联合概率计算，又不至于彻底忽略比较强的属性依赖关系，“独依赖估计”（ODE）是最常用的一种策略，即每个属性在类别之外最多依赖于一个其他属性。以下又提供了两种半朴素贝叶斯分类器所考虑的属性依赖关系。

例如，样本集合如下

我们规定：脐部的依赖属性为敲声，且属性取值为凹陷时，依赖敲声为浊响

并且需要对概率进行修正

7.6 EM算法

• 极大似然估计存在着问题是:

  • 对于许多具体问题不能构造似然函数解析表达式
  • 似然函数的表达式过于复杂而导致求解方程组非常困难

• 所以提出了EM算法。EM算法主要用于非完全数据参数估计，它是通过假设隐变量的存在，极大化地简化了似然函数方程，从而解决了方程求解问题。

• 观测数据:观测到的随机变量y的独立同分布(IID)样本
  $$Y=(Y_1,Y_2,\dots ,Y_n)$$

• 缺失数据:未观测到的随机变量Z的值
  $$Z=(Z_1,Z_2.\dots ,Z_n)$$

• 完整数据:包含观测到的随机变量Y和未观测到的随机变量Z的数据,Z=(X,Y)
  $$Z=((X_1,Y_1),..,(X_n,Y_n))$$

• EM算法是一种迭代算法，用于含有隐变量的概率模型参数的极大似然估计，或极大后验概率估计。

• EM算法的每次迭代由两步组成

  • E步，求期望
  • M步，求极大

• 三硬币模型

  • 假设有3枚硬币，分别记作A,B,C。这些硬币正面出现的概率概率分别是a,b,c。进行如下掷硬币试验:先掷A，根据其结果选出硬币B或C，正面选择硬币B，反面选硬币C;然后掷选出的硬币，掷硬币的结果，出现正面记作1，出现反面记作0;独立地重复n次试验(这里n=10),观测结果如下:
    $$1,1,0,1,0,0,1,0,1,1$$
    假设只能观测到掷硬币的结果，不能观测到掷硬币的过程。问如何估计三硬币正面出现的概率，即三硬币模型的参数。

  • 模型可写作

    

    $$\begin{gathered} 随机变量y是观测变量，表示一次试验观测的结果是1或0 \\ 随机变量z是隐变量，表示未观测到的掷硬币A的结果 \\ \theta =(a,b,c)是模型参数 \end{gathered}$$

  • 将观测数据表示为
    $$Y=(y_1,y_2,...,y_n{)}^T$$
    未观测数据表示为
    $$Z=(z_1,z_2,...,z_n{)}^T$$
    则观测数据的似然函数为
    $$P(Y|\theta )=\sum _{Z}P(Z|\theta )P(Y|Z,\theta )$$
    即
    $$P(Y|\theta )=\prod _{j=1}^n[a{b}^{y_j}(1-b{)}^{1-y_j}+(1-a)c^{y_j}(1-c{)}^{1-y_j}]$$
    考虑求模型参数的极大似然估计
    $$\hat{\theta }=arg\underset{\theta }{\max }\log P(Y|\theta )$$
    只能通过迭代的方式进行求解

  • EM算法首先选取参数的初值
    $${\theta }^{(0)}=(a^{(0)},b^{(0)},c^{(0)})$$
    然后通过以下步骤迭代计算参数的估计值，直至收敛为止。第i次迭代参数的估计值为
    $${\theta }^{(i)}=(a^{(i)},b^{(i)},c^{(i)})$$
    第i+1次迭代如下

    • E步：计算在模型参数
      $$a^{(i)},b^{(i)},c^{(i)}$$

      $$上标为迭代次数$$

      下观测数据
      $$y_j$$

      $$下标为样本编号$$

      来自掷到硬币B的概率
      $$P(z_j=B|y_j,\theta )=\frac{a^{(i)}(b^{(i)}{)}^{y_j}(1-b^{(i)}{)}^{1-y_j}}{a^{(i)}(b^{(i)}{)}^{y_j}(1-b^{(i)}{)}^{1-y_j}+(1-a^{(i)})(c^{(i)}{)}^{y_j}(1-c^{(i)}{)}^{1-y_j}}(1)$$
      并把
      $$\begin{gathered} P(z_j=B|y_j,\theta )记作{\mu }^{(j)} \\ 1-{\mu }^{(j)}即表示为掷到硬币C的概率 \end{gathered}$$
      计算似然函数的期望
      Q ( θ ) = P ( Z ∣ Y , θ ) log ⁡ P ( Y , Z ∣ θ ) = ∑ i = 1 n [ μ i log ⁡ a b y i ( 1 − b ) 1 − y i + ( 1 − μ i ) log ⁡ ( 1 − a ) c y i ( 1 − c ) 1 − y i ] = ∑ i = 1 n { μ i ( log ⁡ a + y i log ⁡ b + ( 1 − y i ) log ⁡ ( 1 − b ) ) + ( 1 − μ i ) ( log ⁡ ( 1 − a ) + y i log ⁡ c + ( 1 − y i ) log ⁡ ( 1 − c ) ) } ( 2 ) Q(\theta)=P(Z|Y,\theta)\log P(Y,Z|\theta) \\=\sum_{i=1}^n[\mu_i\log ab^{y_i}(1-b)^{1-y_i}+(1-\mu_i)\log(1-a)c^{y_i}(1-c)^{1-y_i}]\\ =\sum_{i=1}^n\{\mu_i(\log a +y_i\log b+(1-y_i)\log(1-b))+\\(1-\mu_i)(\log(1-a) +y_i\log c+(1-y_i)\log(1-c))\}\quad(2) Q(θ)=P(Z∣Y,θ)logP(Y,Z∣θ)=i=1∑n​[μi​logabyi​(1−b)1−yi​+(1−μi​)log(1−a)cyi​(1−c)1−yi​]=i=1∑n​{μi​(loga+yi​logb+(1−yi​)log(1−b))+(1−μi​)(log(1−a)+yi​logc+(1−yi​)log(1−c))}(2)

    • M步：求似然函数的极大值

      • 对a求偏导
        $$\begin{gathered} \frac{\partial Q}{\partial a}=\sum _{i=1}^n(\frac{{\mu }_i}{a}-\frac{1-{\mu }_i}{1-a})=\frac{1}{a(1-a)}\sum _{i=1}^n({\mu }_i-a)=0 \\ \frac{1}{a(1-a)}(\sum _{i=1}^n{\mu }_i-na)=0 \\ a=\frac{1}{n}\sum _{i=1}^n{\mu }_i \end{gathered}$$

      • 对b求偏导
        $$\begin{gathered} \frac{\partial Q}{\partial b}=\sum _{i=1}^n{\mu }_i(\frac{y_i}{b}-\frac{1-y_i}{1-b})=\frac{1}{b(1-b)}\sum _{i=1}^n{\mu }_i(y_i-b)=0 \\ \sum _{i=1}^n{\mu }_i{y}_i-b\sum _{i=1}^n{\mu }_i=0 \\ b=\frac{{\sum }_{i=1}^n{\mu }_i{y}_i}{{\sum }_{i=1}^n{\mu }_i} \end{gathered}$$

      • 对c求偏导
        $$\begin{gathered} \frac{\partial Q}{\partial c}=\sum _{i=1}^n(1-{\mu }_i)(\frac{y_i}{c}-\frac{1-y_i}{1-c})=\frac{1}{c(1-c)}\sum _{i=1}^n(1-{\mu }_i)(y_i-c)=0 \\ \sum _{i=1}^n(1-{\mu }_i)y_i-c\sum _{i=1}^n(1-{\mu }_i)=0 \\ c=\frac{{\sum }_{i=1}^n(1-{\mu }_i)y_i}{{\sum }_{i=1}^n(1-{\mu }_i)} \end{gathered}$$

    • 计算完成后，可以计算模型参数的新估计值
      $$a^{(j+1)}=\frac{1}{n}\sum _{i=1}^n{\mu }_i^{(j+1)}(3)$$

      $$b^{(j+1)}=\frac{{\sum }_{i=1}^n{\mu }_i^{(j+1)}{y}_i}{{\sum }_{i=1}^n{\mu }_i^{(j+1)}}(4)$$

      $$c^{(j+1)}=\frac{{\sum }_{i=1}^n(1-{\mu }_i^{(j+1)})y_i}{{\sum }_{i=1}^n(1-{\mu }_i^{(j+1)})}(5)$$

    • 然后进行计算，假设模型参数的初始值为
      $$a^{(0)}=0.5,b^{(0)}=0.5,c^{(0)}=0.5$$
      由式（1），对十次试验进行计算
      $$回顾十次试验结果:1,1,0,1,0,0,1,0,1,1$$

      $${\mu }_1^{(1)}=\frac{0.5\times 0.5^1\times (1-0.5{)}^{1-1}}{0.5\times 0.5^1\times (1-0.5{)}^{1-1}+(1-0.5)\times 0.5^1\times (1-0.5{)}^{1-1}}=0.5$$

      $${\mu }_3^{(1)}=\frac{0.5\times 0.5^0\times (1-0.5{)}^{1-0}}{0.5\times 0.5^0\times (1-0.5{)}^{1-0}+(1-0.5)\times 0.5^0\times (1-0.5{)}^{1-0}}=0.5$$

      经过计算可知：
      $${\mu }_j^{(1)}=0.5,j=1,2,...,10$$

      • 利用（3）（4）（5）式进行第一次迭代

      $$a^{(1)}=\frac{1}{10}\sum _{i=1}^{10}{\mu }_i^{(1)}=0.5$$

      $$\begin{gathered} b^{(1)}=\frac{{\sum }_{i=1}^{10}{\mu }_i^{(1)}{y}_i}{{\sum }_{i=1}^{10}{\mu }_i^{(1)}} \\ =\frac{0.5\times (1+1+0+1+0+0+1+0+1+1)}{0.5\times 10} \\ =0.6 \end{gathered}$$

      $$\begin{gathered} c^{(1)}=\frac{{\sum }_{i=1}^{10}(1-{\mu }_i^{(1)})y_i}{{\sum }_{i=1}^{10}(1-{\mu }_i^{(1)})} \\ =\frac{(1-0.5)\times (1+1+0+1+0+0+1+0+1+1)}{(1-0.5)\times 10} \\ =0.6 \end{gathered}$$

      经过一轮计算后，再由式（1）计算，可得：
      $$\begin{gathered} {\mu }_1^{(2)}=\frac{0.5\times (0.6{)}^1\times (1-0.6{)}^{1-1}}{0.5\times (0.6{)}^1\times (1-0.6{)}^{1-1}+(1-0.5)\times (0.6{)}^1\times (1-0.6{)}^{1-1}} \\ =0.5 \end{gathered}$$

      $$\begin{gathered} {\mu }_3^{(2)}=\frac{0.5\times (0.6{)}^0\times (1-0.6{)}^{1-0}}{0.5\times (0.6{)}^0\times (1-0.6{)}^{1-0}+(1-0.5)\times (0.6{)}^0\times (1-0.6{)}^{1-0}} \\ =0.5 \end{gathered}$$

      经过计算可知：
      $${\mu }_j^{(2)}=0.5,j=1,2,...,10$$

      • 进行第二次迭代

      $$a^{(2)}=\frac{1}{10}\sum _{i=1}^{10}{\mu }_i^{(2)}=0.5$$

      $$\begin{gathered} b^{(2)}=\frac{{\sum }_{i=1}^{10}{\mu }_i^{(2)}{y}_i}{{\sum }_{i=1}^{10}{\mu }_i^{(2)}} \\ =0.6 \end{gathered}$$

      $$\begin{gathered} c^{(2)}=\frac{{\sum }_{i=1}^{10}(1-{\mu }_i^{(2)})y_i}{{\sum }_{i=1}^{10}(1-{\mu }_i^{(2)})} \\ =0.6 \end{gathered}$$

    • 因为
      $$a^{(1)}=a^{(2)},b^{(1)}=b^{(2)},c^{(1)}=c^{(2)}$$
      即达到的平衡，所以
      $$\hat{a}=0.5,\hat{b}=0.6,\hat{c}=0.6$$

7.7 EM算法实现

目标：估算两个正态分布的均值(方差相同)，但现在1000条数据混在了一起

import math
import copy
import numpy as np
import matplotlib.pyplot as plt
isdebug=False

• 创建数据

#指定k个高斯分布参数，这里指定k=2，注意2个高斯分布具有同样均方差Sigma，均值分别为Mu1，Mu2
def int_data(Sigma,Mu1,Mu2,k,N):
    global X
    global Mu
    global Expectations
    X=np.zeros((1,N))
    Mu=np.random.random(2)#随机初始化均值
    Expectations=np.zeros((N,k))
    for i in range(0,N):
        if np.random.random(1)>0.5:
            X[0,i]=np.random.normal()*Sigma+Mu1#normal:从正态分布中抽取随机样本
        else:
            X[0,i]=np.random.normal()*Sigma+Mu2
#最终生成的X中的数据是两个分布混合的数据，我们并不知道数据是来自哪个分布

• 步骤1：E步

$$\begin{gathered} E[i][j]=\frac{P(x=x_i|\mu =\mu =j)}{P(x=x_i)} \\ =\frac{e^{-\frac{1}{2{\sigma }^2}(x_i-{\mu }_j{)}^2}}{{\sum }_{i=1}^k{e}^{-\frac{1}{2{\sigma }^2}(x_i-{\mu }_j{)}^2}} \\ =\frac{Number}{Denom} \end{gathered}$$

#EM算法:步骤1，计算E[zij]
def e_step(Sigma,k,N):
    global X
    global Mu
    global Expectations
    for i in range(0,N):
        Denom=0
        for j in range(0,k):
            Denom+=math.exp((-1/(2*(float(Sigma**2))))*(float(X[0,i]-Mu[j]))**2)
        for j in range(0,k):
            Numer=math.exp((-1/(2*(float(Sigma**2))))*(float(X[0,i]-Mu[j]))**2)
            Expectations[i,j]=Numer/Denom 

• 步骤2：M步

$$\begin{gathered} {\mu }_j=\frac{{\sum }_{i=1}^{n}Expectations[z_{ij}]x_j}{{\sum }_{i=1}^{n}Expectations[z_{ij}]} \\ =\frac{Numer}{Denom} \end{gathered}$$

#EM算法:步骤2，求最大化E[zij]的参数Mu
def m_step(k,N):
    global Expectations
    global X
    for j in range(0,k):
        Numer=0
        Denom=0
        for i in range(0,N):
            Numer+=Expectations[i,j]*X[0,i]
            Denom+=Expectations[i,j]
        Mu[j]=Numer/Denom

• 算法运行

#算法迭代iter_num次，或达到精度Epsilon停止迭代
def run(Sigma,Mu1,Mu2,k,N,iter_num,Epsilon):#Sigma方差,iter_num迭代次数,Mu1和Mu2是生成数据的均值
    #N样本数
    int_data(Sigma,Mu1,Mu2,k,N)
    print(u"初始<u1,u2>:",Mu)
    for i in range(iter_num):
        Old_Mu=copy.deepcopy(Mu)#深拷贝
        e_step(Sigma,k,N)
        m_step(k,N)
        print(i,Mu)
        if sum(abs(Mu-Old_Mu))<Epsilon:
            break

if __name__=='__main__':
    run(6,40,20,2,1000,100,0.0001)
    plt.hist(X[0,:],50)
    plt.show()

• 运行结果

初始<u1,u2>: [0.29065113 0.21786728]
0 [30.60448153 30.32328553]
1 [31.01065788 29.92567547]
2 [32.52251613 28.41379168]
3 [36.79745623 24.13331471]
4 [40.10662972 20.77667049]
5 [40.47241312 20.351192  ]
6 [40.48060384 20.31459203]
7 [40.47597691 20.30740231]
8 [40.4736702  20.30487205]
9 [40.47270634 20.30385998]
10 [40.47231156 20.30344781]
11 [40.47215027 20.30327955]
12 [40.47208441 20.30321084]
13 [40.47205751 20.30318278]