混合高斯分布与 EM 算法

概率论中的 Jensen 不等式

对于 Jensen 不等式，通常情况下是这样的：对于 (f^{prime prime}(x) geq 0) 也就是对于凸函数而言，这个可以用中值定理来证明，同时这个公式还可以用于证明多项式的均值不等式与调和不等式，这里就不给出证明了，在概率的条件下就是：实际上，这就是均值不等式来的啊， E 表示数学期望

[mathrm{E}[f(X)] geq f(mathrm{E} X) ]

EM 算法讲了什么

对于下面的第一行的式子, 本质上这个式子是函数隐变量的, 将这些隐变量显示的写出来就是下面的第二行的式子：

[egin{aligned} ell( heta) &=sum_{i=1}^{N} log p(x ; heta) \ &=sum_{i=1}^{N} log sum_{z} p(x, z ; heta) end{aligned} ]

现在，我们假设对于每一个 $ i $ , (Q_{i}) 表示 $ z $ 的分布，那么我们有 (sum_{z} Q_{i}(z)=1, Q_{i}(z) geq0) 。然后我们使用 Jensen 不等式将上面的式子进行放缩，写成下面的这样形式，

[egin{aligned} ell( heta) =sum_{i} log pleft(x^{(i)} ; heta ight) &=sum_{i} log sum_{z^{(i)}} pleft(x^{(i)}, z^{(i)} ; heta ight) \ &=sum_{i} log sum_{z^{(i)}} Q_{i}left(z^{(i)} ight) frac{pleft(x^{(i)}, z^{(i)} ; heta ight)}{Q_{i}left(z^{(i)} ight)} \ & geq sum_{i} sum_{z^{(i)}} Q_{i}left(z^{(i)} ight) log frac{pleft(x^{(i)}, z^{(i)} ; heta ight)}{Q_{i}left(z^{(i)} ight)} end{aligned} ]

这里的 (Q_{i}) 应该就是 EM 算法不断的迭代的关键。

上面的式子表明 (ell( heta)) 有一个下界，从极大似然的角度考虑，我们就是要最大化这个下界，得到 $ heta $ 的取值，但是其中的 (Q_{i}) 是一个隐变量的分布，我们并不知道这个分布是什么样子的。

上面使用 Jensen 不等式关键的一步是：

[fleft(mathrm{E}_{z^{(i)} sim Q_{i}}left[frac{pleft(x^{(i)}, z^{(i)} ; heta ight)}{Q_{i}left(z^{(i)} ight)} ight] ight) geq mathrm{E}_{z^{(i)} sim Q_{i}}left[fleft(frac{pleft(x^{(i)}, z^{(i)} ; heta ight)}{Q_{i}left(z^{(i)} ight)} ight) ight] ]

这个式子取等式的时候，对于数学期望而言，只有常数可以满足数学期望中的 Jensen 不等式相等。这里不做具体的证明，我们可以考虑从均值不等式来理解这个问题，假设这个常数是 $ c$ :

[frac{pleft(x^{(i)}, z^{(i)} ; heta ight)}{Q_{i}left(z^{(i)} ight)}=c ]

也就是说：

[Q_{i}left(z^{(i)} ight) propto pleft(x^{(i)}, z^{(i)} ; heta ight) ]

我们知道：(sum_{z} Q_{i}left(z^{(i)} ight)=1) ，那么我可以考虑这样的情况：

[1 = sum{Q_ileft( z^{left( i ight)} ight)} propto sum{pleft( x^{left( i ight)},z^{left( i ight)}; heta ight) } = sum_{z} pleft(x^{(i)}, z^{left(i ight)} ; heta ight) ]

这样的话，就有下面的公式：

[egin{aligned} Q_{i}left(z^{(i)} ight) &=frac{pleft(x^{(i)}, z^{(i)} ; heta ight)}{sum_{z} pleft(x^{(i)}, z^{left(i ight)} ; heta ight)} \ &=frac{pleft(x^{(i)}, z^{(i)} ; heta ight)}{pleft(x^{(i)} ; heta ight)} \ &=pleft(z^{(i)} | x^{(i)} ; heta ight) end{aligned} ]

这一步我们称为 E 步骤，可以得到等式成立的条件， (Q_{i}left(z^{(i)} ight) :=pleft(z^{(i)} | x^{(i)} ; heta ight)) ，下界成立的条件，我们将这个条件带入 $ell( heta) $ 得到的结果就是下界，也就是说，这个是最大似然函数的条件之一, 我们需要最大似然的就是：那么我们可以用下面的步骤来计算 $ heta$ ,

[ heta :=arg max _{ heta} sum_{i} sum_{z^{(i)}} Q_{i}left(z^{(i)} ight) log frac{pleft(x^{(i)}, z^{(i)} ; heta ight)}{Q_{i}left(z^{(i)} ight)} ]

那么这个 $ heta$ 是不是最优的呢？接下来我们来证明 EM 算法中一个很关键的问题，那就是上述的 E 步骤与 M 步骤不断放缩与收敛，最终得到最优的 $ heta$ 的值。

所以 EM 算法的步骤可以表示成下面这样：

1. 函数声明 (operatorname{EM}left(p_{Y, C}(y, c ; heta), p_{C | Y}(c | y ; heta), heta^{(0)} ight))
2. for iteration (t in 1,2, ldots) do
3. (Q_{i}^{(t)} leftarrow Pleft(z_i | x_i ; heta^{(t-1)} ight) quad( ext { E-step }))
4. ( heta^{(t)} leftarrow operatorname{argmax}_{ heta} mathbb{E}_{Q_{i}^{(t)}}left[P(y, C ; heta) ight] quadleft(mathrm{M}{-mathrm{Step}} ight))
5. if ( heta^{(t)} approx heta^{(t-1)}) then
6. return $ heta^{(t)} $

EM 算法的收敛问题证明

这个收敛的思想是这样的：我们的 $ heta $ 是直接通过极大似然函数算出的来的。那么 EM 算法迭代的步骤就是，我们不断地最大化极大似然估计函数，也就是说 我们每次都最大化了 (ellleft( heta^{(t)} ight)) 的下界，只有保证 (ellleft( heta^{(t)} ight) leq ellleft( heta^{(t+1)} ight)) , 这样就会不断地逼近最优解, 说明算法是正确的。在 EM 迭代的过程中，我们不断改变是 (Q_{i}) 也就是分布函数的选择(我们本质是更新 ( heta) , 改变 ( heta) 等价于改变 (Q_i) )，

为了更好的说明，假设上一步我们已经得到了一个最优解 (ellleft( heta^{(t)} ight)) 以及 $ heta^{(t)} $ 那么在这一步，我们满足下面的情况：

[Q_{i}^{(t)}left(z^{(i)} ight) :=pleft(z^{(i)} | x^{(i)} ; heta^{(t)} ight) Jensen不等式条件\ 用于计算 ellleft( heta^{(t+1)} ight) ]

[ellleft( heta^{(t)} ight)=sum_{i} sum_{z^{(i)}} Q_{i}^{(t)}left(z^{(i)} ight) log frac{pleft(x^{(i)}, z^{(i)} ; heta^{(t)} ight)}{Q_{i}^{(t)}left(z^{(i)} ight)} ]

而更新的参数 ( heta^{(t+1)}) 就是来自最大化上面右边的式子, 那么我们可以推出下面的不等式：

[egin{aligned} ellleft( heta^{(t+1)} ight) & geq sum_{i} sum_{z^{(i)}} Q_{i}^{(t)}left(z^{(i)} ight) log frac{pleft(x^{(i)}, z^{(i)} ; heta^{(t+1)} ight)}{Q_{i}^{(t)}left(z^{(i)} ight)} \ & geq sum_{i} sum_{z^{(i)}} Q_{i}^{(t)}left(z^{(i)} ight) log frac{pleft(x^{(i)}, z^{(i)} ; heta^{(t)} ight)}{Q_{i}^{(t)}left(z^{(i)} ight)} \ &=ellleft( heta^{(t)} ight) end{aligned} ]

对于上面的两个不等式做出以下解释：

1. 第一个不等式来自 Jensen 不等式，原式是这样的：(ell( heta) geq sum_{i} sum_{z^{(i)}} Q_{i}left(z^{(i)} ight) log frac{pleft(x^{(i)}, z^{(i)} ; heta ight)}{Q_{i}left(z^{(i)} ight)}) 对任意的$ Q_i $ 与( heta) 均成立，这里为什么对任意的 ( heta) 都成立, 是因为 (Q) 本质也是由 ( heta) 决定的.
2. 对于第二个不等式, 就是我们前面 M 步骤的结果, 最大化下面的极大似然函数的时候得到$ heta^{(t+1)}$ 所以第二个式子显然:

[sum_{i} sum_{z^{(i)}} Q_{i}^{(t)}left(z^{(i)} ight) log frac{pleft(x^{(i)}, z^{(i)} ; heta^{(t)} ight)}{Q_{i}^{(t)}left(z^{(i)} ight)} ]

综上证明了 EM 算法的收敛性问题

混合高斯分布

在一般的情况下，对于所得到的样本集，(X=left{x_{1}, dots, x_{N} ight})，我们的目标是最大化似然函数，通过最大化似然函数来获取参数的值。这是似然函数往往取对数表示就是：

[egin{aligned} L( heta | X) &=log left(prod_{i=1}^{N} pleft(x_{i} | heta ight) ight) \ &=sum_{i=1}^{N} log pleft(x_{i} | heta ight) end{aligned} ]

这个过程的参数估计可以描述成：

[hat{ heta}=arg max _{ heta} L( heta | X) ]

这个结果是可以直接计算出来的，那么引入混合高斯分布会是什么样呢？

混合高斯分布：

简单的说就是，非混合的情况下，我们的数据集满足高斯分布，这样用极大似然就直接算出高斯分布中的参数就可以了。那么混合的情况下就是，数据集是由多个满足高斯分布的数据集线性组合起来，这个时候我们理解为：有 (k) 个满足不同参数的高斯分布的原数据集，并且 (sum_{j=1}^{k} pi_{j}=1).

此时:

[p(oldsymbol{x})=sum_{k=1}^{K} pi_{k} mathcal{N}left(oldsymbol{x} | oldsymbol{mu}_{k}, oldsymbol{Sigma}_{k} ight) ]

举个例子说明:

我们假设男生的身高符合正态分布, (mathcal{N_1}left(oldsymbol{mu}_{1}, oldsymbol{Sigma}_{1} ight)) 女生的身高符合正态分布 (mathcal{N_2}left(oldsymbol{mu}_{2}, oldsymbol{Sigma}_{2} ight)) , 现在我们得到一个关于身高的样本数据集, 但是不知道男生与女生的比例, 要拟合这个身高数据集符合哪种分布, 这里就可以使用混合高斯分布,

[p(oldsymbol{x})= pi_{1} mathcal{N}left(oldsymbol{x} | oldsymbol{mu}_{1}, oldsymbol{Sigma}_{1} ight) + pi_{2} mathcal{N}left(oldsymbol{x} | oldsymbol{mu}_{2}, oldsymbol{Sigma}_{2} ight) ]

其中 (pi_1) 和 (pi_2) 就表示男生与女生的比例. 也表示 (p(x)) 对不同高斯分布的可能性,

在上面的问题中, (pi_1) 与 (pi_2) 是未知参数 , 那么这个例子与一个问题很类似, 那就是贝叶斯分布的问题, 如果我们在已知后验概率的情况下, 就可以很简单的计算先验概率了, 也就是说, 如果我们已知男生满足的高斯分布的参数, (mu_1) 与 ({Sigma}_{1}) 和 女生高斯分布的参数 (mu_2) 与 ({Sigma}_{2}) , 这时就可以反过来使用贝叶斯计算出 (pi_1) 与 (pi_2) 的值.

混合高斯分布与贝叶斯

对于一个混合高斯分布来说, (pi_1, pi_2, dots, pi_n) 对这个分布来说也是未知参数, 这里我们从贝叶斯的角度来改写这个未知参数的形式, 这样做的目的是为了更好的使用下面的 EM 算法以及 EM 算法的计算方式, 对于上面的 (N) 维混合高斯分布, 我们取一个向量(N) 维向量 ({z}) 表示取得的, 其中 $ z_k(1 leq k leq K) $, 只能取 0 与 1 两个值, 其中 (p(z_i == 1)) 表示第 i 中高斯分布被选中的概率, 而 (p(z_i = 0)) 表示未被选中的概率, 那么用向量 (z) 表示混合高斯分布就有两种方式, 一种是类似于

[[1,1,1,0,0,dots, 1, 0, 1] ]

的形式, 这样直接使用向量的乘法的到概率,

而另一种是

[[0,1,0, 0dots, dots, 0] ]

的形式, 这里只有一个 (z_i) 为 1, 因为我们要使用贝叶斯分布, 反过来使用全概率公式, 所以我们使用第二种形式, 在这种形式下, (z_i) 之间是相互独立的, 所以 (z) 向量的概率可以写成:

[p(oldsymbol{z})=pleft(z_{1} ight) pleft(z_{2} ight) ldots pleft(z_{K} ight)=prod_{k=1}^{K} pi_{k}^{z_{k}} ]

然后对应于 (z_i) 的那一个分布是服从高斯分布的, 可以将对应的 (z_i) 分布的概率表示如下:

[pleft(oldsymbol{x} | z_{i}=1 ight)=mathcal{N}left(oldsymbol{x} | oldsymbol{mu}_{i}, oldsymbol{Sigma}_{i} ight) ]

进而上式有可以写成如下形式：

[p(oldsymbol{x} | oldsymbol{z})=prod_{i=1}^{N} mathcal{N}left(oldsymbol{x} | oldsymbol{mu}_{i}, oldsymbol{Sigma}_{i} ight)^{z_{i}} ]

我们知道, 在向量 (z) 中, 假设 (z_j = 1) 那么对于所有的 (z_i (i eq j)) 都有 (z_i = 0), 所以上式成立很显然, 我们得到了 (p(z)) 与 (p(x|z)) 就可以得到 (p(x)) 的公式了:

[egin{aligned} p(oldsymbol{x}) &=sum_{oldsymbol{z}} p(oldsymbol{z}) p(oldsymbol{x} | oldsymbol{z}) \ &=sum_{oldsymbol{z}}left(prod_{k=1}^{K} pi_{k}^{z_{k}} mathcal{N}left(oldsymbol{x} | oldsymbol{mu}_{k}, oldsymbol{Sigma}_{k} ight)^{z_{k}} ight) \ &=sum_{k=1}^{K} pi_{k} mathcal{N}left(oldsymbol{x} | oldsymbol{mu}_{k}, oldsymbol{Sigma}_{k} ight) end{aligned} ag{1} ]

公式 (1) 本质上使用的是全概率公式, 因为向量 (z) 的表示方式我们取的是第二种, 上式中 (p(z)p(x|z)) 的结果是 (p(x,z)) 所以要用全概率公式相加, 从上面的式子我们惊奇的发现, 我们 (p(x)) 改写成了后验概率的形式, 因为我们划分了所有的 (z_i) , 在不同的向量 (z) 的基础上得到的概率 (p(x)), 这时反过来就可以求先验概率:

[egin{aligned} gammaleft(z_{i} ight) &=pleft(z_{i}=1 | oldsymbol{x} ight) \ &=frac{pleft(z_{i}=1 ight) pleft(oldsymbol{x} | z_{i}=1 ight)}{pleft(oldsymbol{x}, z_{i}=1 ight)} \ &=frac{pleft(z_{i}=1 ight) pleft(oldsymbol{x} | z_{i}=1 ight)}{sum_{j=1}^{N} pleft(z_{j}=1 ight) pleft(oldsymbol{x} | z_{j}=1 ight)}( ext { 全概率公式 }) \ &=frac{pi_{i} mathcal{N}left(oldsymbol{x} | oldsymbol{mu}_{i}, oldsymbol{Sigma}_{oldsymbol{i}} ight)}{sum_{j=1}^{N} pi_{j} mathcal{N}left(oldsymbol{x} | oldsymbol{mu}_{j}, oldsymbol{Sigma}_{j} ight)}( ext {结合}(1) . end{aligned} ag{2} ]

EM 算法用于混合高斯分布

前面的 EM 算法是一般情况下的 EM 算法求解以及证明的过程, 那么在实际的计算的过程中, EM 算法求解 GMM 问题是什么样的呢? 下面的公式 (3) 是基本的公式,

[p(oldsymbol{x|pi, mu, Sigma})=sum_{k=1}^{K} pi_{i} mathcal{N}left(oldsymbol{x} | oldsymbol{mu}_{i}, oldsymbol{Sigma}_{i} ight) ag{3} ]

这个公式与 EM 算法最终的最大化的公式是一样的(因为 Jensen 不等式的条件), 所以对这个公式求极大似然,

[prod_{i}^{N} p(oldsymbol{x|pi, mu, Sigma}) = prod_{n}^{N} sum_{i}^{N} pi_{i} mathcal{N}left(oldsymbol{x} | oldsymbol{mu}_{i}, oldsymbol{Sigma}_{i} ight) ag{4} ]

然后对这个公式求对数, 以及求偏导, 可以很容易得到下面的公式:

[0=-sum_{n=1}^{N} frac{pi_{i} mathcal{N}left(oldsymbol{x}_{n} | oldsymbol{mu}_{i}, oldsymbol{Sigma}_{i} ight)}{sum_{j} pi_{j} mathcal{N}left(oldsymbol{x}_{n} | oldsymbol{mu}_{j}, oldsymbol{Sigma}_{j} ight)} oldsymbol{Sigma}_{i}^{-1}left(oldsymbol{x}_{n}-oldsymbol{mu}_{i} ight) ag{5} ]

对于公式(5), 前面部分是通常的对数求导, 将 (sum^{K} pi_{i} mathcal{N}left(oldsymbol{x} | oldsymbol{mu}_{i}, oldsymbol{Sigma}_{i} ight)) 放在分母, 然后后面的 (oldsymbol{Sigma}_{i}^{-1}left(oldsymbol{x}_{n}-oldsymbol{mu}_{i} ight)) 就是一维的正态分布的求导公式, 这样就可以求出 (mu_i) 的值了,

[oldsymbol{mu}_{i}=frac{1}{N_{i}} sum_{n=1}^{N} gammaleft(z_{n i} ight) oldsymbol{x}_{n} ag{6.1} ]

其中

[N_i = sum_{n=1}^N gammaleft(z_{n i} ight) ag{6.2} ]

对上面公式的解释:

1. (gammaleft(z_{n i} ight)) 表示对于样本 (x_n) 属于第 (i) 类的概率, 在男生女生问题中就是对于某一个样本属于男生或者女生的概率
2. (N_i) 就是对于所有样本, 属于第 (i) 类的概率的个数总和
3. (mu_i) 就将这个值进行加权平均, 在公式 6.1 中, (x_n) 表示权值, (gammaleft(z_{n i} ight)) 表示概率

同理, 使用极大似然的方法可以求得:

[oldsymbol{Sigma}_{i}=frac{1}{N_{i}} sum_{n=1}^{N} gammaleft(z_{n i} ight)left(x_{n}-oldsymbol{mu}_{i} ight)left(x_{n}-oldsymbol{mu}_{i} ight)^{T} ]

最后是 (pi_1, pi_2, dots, pi_n) , 这个的极大值是条件极大值, 条件限制是 (sum_{1}^{N} pi_i = 1) , 也就是在求偏导的时候需要加入拉格朗日算子, 使用拉格朗日乘子法:

[ln p(oldsymbol{x} | oldsymbol{pi}, oldsymbol{mu}, oldsymbol{Sigma})+lambdaleft(sum_{i=1}^{N} pi_{i}-1 ight) ag{7.1} ]

求上式关于 (pi_i) 的极大似然函数, 得到:

[0=sum_{n=1}^{N} frac{mathcal{N}left(oldsymbol{x}_{n} | oldsymbol{mu}_{i}, oldsymbol{Sigma}_{i} ight)}{sum_{j} pi_{j} mathcal{N}left(oldsymbol{x}_{n} | oldsymbol{mu}_{j}, oldsymbol{Sigma}_{j} ight)}+lambda ag{7.2} ]

公式 7.2 求导的过程就是, 第一部分与前面对 (mu) 求导类似, 这里分子不变很显然, 对于结果两边同乘以 (pi_i) 得到:

[0=sum_{n=1}^{N} frac{pi_{i} mathcal{N}left(oldsymbol{x}_{n} | oldsymbol{mu}_{i}, oldsymbol{Sigma}_{i} ight)}{sum_{j} pi_{j} mathcal{N}left(oldsymbol{x}_{n} | oldsymbol{mu}_{j}, oldsymbol{Sigma}_{j} ight)}+lambda pi_{i} ]

结合公式 (2) 可以得到:

[0 = N_i + lambda pi_i \ 0 = sum_{i = 1}^{N} N_i + lambda sum_{i=1}^{N} pi_i \ 0 = N+k ag{7.3} ]

最终得到就是:

[pi_i = frac{N_i}{N} ]

EM 算法求解的迭代过程

EM 算法的不断迭代收敛过程就是不断地计算 E 步骤与 M 步骤, 直到期望的似然函数收敛,

E Step

根据当前的 (pi_i) , (mu_i) , (Sigma_i) 计算:

[y(z_nk) = frac{pi_{i} mathcal{N}left(oldsymbol{x}_{n} | oldsymbol{mu}_{i}, oldsymbol{Sigma}_{i} ight)}{sum_{j}^{N} pi_{j} mathcal{N}left(oldsymbol{x}_{n} | oldsymbol{mu}_{j}, oldsymbol{Sigma}_{j} ight)} ]

同时需要记录的值为:

[p(oldsymbol{x})=sum_{i=1}^{N} pi_{i} mathcal{N}left(oldsymbol{x} | oldsymbol{mu}_{i}, oldsymbol{Sigma}_{i} ight) ]

M Step

使用上面的公式更新参数:

[oldsymbol{mu}_{i+1}=frac{1}{N_{i}} sum_{n=1}^{N} gammaleft(z_{n i} ight) oldsymbol{x}_{n} \ oldsymbol{Sigma}_{i+1}=frac{1}{N_{i}} sum_{n=1}^{N} gammaleft(z_{n i} ight)left(x_{n}-oldsymbol{mu}_{i} ight)left(x_{n}-oldsymbol{mu}_{i} ight)^{T} \\ pi_{i+1} = frac{N_i}{N} ag{8} ]

其中:

[N_i = sum_{n=1}^N gammaleft(z_{n i} ight) ag{6.2} ]

最后再计算,

[ln p(oldsymbol{x} | oldsymbol{pi}, oldsymbol{mu}, oldsymbol{Sigma})=sum_{n=1}^{N} ln left{sum_{i=1}^{N} pi_{i} mathcal{N}left(oldsymbol{x}_{i} | oldsymbol{mu}_{i}, oldsymbol{Sigma}_{i} ight) ight} ]

判断是否收敛, 不收敛, 返回至 E Step, 然后在 M Step 继续更新, 直到收敛.