sklearn源码阅读，用em算法计算高斯混合模型GMM
代码实现高斯混合模型
参考这篇博客Regularized Gaussian Covariance Estimation非常值得一读，同事这篇博客很深入的讲了协方差怎么求的问题，在前文中我也有提到～但我解释的很low。。
代码直接就看sklearn里面的源码吧～网上很多不靠谱。。。 github源码
类初始化
1 | class GaussianMixture(BaseMixture): |
先看初始化构造函数，参数是真的多。。。
- 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_paramsmethod.
先验权重初始化，对应的就是隐藏变量有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 | def _e_step(self, X): |
那么如何求\(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# 计算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()
- 其中前者是高斯分布概率的对数,根据均值，协方差矩阵的cholesky分解可求得。
1
2
3def _estimate_log_prob(self, X):
return _estimate_log_gaussian_prob(
X, self.means_, self.precisions_cholesky_, self.covariance_type)
这个函数，1
2
3
4
5
6
2. 后者是每一类高斯分布所占的权重，也就是$\phi_j$
```python
def _estimate_log_weights(self):
# 刚开始是初始值，后面随着m step而更新
return np.log(self.weights_)
再计算\(w_j^i\)
1 | def _estimate_log_prob_resp(self, X): |
M setp
1 | def _m_step(self, X, log_resp): |
具体怎么求，就是根据前面推导的公式了。根据前面的公式分别求对应的估计参数： \[\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 | def _estimate_gaussian_covariances_full(resp, X, nk, means, reg_covar): |
然后是正态分布的参数估计\(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
23def _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 | def fit(self, X, y=None): |
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
39def _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 | # Gaussian mixture probability estimators |
计算分子: \(logp(x^{(i)}|z^{(i)}=j;\mu,\Sigma)\)
1 | def _estimate_log_gaussian_prob(X, means, precisions_chol, covariance_type): |
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 | import matplotlib as mpl |
