EM算法 高斯混合模型 方差估计

1.EM算法介绍

E：Expection，期望步，利用估计的参数，来确定未知因变量的概率，并利用其来计算期望值。

M：Maximization，最大化，使用最大似然法更新参数值，使E步中期望值出现的概率最大。

例如网上较多的硬币例子，可以先估算硬币正反面参数A，但是无法获知隐变量B（无法知道某一次实验选择哪一枚硬币），因此可以分别计算每次试验选择了某一枚硬币的概率，也就是说计算了隐变量B的概率。明确了隐变量B的概率后，就可以依此概率，计算出一个期望值E。

通过更新参数A，使期望值为E的概率最大化，也就是M步中“最大化”的含义。

可以参考链接https://zhuanlan.zhihu.com/p/78311644，这里面详细记录了硬币案例的EM迭代过程。

同样，将情况复杂化，例如现在有20位男生和20位女生，数据记录了他们的身高（但是不知道被记录者的性别），需要分别估计男女生身高分布的均值和方差。

这个时候有人就想到我们必须从某一点开始，并用迭代的办法去解决这个问题：我们先设定男生身高和女生身高分布的几个参数（初始值），然后根据这些参数去判断每一个样本（人）是男生还是女生，之后根据标注后的样本，计算出一个期望值，然后反过来重新估计参数，使期望值出现的概率最大化。之后再多次重复这个过程，直至稳定。这个算法也就是EM算法。
https://zhuanlan.zhihu.com/p/56377602

2.EM算法 高斯混合模型 方差估计

参考https://blog.csdn.net/chasdmeng/article/details/38709063，可以了解均值参数是如何估计的，参考https://zhuanlan.zhihu.com/p/411925257，可以了解均值和方差两个参数是如何估计的。
在sklearn中，也有高斯混合模型的实现，但是没有给出方差如何求解。下面自己实现了一个近似的求解，可以模拟前面提到的男女生身高问题。

数据生成：

def generate_data(sigma_boy, sigma_girl, mu_boy, mu_girl, num_boy, num_girl):
    '''
    根据参数生成正态分布的男生身高数据boy，和女生身高数据girl
    将二者合并为all_Student,打散并返回

    :param sigma_boy: 男生均值
    :param sigma_girl: 女生均值
    :param mu_boy: 男生方差
    :param mu_girl: 女生方差
    :param num_boy: 男生数量
    :param num_girl: 女生数量
    :return: 返回打乱的男女生身高排序
    '''

    boy = list(np.random.normal(sigma_boy, mu_boy, num_boy))
    girl = list(np.random.normal(sigma_girl, mu_girl, num_girl))

    all_student = boy + girl
    random.shuffle(all_student)

    return all_student

根据训练好的gmm模型，反过来倒推每一条身高数据的来源，并根据倒推结果，重新计算均值方差，误差较大，但是这种简单的思路可以借鉴：

def get_covariances(gmm, X, sign=[0.01, 1], real_mean=[1.80, 1.80]):
    '''
        根据predict_proba,保留X中差异度大于sign的数据,
        判断保留的数据属于男性还是女性，
        根据判断结果，将身高数据分别存到new_x_0和new_x_1中，
        最后计算new_x_0和new_x_0的均值、方差。

        尝试不同的sign，根据sign计算出对应的均值mean
        记录mean与real_mean差异最小时,对应的均值和方差到result中

    '''

    result = {}

    X_array = np.array(X).reshape(-1, 1)
    
    # 根据训练好的gmm模型，重新对数据X进行一次预测
    predict_result = gmm.predict(np.array(X).reshape(-1, 1))
    predict_proba = gmm._estimate_log_prob(X_array)

    num_0 = len(predict_result[predict_result == 0])
    num_1 = len(predict_result[predict_result == 1])

    min_diff = float("inf")

    if num_0 > num_1:
        predict_proba_diff = predict_proba[:, 0] - predict_proba[:, 1]

    else:
        predict_proba_diff = predict_proba[:, 1] - predict_proba[:, 0]

    for s in range(sign[0], sign[1], 50):

        predict_proba_diff_s = predict_proba_diff[np.where(predict_proba_diff > s)]
        predict_proba_diff_s_sort = sorted(predict_proba_diff_s)

        # 假设人数相等，这里只取前50%的数据
        tol = predict_proba_diff_s_sort[round(0.5 * len(predict_proba_diff_s))]

        # 过滤数据,并将满足条件的数据标为-1和1
        predict_proba_diff[np.where(predict_proba_diff > s) and np.where(predict_proba_diff < tol)] = -1
        predict_proba_diff[predict_proba_diff >= tol] = 1

        # 取出原始的数据
        new_x_0 = X_array[np.where(predict_proba_diff == -1)]
        new_x_1 = X_array[np.where(predict_proba_diff == 1)]

        # 获取平均数据的最大值、最小值
        max_new_x = max(np.mean(new_x_0), np.mean(new_x_1))
        min_new_x = min(np.mean(new_x_0), np.mean(new_x_1))
        max_real_mean = max(real_mean)
        min_real_mean = min(real_mean)

        # 更新均值和方差数据，并记录到result中
        tmp_diff = abs(max_new_x - max_real_mean) + abs(min_new_x - min_real_mean)

        if tmp_diff < min_diff:
            min_diff = tmp_diff
            result["数据组0"] = [np.mean(new_x_0), np.std(new_x_0)]
            result["数据组1"] = [np.mean(new_x_1), np.std(new_x_1)]

    return result

在sklearn中，给出的方差参数好像与一般理解的方差不是同一个数，使用上面的get_covariances函数，能够用筛选、重建数据、重新估计参数的方式，得到一个非常粗略的近似值：

if __name__ == "__main__":
    # 真实数据分布，模拟男女身高
    sigma_boy = 75
    mu_boy = 10

    sigma_girl = 65
    mu_girl = 8

    # 人数，假设均为20000人
    num_boy = 20000
    num_girl = 20000

    # 生成模拟数据，
    X = generate_data(sigma_boy, sigma_girl, mu_boy, mu_girl, num_boy, num_girl)

    # 训练高斯模型
    gmm = GaussianMixture(n_components=2, tol=1e-3, max_iter=3000).fit(np.array(X).reshape(-1, 1))

    print("************gmm模型估计的参数：************")
    print("均值：", list(gmm._get_parameters())[1])
    print("方差：", list(gmm._get_parameters())[2])

    print("************利用gmm模型，重新估算有偏参数************")
    print(get_covariances(gmm, X, sign=[0, 5], real_mean=list(gmm._get_parameters())[1]))

对应的结果为，可以看到估算有偏参数，数量级基本正确

************gmm模型估计的参数：************
均值： [[64.14610149], [77.51365454]]
方差： [[[54.43250651]],[[72.62747836]]]
 
************利用gmm模型，重新估算有偏参数************
{'数据组0': [74.64721634675828, 7.857091825513232], '数据组1': [57.96871415889782, 4.492870796252531]}