Mix Gaussian In EM

发表于 | 分类于

Mix Gaussian In EM

算法原理简述：

概率模型有时即含有观测变量，又含有隐变量或者潜在变量。如果概率模型的变量都是观测变量，那么给定数据，可以直接用最大似然估计的方法或者贝叶斯估计法估计参数模型。当含有隐变量时，我们就不能直接使用这些方法，EM方法就是含有隐含变量的概率模型参数的极大似然估计法。
EM算法主要包括了:

  • 得到 Q函数
  • E步计算
  • Q步计算

其中Q函数是完全数据的对数似然函数在给定观测数据和当前参数下，对未观测数据的条件概率分布的期望。

本文中实现的是利用EM算法对高斯混合模型的进行参数估计.

代码如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
#!/usr/bin/env python
# -*- coding: utf-8 -*-
# @File  : MixGaussianInEM.py
# @Author: mfzhu
# @Date  : 2017/5/3
# @Desc  :
import numpy as np
import copy
class mixGaussianInEM(object):
    def __init__(self, train_data, num_Gaussian):
        self.epsilon = 0.0000001
        # 设定模型停止迭代的阈值
        self.num = num_Gaussian
        # 设定高斯函数的个数
        self.num_sample = len(train_data)
        # 初始化训练数据个数
        self.sample = train_data
        # 初始化训练数据
        self.mean = np.array([15, 35], float)
        # 初始化高斯函数的均值
        self.var = np.array([15, 35], float)
        # 初始化高斯函数的方差
        self.alpha = np.array([0.4, 0.6], float)
        # 初始化高斯函数的系数
        self.gamma = np.zeros([self.num_sample, self.num], float)
        # 初始化每个观测值在每个高斯函数上的分量
    def norm_pdf(self, miu, sigma, value):
        # 计算某个均值和方差下的value对应的概率值
        return (1 / (np.sqrt(2 * np.pi) * sigma)) * np.exp(-1 * (value - miu) ** 2 / (2 * (sigma ** 2)))
    def E_step(self):
        # E步函数更新gamma矩阵中的各个系数
        for i in range(self.num_sample):
            demon = 0
            for j in range(self.num):
                demon += self.alpha[j] * self.norm_pdf(self.mean[j], np.sqrt(self.var[j]), self.sample[i])
                # 计算某个观测值在所有高斯函数上分量的和
            for j in range(self.num):
                self.gamma[i][j] = self.alpha[j] * self.norm_pdf(self.mean[j], np.sqrt(self.var[j]), self.sample[i]) / demon
                # 计算某个高斯函数的分量
    def Q_step(self):
        #  Q步函数更新均值、方差和系数
        for i in range(self.num):
            # 更新方差
            difference = np.square(self.sample - self.mean[i])
            self.var[i] = np.sum(difference * self.gamma[:, i]) / np.sum(self.gamma[:, i])
        for i in range(self.num):
            # 更新均值
            self.mean[i] = (np.sum(self.gamma[:, i] * self.sample)) / np.sum(self.gamma[:, i])
        for i in range(self.num):
            # 更新系数
            self.alpha[i] = np.sum(self.gamma[:, i]) / self.num_sample
    def train(self):
        time = 0
        # 训练次数
        while time < 1000:
            # print(self.mean, self.var, self.alpha, time)
            old_mean = copy.deepcopy(self.mean)
            old_var = copy.deepcopy(self.var)
            # 记录上一次的均值和方差
            self.E_step()
            self.Q_step()
            # 训练
            if np.sum(np.abs(self.mean - old_mean)) < self.epsilon\
                    or np.sum(np.abs(self.var - old_var)) < self.epsilon:
                break
                # 如果两次结果的差值小于阈值则结束训练
            time += 1
    def output(self):
        # 输出训练结果
        print("the estimation of mean:")
        print(self.mean)
        print("the estimation of var:")
        print(self.var)
        print("the estimation of alpha:")
        print(self.alpha)
if __name__ == '__main__':
    test_data = []
    sigma = [20, 40]
    miu = [20, 40]
    for i in range(10000):
        if np.random.random(1) > 0.5:
            test_data.append(np.random.normal() * sigma[0] + miu[0])
        else:
            test_data.append(np.random.normal() * sigma[1] + miu[1])
    # 初始化训练数据
    test_data = np.array(test_data)
    mix = mixGaussianInEM(test_data, 2)
    mix.train()
    mix.output()

结果说明：

在利用两个高斯分布分别为N(20, 20)和N(40, 40)的情况下，按照概率0.5和0.5分别产生10000个数据。通过EM算法的估计结果如下：
img

可以看到对于参数的估计还是相当准确的。需要注意的是在训练模型的初始值的设定时，偏离与待估计参数太大的情况容易导致估计失败，即EM算法对于初始值敏感，并且不能保证收敛到全局最优。所以我们可以多次选取初始值而后对估计结果取平均来保证参加估计的正确性。