机器学习-生成模型到高斯判别分析再到GMM和EM算法
生成模型-高斯判别分析GDA- 高斯混合模型GMM- EM算法
在学习生成模型之前，先学习了解下密度估计和高斯混合模型。
生成学习算法(cs229,Ng)
生成算法和判别算法的区别
举个栗子：
我们要区分elephants(y=1)和dogs(y=0)
- 对判别模型（discriminative），以logistic回归为例：
logistic回归模型：$p(y|x;\theta)，h_{\theta}=g(\theta^Tx)$,对应的模型其中g是sigmoid函数。通过logistic回归，我们找到一条决策边界decision boundary，能够区分elephants和dogs.
这个学习的过程就是找到表征这个决策过程的参数 $\theta$.
- 生成模型（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，再带入求导。。。迭代直到收敛。