EM算法

EM算法两步求解：
E步：求期望（expectation）
M步：求极大（maximization)
EM 算法是一种求解最大似然估计的方法，通过观测样本，来找出样本的模型参数。
最大似然可以理解为最大可能性，最大似然估计是通过已知结果，估计参数的方法。
如果使用基于最大似然估计的模型，模型中存在隐变量的时候，就要用到EM算法去做估计。
举个例子，扔两枚硬币A、B，每次投掷一枚，分别做5组实验，每组投10次，统计正反面结果为：

若要计算硬币A和B正面向上的概率，该怎么计算？这就无法计算了，因为不知道结果对应的硬币到底是A还是B，而实际上这里A和B就是隐含信息，即隐变量。假设投掷结果如下图：

则硬币A和B正面向上的概率即为：

这里的第二列，就是隐含的数据，而A和B就是隐变量。实际中我们是不知道这一列的，就是开始给你的只有实验组数和正面的次数，那么你该怎么办呢？也就是说，我们如果不知道每一组扔的是A还是B，那么我们就无法去估计θA和θB，而如果想知道每一组扔的是A还是B，我们就必须先知道A和B正面朝上的概率θA和θB，然后利用极大似然的思想，根据每一组实验正面朝上的次数去估计出这一轮究竟用的A还是B。这里就采用了EM算法的思想。
先随机初始化一个θA和θB，有了这两个参数，我们就能按照极大似然估计出每一组用的是A还是B，然后基于每一组用的是A还是B，我们又能按照极大似然反过来计算出θA和θB，然后又能去估计新的用的是A还是B，然后又能计算新的θA和θB，这样一轮轮的下去，当计算出的新的θA和θB与我们前一轮θA和θB一样的时候，说明这个θA和θB有可能就是真实的值了。这个就是EM初级版。

EM 算法中的 E 步骤就是通过旧的参数来计算隐藏变量。然后在 M 步骤中，通过得到的隐藏变量的结果来重新估计参数。直到参数不再发生变化，得到我们想要的结果。
K-means算法与EM算法关系：k-means是最简单的EM算法
K-means算法先假定k个中心，然后进行最短距离聚类，之后根据聚类结果重新计算各个聚类的中心点，一次迭代，是不是很像，而且K-means也是初始值敏感，因此其实K-means算法也包含了EM算法思想，只是这边EM算法中用P概率计算，而K-means直接用最短距离计算。
原理：


KMeans解决不了的问题，EM聚类是可以解决的。在 EM 框架中，我们将潜在类别当做隐藏变量，样本看做观察值，把聚类问题转化为参数估计问题，最终把样本进行聚类。EM 算法相当于一个框架，你可以采用不同的模型来进行聚类，比如 GMM（高斯混合模型），或者 HMM（隐马尔科夫模型）来进行聚类。
GMM 是通过概率密度来进行聚类，聚成的类符合高斯分布（正态分布）。
HMM 用到了马尔可夫过程，在这个过程中，我们通过状态转移矩阵来计算状态转移的概率。HMM 在自然语言处理和语音识别领域中有广泛的应用。
GMM模型：
MM模型全称为Gaussian Mixture Model，即高斯混合模型。其主要是针对普通的单个高斯模型提出来的。我们知道，普通高斯模型对实际数据拟合效果还不错，但是其有一个致命的缺陷，就是其为单峰函数，如果数据的真实分布为复杂的多峰分布，那么单峰高斯的拟合效果就不够好了。
与单峰高斯模型不同，GMM模型是多个高斯模型的加权和，

这是一个多峰分布，理论上，只要k足够大，GMM模型能拟合任何分布。

EM-E步代码：

import numpy as np
import random
import math
import time
'''
数据集：伪造数据集（两个高斯分布混合）
数据集长度：1000
------------------------------
运行结果：
----------------------------
the Parameters set is:
alpha0:0.3, mu0:0.7, sigmod0:-2.0, alpha1:0.5, mu1:0.5, sigmod1:1.0
----------------------------
the Parameters predict is:
alpha0:0.4, mu0:0.6, sigmod0:-1.7, alpha1:0.7, mu1:0.7, sigmod1:0.9
----------------------------
'''

def loadData(mu0, sigma0, mu1, sigma1, alpha0, alpha1):
    '''
    初始化数据集
    这里通过服从高斯分布的随机函数来伪造数据集
    :param mu0: 高斯0的均值
    :param sigma0: 高斯0的方差
    :param mu1: 高斯1的均值
    :param sigma1: 高斯1的方差
    :param alpha0: 高斯0的系数
    :param alpha1: 高斯1的系数
    :return: 混合了两个高斯分布的数据
    '''
    # 定义数据集长度为1000
    length = 1000

    # 初始化第一个高斯分布，生成数据，数据长度为length * alpha系数，以此来
    # 满足alpha的作用
    data0 = np.random.normal(mu0, sigma0, int(length * alpha0))
    # 第二个高斯分布的数据
    data1 = np.random.normal(mu1, sigma1, int(length * alpha1))

    # 初始化总数据集
    # 两个高斯分布的数据混合后会放在该数据集中返回
    dataSet = []
    # 将第一个数据集的内容添加进去
    dataSet.extend(data0)
    # 添加第二个数据集的数据
    dataSet.extend(data1)
    # 对总的数据集进行打乱（其实不打乱也没事，只不过打乱一下直观上让人感觉已经混合了
    # 读者可以将下面这句话屏蔽以后看看效果是否有差别）
    random.shuffle(dataSet)

    #返回伪造好的数据集
    return dataSet

    # 高斯分布公式，没有什么特殊的
def calcGauss(dataSetArr, mu, sigmod):
    '''
    根据高斯密度函数计算值
    依据：“9.3.1 高斯混合模型” 式9.25
    注：在公式中y是一个实数，但是在EM算法中(见算法9.2的E步)，需要对每个j
    都求一次yjk，在本实例中有1000个可观测数据，因此需要计算1000次。考虑到
    在E步时进行1000次高斯计算，程序上比较不简洁，因此这里的y是向量，在numpy
    的exp中如果exp内部值为向量，则对向量中每个值进行exp，输出仍是向量的形式。
    所以使用向量的形式1次计算即可将所有计算结果得出，程序上较为简洁
    
    :param dataSetArr: 可观测数据集
    :param mu: 均值
    :param sigmod: 方差
    :return: 整个可观测数据集的高斯分布密度（向量形式）
    '''
    # 计算过程就是依据式9.25写的，没有别的花样
    result = (1 / (math.sqrt(2*math.pi)*sigmod**2)) * np.exp(-1 * (dataSetArr-mu) * (dataSetArr-mu) / (2*sigmod**2))
    # 返回结果
    return result


def E_step(dataSetArr, alpha0, mu0, sigmod0, alpha1, mu1, sigmod1):
    '''
    EM算法中的E步
    依据当前模型参数，计算分模型k对观数据y的响应度
    :param dataSetArr: 可观测数据y
    :param alpha0: 高斯模型0的系数
    :param mu0: 高斯模型0的均值
    :param sigmod0: 高斯模型0的方差
    :param alpha1: 高斯模型1的系数
    :param mu1: 高斯模型1的均值
    :param sigmod1: 高斯模型1的方差
    :return: 两个模型各自的响应度
    '''
    # 计算y0的响应度
    # 先计算模型0的响应度的分子
    gamma0 = alpha0 * calcGauss(dataSetArr, mu0, sigmod0)
    #print("gamma0=",gamma0.shape) # 1000, 维向量
    # 模型1响应度的分子
    gamma1 = alpha1 * calcGauss(dataSetArr, mu1, sigmod1)

    # 两者相加为E步中的分布
    sum = gamma0 + gamma1
    # 各自相除，得到两个模型的响应度
    gamma0 = gamma0 / sum
    gamma1 = gamma1 / sum

    # 返回两个模型响应度
    return gamma0, gamma1

def M_step