🛂 🏇🏽 🙌🏻 我们将EM算法分解成小块 🧔🏼 👩🏿‍🤝‍👨🏼 🏒

正如您可能已经猜到的那样，在本文中，我们将讨论设备EM算法。首先，这篇文章可能对那些慢慢加入数据sainists社区的人们感兴趣。本文中介绍的材料对于最近刚开始参加MIPT和Yandex的“机器学习和数据分析”专业的第三门课程“搜索数据结构”的人来说将更为有用。

从某种意义上说，本文中介绍的材料是上述课程的第一周培训的补充，即，它使您可以回答有关EM算法的操作原理的一些重要问题。为了更好地理解该材料，建议我们的受尊敬的读者能够使用矩阵执行操作（矩阵乘法，找到矩阵和逆矩阵的行列式），了解概率论和matstat的基础，当然，至少要对基本聚类算法有基本的了解，并了解什么聚类发生在机器学习中。虽然，当然，如果没有此知识，您可以阅读本文，但可以肯定，这很清楚:)

同样，根据旧的传统，本文将不包含深入的理论研究，而将充满简单易懂的示例。每个后续示例将比前一个示例更深入地说明EM算法的操作，最终使我们直接了解算法本身。对于每个示例，将编写代码。所有代码都是用python 2.7编写的，为此我深表歉意。原来，我现在使用的是此特定版本，但是切换到python 3之后，我将尝试更改本文中的代码。

让我们熟悉本文的概述： 1）让我们从总体上考虑EM算法的工作原理。

2）我们重复贝叶斯定理（公式）的要点：

$P_{x_i}(B_j) = \frac{P(B_j) \centerdot P_{B_j}(x_i)}{P(x_i)}$

3）考虑贝叶斯定理的第一个示例，其中公式的所有元素（或更确切地说是概率）

$P(B_j), P_{B_j}(x_i)$ 右边等号之后的

）是已知的。在这种情况下，我们只需要正确了解要替换的数字，然后执行基本操作即可。

4）考虑贝叶斯定理的第二个例子。为了解决该问题，我们将需要另外计算某个事件的概率密度（

$x_i$ ）受假设（

$B_j$ ）-

$\rho_{B_j}(x_i)$ 。为了进行计算，我们将获得随机变量的参数，即数学期望-

$\mu$ 和标准偏差-

$\sigma$ 。概率密度计算将根据以下公式进行：

$\rho_{B_j}(x_i) = \frac{1}{\sigma_j \centerdot \sqrt{2 \pi}} \centerdot \exp(- \frac{(x_i-\mu_j)^2}{2 \sigma_j^2})$

使用以上公式假定随机变量具有一维维度并且呈正态分布（换句话说，高斯或高斯-拉普拉斯分布是概率分布）。

5）考虑第三个示例，该示例与前一个示例的不同之处仅在于它考虑了随机变量的多维（在我们的二维版本中）正态分布的情况。在这种情况下，密度计算根据以下公式进行：

$\rho_{B_j}(x_i) = \frac{1}{(2\pi)^{k/2} \centerdot \mid \Sigma_j\mid ^{1/2}} \centerdot \exp(- \frac{1}{2} \centerdot (x_i-\mu_j)^T \Sigma_j^{-1} (x_i-\mu_j))$

6）最后，我们修改第三个示例，以清楚地说明EM算法的操作。

7）最后，我们将推论得出的算法质量与EM算法的质量进行比较，该算法内置于sklearn库（类sklearn.mixture.GaussianMixture）中。

如果您看到的是日语字符而不是本文概述中指定的公式，请不要害怕-根据贝叶斯定理（本文概述的第二段），我们将分析这样的示例，您可能至少会直观地理解这些示例。好吧，走吧！

一般的EM算法

在该课程的基础上，根据文章的写作，提出了EM算法作为聚类方法之一。换句话说，当我们事先不知道真正的答案时，这是一种没有老师的机器学习方法。在本文中，我们还将算法视为具有正态分布的随机变量的聚类方法之一。但是，应注意该算法具有更广泛的应用。 EM算法（英语期望最大化）是一种在具有隐藏变量的模型中查找似然函数估计值的通用方法，该方法可以通过混合分布来构建（近似）复杂概率分布。

聚类问题中的EM算法用作迭代算法，该算法在每次迭代时执行两个步骤：

E步骤。在第一个E步中，以某种方式（例如随机地），我们选择隐藏变量，在这种情况下，这是数学上的期望-

$\mu$ 和标准偏差-

$\sigma$ 。使用选择的变量，我们计算将每个对象分配给特定群集的概率。后续的E步骤使用M步骤中定义的隐藏变量。

M步。在M步，根据将每个对象分配给在E步获得的特定聚类的概率，重新计算隐藏变量

$\mu$ 和

$\sigma$

重复迭代，直到发生收敛为止。我们可以假设当隐藏变量的值没有显着变化时发生收敛，例如在给定常数内。在本文中考虑的最后一个示例中，我们将不会设置常量，因此不会确定每次迭代中隐藏变量的值的变化。让我们做得更简单-我们将算法限制为固定数量的步骤。

M步骤很容易理解，我们将在最后一个示例中看到。在本文中，我们将重点介绍E步骤。

E步基于贝叶斯定理，在此基础上我们确定将每个对象分配给一个或另一个聚类的概率。让我们记住这个定理是关于什么的。

回顾贝叶斯公式

$P_{x_i}(B_j) = \frac{P(B_j) \centerdot P_{B_j}(x_i)}{P(x_i)}$

哪里

$P_{x_i}(B_j)$ -事件发生的概率

$x_i$ 一个假设

$B_j$

$P(B_j)$ -假设的可能性

$B_j$

$P_{B_j}(x_i)$ -事件发生的概率

$x_i$ 服从假设

$B_j$

$P(x_i)$ -事件发生的概率

$x_i$ 假设每个假设均成立（由该事件的总概率公式计算）

假设该事件

$x_i$ 已经发生，满足假设的概率

$B_j$ 使用上述公式对

重估。我们可以这样说

$P(B_j)$ 是假设的先验概率（在检验之前），并且

$P_{x_i}(B_j)$ 事件发生后是否评估了相同假设的后验概率

$x_i$ 也就是说，鉴于该事件

$x_i$ 可靠地发生了。

贝叶斯定理的例子一

想象一下，我们收到的零件是在两台不同的机器上生产的。仓库收到的产品的以下特征是已知的：

总计，零件-10,000个（

$N = 10000$ ），其中第一台机器上生产的零件-6000个。（

$N1 = 6000$ ），第二个-4000个。（

$N2 = 4000$ ）

在第一台机器上制造的标准（即无缺陷）产品的比例为0.9，在第二台机器上制造的标准产品的份额为0.8。

我们从收到的零件中随机抽取了一个零件，结果证明是标准的。

我们需要确定以下几率：

1）在第一台机器上生产零件；

2）该项目是在第二台机器上生产的。

决定

1）事件

$x$ 在我们的示例中，提取标准产品。

2）现在，让我们确定假设

$B_j$ 。我们有两个部分，这意味着两个假设：

假设

$B_1$ ：在第1号机器上生产了意外取出的产品。

假设

$B_2$ ：在2号机器上生产了被意外取出的产品。

因此，提取在第一台机器上制造的产品的概率为

$P(B_1)$ 弥补

$N1/N = 6000/10000 = 0.6$ 。

提取第二台机器上制造的产品的概率为

$P(B_2)$ 弥补

$N2/N = 0.4$ 。此外，我们不在乎标准产品或有缺陷的产品，重要的是批次是多少。

3）从1号机器生产的产品（即第一个假设）中提取标准产品的概率对应于1号机器生产的标准产品的比例，为0.9，

$P_{B_1}(x_i) = 0.9$ 。

提取产品的概率是受假设2约束的标准产品-

$P_{B_2}(x_i) = 0.8$

4）确定从整套产品中提取标准产品的可能性-

$P(x_i)$ 。按照公式计算事件的总概率

$P(x_i) = P(B_1) \centerdot P_{B_1}(x_i) + P(B_2) \centerdot P_{B_2}(x_i) = 0.6 \centerdot 0.9 + 0.4 \centerdot 0.8 = 0.54 + 0.32 = 0.86$ 。在这里我们知道

$0.14$ -这是从到达仓库的所有零件中提取有缺陷的零件的概率，总计为1

$(0.86 + 0.14 = 1)$ 。

5）现在我们拥有所有数据以解决问题。

5.1）确定在机器上生产随机提取的标准零件的可能性

№

$№1$ ：

$P_{x_i}(B_1) = \frac{P(B_1) \centerdot P_{B_1}(x_i)}{P(x_i)} = \frac{0.6 \centerdot 0.9}{0.86} \approx 0.63$

5.2）我们确定在机器上生产随机提取的标准零件的可能性

№

$№2$ ：

$P_{x_i}(B_2) = \frac{P(B_2) \centerdot P_{B_2}(x_i)}{P(x_i)} = \frac{0.4 \centerdot 0.8}{0.86} \approx 0.37$

因此，我们重新评估了假设

$B_1$ 和

$B_2$ 。重新评估之后，假设也构成了完整的事件组：

$0.63 + 0.37 = 1$ 。

好了，现在是时候继续第二个例子了。

使用随机变量正态分布参数的贝叶斯定理的第二个例子：数学期望 $\mu$ 和标准偏差 $\sigma$

让我们稍微修改前面示例的条件。我们假设我们不知道在1号和2号机器上制造的标准产品的比例，而是假设，我们知道产品直径的平均尺寸和每台机器上产品直径的标准偏差。

1号机生产的零件直径为64毫米，标准偏差为4毫米。

2号机生产的零件直径为52毫米，标准偏差为2毫米。

一个重要条件是，整个产品集都由正态分布或高斯分布来描述。

否则，条件是相同的，我们将它们编写为：

$N1 = 6000, \mu_1 = 64, \sigma_1 = 4$

$N2 = 4000, \mu_2 = 52, \sigma_2 = 2$

我们假设在产品验收过程中发生的事件很小，因此所有产品混合在一起。

我们的任务是整理细节，并确定每个细节在1号机或2号机上生产的可能性。我们还将假定零件是在机器上生产的，其机率更高。

解决方案

首先，我们将分析解决方案算法。我们可以轻松计算出每个假设的概率，但是，我们已经从过去的例子中知道了它们的值：

$P(B_1)=0.6$ ，

$P(B_2)=0.4$ 。但是，我们如何计算每个可用假设分别扣押标准产品的概率？一切都很简单！我们将不考虑概率，相反，我们将为每个假设分别确定要除去的项目的概率密度及其直径值。

为此，我们使用众所周知的公式：

$\rho_{B_j}(x_i) = \frac{1}{\sigma_j \centerdot \sqrt{2 \pi}} \centerdot \exp(- \frac{(x_i-\mu_j)^2}{2 \sigma_j^2})$

哪里

$x_i$ -随机值（在我们的情况下，是产品直径的实际大小），

$\mu_j$ -随机变量的数学期望

$j$ 假设（在我们的案例中，是指根据

$j$ 机床）

$\sigma_j$ -随机变量的标准差

$j$ 假设（在我们的案例中，是指从

$j$ 机床）。

因此，我们在贝叶斯公式的分子中巧妙地将概率替换为概率密度，在分母中我们将进行类似的替换。

我们专门针对我们的情况调整公式。

具有直径尺寸的零件的可能性

$x_i$ 在1号机器上生产的产品由以下公式确定：

$P_{x_i}(B_1) = \frac{w_1 \centerdot \rho_{B_1}(x_i)}{w_1 \centerdot \rho_{B_1}(x_i) + w_2 \centerdot \rho_{B_1}(x_i)}$

具有直径尺寸的零件的可能性

$x_i$ 在2号机器上生产

$P_{x_i}(B_2) = \frac{w_2 \centerdot \rho_{B_2}(x_i)}{w_1 \centerdot \rho_{B_1}(x_i) + w_2 \centerdot \rho_{B_1}(x_i)}$

请注意，我们将假设概率的表示法替换为

$P(B_j)$ 在

$w_j$ 。这是由于先前指定的课程有相应的指定。此外，我们将仅使用这种表示法。

现在我们已100％完成并准备解决问题。

我们看一下代码

导入库

#  , 
import numpy as np
import random
import math
import matplotlib.pyplot as plt
%matplotlib inline
from sklearn.metrics import accuracy_score

第二个示例中使用的函数

#      ,    
#      : .,  . 
def gaus_func_01(mu,sigma,x):
    return math.e**(-(x-mu)**2/(2*sigma**2)) / (sigma*(2*math.pi)**0.5)

#        
def proba_x(X, w1, w2, mu_1, mu_2, sigma_1, sigma_2):
    for i in X:
        P1_x = gaus_func_01(mu_1,sigma_1,i)
        P2_x = gaus_func_01(mu_2,sigma_2,i)
        P_x = w1*P1_x + w2*P2_x
        P_x_1 = (w1*P1_x)/P_x
        P_x_2 = (w2*P2_x)/P_x
        proba_temp = []
        proba_temp.append(P_x_1)
        proba_temp.append(P_x_2)
        proba_X.append(proba_temp)
    return proba_X

#        
def pred_x(proba_X, limit_proba):
    pred_X = []
    for x in proba_X:
        if x[0] >= limit_proba:
            pred_X.append(1)
        else:
            pred_X.append(2)
    return np.array(pred_X)

#    
def graph_01(X, pred_X, mu_1, sigma_1, mu_2, sigma_2):
    true_pred = []
    false_pred_1 = []
    false_pred_2 = []
    for i in range(X.shape[0]):
        if pred_X[i] == y[i]:
            true_pred.append([X[i], -0.025])
        else:
            if y[i] == 1:
                false_pred_1.append([X[i], -0.0075])
            else:
                false_pred_2.append([X[i], -0.015])

    false_pred_1 = np.array(false_pred_1)            
    false_pred_2 = np.array(false_pred_2)
    true_pred = np.array(true_pred)

    x_theory = np.linspace(42, 85, 20000)
    y_theory_1 = []
    for x in x_theory:
        y_theory_1.append(gaus_func_01(mu_1,sigma_1,x))
    y_theory_2 = []
    for x in x_theory:
        y_theory_2.append(gaus_func_01(mu_2,sigma_2,x))

    plt.figure(figsize=(18, 8))    
    plt.plot(
        x_theory, y_theory_1, color = 'green', lw = 2, label = 'Theoretical probability density for machine 1')
    plt.plot(
        x_theory, y_theory_2, color = 'firebrick', lw = 2, label = 'Theoretical probability density for machine 2')
    plt.hist(
        X[:N1], bins = 'auto', color='#539caf', normed = True, alpha = 0.35, label = 'machine tool products 1')
    plt.hist(
        X[N1:N], bins = 'auto', color='sandybrown', normed = True, alpha = 0.75, label = 'machine tool products 2')
    plt.plot(mu_1, 0, 'o', markersize = 11, color = 'blue', label = 'Mu 1')
    plt.plot(mu_2, 0, 'o', markersize = 11, color = 'red', label = 'Mu 2')

    plt.plot([mu_1 - sigma_1, mu_1 - sigma_1], [0,0.85*np.max(y_theory_1)],
             ':', lw = 3, color = 'blue', alpha = 0.55, label = 'Mu1 - sigma1')
    plt.plot([mu_1 + sigma_1, mu_1 + sigma_1], [0,0.85*np.max(y_theory_1)],
             ':', lw = 3, color = 'blue', alpha = 0.55, label = 'Mu1 + sigma1')
    plt.plot([mu_2 - sigma_2, mu_2 - sigma_2], [0,0.85*np.max(y_theory_2)],
             ':', lw = 3, color = 'red', alpha = 0.55, label = 'Mu2 - sigma2')
    plt.plot([mu_2 + sigma_2, mu_2 + sigma_2], [0,0.85*np.max(y_theory_2)],
             ':', lw = 3, color = 'red', alpha = 0.55, label = 'Mu2 + sigma2')

    plt.plot([mu_1 - 2 * sigma_1, mu_1 - 2 * sigma_1], [0, 0.9*0.5 * np.max(y_theory_1)],
             ':', lw = 2.5, color = 'blue', alpha = 0.35, label = 'Mu1 - 2*sigma1')
    plt.plot([mu_1 + 2 * sigma_1, mu_1 + 2 * sigma_1], [0, 0.9*0.5 * np.max(y_theory_1)],
             ':', lw = 2.5, color = 'blue', alpha = 0.35, label = 'Mu1 + 2*sigma1')
    plt.plot([mu_2 - 2 * sigma_2, mu_2 - 2 * sigma_2], [0, 0.9*0.5 * np.max(y_theory_2)],
             ':', lw = 2.5, color = 'red', alpha = 0.35, label = 'Mu2 - 2*sigma2')
    plt.plot([mu_2 + 2 * sigma_2, mu_2 + 2 * sigma_2], [0, 0.9*0.5 * np.max(y_theory_2)],
             ':', lw = 2.5, color = 'red', alpha = 0.35, label = 'Mu2 + 2*sigma2')

    plt.plot(false_pred_1[:,0], false_pred_1[:,1], 'o', markersize = 2.5, color = 'blue', alpha = 0.2, label = 'errors1')
    plt.plot(false_pred_2[:,0], false_pred_2[:,1], 'o', markersize = 2.5, color = 'red', alpha = 0.3, label = 'errors2')
    plt.plot(true_pred[:,0], true_pred[:,1], 'o', markersize = 3, color = 'green', alpha = 0.2, label = 'right answers')

    plt.xlabel('Caliber')
    plt.ylabel('Probability density')
    plt.legend()
    plt.show()

指定参数

#    
#      №1
N1 = 6000
#      №2
N2 = 4000
#      
N = N1+N2

#    №1
mu_1 = 64.
#        №1
sigma_1 = 3.5

#    №2
mu_2 = 52
#        №2
sigma_2 = 2.

我们将计算

X = np.zeros((N))
np.random.seed(seed=42)
#    ,   №1
X[:N1] = np.random.normal(loc=mu_1, scale=sigma_1, size=N1)
#  ,   №2
X[N1:N] = np.random.normal(loc=mu_2, scale=sigma_2, size=N2)

#   
y = np.zeros((N))
y[:N1] = np.array((1))
y[N1:N] = np.array((2))

#     ,    №1
w1 = float(N1)/N
#     ,    №2
w2 = float(N2)/N

#           
proba_X = []
proba_X = proba_x(X, w1, w2, mu_1, mu_2, sigma_1, sigma_2)

#   ,   ,        
limit_proba = 0.5

#     
pred_X = []
pred_X = pred_x(proba_X, limit_proba)

#    
print '   :', round(accuracy_score(y, pred_X),3)
print

print ' №1'
graph_01(X, pred_X, mu_1, sigma_1, mu_2, sigma_2)

图表1

我们刚刚做了什么？我们以伪随机方式生成了由正态分布描述的10,000个值，其中6,000个值的平均值为64毫米，标准偏差为4毫米，有4,000个值的平均值为52毫米，标准偏差为2毫米。另外，根据上述算法，对于每个零件，确定了其在1号或2号机器上生产的可能性。之后，我们根据最有可能的假设选择关于在一台或另一台机器上生产产品的假设。最后，我们将算法的结果与真实答案进行了比较。

在输出中，我们得到了正确答案-0.986。总体而言，这非常好，因为我们在进行培训时并未使用真实答案。

让我们看一下图表。建议注意什么？查看未由算法正确识别的产品位于何处。

首先，我们仅在使用不同机器制造的产品具有相同直径的区域中看到算法错误。看起来很合逻辑。

其次，我们看到该算法仅在那些与对象的真实平均值相距最远且同时又非常接近错误中心的对象中才被误认为。

但最有趣的是，问题主要在超过近似相等的值后开始

$\mu +- 2\sigma$ ，特别是在交叉路口

$\mu_1 - 2\sigma_1$ 和

$\mu_2 + 2\sigma_2$

那些希望的人可以“扭曲”参数，并查看算法质量如何变化。这将有助于更好地吸收材料。

好吧，我们继续下一个例子

第三个例子。二维情况

这次我们将考虑这样一种情况，即我们不仅具有有关产品直径的数据，而且还具有例如每种产品的重量的数据。假设第一台机器生产的零件重量为12 g，标准偏差为1 g，第二台机器生产的产品重量为10 g，标准偏差为0.8 mm

，其余的，我们将使用上一个示例中的数据，将其记录下来。

$N1 = 6000, \mu_{1 1} = 64, \sigma_{1 1} = 4, \mu_{1 2} = 14, \sigma_{1 2} = 1$

$N2 = 4000, \mu_{2 1} = 52, \sigma_{2 1} = 2, \mu_{2 2} = 10, \sigma_{2 2} = 0.9$

情况是一样的-所有细节都混合在一起，我们需要对其进行分类，并确定每个零件在1号机和2号机上生产的可能性。

解决方案

通常，除了用于确定随机变量的概率密度的公式之外，对上一示例的解决方案方法不会更改。对于多维情况，使用以下公式：

$p_{B_j}(x_i) = \frac{1}{(2\pi)^{k/2} \centerdot \mid \Sigma_j\mid ^{1/2}} \centerdot \exp(- \frac{1}{2} \centerdot (x_i-\mu_j)^T \Sigma_j^{-1} (x_i-\mu_j))$

哪里

$\Sigma_j$ 是随机变量协方差的矩阵

$j$ 假设（在我们的案例中，是在以下条件下产生的产品的参数协方差矩阵

$j$ 机床）

$\mid \Sigma_j\mid$ 是随机变量协方差矩阵的行列式

$j$ 假设

$\mu_j$ -随机变量的数学期望

$j$ 假设（在我们的案例中，是根据

$j$ 机床）

$x_i$ -随机变量（在本例中为参数

$i$ 产品）

这一次，要理解代码，受尊敬的读者将需要有关矩阵运算的知识

第三个示例中使用的函数

#      ,    
#      : .,  . 
def gaus_func_02(k, m, x, mu, sigma):
    pj_xi = []
    for j in range(k):
        det_sigma_j = np.linalg.det(sigma[j])
        factor_1 = 1 / (((2 * math.pi)**(k/2)) * ((det_sigma_j)**0.5))
        factor_2 = []
        for i in x:
#    ,            
#          
#   .    ,      
#     
            factor_2.append(math.e**float(
            -0.5 * np.matrix(i - mu[j]) * np.matrix(np.linalg.inv(sigma[j])) * np.matrix(i - mu[j]).T))
        pj_xi.append(factor_1 * np.array(factor_2))
    return np.array(pj_xi)

#     ,       №1  №2
def proba_func_02(pjxi, w, k):
    #          
    P_X = []
    for j in range(k):
        P_X.append(w[j] * pjxi[j])
    P_X = np.sum(np.array(P_X), axis = 0)
    #    ,       №1  №2
    P_J_X = []
    for j in range(k):
        P_J_X.append(w[j] * pjxi[j] / P_X)
    return np.array(P_J_X)

#        
def pred_x_02(proba_X, limit_proba):
    pred_X = []
    for x in proba_X[0]:
        if x >= limit_proba:
            pred_X.append(1)
        else:
            pred_X.append(2)
    return np.array(pred_X)

#             
def graph_02_algorithm(pred_X, mu):
    #   
    pred_X = np.array(pred_X)

    #   ,         
    answers_1 = []
    answers_2 = []

    for i in range(pred_X.shape[0]):
        if pred_X[i] == 1:
            answers_1.append(X[i])
        else:
            answers_2.append(X[i])
    
    print ' "     "'
    plt.figure(figsize=(16, 6))  
    plt.plot(
        np.array(answers_1)[:,0], np.array(answers_1)[:,1], 'o', alpha = 0.7, color='sandybrown', label = 'Produced on machine #1')
    plt.plot(
        np.array(answers_2)[:,0], np.array(answers_2)[:,1], 'o', alpha = 0.45, color = 'darkblue', label = 'Produced on machine #2')
    plt.plot(mu[0][0], mu[0][1], 'o', markersize = 16, color = 'red', label = 'Mu 1')
    plt.plot(mu[1][0], mu[1][1], 'o',  markersize = 16, color = 'slateblue', label = 'Mu 2')
    plt.xlabel('Diameter')
    plt.ylabel('Weight')
    plt.legend()
    plt.show()
    
#          
def graph_02_true(X, mu):
    print ' "  "'
    plt.figure(figsize=(16, 6))  
    plt.plot(
        X[0:N1,0], X[0:N1,1], 'o', alpha = 0.7, color='sandybrown', label = 'Produced on machine #1')
    plt.plot(
        X[N1:N,0], X[N1:N,1], 'o', alpha = 0.45, color = 'darkblue', label = 'Produced on machine #2')
    plt.plot(mu[0][0], mu[0][1], 'o', markersize = 16, color = 'red', label = 'Mu 1')
    plt.plot(mu[1][0], mu[1][1], 'o',  markersize = 16, color = 'slateblue', label = 'Mu 2')
    plt.xlabel('Diameter')
    plt.ylabel('Weight')
    plt.legend()
    plt.show()

指定参数

#  
k = 2

#      №1
N1 = 6000
#      №2
N2 = 4000
N = N1+N2

#   (  )
m = 2

#    №1
mu_1_1 = 64.
#    №1
mu_1_2 = 14.
#  
sigma_1_1 = 3.5
sigma_1_2 = 1.

#    №2
mu_2_1 = 52.
#    №2
mu_2_2 = 9.5
#  
sigma_2_1 = 2.
sigma_2_2 = 0.7

我们将计算

X = np.zeros((N, m))
np.random.seed(seed=42)
#    ,   №1
X[:N1, 0] = np.random.normal(loc=mu_1_1, scale=sigma_1_1, size=N1)
X[:N1, 1] = np.random.normal(loc=mu_1_2, scale=sigma_1_2, size=N1)
#  ,   №2
X[N1:N, 0] = np.random.normal(loc=mu_2_1, scale=sigma_2_1, size=N2)
X[N1:N, 1] = np.random.normal(loc=mu_2_2, scale=sigma_2_2, size=N2)

#    (   ,    )
y = np.zeros((N))
y[:N1] = np.array((1))
y[N1:N] = np.array((2))

#           (  )
mu  = np.array(([mu_1_1, mu_1_2], [mu_2_1, mu_2_2]))

#        (  )
sigma = np.array(([sigma_1_1, 0.],[0., sigma_1_2], [sigma_2_1, 0.],[0., sigma_2_2]))
sigma = sigma.reshape(k, m, m)

#     ,    №1  №2
w = np.array([float(1./k), float(1./k)])

#   
pj_xi = gaus_func_02(k, m, X, mu, sigma)
proba_X = proba_func_02(pj_xi, w, m)

#   ,   ,        
limit_proba = 0.5

pred_X = pred_x_02(proba_X, limit_proba)
        
#    
print '   :', round(accuracy_score(y, pred_X),3)
print

graph_02_algorithm(pred_X, mu)

graph_02_true(X, mu)

图2.1“根据算法分配产品”

图2.2“产品真实分配”

与上一个示例类似，我们根据上述参数生成了10,000个值

$\mu$ 和

$\sigma$ ，为算法的运行编写了几个函数并启动了它。该代码与上一个示例的代码之间的根本区别在于，这次我们使用矩阵表达式进行计算。

结果，我们的算法显示正确答案的比例等于0.998，这实际上是相当不错的。

正如我们所看到的，这些问题都是相同的-在不同机器上产生的细节错误同时具有相似的尺寸和重量。

编写代码，以便读者可以替换任何参数值

$\mu$ 和

$\sigma$ 并查看算法的工作原理：在这种情况下，质量会下降，在质量上会有所提高。最主要的是不要忘记学习时间表。好吧，我们在解析EM算法（实际上是EM算法本身）时已走到最后一步。

符合EM算法

我们继续以到达仓库的零件为例。但是这一次我们只会知道产品是在两台不同的机器上制造的，有10,000台，每个零件都有直径和大小，我们一无所知。但是任务并没有改变-与以前一样，我们需要从整个大堆随机混合的物品中确定该零件属于哪台机器。

乍看之下，这听起来似乎不切实际，但实际上，我们手中有一个功能强大的工具包：贝叶斯公式和随机变量的概率密度公式。让我们利用所有这些好处。

解决方案

我们将做什么？正如在EM算法中应有的那样，我们首先初始化参数：

假设的可能性是提取在1号机器上生产的零件-

$w_1$ 我们确定假设的可能性，以提取第2号机器上生产的零件

$w_2$ 。只有两个假设，这意味着第一步中的每个假设都等于0.5。

随机变量的数学期望

$\mu$ 定义如下。我们使用随机函数混合所有产品，将集合平均分为两部分，对于每个部分的每个参数（直径，重量），我们确定平均值。

标准偏差

$\sigma$ 从天花板上拿起所谓的东西-在所有方面都将其设置为等于统一。我们以协方差矩阵的格式编写。

我们准备进行算法的第一个E步。使用随机变量的初始化参数，我们确定将每个零件分配给1号机或2号机的概率。

实际上，通过这种方式，我们迈出了第一步。

现在由M步骤决定。这里的一切都很简单。确定了在特定机器上生产每个零件的概率后，我们可以重新计算每个假设的概率-

$w_1$ ，

$w_2$ 和

$\mu$ 和

$\sigma$ 。

我们将进行15个这样的迭代，每个迭代分为两个步骤，

我们看一下代码。

EM算法中使用的函数

#   E-
def e_step(x, k, m, n, w, mu, sigma):
    #      i-     j-  
    pj_xi = []
    for j in range(k):
        det_sigma_j = np.linalg.det(sigma[j])
        factor_1 = 1 / (((2 * math.pi)**(k/2)) * ((det_sigma_j)**0.5))
        factor_2 = []
        for i in x:
            factor_2.append(math.e**float(
                -0.5 * np.matrix(i - mu[j]) * np.matrix(np.linalg.inv(sigma[j])) * np.matrix(i - mu[j]).T))
        pj_xi.append(factor_1 * np.array(factor_2))
    pj_xi = np.array(pj_xi)
    
    #     ,  i-    j- 
    pj_xi_w = []
    for j in range(k):
        pj_xi_w.append(pj_xi[j] * w[j])
    pj_xi_w = np.array(pj_xi_w)
    
    #     i-      
    sum_pj_xi_w = np.sum(pj_xi_w, axis = 0)
    
    #         
    proba_xi = []
    for j in range(k):
        proba_xi.append(pj_xi_w[j] / sum_pj_xi_w)
    
    return np.array(proba_xi)

#  ,    ,            ,
#       
def x_new(proba_xi):
    X1_new_ind = []
    X2_new_ind = []
    X_answers = []

    count = 0
    for x in proba_xi[0]:
        if x >= 0.5:
            X1_new_ind.append(count)
            X_answers.append(1)
        else:
            X2_new_ind.append(count)
            X_answers.append(2)
        count += 1
    #      ,    №1, №2   
    return X1_new_ind, X2_new_ind, X_answers


#   M-
def m_step(x, proba_xi,N):
    w_new = np.sum(proba_xi, axis = 1) / N
    
    #   
    mu_new = (np.array((np.matrix(proba_xi) * np.matrix(X))).T / np.sum(proba_xi, axis = 1)).T
    
    #  
    cov_new = []
    for mu in range(mu_new.shape[0]):
        X_cd = []
        X_cd_proba = []
        count = 0
        for x_i in x:
            cd = np.array(x_i - mu_new[mu])
            X_cd.append(cd)
            X_cd_proba.append(cd * proba_xi[mu][count])
            count += 1
        X_cd = np.array(X_cd)
        X_cd = X_cd.reshape(N, m)
        X_cd_proba = np.array(X_cd_proba)
        X_cd_proba = X_cd_proba.reshape(N, m)

        cov_new.append(np.matrix(X_cd.T) * np.matrix(X_cd_proba))
    cov_new = np.array((np.array(cov_new) / (np.sum(proba_xi, axis = 1)-1)))
    #             ,    
    #       ,        
    if cov_new[0][0][1] < 0:
        cov_new[0][0][1] = 0
    if cov_new[0][1][0] < 0:
        cov_new[0][1][0] = 0
    
    if cov_new[1][0][1] < 0:
        cov_new[1][0][1] = 0
    if cov_new[1][1][0] < 0:
        cov_new[1][1][0] = 0
    
    #   
    sigma_new = cov_new**0.5
    return w_new, mu_new, sigma_new

指定参数

#   
#      №1 ( 1)
N1 = 6000
#      №2 ( 2)
N2 = 4000
#       
N = N1 + N2

#  
k = 2

#    №1
mu_samples_1_1 = 64.
#    №1
mu_samples_1_2 = 14.

#    №2
mu_samples_2_1 = 52.
#    №2
mu_samples_2_2 = 9.5

#      №1
sigma_samples_1_1 = 3.5
#      №1
sigma_samples_1_2 = 1.

#      №2
sigma_samples_2_1 = 2.
#      №2
sigma_samples_2_2 = 0.7

我们将计算

X = np.zeros((N, 2))

np.random.seed(seed=42)
#    ,    №1
X[:N1, 0] = np.random.normal(loc=mu_samples_1_1, scale=sigma_samples_1_1, size=N1)
X[:N1, 1] = np.random.normal(loc=mu_samples_1_2, scale=sigma_samples_1_2, size=N1)
#    ,    №2
X[N1:N, 0] = np.random.normal(loc=mu_samples_2_1, scale=sigma_samples_2_1, size=N2)
X[N1:N, 1] = np.random.normal(loc=mu_samples_2_2, scale=sigma_samples_2_2, size=N2)

#   
m = X.shape[1]

#   
n = X.shape[0]

#        (   )
y = np.zeros((N))
y[:N1] = np.array((1))
y[N1:N] = np.array((2))

#     ,    №1  №2
w = np.array([float(1./k), float(1./k)])

np.random.seed(seed = None)
#        (   )
mu  = np.array(
    (np.mean(X[np.random.choice(n, n/k)], axis = 0), np.mean(X[np.random.choice(n, n/k)], axis = 0)))
# mu = np.array(([mu_samples_1_1, mu_samples_1_2],[mu_samples_2_1, mu_samples_2_2]))

#         (    )
sigma = np.array(([1., 0.],[0., 1.], [1., 0.],[0., 1.]))
# sigma = np.array(([sigma_samples_1_1, 0.],[0., sigma_samples_1_2], [sigma_samples_2_1, 0.],[0., sigma_samples_2_2]))
sigma = sigma.reshape(k, m, m)

#    EM-
steps = 15
#   EM-
for i in range(steps):
    proba_xi = e_step(X, k, m, n, w, mu, sigma)
    w, mu, sigma = m_step(X, proba_xi,N)
    X1_new_ind, X2_new_ind, X_answers = x_new(proba_xi)
    print ' №', i+1
    print
    print '   '
    print mu
    print
    print '   '
    print sigma
    print
    print '   '
    print round(accuracy_score(y, X_answers),3)
    
    plt.figure(figsize=(16, 6))  
    plt.plot(
        X[X1_new_ind,0], X[X1_new_ind,1], 'o', alpha = 0.7, color='sandybrown', label = 'Produced on machine #1')
    plt.plot(
        X[X2_new_ind,0], X[X2_new_ind,1], 'o', alpha = 0.45, color = 'darkblue', label = 'Produced on machine #2')
    plt.plot(mu[0][0], mu[0][1], 'o', markersize = 16, color = 'red', label = 'Mu 1')
    plt.plot(mu[1][0], mu[1][1], 'o',  markersize = 16, color = 'purple', label = 'Mu 2')
    plt.xlabel('Diameter')
    plt.ylabel('Weight')
    plt.legend()
    plt.show()

图表3.1“ EM算法每次迭代中的零件分布”

所有迭代

结局是什么？结果表明，没有关于提取1号机器上生产的零件的假设可能性的初步信息-

$w_1$ 或机器编号2

$w_2$ ，随机变量的数学期望和标准偏差-

$\mu$ 和

$\sigma$ ，我们达到了与上述示例相同的质量，在该示例中，上述所有参数都是已知的。

比较生成产品集时我们设置的参数，所谓的真实参数和EM算法为我们确定的参数：

真实参数

$\mu_{11} = 64, \mu_{12} = 14, \sigma_{11} = 3.5, \sigma_{12} = 1.$

$\mu_{21} = 52, \mu_{12} = 9.5, \sigma_{21} = 2.0, \sigma_{22} = 0.7$

具体参数

$\mu_{11} \approx 64.03, \mu_{12} \approx 13.99, \sigma_{11} \approx 3.44, \sigma_{12} \approx 1.23$

$\mu_{21} \approx 52.06, \mu_{12} \approx 9.54, \sigma_{21} \approx 1.65, \sigma_{22} \approx 0.77$

关于本文的结尾。最后，我要补充一点，出于数据分析的目的，重新发明轮子并自行编写EM算法是没有意义的，您只需使用提供的sklearn库函数即可。

让我们看看sklearn GaussianMixture库如何完成我们的示例。

导入模块

from sklearn.mixture import GaussianMixture as GMM

指定参数

#   
#      №1 ( 1)
N1 = 6000
#      №2 ( 2)
N2 = 4000
#       
N = N1 + N2

#    №1
mu_samples_1_1 = 64.
#    №1
mu_samples_1_2 = 14.

#    №2
mu_samples_2_1 = 52.
#    №2
mu_samples_2_2 = 9.5

#      №1
sigma_samples_1_1 = 3.5
#      №1
sigma_samples_1_2 = 1.

#      №2
sigma_samples_2_1 = 2.
#      №2
sigma_samples_2_2 = 0.7


#  
k = 2

X = np.zeros((N, 2))

np.random.seed(seed=42)
#    ,    №1
X[:N1, 0] = np.random.normal(loc=mu_samples_1_1, scale=sigma_samples_1_1, size=N1)
X[:N1, 1] = np.random.normal(loc=mu_samples_1_2, scale=sigma_samples_1_2, size=N1)
#    ,    №2
X[N1:N, 0] = np.random.normal(loc=mu_samples_2_1, scale=sigma_samples_2_1, size=N2)
X[N1:N, 1] = np.random.normal(loc=mu_samples_2_2, scale=sigma_samples_2_2, size=N2)

#        (   )
y = np.zeros((N))
y[:N1] = np.array((1))
y[N1:N] = np.array((2))

运行EM算法

np.random.seed(1)
model = GMM(n_components=k, covariance_type='full')
model.fit(X)


temp_predict_X = model.predict(X)
X_answers = []
for i in range(X.shape[0]):
    if temp_predict_X[i] == 0:
        X_answers.append(1)
    else:
        X_answers.append(2)
        

print '   '
print round(accuracy_score(y, X_answers),3)

一起工作了！没错！如您所料，skussian的GaussianMixture库比我们的算法更快，更好。但是，正如我们回想的那样，本文的目的是了解EM算法的工作原理。我希望这个目标已经实现：)

信息来源

文献

1）应用统计：分类和降维，S.A。阿瓦兹扬（弗吉尼亚州）布赫斯塔伯（美国）恩尤科夫（美国） Meshalkin，莫斯科，金融和统计，1989

互联网资源

1）国际计算机科学学院，Jeff A. Bilmes，1998年

2）“ EM算法的温和教程及其在高斯混合和隐马尔可夫模型的参数估计中的应用”。纽约大学戴维·桑塔格（David Sontag

）的“混合模型与EM算法”演示 3）课程“搜索数据结构”

4）完全概率公式和贝叶斯公式

5）函数“ sklearn.mixture.GaussianMixture”

6）正态分布，维基百科

7）多维正态分布，维基百科

8）github.com/diefimov/MTH594_MachineLearning/blob/master/ipython/Lecture10.ipynb

9）关于EM算法的文章

10）关于EM算法的另一篇文章

11）Github github.com/AlexanderPetrenko83/Articles

我们将EM算法分解成小块

一般的EM算法

回顾贝叶斯公式

贝叶斯定理的例子一

使用随机变量正态分布参数的贝叶斯定理的第二个例子：数学期望 μ\mu 和标准偏差 σ\sigma

第三个例子。二维情况

符合EM算法

信息来源

More articles:

使用随机变量正态分布参数的贝叶斯定理的第二个例子：数学期望 $\mu$ 和标准偏差 $\sigma$