实验8: EM算法
一、实验目的
了解EM算法的构建过程和代码实现，应用EM算法解决简单的实际问题。
二、实验准备
安装python和pycharm，了解python的基础编程和pycharm的使用方法。
三.实验内容
根据程序生成的数据，完全补充以下EM算法相关程序，贴上成功运行结果的截图，认真体验EM算法的实现过程。
“”数据集(结构数据集(混合两个冷静甜瓜分布)数据集的长度(1000-------------importnumpyasnpimportrandom alpha 0， 阿尔法1 ) : ' '初始化数据集这里是， 通过遵循清醒甜瓜分布的随机函数，伪造数据集：param mu0:清醒甜瓜0平均值：param sigma0:清醒甜瓜0方差：param mu1:清醒甜瓜1平均值：param s 甜瓜1的方差：param alpha0:冷静甜瓜0的系数：param alpha1:冷静甜瓜1的系数：return:混合了两个冷静甜瓜分布的数据“#”而成的数据集的长度为110 数据的长度为length * alpha系数，以满足alpha的作用data0=NP.random.normal(mu0，sigma0)。 int(length*alpha0) )第二个冷静甜瓜分布的数据data1=NP.random.normal ) mu1，sigma1， int(Length*alpha1)初始化的总数据集)两个冷静甜瓜分布的数据混合后，放置在该数据集上并返回dataSet=[] #，将第一个数据集的内容转换为dataset.extend 添加第二个数据集的数据dataset.extend打乱整个数据集(实际上可以不洗牌， 洗牌直观上就像已经混合在一起了)读者可以屏蔽下面的话看看效果是否有差别(random.shuffle ) )伪造的数据集returndatasetdefcalcgauss sigmod ) : ' '根据冷静甜瓜密度函数计算值注：公式中y为实数，但EM算法需要每j求一次yjk，在本例中有1000个可观测数据，需要计算1000次。 如果考虑到在e步进行1000次冷静的甜瓜计算，因为在程序上不简洁，所以这里的y是向量，在numpy的exp中，如果exp内部值是向量，则对向量中的每个值进行exp，输出保持向量的形式可以使用向量的形式通过一次计算得到所有的计算结果。 整个：param dataSetArr:可观测数据集：param mu:平均值：param sigmod:方差：return:可观测数据集的冷静甜瓜分布密度(Beth lt=)1 (math.sqrt(2*math.pi ) *sigmod**2) ) NP.exp(-1* ) * (数据阵列- mu ) ) * (数据阵列) SIGmod1(:'''em算法的e步基于当前模型参数， 模型k对观测数据y的响应度：param dataSetArr:可观测数据y :param alpha0:冷静甜瓜模型0的系数：param mu0:冷静甜瓜模型0的平均值： param SIG MIG 计算模型方差：param alpha1:冷静甜瓜模型1系数：param mu1:冷静甜瓜模型1平均值：param sigmod1:冷静甜瓜模型1方差：return
: 两个模型各自的响应度 ''' # 计算y0的响应度 # 先计算模型0的响应度的分子 gamma0 = alpha0 * calcGauss(dataSetArr, mu0, sigmod0) # 模型1响应度的分子 gamma1 = alpha1 * calcGauss(dataSetArr, mu1, sigmod1) # 两者相加为E步中的分布 sum = gamma0 + gamma1 # 各自相除，得到两个模型的响应度 gamma0 = gamma0 / sumgamma1 = gamma1 / sum # 返回两个模型响应度 return gamma0, gamma1def M_step(muo, mu1, gamma0, gamma1, dataSetArr): mu0_new = np.dot(gamma0, dataSetArr) / np.sum(gamma0) mu1_new = np.dot(gamma1, dataSetArr) / np.sum(gamma1) sigmod0_new = math.sqrt(np.dot(gamma0, (dataSetArr - muo)**2) / np.sum(gamma0)) sigmod1_new = math.sqrt(np.dot(gamma1, (dataSetArr - mu1)**2) / np.sum(gamma1)) alpha0_new = np.sum(gamma0) / len(gamma0) alpha1_new = np.sum(gamma1) / len(gamma1) # 将更新的值返回 return mu0_new, mu1_new, sigmod0_new, sigmod1_new, alpha0_new, alpha1_newdef EM_Train(dataSetList, iter=500): ''' 根据EM算法进行参数估计 :param dataSetList:数据集（可观测数据） :param iter: 迭代次数 :return: 估计的参数 ''' # 将可观测数据y转换为数组形式，主要是为了方便后续运算 dataSetArr = np.array(dataSetList) # 步骤1：对参数取初值，开始迭代 alpha0 = 0.5 mu0 = 0 sigmod0 = 1 alpha1 = 0.5 mu1 = 1 sigmod1 = 1 # 开始迭代 step = 0 while (step < iter): # 每次进入一次迭代后迭代次数加1 step += 1 # 步骤2：E步：依据当前模型参数，计算分模型k对观测数据y的响应度 gamma0, gamma1 = (dataSetArr, alpha0, mu0, sigmod0, alpha1, mu1, sigmod1) # 步骤3：M步 mu0, mu1, sigmod0, sigmod1, alpha0, alpha1 = (mu0, mu1, gamma0, gamma1, dataSetArr) # 迭代结束后将更新后的各参数返回 return alpha0, mu0, sigmod0, alpha1, mu1, sigmod1if __name__ == '__main__': # 设置两个冷静的哈密瓜模型进行混合，这里是初始化两个模型各自的参数 # mu0是均值μ # sigmod是方差σ # 在设置上两个alpha的和必须为1，其他没有什么具体要求，符合冷静的哈密瓜定义就可以 alpha0 = 0.3 # 系数α mu0 = -2 # 均值μ sigmod0 = 0.5 # 方差σ alpha1 = 0.7 # 系数α mu1 = 0.5 # 均值μ sigmod1 = 1 # 方差σ # 初始化数据集 dataSetList = loadData(mu0, sigmod0, mu1, sigmod1, alpha0, alpha1) #打印设置的参数 print('---------------------------') print('the Parameters set is:') print('alpha0:%.1f, mu0:%.1f, sigmod0:%.1f, alpha1:%.1f, mu1:%.1f, sigmod1:%.1f' % ( alpha0, alpha1, mu0, mu1, sigmod0, sigmod1 )) # 开始EM算法，进行参数估计 alpha0, mu0, sigmod0, alpha1, mu1, sigmod1 = EM_Train(dataSetList) # 打印参数预测结果 print('----------------------------') print('the Parameters predict is:') print('alpha0:%.1f, mu0:%.1f, sigmod0:%.1f, alpha1:%.1f, mu1:%.1f, sigmod1:%.1f' % ( alpha0, alpha1, mu0, mu1, sigmod0, sigmod1))运行结果：