• 有限状态机 FSA

• 确定性的识别器 DFSA

• 确定性的识别器 NFSA：深度优先搜索 or 广度优先搜索

• 形式语言

• 各种正则表达式，析取，组合和优先关系

• 文本标准化：各种预处理方法

• 编辑距离：动态规划

chapter6: 隐马尔科夫模型 standford tutorial, 可以说是看过的关于隐马尔科夫最好的教程了。要是看英文原版能和看中文一样easy该多好，可惜文化差异的情况下，即使是单词都认识，理解起来也会有点困难。对此，只能自己重新总结一下吧～

可以说，斯坦福从未让人失望过，太赞了！

也是无意中在google上看到这篇文章，才发现了这么好的一本书, Speech and language processing.

sklearn源码阅读，用em算法计算高斯混合模型GMM

代码实现高斯混合模型

参考这篇博客Regularized Gaussian Covariance Estimation非常值得一读，同事这篇博客很深入的讲了协方差怎么求的问题，在前文中我也有提到～但我解释的很low。。

代码直接就看sklearn里面的源码吧～网上很多不靠谱。。。

github源码

类初始化

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

class GaussianMixture(BaseMixture):

  """

  Representation of a Gaussian mixture model probability distribution.

  This class allows to estimate the parameters of a Gaussian mixture

  distribution. 对混合的高斯分布进行参数估计～

  """

  def __init__(self, n_components=1, covariance_type='full', tol=1e-3,

             reg_covar=1e-6, max_iter=100, n_init=1, init_params='kmeans',

             weights_init=None, means_init=None, precisions_init=None,

             random_state=None, warm_start=False,

             verbose=0, verbose_interval=10):

    super(GaussianMixture, self).__init__(

        n_components=n_components, tol=tol, reg_covar=reg_covar,

        max_iter=max_iter, n_init=n_init, init_params=init_params,

        random_state=random_state, warm_start=warm_start,

        verbose=verbose, verbose_interval=verbose_interval)



    # 主要是针对3个要学习的参数的初始化

    self.covariance_type = covariance_type # 协方差矩阵形式

    self.weights_init = weights_init # 多项式分布，每一类的概率

    self.means_init = means_init  # 均值 (n_components, n_features)

    self.precisions_init = precisions_init # 协方差

先看初始化构造函数，参数是真的多。。。

• n_components=1: The number of mixture components.表示混合类别的个数，也就是混合高斯分布的个数

• covariance_type=’full’: 协方差矩阵的类型。{‘full’, ‘tied’, ‘diag’, ‘spherical’} 分别对应完全协方差矩阵（元素都不为零），相同的完全协方差矩阵（HMM会用到），对角协方差矩阵（非对角为零，对角不为零），球面协方差矩阵（非对角为零，对角完全相同，球面特性），默认‘full’ 完全协方差矩阵

• tol=1e-3 收敛阈值，EM iterations will stop when the lower bound average gain is

below this threshold.也就是当下界的平均增益小于阈值时，em迭代就停止。这里的下界指的是公式

（3）中的下界凸函数。我们知道em算法分两步，e step是期望，也就是不等式相等，m setp是最大化，

也就是下界凸函数最大化。这里的阈值平均增益就是指凸函数的最大化过程中的增益。

• reg_covar=1e-6： Non-negative regularization added to the diagonal of

covariance.Allows to assure that the covariance matrices are all positive.

非负正则化添加到协方差矩阵对角线上，保证协方差矩阵都是正定的。

• max_iter=100: em算法的最大迭代次数

• n_init: int, defaults to 1.初始化的次数

• init_params: {‘kmeans’, ‘random’}, defaults to ‘kmeans’.

The method used to initialize the weights, the means and the precisionsself.

Must be one of::

-  'kmeans' : responsibilities are initialized using kmeans.

- 'random' : responsibilities are initialized randomly.

- 这里对应的初始化，是指的隐藏变量z的分类所占比例，也就是weight_init，kmeans表示“hard”guess， {0, 1} or {1, . . . , k})

random应该就是”soft”guess吧。

• weights_init : shape (n_components, ), optional The user-provided initial weights, defaults to None. If it None, weights are initialized using the init_params method.

先验权重初始化，对应的就是隐藏变量有n_components类，而每一类所占的比例，也就是多项式分布的初始化～对应$\phi_i$

• means_init : array-like, shape (n_components, n_features), optional. The user-provided initial means, defaults to None, If it None, means are initialized using the init_params method.混合高斯分布的均值初始化，注意shape=(n_components, n_features),有n_components这样的多维高斯分布，每个高斯分布有n_features维度

• precisions_init : The user-provided initial precisions (inverse of the covariance matrices), defaults to None. If it None, precisions are initialized using the ‘init_params’ method.The shape depends on ‘covariance_type’::

  • (n_components,) if ‘spherical’,

  • (n_features, n_features) if ‘tied’,

  • (n_components, n_features) if ‘diag’,

  • (n_components, n_features, n_features) if ‘full’

  • 用来初始化高斯分布中的协方差矩阵，协方差矩阵代表的是n_features维向量中每一维特征与其他维度特征的关系，对于一个高斯分布来说是n_featuresn_features，n_components个混合也就是’full’。其中要学习的参数个数是(n_features+1) n_features/2.具体关于协方差矩阵参考前面那篇博客

• random_state : int, RandomState instance or None, optional (default=None) 随机数生成器

• warm_start : bool, default to False.If ‘warm_start’ is True, the solution of the last fitting is used as initialization for the next call of fit(). This can speed up convergence when fit is called several times on similar problems. 若为True，则fit（）调用会以上一次fit（）的结果作为初始化参数，适合相同问题多次fit的情况，能加速收敛，默认为False。

• verbose : int, default to 0. Enable verbose output. If 1 then it prints the current initialization and each iteration step. If greater than 1 then it prints also the log probability and the time needed for each step. 使能迭代信息显示，默认为0，可以为1或者大于1（显示的信息不同）

• verbose_interval: 与13挂钩，若使能迭代信息显示，设置多少次迭代后显示信息，默认10次。

E step

就是求$w_j^i$

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

def _e_step(self, X):

    """E step.

    Parameters

    ----------

    X : array-like, shape (n_samples, n_features)

    Returns

    -------

    log_prob_norm :

    log_responsibility : 后验概率，样本i是j类的概率w_j^{i}

    """

    log_prob_norm, log_resp = self._estimate_log_prob_resp(X)

    return np.mean(log_prob_norm), log_resp

那么如何求$w_j^{i}$呢？

$$w_j^{(i)}:=p(z^{(i)}=j|x^{(i)};\phi,\mu,\Sigma)=\dfrac{p(x^{(i)}|z^{(i)}=j;\mu,\Sigma)p(z^{(i)}=j;\phi)}{\sum_{l=1}^kp(x^{(i)}|z^{(i)}=l;\mu,\Sigma)p(z^{(i)}=l;\phi)}$$

要注意的是，为了计算方便，有以下几点：

• 因为分子分母中设计到正态分布，即指数形式，故先计算其log形式。然后带入到M step中取回指数形式即可。

• 对于协方差矩阵，如果n_features很大的话，计算其逆矩阵和行列式就很复杂，因此可以先计算其precision矩阵，然后进行cholesky分解，以便优化计算。

先计算分子对数形式，两个对数相加：

$$logp(x^{(i)}|z^{(i)}=j;\mu,\Sigma)+logp(z^{(i)}=j;\phi)$$

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23

# 计算P(x|z)p(z)的对数形式

def _estimate_weighted_log_prob(self, X):

        """Estimate the weighted log-probabilities, log P(X | Z) + log weights.

        Parameters

        ----------

        X : array-like, shape (n_samples, n_features)

        Returns

        -------

        weighted_log_prob : array, shape (n_samples, n_component)

        """

        return self._estimate_log_prob(X) + self._estimate_log_weights()

1. 其中前者是高斯分布概率的对数,根据均值，协方差矩阵的cholesky分解可求得。

1
2
3
4
5
6
7

def _estimate_log_prob(self, X):

    return _estimate_log_gaussian_prob(

            X, self.means_, self.precisions_cholesky_, self.covariance_type)

这个函数，_estimate_log_gaussian_prob根据高斯分布的参数计算概率，涉及到协方差矩阵，要优化计算，很复杂，放在最后说。先把整个流程走完。

2. 后者是每一类高斯分布所占的权重，也就是$\phi_j$

1
2
3
4
5
6
7

def _estimate_log_weights(self):

    # 刚开始是初始值，后面随着m step而更新

    return np.log(self.weights_)

再计算$w_j^i$

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

def _estimate_log_prob_resp(self, X):

    """Estimate log probabilities and responsibilities for each sample.

    Compute the log probabilities, weighted log probabilities per

    component and responsibilities for each sample in X with respect to

    the current state of the model.

    Parameters

    ----------

    X : array-like, shape (n_samples, n_features)

    Returns

    -------

    log_prob_norm : array, shape (n_samples,)

        log p(X)

    log_responsibilities : array, shape (n_samples, n_components)

        logarithm of the responsibilities

    """

    # 计算分子log P(X | Z) + log weights.

    weighted_log_prob = self._estimate_weighted_log_prob(X)

    # 计算分母log P(x)

    log_prob_norm = logsumexp(weighted_log_prob, axis=1)

    with np.errstate(under='ignore'):

        # 忽略下溢，计算log(w_J^j)，也就是两个对数相减

        log_resp = weighted_log_prob - log_prob_norm[:, np.newaxis]

    return log_prob_norm, log_resp

M setp

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

def _m_step(self, X, log_resp):

    """M step.

    Parameters

    ----------

    X : array-like, shape (n_samples, n_features)

    log_resp : array-like, shape (n_samples, n_components)

        Logarithm of the posterior probabilities (or responsibilities) of

        the point of each sample in X.

    """

    n_samples, _ = X.shape

    # 根据E step中求得的log_resp,更新权重，均值和协方差

    self.weights_, self.means_, self.covariances_ = (

        _estimate_gaussian_parameters(X, np.exp(log_resp), self.reg_covar,

                                      self.covariance_type))

    # 更新类别权重phi_j

    self.weights_ /= n_samples

    # 更新协方差矩阵的精度矩阵

    self.precisions_cholesky_ = _compute_precision_cholesky(

        self.covariances_, self.covariance_type)

具体怎么求，就是根据前面推导的公式了。根据前面的公式分别求对应的估计参数：

$$\Sigma_j:=\dfrac{\sum_{i=1}^mw_j^{(i)}(x^{(i)}-\mu_j)(x^{(i)}-\mu_j)^T}{\sum_{i=1}^mw_j^{(i)}}$$

协方差矩阵：以‘full’为例

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

def _estimate_gaussian_covariances_full(resp, X, nk, means, reg_covar):

    """Estimate the full covariance matrices.

    resp:表示E step中猜测的w_j^{i}

    """

    n_components, n_features = means.shape

    # 协方差矩阵

    covariances = np.empty((n_components, n_features, n_features))

    for k in range(n_components):

        diff = X - means[k]

        covariances[k] = np.dot(resp[:, k] * diff.T, diff) / nk[k]

        # 正则化，flat表示展开成一维，然后每隔n_features取一个元素，单个协方差矩阵shape是

        # [n_features, n_features],所以就是对角线元素加上reg_covar

        covariances[k].flat[::n_features + 1] += reg_covar

    return covariances

然后是正态分布的参数估计$u_j, \phi_j$

$$\phi_j:=\frac{1}{m}\sum_{i=1}^mw_j^{(i)}$$

$$\mu_j:=\dfrac{\sum_{i=1}^mw_j^{(i)}x^{(i)}}{\sum_{i=1}^mw_j^{(i)}}$$

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

def _estimate_gaussian_parameters(X, resp, reg_covar, covariance_type):

    """Estimate the Gaussian distribution parameters.

    Parameters

    ----------

    X : 样本数据 (n_samples, n_features)

    resp : Estep猜测的样本i是j类的概率w_i^{j}, shape (n_samples, n_components)

    reg_covar : 对角线正则化项

    covariance_type : {'full', 'tied', 'diag', 'spherical'}

    Returns

    -------

    nk : 当前类别下的样本和 (n_components,) 也就是\sum_i^{m}(w_j^{i})

    means : k个n维正态分布的均值, shape (n_components, n_features)

    covariances : 协方差矩阵

    """

    # 因为要做分母，避免为0

    nk = resp.sum(axis=0) + 10 * np.finfo(resp.dtype).eps

    means = np.dot(resp.T, X) / nk[:, np.newaxis]

    covariances = {"full": _estimate_gaussian_covariances_full,

                   "tied": _estimate_gaussian_covariances_tied,

                   "diag": _estimate_gaussian_covariances_diag,

                   "spherical": _estimate_gaussian_covariances_spherical

                   }[covariance_type](resp, X, nk, means, reg_covar)

    return nk, means, covariances

迭代收敛，重复以上过程

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
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161

def fit(self, X, y=None):

    """Estimate model parameters with the EM algorithm.



    The method fit the model `n_init` times and set the parameters with

    which the model has the largest likelihood or lower bound. Within each

    trial, the method iterates between E-step and M-step for `max_iter`

    times until the change of likelihood or lower bound is less than

    `tol`, otherwise, a `ConvergenceWarning` is raised.



    迭代终止条件： 迭代次数ｎ_init，极大似然函数或下界函数的增益小于`tol`



    Parameters

    ----------

    X : array-like, shape (n_samples, n_features)



    Returns

    -------

    self

    """

    X = _check_X(X, self.n_components)

    self._check_initial_parameters(X)



    # if we enable warm_start, we will have a unique initialisation

    do_init = not(self.warm_start and hasattr(self, 'converged_'))

    n_init = self.n_init if do_init else 1



    max_lower_bound = -np.infty

    self.converged_ = False



    random_state = check_random_state(self.random_state)



    n_samples, _ = X.shape

    # 初始化次数

    for init in range(n_init):

        self._print_verbose_msg_init_beg(init)



        # 先初始化参数

        if do_init:

            self._initialize_parameters(X, random_state)

            self.lower_bound_ = -np.infty



        # 迭代次数

        for n_iter in range(self.max_iter):

            prev_lower_bound = self.lower_bound_



            # E step求出后验概率w_j^i或是Q分布

            log_prob_norm, log_resp = self._e_step(X)

            # Ｍ step更新参数

            self._m_step(X, log_resp)

            # 求出下界函数的最大值

            self.lower_bound_ = self._compute_lower_bound(

                log_resp, log_prob_norm)

            # 下界函数的增益

            change = self.lower_bound_ - prev_lower_bound

            self._print_verbose_msg_iter_end(n_iter, change)



            # 比较下界函数增益与ｔｏｌ

            if abs(change) < self.tol:

                self.converged_ = True

                break



        self._print_verbose_msg_init_end(self.lower_bound_)



        #

        if self.lower_bound_ > max_lower_bound:

            max_lower_bound = self.lower_bound_

            best_params = self._get_parameters()

            best_n_iter = n_iter



    if not self.converged_:

        warnings.warn('Initialization %d did not converge. '

                      'Try different init parameters, '

                      'or increase max_iter, tol '

                      'or check for degenerate data.'

                      % (init + 1), ConvergenceWarning)



    self._set_parameters(best_params)

    self.n_iter_ = best_n_iter



    return self

E step中p(x|z=j)

根据高斯分布的参数计算概率,优化的计算方法。

先计算协方差矩阵的precision矩阵，并进行cholesky分解

Precision matrix 协方差矩阵的逆矩阵：https://www.statlect.com/glossary/precision-matrix

然后根据精度矩阵的cholesky分解形式,这样可以优化矩阵运算

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

def _compute_precision_cholesky(covariances, covariance_type):

    """Compute the Cholesky decomposition of the precisions.

    Parameters

    ----------

    covariances : array-like

        The covariance matrix of the current components.

        The shape depends of the covariance_type.

    covariance_type : {'full', 'tied', 'diag', 'spherical'}

        The type of precision matrices.

    Returns

    -------

    precisions_cholesky : array-like

        The cholesky decomposition of sample precisions of the current

        components. The shape depends of the covariance_type.

    """

    if covariance_type in 'full':

        n_components, n_features, _ = covariances.shape

        precisions_chol = np.empty((n_components, n_features, n_features))

        for k, covariance in enumerate(covariances):

            try:

                cov_chol = linalg.cholesky(covariance, lower=True)

            except linalg.LinAlgError:

                raise ValueError(estimate_precision_error_message)

            precisions_chol[k] = linalg.solve_triangular(cov_chol,

                                                         np.eye(n_features),

                                                         lower=True).T

    elif covariance_type == 'tied':

        _, n_features = covariances.shape

        try:

            cov_chol = linalg.cholesky(covariances, lower=True)

        except linalg.LinAlgError:

            raise ValueError(estimate_precision_error_message)

        precisions_chol = linalg.solve_triangular(cov_chol, np.eye(n_features),

                                                  lower=True).T

    else:

        if np.any(np.less_equal(covariances, 0.0)):

            raise ValueError(estimate_precision_error_message)

        precisions_chol = 1. / np.sqrt(covariances)

    return precisions_chol

计算cholesky分解的行列式

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

# Gaussian mixture probability estimators

# 根据cholesky分解计算行列式的log

def _compute_log_det_cholesky(matrix_chol, covariance_type, n_features):

    """Compute the log-det of the cholesky decomposition of matrices.

    Parameters

    ----------

    matrix_chol : 协方差矩阵的cholesky分解

    covariance_type : {'full', 'tied', 'diag', 'spherical'}

    n_features : int

        Number of features.

    Returns

    -------

    log_det_precision_chol : array-like, shape (n_components,)

        The determinant of the precision matrix for each component.

    """

    if covariance_type == 'full':

        n_components, _, _ = matrix_chol.shape

        log_det_chol = (np.sum(np.log(

            matrix_chol.reshape(

                n_components, -1)[:, ::n_features + 1]), 1))



    elif covariance_type == 'tied':

        log_det_chol = (np.sum(np.log(np.diag(matrix_chol))))



    elif covariance_type == 'diag':

        log_det_chol = (np.sum(np.log(matrix_chol), axis=1))



    else:

        log_det_chol = n_features * (np.log(matrix_chol))



    return log_det_chol

计算分子: $logp(x^{(i)}|z^{(i)}=j;\mu,\Sigma)$

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

def _estimate_log_gaussian_prob(X, means, precisions_chol, covariance_type):

    """Estimate the log Gaussian probability.

    Parameters

    ----------

    X : 样本数据(n_samples, n_features)

    means : k个n维正态分布的均值(n_components, n_features)

    precisions_chol : 精度矩阵的Cholesky分解

    covariance_type : {'full', 'tied', 'diag', 'spherical'}

    Returns

    -------

    log_prob : (n_samples, n_components)

    """

    n_samples, n_features = X.shape

    n_components, _ = means.shape

    # det(precision_chol) is half of det(precision)

    log_det = _compute_log_det_cholesky(

        precisions_chol, covariance_type, n_features)



    if covariance_type == 'full':

        log_prob = np.empty((n_samples, n_components))

        for k, (mu, prec_chol) in enumerate(zip(means, precisions_chol)):

            y = np.dot(X, prec_chol) - np.dot(mu, prec_chol)

            log_prob[:, k] = np.sum(np.square(y), axis=1)



    elif covariance_type == 'tied':

        pass



    elif covariance_type == 'diag':

        pass



    elif covariance_type == 'spherical':

        pass



    return -.5 * (n_features * np.log(2 * np.pi) + log_prob) + log_det

sklearn中实例

Although GMM are often used for clustering, we can compare the obtained clusters with the actual classes from the dataset. We initialize the means of the Gaussians with the means of the classes from the training set to make this comparison valid.

We plot predicted labels on both training and held out test data using a variety of GMM covariance types on the iris dataset. We compare GMMs with spherical, diagonal, full, and tied covariance matrices in increasing order of performance. Although one would expect full covariance to perform best in general, it is prone to overfitting on small datasets and does not generalize well to held out test data.

On the plots, train data is shown as dots, while test data is shown as crosses. The iris dataset is four-dimensional. Only the first two dimensions are shown here, and thus some points are separated in other dimensions.

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
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231

import matplotlib as mpl

import matplotlib.pyplot as plt

import numpy as np

from sklearn import datasets

from sklearn.mixture import GaussianMixture

from sklearn.model_selection import StratifiedKFold



print(__doc__)



colors = ['navy', 'turquoise', 'darkorange']



iris = datasets.load_iris()  # data, target



# Break up the dataset into non-overlapping training (75%) and testing

# (25%) sets.

skf = StratifiedKFold(n_splits=4)

# Only take the first fold.

train_index, test_index = next(iter(skf.split(iris.data, iris.target)))



X_train = iris.data[train_index]     # (111, 4)

y_train = iris.target[train_index]   # (111,)

X_test = iris.data[test_index]       # (39, 4)

y_test = iris.target[test_index]     # (39,)



n_classes = len(np.unique(y_train))



# Try GMMs using different types of covariances. 根据协方差矩阵，有4中不同的GMM模型

estimators = dict((cov_type, GaussianMixture(n_components=n_classes,

                   covariance_type=cov_type, max_iter=20, random_state=0))

                  for cov_type in ['spherical', 'diag', 'tied', 'full'])



n_estimators = len(estimators)  # 4



# figsize表示图像的尺寸（width, height in inches）

plt.figure(figsize=(3 * 5 // 2, 6))

# 图像之间的间距

plt.subplots_adjust(bottom=.01, top=0.95, hspace=.15, wspace=.05,

                    left=.01, right=.99)



# 椭圆

def make_ellipses(gmm, ax):

    for n, color in enumerate(colors):

        if gmm.covariance_type == 'full':

            covariances = gmm.covariances_[n][:2, :2]

        elif gmm.covariance_type == 'tied':

            covariances = gmm.covariances_[:2, :2]

        elif gmm.covariance_type == 'diag':

            covariances = np.diag(gmm.covariances_[n][:2])

        elif gmm.covariance_type == 'spherical':

            covariances = np.eye(gmm.means_.shape[1]) * gmm.covariances_[n]



        # 参数估计得到的混合二维高斯分布，将其用椭圆表示出来～

        v, w = np.linalg.eigh(covariances) # 返回协方差矩阵的特征值和列向量由特征矩阵构成的矩阵

        u = w[0] / np.linalg.norm(w[0]) # order=None 表示 Frobenius norm，2-norm

        angle = np.arctan2(u[1], u[0])

        angle = 180 * angle / np.pi  # 转换为角度

        v = 2. * np.sqrt(2.) * np.sqrt(v)

        ell = mpl.patches.Ellipse(gmm.means_[n, :2], v[0], v[1],

                                  180 + angle, color=color)

        ell.set_clip_box(ax.bbox)

        ell.set_alpha(0.5)

        ax.add_artist(ell)  # 增加文字



for index, (name, estimator) in enumerate(estimators.items()):

    # Since we have class labels for the training data, we can

    # initialize the GMM parameters in a supervised manner.

    # 这里因为有类标签，所以直接用真实均值来初始化GMM的均值。在无标签或者标签较少的情况下，则需要随机初始化

    estimator.means_init = np.array([X_train[y_train == i].mean(axis=0)

                                    for i in range(n_classes)])



    # Train the other parameters using the EM algorithm.

    # 用em算法来估计其他参数

    estimator.fit(X_train)



    # 画椭圆

    h = plt.subplot(2, n_estimators // 2, index + 1)

    make_ellipses(estimator, h)







    for n, color in enumerate(colors):

        data = iris.data[iris.target == n]

        # 不同的种类数据用不同的点表示

        plt.scatter(data[:, 0], data[:, 1], s=0.8, color=color,

                    label=iris.target_names[n])

    # 用xx表示测试集

    for n, color in enumerate(colors):

        data = X_test[y_test == n]

        plt.scatter(data[:, 0], data[:, 1], marker='x', color=color)



    # 训练集的准确率

    y_train_pred = estimator.predict(X_train) # 预测是选取概率最大的一类



    # 当无标签时是没有办法计算准确率的。但是这里有标签，

    # y_train_pred返回的是概率最大的索引，　y_train的元素是[0,1,2,3]中的一个

    # 因此可以求得准确率

    print(y_train_pred[:5])

    print(y_train[:5])

    train_accuracy = np.mean(y_train_pred.ravel() == y_train.ravel()) * 100

    plt.text(0.05, 0.9, 'Train accuracy: %.1f' % train_accuracy,

             transform=h.transAxes)



    y_test_pred = estimator.predict(X_test)

    test_accuracy = np.mean(y_test_pred.ravel() == y_test.ravel()) * 100

    plt.text(0.05, 0.8, 'Test accuracy: %.1f' % test_accuracy,

             transform=h.transAxes)



    plt.xticks(())

    plt.yticks(())

    plt.title(name)



plt.legend(scatterpoints=1, loc='lower right', prop=dict(size=12))





plt.show()

生成模型-高斯判别分析GDA- 高斯混合模型GMM- EM算法

在学习生成模型之前，先学习了解下密度估计和高斯混合模型。

生成学习算法(cs229,Ng)

生成算法和判别算法的区别

举个栗子：

我们要区分elephants(y=1)和dogs(y=0)

1. 对判别模型（discriminative），以logistic回归为例：

• logistic回归模型：$p(y|x;\theta)，h_{\theta}=g(\theta^Tx)$,对应的模型其中g是sigmoid函数。通过logistic回归，我们找到一条决策边界decision boundary，能够区分elephants和dogs.

• 这个学习的过程就是找到表征这个决策过程的参数 $\theta$.

2. 生成模型（generative）：

同样的我们也是要通过给定的特征x来判别其对应的类别y。但我们换个思路，就是先求p(x|y),也就是通过y来分析对应x满足的一个概率模型p(x|y)。然后在反过来看特征x，以二分类为例，p(x|y=0)和p(x|y=1)哪个概率大，那么x就属于哪一类。

• 模型：p(x|y)，在给定了样本所属的类的条件下，对样本特征建立概率模型。

• p(x|y=1)是elephants的分类特征模型

• p(x|y=0)是dogs的分类特征模型

然后通过p(x|y)来判断特征x所属的类别，根据贝叶斯公式：

$$p(y=1|x) = \dfrac{p(x|y=1)p(x)}{p(x)}$$

在给定了x的情况下p(x)是个定值，p(y)是先验分布，那么计算方法如下：

$$arg\max_yp(y|x) = arg\max_{y}\dfrac{p(x|y)p(y)}{p(x)}= arg\max_{y}p(x|y)p(y)$$

总结下就是：

• 生成模型：一般是学习一个代表目标的模型，然后通过它去搜索图像区域，然后最小化重构误差。类似于生成模型描述一个目标，然后就是模式匹配了，在图像中找到和这个模型最匹配的区域，就是目标了。

• 判别模型：以分类问题为例，然后找到目标和背景的决策边界。它不管目标是怎么描述的，那只要知道目标和背景的差别在哪，然后你给一个图像，它看它处于边界的那一边，就归为哪一类。

• 由生成模型可以得到判别模型，但由判别模型得不到生成模型。

然鹅，生成模型p(x|y)怎么得到呢？不慌，我们先了解下多维正态分布～

多维正态分布(the multivariate nirmal distribution)

关于一维正态分布怎么推导出多维正态分布的概率密度函数，可参考知乎:多维高斯分布是如何由一维发展而来的？

首先一维正态分布:

$p(x) = \dfrac{1}{\sqrt{2\pi}}exp(\dfrac{-x^2}{2})$

二维标准正态分布，就是两个独立的一维标准正态分布随机变量的联合分布：

$p(x,y) = p(x)p(y)=\dfrac{1}{2\pi}exp(-\dfrac{x^2+y^2}{2})$

把两个随机变量组合成一个随机向量：$v=[x\quad y]^T$

$p(v)=\dfrac{1}{2\pi}exp(-\dfrac{1}{2}v^Tv)\quad$ 显然x,y相互独立的话，就是上面的二维标准正态分布公式～

然后从标准正态分布推广到一般正态分布，通过一个线性变化：$v=A(x-\mu)$

$p(x)=\dfrac{|A|}{2\pi}exp[-\dfrac{1}{2}(x-\mu)^TA^TA(x-\mu)]$

注意前面的系数多了一个|A|（A的行列式）。

可以证明这个分布的均值为$\mu$，协方差为$(A^TA)^{-1}$。记$\Sigma = (A^TA)^{-1}$，那就有

$$p(\mathbf{x}) = \frac{1}{2\pi|\Sigma|^{1/2}} \exp \left[ -\frac{1}{2} (\mathbf{x} - \mu) ^T \Sigma^{-1} (\mathbf{x} - \mu) \right]$$

推广到n维：

$$p(\mathbf{x}) = \frac{1}{(2\pi)^{n/2}|\Sigma|^{1/2}} \exp \left[ -\frac{1}{2} (\mathbf{x} - \mu) ^T \Sigma^{-1} (\mathbf{x} - \mu) \right]$$

需要注意的是：这里的二维、n维到底指的是什么？

• 以飞机检测的数据点为例，假设它由heat和time决定，那么这就是个二维正态分布，数据点的生成所处的位置由其概率决定，也就是$p(\mathbf{x})$

• 如果这个数据有n个特征，那么其分布就是n维正态分布。

• 之前一直理解的是，n维正态分布是两个向量巴拉巴拉。。好像一直没搞懂。。

再顺便了解下协方差矩阵吧～

关于协方差矩阵，参考blog

对多维随机变量$X=[X_1,X_2,…,X_n]^T$，我们往往需要计算各维度之间的协方差，这样协方差就组成了一个n×n的矩阵，称为协方差矩阵。协方差矩阵是一个对角矩阵，对角线上的元素是各维度上随机变量的方差,非对角线元素是维度之间的协方差。 我们定义协方差为$\Sigma$, 矩阵内的元素$\Sigma_{ij}$为:

$$\Sigma_{ij} = cov(X_i,X_j) = E[(X_i - E(X_i)) (X_j - E(X_j))]$$

则协方差矩阵为:

$$
\Sigma = E[(X-E(X)) (X-E(X))^T] = \left[

\begin{array}{cccc}

cov(X_1,X_1) & cov(X_1,X_2) & \cdots & cov(X_1,X_n) \

cov(X_2,X_1) & cov(X_2,X_2) & \cdots &cov(X_2,X_n) \

\vdots & \vdots& \vdots & \vdots \

cov(X_n,X_1) & cov(X_n,X_2,)&\cdots& cov(X_n,X_n)

\end{array}

\right]
$$
如果X~$N(\mu,\Sigma)$,则$Cov(X)=\Sigma$

可以这么理解协方差，对于n维随机变量X，第一维是体重$X_1$，第二维是颜值$X_2$，显然这两个维度是有一定联系的，就用$cov(X_1,X_2)$来表征，这个值越小，代表他们越相似。协方差怎么求，假设有m个样本，那么所有的样本的第一维就构成$X_1$…不要把$X_1$和样本搞混淆了。

了解了多维正态分布和协方差，我们再回到生成模型p(x|y)。。其实我们就是假设对于n维特征，p(x|y)是n维正态分布～怎么理解呢，下面就说！

高斯判别分析模型The Gaussian Discriminant Analysis model

高斯判别模型就是：假设p(x|y)是一个多维正态分布，为什么可以这么假设呢？因为对于给定y的条件下对应的特征x都是用来描述这一类y的，比如特征是n维的，第一维描述身高，一般都是满足正态分布的吧，第二维描述体重，也可认为是正态分布吧～

则生成模型：

y ~ Bernoulli($\phi)$ 伯努利分布，又称两点分布，0-1分布

x|y=0 ~ $N(u_0,\Sigma)$

x|y=1 ~ $N(u_1,\Sigma)$

• 这里可以看作是一个二分类，y=0和y=1,可以看作是伯努利分布，则$p(y)=\phi^y(1-\phi)^{1-y}$，要学的参数之一: $\phi=p(y=1)$，试想如果是多分类呢，那么要学习的参数就有$\phi_1,\phi_2,….\phi_k$

• 其中类别对应的特征x|y=0,x|y=1服从正态分布。怎么理解呢？就是既然你们都是一类人，那么你们的身高啊，体重啊等等应该满足正态分布。。有几维特征就满足几维正态分布

• 这里x是n维特征，身高，体重，颜值…balabala，所以x|y=0满足n维正态分布～x|y=1也是啦，只不过对于不同的类，对应n维特征的均值不一样，奇怪为什么协方差矩阵是一样的？？这里是将它特殊化了，后面会讲的一般性的em算法就不是这样的了

• 每个分类对应的n维特征的分布显然不是独立的，比如体重和颜值还是有关系的吧～他们的协方差，方差就统统都在$\Sigma$协方差矩阵里面了

$$p(x|y=0) = \frac{1}{(2\pi)^{n/2}|\Sigma|^{1/2}} \exp \left[ -\frac{1}{2} (x - \mu_0) ^T \Sigma^{-1} (x - \mu_0) \right]$$

$$p(x|y=1) = \frac{1}{(2\pi)^{n/2}|\Sigma|^{1/2}} \exp \left[ -\frac{1}{2} (x - \mu_1) ^T \Sigma^{-1} (x - \mu_1) \right]$$

这样，模型中我们要学习的参数有$\phi,\Sigma, \mu_0,\mu_1$，对于训练数据，就是观测到的数据x,y，既然他们出现了，那么他们的联合概率，也就是似然函数$\prod_{i=1}^mp(x,y)$就要最大～其对数似然log-likelihood：

$$\begin{equation}\begin{aligned}

L(\phi,\Sigma, \mu_0,\mu_1) &= log\prod_{i=1}^mp(x^{(i)},y^{(i)};\phi,\Sigma, \mu_0,\mu_1$) \

&= log\prod_{i=1}^mp(x^{(i)}|y^{(i)};\phi,\Sigma, \mu_0,\mu_1$) p(y^{(i)};\phi)\

\end{aligned}\end{equation}\label{eq2}$$

其中$p(y^{(i)};\phi)$是已知的，也就是先验概率(class priors)，$p(x^{(i)}|y^{(i)})$就是上面推导的～代入后，分别对参数求导即可：

在回过头来看这些公式，

• $\phi$很好理解，就是样本中正分类的概率。

• $\mu_0$就是负分类中x对应的均值

• $\mu_1$就是正分类中x对应的均值

• $\Sigma$就是$(x-\mu_1)$和$x-\mu_2$的协方差矩阵

然后通过p(x|y=0),p(x|y=1)即可对需要预测的x求出对应的概率，然后做出判别了。这样看来，如果直接对x|y=1,和x|y=0做出了正态分布的猜测，就可以直接写出来了。只不过，我们用极大似然估计重新推导了一遍。

高斯混合模型GMM

GMM

前面GDA是有标签的，也算是有监督学习。而在没有标签的情况下呢，就是无监督学习了，虽然我们无法给出x所属的类叫啥，但是我们可以判断出哪些x是同一类，以及样本中总共有多少类（虽然这个类数嘛。。类似于k-means的类数，可根据交叉验证选择）。

其实和GDA非常相似，不过这里没有了类标签，只有一堆样本特征，${x^{(1)},x^{(2)},…,x^{(m)}}$,

我们不知道这些样本属于几个类别，也不知道有哪些类了。但虽然不知道，我们确定他们是存在的，只是看不见而已。我们可以假设存在k类，${z^{(1)},z^{(2)},…,z^{(k)}}$,看不见的，我们就叫它们隐藏随机变量(latent random variable)，

这样一来，就训练样本就可以用这样的联合分概率模型表示了，$p(x^{(i)},z^{(i)})=p(x^{(i)}|z^{(i)})p(z^{(i)})$

• 同GDA不一样的是，这里是多分类，可假定$z^{(i)}\sim Multinomial(\phi)$，多项式分布（二项分布的拓展～），那么$p(z^{(i)})=\phi_j$

• 同GDA相同的是，对于每一个类别，其对应的样本满足n维正态分布，也就是：$x^{(i)}|z^{(i)}=j\sim N(\mu_j,\Sigma_j)$,但注意哦，这里每个高斯分布使用了不同的协方差矩阵$\Sigma_j$

$$p(x|z^{(1)}) = \frac{1}{(2\pi)^{n/2}|\Sigma_0|^{1/2}} \exp \left[ -\frac{1}{2} (x - \mu_0) ^T \Sigma_0^{-1} (x - \mu_0) \right]$$

$$p(x|z^{(2)}) = \frac{1}{(2\pi)^{n/2}|\Sigma_1|^{1/2}} \exp \left[ -\frac{1}{2} (x - \mu_1) ^T \Sigma_1^{-1} (x - \mu_1) \right]$$

$$….$$

$$p(x|z^{(k)}) = \frac{1}{(2\pi)^{n/2}|\Sigma_k|^{1/2}} \exp \left[ -\frac{1}{2} (x - \mu_k) ^T \Sigma_k^{-1} (x - \mu_k) \right]$$

然后带入到训练样本的对数似然（log-likelihood）:

$$L(\phi,\mu,\Sigma)=\sum_{i=1}^{m}logp(x^{(i)};\phi,\mu,\Sigma)$$

$$L(\phi,\mu,\Sigma)=\sum_{i=1}^{m}log\sum_{z^{(i)}=1}^{k}p(x^{(i)}|z^{(i)};\mu,\Sigma) p(z^{(i)};\phi)\$$

这里需要注意下标：对于类别有k类，第一个求和符号是对第i个样本在k个类别上的联合概率，第二个求和符号是m个样本的联合概率。

我们可以注意到，如果我们知道$z^{(i)}$,那么这个似然函数求极大值就很容易了，类似于高斯判别分析，这里的$z^{(i)}$相当于标签，分别对参数求导可得：

其中的参数:

• $1{z^{(i)}=j}$表示第i个样本为j类时，这个值就为１，那么$\phi_j=\frac{1}{m}\sum_{i=1}^m1{z^{(i)}=j}$表示样本中类别为j的概率

• 其中$p(z^{(i)};\phi)$是根据伯努利分布得到的，在GDA中$p(y|\phi)$是已知的频率概率。

So $z^{(i)}$ 到底有多少个分类？每个类别的概率是多少？譬如上式中 $\sum_{i=1}^{m}1{z^{(i)}=j}$ 这个没法求对吧～它是隐藏变量！所以还是按照这个方法是求不出来的～

这个时候EM算法就登场了～～～

用EM算法求解GMM模型

上面也提到了，如果$z^({i})$是已知的话，那么$\phi_j=\frac{1}{m}\sum_{i=1}^m1{z^{(i)}=j}$表示类别j的概率$p(z^{(i)}=j)$也就已知了，但是呢？我们不知道。。所以我们要猜测$p(z^{(i)}=j)$这个值，也就是EM算法的第一步：

Repeat until convergence 迭代直到收敛:{

(E-step):for each i,j,set:

$w_j^{(i)}:=p(z^{(i)}=j|x^{(i)};\phi,\mu,\Sigma)$

$w_j^{(i)}$什么意思呢?就是对于i样本,它是j类的后验概率。在GDA里面，x_i的类别是确定的，在GMM里面呢？不知道它的类别，所以只能假设k类都有可能，它是j类别的概率就是$w_j^{(i)}$，它仅仅取决于$\phi_j$,而在GMM里面，它取决于$\phi_j,\mu_j,\Sigma_j$，实际上$w_j^{(i)}$的值，就包含了两个我们在GMM所做的假设，多项式分布和正态分布。

The values $w_j$ calculated in the E-step represent our “soft” guesses for

the values of $z^{(i)}$ .

The term “soft” refers to our guesses being probabilities and taking values in [0, 1]; in

contrast, a “hard” guess is one that represents a single best guess (such as taking values

in {0, 1} or {1, . . . , k}).

硬猜测是k均值聚类，GMM是软猜测。

这样一来，参数更新就可以这样写了，也就是EM算法的第二步：

(M-step) Updata the parameters:

然后对似然函数求导，后面会详细介绍

$$\phi_j:=\frac{1}{m}\sum_{i=1}^mw_j^{(i)}$$

$$\mu_j:=\dfrac{\sum_{i=1}^mw_j^{(i)}x^{(i)}}{\sum_{i=1}^mw_j^{(i)}}$$

$$\Sigma_j:=\dfrac{\sum_{i=1}^mw_j^{(i)}(x^{(i)}-\mu_j)(x^{(i)}-\mu_j)^T}{\sum_{i=1}^mw_j^{(i)}}$$

｝

训练过程的理解可参考blog

$w_j^{(i)}表示第i个样本为ｊ类别的概率，而\phi_j$表示m个样本中j类别的概率，$\mu_j,\Sigma_j$分别表示j类别对应的n维高斯分布的期望和协方差矩阵

所以，求出$w_j^{(i)}$，一切就都解决了吧？对于后验概率$p(z^{(i)}=j|x^{(i)})$可以根据Bayes公式：

$$p(z^{(i)}=j|x^{(i)};\phi,\mu,\Sigma)=\dfrac{p(x^{(i)}|z^{(i)}=j;\mu,\Sigma)p(z^{(i)}=j;\phi)}{\sum_{l=1}^kp(x^{(i)}|z^{(i)}=l;\mu,\Sigma)p(z^{(i)}=l;\phi)}$$

• 其中先验概率$p(x^{(i)}|z^{(i)}=j;\mu,\Sigma)$可以根据高斯分布的密度函数来求，$z^{(i)}=j\sim N(\mu_j,\Sigma_j)$

• 类先验(class priors)$p(z^{(i)}=j;\phi)$可以取决于多项式分布中j类的概率$\phi_j$

The EM-algorithm is also reminiscent of the K-means clustering algorithm, except that instead of the “hard” cluster assignments $c^(i)$, we instead have the “soft” assignments $w_j$ . Similar to K-means, it is also susceptible to local optima, so reinitializing at several different initial parameters may be a good idea.

EM算法使我们联想起了k-means,区别在于k-means的聚类是通过欧氏距离c(i)来定义的，而EM是通过$w_j$probabilities来分类的。同k-means一样，这里的EM算法也是局部优化，因此最好采用不同的方式初始化～

convergence?

我们知道k-means一定是收敛的，虽然结果不一定是全局最优解，但它总能达到一个最优解。但是EM算法呢，也是收敛的。

The EM algorithm

前面我们讲的是基于高斯混合模型的EM算法，但一定所有的类别都是高斯分布吗？还有卡方分布，泊松分布等等呢，接下来我们就将讨论EM算法的一般性。

在学习一般性的EM算法前，先了解下Jensen’s inequality

Jensen’s inequality

如果函数$f$，其二阶导恒大与等于0 $(f^{‘’}\ge 0)$，则它是凸函数f(convec function)。

如果凸函数的输入是向量vector-valued inputs，那么它的海森矩阵(hessian)H是半正定的。Jensen’s 不等式：

Let f be a convex function, and let X be a random variable.

Then:

$$E[f (X)] ≥ f (EX).$$

Moreover, if f is strictly convex, then $E[f (X)] = f (EX)$ holds true if and

only if $X = E[X]$ with probability 1 (i.e., if X is a constant).

举个栗子来解释jensen不等式：

假设输入随机变量X是一维的哈，然后Ｘ取a,b的概率都是0.5,那么

$$EX=(a+b)/2,f(EX)=f(\dfrac{a+b}{2})$$,$$E[f(X)]=\dfrac{f(a)+f(b)}{2}$$

因为是凸函数，所以 $f(EX)\le E[f(X)]$

同理，如果是凹函数(concave function),那么不等式方向相反$f(EX)\ge E[f(X)]$。后面EM算法里面就要用到log(X)，log(x)就是个典型的凹函数～

The EM algorithm

首先，问题是：我们要基于给定的m个训练样本${x^{(1)},x^{(2)},…,x^{(m)}}$来进行密度估计～

像前面一样，创建一个参数模型p(x,z)来最大化训练样本的对数似然：

$$L(\theta)=\sum_{i=1}^mlogp(x;\theta)$$

$$L(\theta)=\sum_{i=1}^mlog\sum_zp(x,z;\theta)$$

一般性就是把前面特殊化的假设去掉，没有了正态分布和多项式分布。

可以看到，$z^{(i)}$是隐藏的随机变量(latent random variable),关于参数$\theta$的最大似然估计就很难计算了。

解释下公式中的推导：

• 这里是针对样本i来说，对于样本i，它可能是$z^1,z^2,…,z^k$都有可能，但他们的probability之和为１，也就是

$\sum_zQ_i(z)=1$

• (2)到(3)的推导：可以将

$\dfrac{p(x^{(i)},z^{(i)};\theta)}{Q_i(z^{(i)})}$

看做随机变量Ｘ,那么（２）式中的后半部分 $log\sum_{z^{(i)})}[\dfrac{p(x^{(i)},z^{(i)};\theta)}{Q_i(z^{(i)})}]$就是log(EX)了，$logx$是一个凹函数，则其大于$E[log(x)]$

EM迭代过程(重点):

• （1）根据上式可以看做$L(\theta)\ge J(Q,\theta)$.两边都是关于$\theta$的函数，那么将$\theta$固定，调整Q在一定条件下能使等式成立。

• （2）然后固定Q,调整$\theta^t$到$\theta^{t+1}$找到下界函数的最大值$J(Q,\theta^{t+1})$.显然在当前Q的条件下，$L(\theta^{t+1})\ne J(Q,\theta^{t+1})$,那么根据Jensen不等式，$L(\theta_{t+1})>J(Q,\theta^{t+1})=L(\theta^{t})$,也就是说找到了使得对数似然L更大的$\theta$.这不就是我们的目的吗？！

• 然后迭代循环(1)(2)步骤，直到在调整$\theta$时，下界函数$J(Q,\theta)$不在增加，即小于某个阈值。

看下Ng画的图：

任意初始化$\theta$和Q,然后找下界函数和$l(\theta)$交接的点，这就是EM算法的第一步：

我们要让不等式相等,即Jensen’s inequality中的随机变量取值是一个常量，看(2)式：

$$\dfrac{p(x^{(i)},z^{(i)};\theta)}{Q_i(z^{(i)})}=c$$

对左边分子分母同时对z求类和：

$$\dfrac{\sum_zp(x^{(i)},z^{(i)};\theta)}{\sum_zQ_i(z^{(i)})}=c$$

根据$\sum_zQ_i(z)=1$：

$$\sum_zp(x^{(i)},z^{(i)};\theta)=c$$

带回去可得：

$$Q_i(z^{(i)})=\dfrac{p(x^{(i)},z^{(i)};\theta)}{\sum_zp(x^{(i)},z^{(i)};\theta)}$$

$$Q_i(z^{(i)})=\dfrac{p(x^{(i)},z^{(i)};\theta)}{p(x^{(i)};\theta)}$$

$$Q_i(z^{(i)})=p(z^{(i)}|x^{(i)};\theta)$$

EM总结下来：

Repeat until convergence {

(E-step)

For each i,找到下界函数, set:

$$Q_i(z^{(i)}):=p(z^{(i)}|x^{(i)};\theta)$$

(M-step)找到下界凹函数的最大值,也就是(3)式 Set:

$$\theta:=arg\max_{\theta}\sum_i^m\sum_{z^{(i)}}^kQ_i(z^{(i)})log\dfrac{p(x^{(i)},z^{(i)};\theta)}{Q_i(z^{(i)})}$$

}

要理解的是：

EM算法只是一种计算方式，对于上式中的$Q_i(z^{(i)})=p(z^{(i)}|x^{(i)};\theta)$我们还是要根据假设来求得，比如GMM中的多类高斯分布。然后带回到对数似然中，通过求导得到参数估计。我们费尽心机证明EM算法收敛，只是为了去证明这样去求似然函数的极大值是可行的，然后应用到类似于GMM，HMM中。

training and will converge?

首先说是否收敛，答案是肯定收敛的。。懒得输公式了。。直接贴图吧，这个比较好理解：

上面写这么多，其实就是证明$L(\theta_{t+1})>L(\theta_t)$.

Mixture of Gaussians revisited

我们知道了em算法是一种计算方式，用来解决含有隐变量似然对数很难求的问题，那么我们把它运用到GMM中。

E step:

$w_j^{(i)}:=p(z^{(i)}=j|x^{(i)};\phi,\mu,\Sigma)$

$$w_j^{(i)}:=p(z^{(i)}=j|x^{(i)};\phi,\mu,\Sigma)=\dfrac{p(x^{(i)}|z^{(i)}=j;\mu,\Sigma)p(z^{(i)}=j;\phi)}{\sum_{l=1}^kp(x^{(i)}|z^{(i)}=l;\mu,\Sigma)p(z^{(i)}=l;\phi)}$$

• 其中先验概率$p(x^{(i)}|z^{(i)}=j;\mu,\Sigma)$可以根据高斯分布的密度函数来求，$z^{(i)}=j\sim N(\mu_j,\Sigma_j)$

• 类先验(class priors)$p(z^{(i)}=j;\phi)$可以取决于多项式分布中j类的概率$\phi_j$

这样我们就完成了对$w_j^{(i)}$的soft ‘guess’，也就是E step.

M step:

然后对参数求导：

详细推导过程，参考cs229-notes8

我们在整体回顾一下整个过程，所谓的E step就是找到$Q_i(z^{j}),w_i^j$（在一定假设下是可以通过bayes公式求得的），使得下界函数与log函数相等，也就是Jensen取等号时。然后是M step就是在Q的条件下找到下界函数最大值，也就是对参数求导，导数为0的地方。

然后在根据求得的参数，再求Q，再带入求导。。。迭代直到收敛。