机器学习-不均衡数据问题

对于机器学习中不均衡数据的问题是很常见的,之前在三星实习也遇到过,但当时对指标性能要求并不高,差不多也就行了。。但这次参加 kaggle 比赛,0.1% 就有很大差距了,所以还是要好好搞搞 imbalanced data problem.

很多时候,数据不均衡不仅仅数据收集和整理的问题,而是数据本身预期的就是这样,比如表征欺诈性交易(characterize fraudulent transactions),大多数交易都是不具有欺诈性的,只有极少数的 具有欺诈性的(Fraud), 或者 kaggle 比赛中的 insincere-questions, 大多数问题都是正常的,只有极少数的 insincere-questions.

More data

显然,这是最直接的。在实际中,也是最可行的。但打比赛时,数据是给定的,就需要用到 resampling.

Metrics

改变评价指标来选择更优的模型。

  • accurayc 显然这个对不均衡数据是不合适的

  • Precision: (tp/(tp+fp))

  • Recall: (tp/(tp+fn))

  • F1 Score (or F-score): precision 和 recall 的一种权衡

  • Kappa: Cohen’s kappa

  • ROC 曲线:这些在之前的笔记中都有介绍 机器学习-常用指标总结

  • PRAUC 损失函数? ai challenger 比赛中看到的

Resampling

  • 过采样 over-sampling

  • 欠采样 under-sampling

实践中可以两种都尝试,并且能做集成增强。

一些采样策略:

  • 数据量很大(上万或者十万?),建议欠采样,数据量较小,建议过采样

  • 可以尝试随机采样,也可以尝试分层采样(划分好)

  • 尝试不同的比例,并不一定要最后是 1:1

Generate Synthetic Samples

SMOTE(Synthetic Minority Over-sampling Technique): 不是简单的 copy, 而是选取两个或更多的相似的样本(根据距离),然后随机扰动一个样本的属性(某一维特征吧),扰动值在它与相邻样本的差异之间。

paper: SMOTE: Synthetic Minority Over-sampling Technique

python tool: UnbalancedDataset

Try Penalized Models

在训练过程中,给类别较少的一类增加惩罚项,通常是正则化.这有利于模型能注重 minority class.

对类和权重的惩罚对不同的算法不太一样,所有有专门的版本的算法:penalized-SVM 和 penalized-LDA.

也有通用的惩罚模型,适用于不同的分类器 CostSensitiveClassifier

惩罚模型提供了另外一种方法来 “balance” 模型,但设置惩罚矩阵是很复杂的,需要很多尝试。

Try a Different Perspective

对于之前说过的欺诈性检测和kaggle insincere 问题的发现,从另外一个角度看,也能看做是异常检测(anomaly detection)和 变异检测(change detection).

类似于 open set recognition 或 out-of-distribution 问题,样本数极少的那个类别 minior class 可以看做是 outliers class.

Anomaly detection is the detection of rare events. This might be a machine malfunction indicated through its vibrations or a malicious activity by a program indicated by it’s sequence of system calls. The events are rare and when compared to normal operation.

Change detection is similar to anomaly detection except rather than looking for an anomaly it is looking for a change or difference. This might be a change in behavior of a user as observed by usage patterns or bank transactions.

机器学习-常用指标总结

参考链接

精确率 Precision, 召回率 Recall 和 F1 值

对于数据不均衡时,使用accuracy是不准确的。

举个栗子:

a million tweet: 999,900条不是关于pie的,只有100条是关于pie的

对于一个stupid分类器,他认为所有的tweet都是跟pie无关的,那么它的准确率是99.99%!但这个分类器显然不是我们想要的,因而accuracy不是一个好的metric,当它目标是rare,或是complete unbalanced.

引入另外两个指标:

  • 精度 precision: 是检索出来的条目(比如:文档、网页等)有多少是准确的。精度是检索出相关文档数与检索出的文档总数的比率,衡量的是检索系统的查准率。
  • 召回 call: 召回率是指检索出的相关文档数和文档库中所有的相关文档数的比率,衡量的是检索系统的查全率.

需要先确定一个正分类: 这里需要检索出来的是与 pie 无关的,也就是与 pie 无关的是正分类。那么精度就是分类器检测出来的正分类中 gold pos 所占的比例,那么 precision = 0/(0+0)

这里需要检索的是与 pie 无关的,检索出来的 true pos 占样本中 gold pos 的比例,那么 recall = 0/(100+0) = 0.

总结下来,precision 和 recall 都是以 true pos 作为分子,precision 是以分类器预测出来的 pos(true pos + false neg) 作为分母,所以是差准率. recall 则是以总的 gold pos(true pos + false neg) 作为分母,所以是查全率。

当然希望检索结果Precision越高越好,同时Recall也越高越好,但事实上这两者在某些情况下有矛盾的。比如极端情况下,我们只搜索出了一个结果,且是准确的,那么Precision就是100%,但是Recall就很低;而如果我们把所有结果都返回,那么比如Recall是100%,但是Precision就会很低。因此在不同的场合中需要自己判断希望Precision比较高或是Recall比较高。如果是做实验研究,可以绘制Precision-Recall曲线来帮助分析。

F-measure

$$F_{\beta}=\dfrac{(\beta^2+1)PR}{\beta^2P+R}$$

当 $\beta>1$时,Recall的比重更大;当 $\beta<1$时,precision的比重更大。使用最多的是 $\beta=1$,也就是 $F_{\beta=1}, F_1$. $\beta$ 的的取值取决于实际应用。

$$F_1 = \dfrac{2PR}{P+R}$$

F-measure 是 precision 和 recall 的 **加权调和平均值(weighted harmonic mean)**。

调和平均值是倒数的算术平均值的倒数。

为什么要使用调和平均值呢?因为它是一个更 保守的度量(conservative metric). 相比直接计算 P 和 R 的平均值, F-measure的值更看重两者中的较小值。

ROC曲线和AUC

考虑一个二分问题,即将实例分成正类(positive)或负类(negative)。

  • TPR, 真正类率(true positive rate ,TPR),: 如果一个实例是正类并且也被 预测成正类,即为真正类(True positive),真正类率是指分类器所识别出来的 正实例占所有正实例的比例。就是 Recall 吧~ TPR = TP / (TP + FN)

  • FPR, 负正类率: 分类器错认为正类的负实例占所有负实例的比例,FPR = FP / (FP + TN)

  • TNR, 真负类率: 分类器认为负类的负实例占所有负实例的比例,也就是负类的 Recall 吧~ TNR = TN /(FP + TN) = 1 - FPR

为什么要引入 ROC 曲线

  • Motivation1:在一个二分类模型中,对于所得到的连续结果,假设已确定一个阀值,比如说 0.6,大于这个值的实例划归为正类,小于这个值则划到负类中。如果减小阀值,减到0.5,固然能识别出更多的正类,也就是提高了识别出的正例占所有正例 的比类,即TPR,但同时也将更多的负实例当作了正实例,即提高了FPR。为了形象化这一变化,引入ROC,ROC曲线可以用于评价一个分类器。
  • Motivation2:在类不平衡的情况下,如正样本90个,负样本10个,直接把所有样本分类为正样本,得到识别率为90%。但这显然是没有意义的。单纯根据Precision和Recall来衡量算法的优劣已经不能表征这种病态问题。

ROC 曲线

前面已经说道,对于数据不均衡的情况下,precision 和 recall 不足以表征这类问题,就比如上面的例子中,把找出不关于 pie 的 tweet看作是正类,那么它的 精度和召回率 都很高。所以引入 ROC

因为我们要关注的是正类,所以关注指标是 真正类率 TPR 和 负正类率 FPR,真正类率越高越好,负正类率越低越好。但显然这两者之间是矛盾的,与分类器的阈值有关,ROC 就是用来表征阈值与 TPR 和 FPR 之间的关系曲线。

ROC(Receiver Operating Characteristic)翻译为“接受者操作曲线”。曲线由两个变量1-specificity 和 Sensitivity绘制. 1-specificity=FPR,即负正类率。Sensitivity即是真正类率,TPR(True positive rate),反映了正类覆盖程度。

为了更好地理解ROC曲线,我们使用具体的实例来说明:

如在医学诊断中,判断有病的样本。那么尽量把有病的揪出来是主要任务,也就是第一个指标TPR,要越高越好。而把没病的样本误诊为有病的,也就是第二个指标FPR,要越低越好。

不难发现,这两个指标之间是相互制约的。如果某个医生对于有病的症状比较敏感,稍微的小症状都判断为有病,那么他的第一个指标应该会很高,但是第二个指标也就相应地变高。最极端的情况下,他把所有的样本都看做有病,那么第一个指标达到1,第二个指标也为1。

我们以FPR为横轴,TPR为纵轴,得到如下ROC空间。

参考链接中 ROC 曲线的解释: http://www.cnblogs.com/maybe2030/p/5375175.html#_label0

上面的图表示,分类器对于 pos 和 neg 的分别是有个阈值的,阈值很高,也就是 A 处,显然它的 TPR 不会很高,阈值越低,TPR 越高。 但是阈值很低的话,预测为正类的负实例也就越多, FPR 也会越高。

曲线距离左上角越近,证明分类器效果越好。我们用一个标量 AUC 来量化这个分类效果。

怎么得到 ROC 曲线

我们知道对一个二值分类器,其预测出来的正类的 score 是一个概率。当一个样本为正类的概率大于 threshold时,我们判别它为正类。所以不同的 threshold,其对应的 TPR 和 FPR 的值也就不一样,这样就得到了 ROC 曲线。

详细可参考:ROC和AUC介绍以及如何计算AUC 非常清楚!!!

AUC

AUC值为ROC曲线所覆盖的区域面积,显然,AUC越大,分类器分类效果越好。

  • AUC = 1,是完美分类器,采用这个预测模型时,不管设定什么阈值都能得出完美预测。绝大多数预测的场合,不存在完美分类器。

  • 0.5 < AUC < 1,优于随机猜测。这个分类器(模型)妥善设定阈值的话,能有预测价值。

  • AUC = 0.5,跟随机猜测一样(例:丢铜板),模型没有预测价值。

  • AUC < 0.5,比随机猜测还差;但只要总是反预测而行,就优于随机猜测。

AUC的物理意义:假设分类器的输出是样本属于正类的socre(置信度),则AUC的物理意义为,任取一对(正、负)样本,正样本的score大于负样本的score的概率。

怎么计算 AUC

AUC(Area Under Curve)被定义为ROC曲线下的面积,显然这个面积的数值不会大于1。又由于ROC曲线一般都处于y=x这条直线的上方,所以AUC的取值范围在0.5和1之间。使用AUC值作为评价标准是因为很多时候ROC曲线并不能清晰的说明哪个分类器的效果更好,而作为一个数值,对应AUC更大的分类器效果更好。

代码实现: sklearn.metrics.roc_auc_score

置信度和置信区间

这里好像跟置信度和置信区间没啥关系。。。但是了解下也没事

置信区间

再来理解置信度,首先要理解 95%置信区间 和 置信度。

知乎:95%置信区间 我的理解就是,首先总体均值是固定的,但只有上帝知道。我们就要用样本去估计总体均值。一次次抽样会得到很多样本均值,但是我们无法判断哪个均值最好,最接近总体均值。于是,我们构造区间来看这个区间是否含有总体均值。

怎么构造区间:

通过一次次抽样得到的样本均值:

$$M=\dfrac{X_1 + X_2+…+X_n}{n}$$

根据大数定律:

$$M\sim N(\mu,\dfrac{\sigma^2}{n})$$

通过查表 标准正太分布表可知,这样可以计算出置信区间,(如果方差为止,则用样本方差代替)

我们以 $1.96\dfrac{\sigma }{\sqrt{n}}$ 为半径做区间,就构造出了 $95%$ 置信区间。按这样去构造的100个区间,其中大约会有95个会包含 $\mu$ :

那么,只有一个问题了,我们不知道、并且永远都不会知道真实的 $\mu$ 是多少:

我们就只有用 $\hat{\mu }$ 来代替 $\mu$ :

$$P(\hat \mu-1.96\dfrac{\sigma}{n}\le M\le\hat \mu+1.96\dfrac{\sigma}{n}) = 0.95$$

这样可以得到置信区间了。如果抽样100次,对应也就有100个置信区间,那么其中含有总体均值的概率约为 95%.

置信度

样本数目不变的情况下,做一百次试验,有95个置信区间包含了总体真值。置信度为95%

代码实现高斯混合模型

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根据高斯分布的参数计算概率,涉及到协方差矩阵,要优化计算,很复杂,放在最后说。先把整个流程走完。

  1. 后者是每一类高斯分布所占的权重,也就是$\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.



迭代终止条件: 迭代次数n_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)

# M 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)



# 比较下界函数增益与tol

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()

机器学习-生成模型到高斯判别分析再到GMM和EM算法

生成模型-高斯判别分析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$.

  1. 生成模型(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类时,这个值就为1,那么$\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个样本为j类别的概率,而\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是一维的哈,然后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之和为1,也就是

$\sum_zQ_i(z)=1$

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

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

看做随机变量X,那么(2)式中的后半部分 $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画的图:

em1.png

任意初始化$\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?

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

em3.png

上面写这么多,其实就是证明$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,再带入求导。。。迭代直到收敛。