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]}