实验十:EM算法_em算法的报告-CSDN博客

本实验是在Anaconda下的jupyter notebook上进行的，使用的代码语言为python。在开始实验前，我们首先需要导入所需要的库与包或者模块。本实验是一个EM算法的实验，需要处理大量的实验数据，需要处理多维数组对象，以及可能还需要画图进行可视化处理，还有一些数学公式的运用，所以我们需要导入的包为numpy、pandas、math以及matplotlib.pyplot。同时我们可能还需要使用random函数来随机选取标签，所以我们还需要导入random。

实验要求：

代码实现：

from collections import Counter
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import math
import random

E步

导入数据集

实验要求：

代码实现：

# 将数据集data.csv载入
D = pd.read_csv('data.csv')

# 将数据集D转换为列表
D=np.array(D)

print(D)

实验结果分析：如图所示，我们首先使用pd.read_csv将数据集data.csv进行导入，导入的数据格式为dataframe，由于在后面的实验中使用该格式的数据不太方便，所有在这里我们需要使用函数np.array()将dataframe的数据转变为列表形式的数据。最后我们再将数据进行打印出来。

初始化t=0时的各个参数

实验要求：

代码实现：

# 将t=0时的所有参数进行初始化(随意初始化即可，但是为了方便进行测试，还是初始化为这个数据集正确的参数)
alpha1=0.625
alpha2=0.375
mu1=175
mu2=165
sigma1=4
sigma2=6

# 将所有的参数组合成一个parameter变量
parameter=[alpha1,alpha2,mu1,mu2,sigma1,sigma2]

print(parameter)

实验结果分析：如图所示，我们首先定义6个参数，分别为alpha1、alpha2、mu1、mu2、sigma1和sigma2，然后将它们分别初始化为这个数据集原本的参数大小。最后我们再将所有的参数组合成为一个parameter变量，同时将这个变量进行打印出来。

编写函数计算每个样本的P（xi,zi|theta）

实验要求：

代码实现：

#parameter=[alpha1,alpha2,mu1,mu2,sigma1,sigma2]
# 编写函数计算每一个样本xi的P_xi
def P(x,parameter,z):
    # 获取alpha1、mu1、sigma1
    Theta1=parameter[::2]
    # 获取alpha2、mu2、sigma2
    Theta2=parameter[1:7:2]
    
    # 计算概率密度函数f_xi
    f1_xi=(1/(math.sqrt(2*math.pi)*Theta1[2]))*(math.e**(-(x-Theta1[1])*(x-Theta1[1])/(2*Theta1[2]*Theta1[2])))
    f2_xi=(1/(math.sqrt(2*math.pi)*Theta2[2]))*(math.e**(-(x-Theta2[1])*(x-Theta2[1])/(2*Theta2[2]*Theta2[2])))
    
    # 当𝑧𝑖=1时
    if z==1:
        P_xi=Theta1[0]*f1_xi
    
    # 当𝑧𝑖=2时
    elif z==2:
        P_xi=Theta2[0]*f2_xi
    
    # 当𝑧𝑖为其它值时
    else:
        P_xi=0
        
    return P_xi

print("计算得到每个样本 𝑥𝑖 的 𝑃(𝑥𝑖,𝑧𝑖|𝜃)：")
print("当𝑧𝑖=1时：")
print(P(D,parameter,1))
print("\n当𝑧𝑖=2时：")
print(P(D,parameter,2))

实验结果分析：如图所示，我们编写一个函数P用于计算每一个样本xi的P_xi，我们首先从parameter中获取Theta1和Theta2的值，然后再通过公式编写计算得到该模型数据的概率密度，随后再对zi进行分类讨论。当zi为1时，我们代入公式使用alpha1乘上前面计算得到的f1_xi；当zi为2时，我们代入公式使用alpha2乘上前面计算得到的f2_xi；而当zi为其它值时，我们将计算的结果返回为0。最后我们再将计算结果作为返回参数返回即可。

随后我们开始进行测试，使用前面初始化的参数，分别当zi为1、zi为2时调用前面编写的函数进行计算，然后将结果打印出来。计算的结果如图所示，我们与其它人或者与这一个数据集本来的数据进行比对，发现这一个计算结果是正确无误的，所以我们可以知道我们的这一个部分的代码编写是没有什么问题的。

编写函数计算样本的y1,i=P((zi=1)|xi,theta)和y2,i=P((zi=2)|xi,theta)

实验要求：

代码实现：

#parameter=[alpha1,alpha2,mu1,mu2,sigma1,sigma2]
# 编写函数计算每一个样本xi的y_zi
def Y(x,parameter,z):
    # 通过以上的P函数计算得到每一个样本xi的P_xi
    P_xi_z1=P(x,parameter,1)
    P_xi_z2=P(x,parameter,2)
    
    # 当𝑧𝑖=1时
    if z==1:
        y_zi=P_xi_z1/(P_xi_z1+P_xi_z2)
    
    # 当𝑧𝑖=2时
    elif z==2:
        y_zi=P_xi_z2/(P_xi_z1+P_xi_z2)
    
    # 当𝑧𝑖为其它值时
    else:
        y_zi=0
        
    return y_zi

print("计算每个样本 𝑥𝑖 的 𝑦1,𝑖=𝑃((𝑧𝑖=1)|𝑥𝑖,𝜃) 和 𝑦2,𝑖=𝑃((𝑧𝑖=2)|𝑥𝑖,𝜃)：")
print("当𝑧𝑖=1时：")
print(Y(D,parameter,1))
print("\n当𝑧𝑖=2时：")
print(Y(D,parameter,2))

实验结果分析：如图所示，我们编写一个函数用于计算每一个样本xi的y_zi，我们首先通过上面定义好的函数P来计算每一个样本xi的P_xi，然后对zi进行分类讨论。当zi为1时，我们代入公式进行计算得到y_zi，然后当zi为2时，我们一样是通过公式进行代入计算得到y_zi，最后是当zi为其它值时，我们将函数的返回值直接设为0。最后我们再将计算得到的结果作为函数的返回值进行返回即可。

编写函数计算对数似然函数在该分布下的期望值Q(theta)

实验要求：

代码实现：

#parameter=[alpha1,alpha2,mu1,mu2,sigma1,sigma2]
# 编写一个函数计算每一个样本的期望值Q_theta
def Q(x,parameter):
    # 使用Y函数获取每一个样本的y_zi
    y_zi_1=Y(x,parameter,1)
    y_zi_2=Y(x,parameter,2)
    
    # 使用P函数获取每一个样本的P_xi_z1和P_xi_z2
    P_xi_z1=P(x,parameter,1)
    P_xi_z2=P(x,parameter,2)
    
    # 计算对数似然函数在该分布下的期望值
    Q_theta=np.sum((y_zi_1*np.log(P_xi_z1)+y_zi_2*np.log(P_xi_z2)))
    
    return Q_theta
   
print("计算对数似然函数在该分布和基于 𝜃(0) 下的期望值 𝑄(𝜃)为：")
print(Q(D,parameter))

实验结果分析：如图所示，我们编写一个函数用于计算每一个样本的期望值Q_theta。我们首先使用前面定义好的函数Y获取每一个样本的y_zi_1和y_zi_2，然后使用P函数用于获取每一个样本的P_xi_z1和P_xi_z2。最后我们便可以使用numpy的内置函数np.sum以及要求中的公式来计算对数似然函数在该分布下的期望值，在这里计算期望值，我搜索的资料是计算所有样本的期望值的平均值或者总值均可，但是其实，无论是平均值还是总值，对于后面的参数更新都是没有任何的影响的，为了方便起见，我们在这里选用的是计算所有样本的期望值的总值。最后我们再将这一个期望值作为函数的返回值进行返回即可。

接下来我们便开始进行测试，我们直接调用该函数Q，然后将返回的结果进行打印即可。结果如图所示，我们将该结果与其他同学以及该数据集的原本数据进行比较，发现我们的结果是正确无误的，所以我们可以知道我们所编写的代码是没有问题的。

M步

编写函数更新第(t+1)轮的alpha

实验要求：

代码实现：

#parameter=[alpha1,alpha2,mu1,mu2,sigma1,sigma2]
# 编写一个函数用来计算更新第t+1轮的alpha1和alpha2的值
def alpha_expection(D,parameter):
    # 初始化两个变量用于存储所有样本的总y1_i和y2_i
    y_zi_1_total=0
    y_zi_2_total=0
    
    # 开始遍历所有的样本
    for x in D:
        # 使用Y函数计算每一个样本的y_zi_1和y_zi_2
        y_zi_1=Y(x,parameter,1)
        y_zi_2=Y(x,parameter,2)
        
        # 将其累加到样本的总y1_i和y2_i中
        y_zi_1_total+=y_zi_1
        y_zi_2_total+=y_zi_2
        
    # 更新alpha1和alpha2
    alpha1_t_add_1=(y_zi_1_total/len(D))
    alpha2_t_add_1=(y_zi_2_total/len(D))
    
    return alpha1_t_add_1,alpha2_t_add_1

print("计算第 1 轮的 (𝛼1(1),𝛼2(1)) 的更新值为：")
print(alpha_expection(D,parameter))

实验结果分析：如图所示，我们编写一个函数用于计算更新第(t+1)轮的alpha1和alpha2的参数值。我们首先需要初始化两个变量用于存储所有样本的总y1_i和y2_i，我们将它们初始化为0。随后我们便可以开始遍历所有的样本，首先使用Y函数计算每一个样本的y_zi_1和y_zi_2，然后将其进行累加到前面定义的y_zi_1_total和y_zi_2_total中。在遍历完所有的样本后，我们便可以更新我们的alpha1和alpha2，直接代入公式即可。最后我们将更新后的参数作为函数的返回值直接返回即可。

接下来我们便可以开始我们的测试，我们直接调用该函数，然后将这个函数的返回值打印出来即可。随后我们将这个结果与其他同学的进行比较，结果发现我们的结果是正确无误的，所以我们可以知道我们的代码是正确无误的。

编写函数更新第(t+1)轮的mu

实验要求：

代码实现：

#parameter=[alpha1,alpha2,mu1,mu2,sigma1,sigma2]
# 编写一个函数用于计算更新第t+1轮的mu1和mu2的值
def mu_expection(D,parameter):
    # 初始化变量用于存储所有样本的总y_zi_1_xi，y_zi_2_xi和y_zi_1,y_zi_2
    y_zi_1_xi_total=0
    y_zi_2_xi_total=0
    y_zi_1_total=0
    y_zi_2_total=0
    
    # 开始遍历每一个样本
    for x in D:
        # 使用Y函数计算每一个样本的y_zi_1和y_zi_2
        y_zi_1=Y(x,parameter,1)
        y_zi_2=Y(x,parameter,2)
        
        # 进行累加
        y_zi_1_xi_total+=(y_zi_1*x)
        y_zi_2_xi_total+=(y_zi_2*x)
        y_zi_1_total+=(y_zi_1)
        y_zi_2_total+=(y_zi_2)
        
    # 更新mu1和mu2
    mu1_add_1=y_zi_1_xi_total/y_zi_1_total
    mu2_add_1=y_zi_2_xi_total/y_zi_2_total
    
    return mu1_add_1,mu2_add_1

print("计算第 1 轮的 (𝜇1(1),𝜇2(1)) 的更新值：")
print(mu_expection(D,parameter))

实验结果分析：如图所示，我们编写一个函数用于计算更新第(t+1)轮的mu1和mu2的值。我们首先需要初始化变量用于存储所有样本的总y_zi_1、y_zi_2和y_zi_1_xi、y_zi_2_xi，我们将它们全部初始化为0。然后我们便开始遍历每一个样本，我们需要使用Y函数来计算每一个样本的y_zi_1和y_zi_2，然后将其按照公式进行累加即可。在遍历完了所有的样本后，我们便可以更新mu1和mu2，调用公式进行更新即可，最后我们再将更新后的参数进行返回。

接下来我们便可以开始进行测试，我们直接调用该函数即可，然后我们再将函数的返回值进行打印出来，结果如图所示。随后我们将我们返回的结果与其它同学的进行比较，发现我们的结果是正确无误的，所以我们可以知道我们的代码编写是没有问题的。

编写函数更新第(t+1)轮的sigma

实验要求：

代码实现：

#parameter=[alpha1,alpha2,mu1,mu2,sigma1,sigma2]
# 编写一个函数用来计算更新第t+1轮的sigma1和sigma2的值
def sigma_expection(D,parameter,mu_next_1,mu_next_2):
    # 初始化变量用于存储所有样本的总的y_zi_1,y_zi_2和y_zi_1_xi_mu1_add1,y_zi_2_xi_mu2_add1
    y_zi_1_total=0
    y_zi_2_total=0
    y_zi_1_xi_mu1_add1_total=0
    y_zi_2_xi_mu2_add1_total=0
    
    # 开始遍历每一个变量
    for x in D:
        # 使用Y函数计算每一个样本的y_zi_1和y_zi_2
        y_zi_1=Y(x,parameter,1)
        y_zi_2=Y(x,parameter,2)
        
        # 开始累加
        y_zi_1_total+=y_zi_1
        y_zi_2_total+=y_zi_2
        y_zi_1_xi_mu1_add1_total+=(y_zi_1*(x-mu_next_1)*(x-mu_next_1))
        y_zi_2_xi_mu2_add1_total+=(y_zi_2*(x-mu_next_2)*(x-mu_next_2))
        
    # 开始更新sigma1和sigma2
    sigma1_add_1=np.sqrt(y_zi_1_xi_mu1_add1_total/y_zi_1_total)
    sigma2_add_1=np.sqrt(y_zi_2_xi_mu2_add1_total/y_zi_2_total)
    
    return sigma1_add_1,sigma2_add_1

print("计算第 1 轮的 (𝜎1(1),𝜎2(1)) 的更新值：")
mu1,mu2=mu_expection(D,parameter)
print(sigma_expection(D,parameter,mu1,mu2))

实验结分析：如图所示，我们编写一个函数用于计算更新第(t+1)轮的sigma1和sigma2的值。我们首先需要初始化变量用于存储所有样本的总的y_zi_1、y_zi_2和y_zi_1_xi_mu1_add1、y_zi_2_xi_mu2_add2。接下来我们便可以开始遍历每一个样本，我们首先需要使用Y函数计算每一个样本的y_zi_1和y_zi_2，然后便可以开始累加，将其累加到我们前面初始化的变量中。在遍历完了所有的样本数据后，我们便可以开始更新sigma1和sigma2，我们直接代入公式更新即可。最后我们将计算更新后的结果作为函数的返回值进行返回即可。

E步和M步的迭代过程

实验要求：

代码实现：

# 定义一个函数用于处理EM算法的迭代过程
def EM_algorithm(D,parameter,max_iterations,threshold):
    # 初始化最大期望值为最小值
    best_Q = float('-inf')
    
    # 开始迭代
    for i in range(max_iterations):
        # E步
        Q_new=Q(D,parameter)
        
        # M步
        # 如果每轮参数更新的差值小于阈值
        if abs(Q_new-best_Q)<threshold:
            break
        
        # 获取更新后的参数
        alpha1,alpha2=alpha_expection(D,parameter)
        mu1,mu2=mu_expection(D,parameter)
        sigma1,sigma2=sigma_expection(D,parameter,mu1,mu2)
        
        if Q_new>best_Q:
            # 更新最大化似然函数期望值
            best_Q=Q_new
        
            # 更新parameter
            parameter=[alpha1,alpha2,mu1,mu2,sigma1,sigma2]
        
    return best_Q,parameter


# 初始化参数、迭代数，阈值
max_iterations=100
threshold=0.001

# 使用数据集正确的参数进行训练作为对比
print("正确参数进行训练")
parameter=[0.625, 0.375, 175, 165, 4, 6]
print("𝜃初始值设为：",
      f"[𝛼1: {parameter[0]}, 𝛼2: {parameter[1]}, 𝜇1: {parameter[2]}, 𝜇2: {parameter[3]}, 𝜎1: {parameter[4]}, 𝜎2: {parameter[5]}]")
print("结果如下（ 𝑄(𝜃),𝜃 ）：")
print(EM_algorithm(D, parameter, max_iterations, threshold))
print('\n')

# 开始第一次迭代训练
print("训练1")
parameter=[0.8, 0.2, 175, 165, 4, 6]
print("𝜃初始值设为：",
      f"[𝛼1: {parameter[0]}, 𝛼2: {parameter[1]}, 𝜇1: {parameter[2]}, 𝜇2: {parameter[3]}, 𝜎1: {parameter[4]}, 𝜎2: {parameter[5]}]")
print("结果如下（ 𝑄(𝜃),𝜃 ）：")
print(EM_algorithm(D, parameter, max_iterations, threshold))

# 开始第二次训练
print("\n训练2")
parameter=[0.6, 0.4, 180, 160, 3, 7]
print("𝜃初始值设为：",
      f"[𝛼1: {parameter[0]}, 𝛼2: {parameter[1]}, 𝜇1: {parameter[2]}, 𝜇2: {parameter[3]}, 𝜎1: {parameter[4]}, 𝜎2: {parameter[5]}]")
print("结果如下（ 𝑄(𝜃),𝜃 ）：")
print(EM_algorithm(D, parameter, max_iterations, threshold))

# 开始第三次训练
print("\n训练3")
parameter = [0.4, 0.6, 172, 168, 6, 4]
print("𝜃初始值设为：",
      f"[𝛼1: {parameter[0]}, 𝛼2: {parameter[1]}, 𝜇1: {parameter[2]}, 𝜇2: {parameter[3]}, 𝜎1: {parameter[4]}, 𝜎2: {parameter[5]}]")
print("结果如下（ 𝑄(𝜃),𝜃 ）：")
print(EM_algorithm(D, parameter, max_iterations, threshold))

实验结果分析：如图所示，我们需要定义一个函数用于处理EM算法的迭代过程。我们首先需要初始化最大似然函数的期望值为最小值，使用了float(‘-inf’)函数，随后我们便可以开始迭代。通过一个for循环，循环次数为最大的迭代数，然后执行E步，开始计算这一个似然函数对应该分布下的期望值。然后开始执行M步，首先判断似然函数的期望值的更新幅度是否小于最小的阈值，如果小于阈值，那么直接返回，退出循环。接下来继续执行下面的语句，首先通过上面定义的函数分别获取各个更新后的参数。随后如果计算得到的似然函数的期望值大于原来的期望值，那么我们更新期望值，同时更新参数值。然后在迭代结束后，我们将最大的似然函数的期望值和参数作为结果返回即可。

接下来我们便可以开始进行训练，我们首先需要进行初始化，初始化最大的迭代次数为100，最小的更新阈值为0.001。随后我们便可以开始我们的数据训练了。我们首先需要使用该数据集的正确标准参数来进行一次训练，作为标准输出方便对比，标准的参数分别为0.625、0.375、175、165、4、6。而我们在接下来进行三次训练，训练的参数初始值选取为只改变前面的alpha1和alpha2，而不改变后面的四个参数。第一次训练的alpha1为0.8，alpha2为0.2；第二次训练的alpha1为0.6，alpha2为0.4；第三次训练的alpha1为0.4，alpha2为0.6。

实验结果如图所示，我们的标准对比训练的最大似然函数的期望值为-7189.9649688552，而第一次训练的最大似然函数的期望值为-7069.713741887326，第二次训练的最大似然函数的期望值为-7212.587632921946，第三次训练的最大似然函数的期望值为-7246.11947380807。我们通过对比发现，第二次训练的实验结果的最大似然函数的值与标准的距离最近，所以我们认为这三个训练的初始数据中，第二次训练的初始参数的设置是最好的，因为它的最大似然函数的期望值与标准的最为接近。而当我们如果没有标准的最大似然函数的期望值进行对比时，我们一般是选择最大似然函数的期望值较大的对应的参数作为最好的初始参数进行设定。但是由于本次实验中，可以求出该数据集的标准最大似然函数的期望值，所以我们选取最为接近的第二次训练的初始参数作为这三次训练中最好的初始参数。

实验总结

EM算法的概念和基本原理：

EM（Expectation-Maximization）算法是一种迭代优化算法，用于在存在隐变量的概率模型中估计参数。它是一种通过迭代的方式逐步提高参数估计的方法。
EM算法的核心思想是通过两个步骤交替进行：E步（Expectation）和M步（Maximization）。在E步，根据当前的参数估计计算隐变量的后验概率；在M步，最大化完整数据的对数似然函数来更新参数估计。通过交替执行E步和M步，EM算法逐渐收敛到参数的局部最优解。

EM算法实现步骤：

初始化参数：选择合适的初始参数。
E步（Expectation）：根据当前参数估计计算隐变量的后验概率，即计算完整数据的条件概率分布。
M步（Maximization）：最大化完整数据的对数似然函数，通过对隐变量的边缘概率进行求和，更新参数估计。
重复执行E步和M步，直到参数收敛或达到最大迭代次数。

EM算法的特点：

EM算法适用于存在隐变量的概率模型参数估计。
EM算法是一种迭代算法，通过交替的E步和M步逐渐提高参数估计的准确性。
EM算法通常用于无监督学习问题，其中只有观测数据而没有标签信息。
EM算法对初始参数的选择较为敏感，不同的初始参数可能导致不同的局部最优解。

EM算法需要注意的方面：

EM算法只能保证收敛到局部最优解，而不是全局最优解。结果的好坏取决于初始参数的选择和算法的收敛性。
对于高维数据或复杂模型，EM算法可能收敛较慢或陷入局部极小值。可以通过多次随机初始化和增加迭代次数来改善结果。
在实际应用中，需要根据具体问题选择合适的概率模型和对应的EM算法。

EM算法的优势：

EM算法是一种经典的参数估计方法，广泛应用于概率模型的参数估计和无监督学习任务。
EM算法在数学理论上具有良好的收敛性和一致性。

EM算法的局限性：

EM算法对初始参数的选择较为敏感，可能收敛到局部最优解。
EM算法的收敛速度可能较慢，特别是对于高维数据或复杂模型。
EM算法要求能够明确建立观测数据和隐变量之间的联合分布，且能够计算后验概率。

EM算法的收敛性：

EM算法的收敛性是基于两个重要性质：似然函数的单调性和隐变量后验概率的增加性。
1. 似然函数的单调性：在每次迭代的M步中，通过最大化完整数据的对数似然函数来更新参数。由于每次迭代都会增加似然函数的值，因此似然函数是单调增加的。
2. 隐变量后验概率的增加性：在每次迭代的E步中，计算隐变量的后验概率。根据贝叶斯定理，后验概率是由似然函数和先验概率相乘得到的。由于似然函数是单调增加的，因此后验概率也是单调增加的。
基于以上性质，EM算法在每次迭代中都会增加似然函数的值，并且隐变量的后验概率也会增加。因此，EM算法保证了似然函数的收敛性，并最终收敛到局部最优解。

EM算法的扩展和改进：

Mixture Models的EM算法：EM算法在混合模型参数估计中得到了广泛应用。混合模型是由多个概率分布组合而成的模型，其中每个分布对应一个潜在的子群体。EM算法通过估计每个子群体的概率分布和权重来进行参数估计。
高斯混合模型的EM算法：高斯混合模型是一种常见的混合模型，其中每个子群体的概率分布是高斯分布。EM算法在高斯混合模型中用于估计每个高斯分布的均值、方差和权重。
快速EM算法：为了加快EM算法的收敛速度，有一些改进的方法被提出。其中一种方法是使用一阶和二阶导数信息来更新参数，以减少迭代次数。
增量式EM算法：当数据集非常大或不断增长时，传统的EM算法可能不适用。增量式EM算法通过逐步处理数据，每次只使用部分数据来更新参数，以适应大规模数据的情况。

EM算法与其他算法的关系：

EM算法与最大似然估计：EM算法可以看作是最大似然估计在存在隐变量的情况下的一种特定实现。通过通过迭代的方式逐步提高对数似然函数的值，EM算法逼近最大似然估计。
EM算法与梯度下降算法：在M步中，EM算法通过最大化对数似然函数来更新参数。这类似于梯度下降算法中的更新步骤，其中参数的更新方向是对数似然函数的梯度方向。
EM算法与K-means聚类算法：K-means算法可以看作是EM算法在特定情况下的一种特例，即高斯混合模型中方差趋于零的情况。K-means算法可以通过EM算法的一次迭代来实现。