混合高斯
单一高斯模型无法应对如老忠实间歇喷泉这些实际的问题，而高斯混合模型提供了一类比单独的高斯分布更强大的概率模型。我们将高斯混合模型看成高斯分量的简单线性叠加，其公式为[注0]：
\[p(\mathbf x) = \sum_{k=1}^{K} \pi_{k} \mathcal N(\mathbf x|\mu_k, \Sigma_k) \tag {9.7} \]
引入一个K维的二值随机变量\(\mathbf z=(z_1, z_2,\cdots,z_K)\),这个变量采取"1-of-K"表示方法，其中一个特定元素\(z_k\)等于1，其余所有的元素等于0.于是\(z_k\)的值满足\(z_k \in {0, 1} 且 \sum_k z_k = 1\)。因此我们可以根据哪个元素非零，判断向量\(z\)有K个可能的状态。现在做出如下假设：向量\(\mathbf z\)的概率分布根据混合系数 \(\pi_k\)
\[p(z_k=1)=\pi_k \]
其中系数 \(\{\pi_k\}\) 满足 \(\pi_k \in [0,1]且\sum_{k}^K z_k = 1\)，使得上述概率是合法的。由于向量\(\mathbf z\)采用了“1-of-K”表示方法，因此可以将这个概率分布写成
\[p(\mathbf z) = \prod_k^{K} \pi_{k}^{z_k} \tag {9.10} \]
类似地，给定一个特定的值，x的条件概率分布是一个高斯分布
\[p(\mathbf x|z_k=1)=\mathcal N(x|\mu_k, \Sigma_k) \]
于是
\[p(\mathbf x|\mathbf z)=\prod_{k=1}^{K} \mathcal N(x|\mu_k, \Sigma_k)^{z_k} \tag{9.11} \]
因此向量\(\mathbf x,\mathbf z\)的联合概率分布为\(p(\mathbf x,\mathbf z)=p(\mathbf z)p(\mathbf x|\mathbf z)\)， 而\(\mathbf x\)的边缘概率分布可以通过将联合概率分布对所有可能的\(\mathbf z\)求和得到，即[注1]
\[p(\mathbf x)=\sum_{\mathbf z} p(\mathbf z)p(\mathbf x|\mathbf z)=\sum_k^{K}\pi_k \mathcal N(x|\mu_k, \Sigma_k) \tag {9.12} \]
\(\mathbf x\)的分布，然后通过引入一个新的变量\(\mathbf z\)，通过这两个变量联合分布对新变量\(\mathbf z\)所有可能求和，又得到这个混合高斯分布。转了一圈又回来了。
- 如果可以通过对联合概率关于某个变量所有可能求和可以表示边缘概率分布，即存在\(p(\mathbf x)=\sum_z p(\mathbf x,\mathbf z)\) 的表示形式，那么对于每个观测数据点\(x_n\)，存在一个与之对应的潜在变量\(z_n\)
- 潜在变量可以显式写出？
- 我们能够对联合概率分布\(p(\mathbf x,\mathbf z)\)操作，而不是对边缘概率分布\(p(\mathbf x)\)操作，这会极大地简化计算
- 给定x的条件下，z的条件概率可以求出，使用贝叶斯定理得到
\[\begin{aligned} p(z_k=1|\mathbf x) & = \frac{p(z_k=1)p(\mathbf x|z_k=1)}{\sum_{j=1}^{K}p(z_j=1)p(\mathbf x|z_j=1)}\\ & = \frac{\pi_k \mathcal N(\mathbf x|\mu_k, \Sigma_k)}{\sum_{j=1}^{K}\pi_j \mathcal N(\mathbf x|\mu_j, \Sigma_j)} \end{aligned} \tag {9.13} \]
我们用 \(\gamma(z_k)\) 表示 \(p(z_k=1|\mathbf x)\)。如果将 \(\pi_k\) 看成 \(z_k=1\) 的先验概率，则 \(\gamma(z_k)\) 就是观测到x之后，对应的后验概率。对应混合模型，\(\gamma(z_k)\) 也可以被看着第k个分量对应"解释"观测值x的"责任"[注2]。（\(z_k=1\)等价于变量z的第k分量是1，其余分量是0.即\(z=(0,0,\cdots,0,1,0,\cdots,0)\)）
最大似然
\({x_1,...,x_n}\), 我们希望使用混合高斯模型来对数据进行建模。这个数据可以表示为一个 \(N \times D\) 的矩阵\(X\)，其中第n行为\(x_n^T\)。类似地对应的潜在变量被表示为一个\(N \times K\) 的矩阵\(Z\)，它的行为\(z_n^T\)。我们假定数据点独立地从概率分布中取样。根据公式（9.7），对数似然函数为
\[ln\{p(X|\pi,\mu,\Sigma)\} = \sum_{n=1}^N ln \{ \sum_{k=1}^K \pi_k \mathcal N(x_n|\mu_k,\Sigma_k) \} \tag {9.14} \]
求解最大似然函数面临的问题：
（1）奇异性的存在。假设混合高斯模型，它的分量的协方差矩阵为 \(\Sigma_k=\sigma_k^2I\)，混合模型的第k分量均值和某个数据点完全相同，即对于某个n值，\(u_k=x_n\)，那么这个数据点为似然函数贡献一项，形式为
\[\mathcal N(x_n|\mu_k,\Sigma_k)= \frac{1}{(2\pi)^{0.5}} \frac{1}{\sigma_k} \tag{9.15} \]
如果考虑这个方差\(\sigma_k\)很小，趋于0，那么对数似然函数也会趋于无穷大。因此对数似然函数的最大化不是一个良好定义的问题[注3]。
（2）最大化高斯混合模型的对数似然函数，比单一的高斯分布更加复杂。困难源自在公式（9.14）中，对k的求和出现在对数计算的内部，从而对数函数无法之间作用于高斯分布，也就无法简化计算。我们令对数似然函数的导数等于零，也不会得到一个解析解。
用于GM的EM
有一种优雅的且强大的求解带有潜在变量的模型的最大似然解的方法被称为期望最大化算法，expectation-maximization algorithm，或者简称EM算法。
首先，令公式（9.14）中关于高斯分量的均值\(\mu_k\)导数等于零[注4]，得到
\[\sum_{n=1}^{N} \frac{\pi_k \mathcal N(x_n|\mu_k,\Sigma_k)}{\sum_{j=1}^{K} \pi_j \mathcal N(x_n|\mu_j,\Sigma_j)} \Sigma_k^{-1}(x_n-\mu_k) = 0 \tag{9.16} \]
用上文提到的后验概率（9.13）代替简化公式得
\[\sum_{n=1}^{N} \gamma(z_k) \Sigma_k^{-1}(x_n-\mu_k) = 0 \]
再对上述公式两侧同时乘以\(\Sigma_k\)
\[\sum_{n=1}^{N} \gamma(z_k) (x_n-\mu_k) = 0 \]
\[\sum_{n=1}^{N} \gamma(z_{nk}) x_n - \sum_{n=1}^{N} \gamma(z_{nk}) \mu_k = 0 \]
\[\mu_k = \frac{\sum_{n=1}^{N} \gamma(z_{nk}) x_n}{\sum_{n=1}^{N} \gamma(z_{nk})} \tag{9.17} \]
可以看到第k个高斯分量的均值\(\mu_k\)可以通过对数据集中所有的数据点\(x_n\)加权平均的方式得到，其中数据点\(x_n\)的权重因子由后验概率\(\gamma(z_{nk})\)给出。
令公式(9.14)中关于高斯分布的协方差矩阵\(\Sigma_k\)的导数等于零[注5]，得到
\[\Sigma_k = \frac{\sum_{n=1}^{N} \gamma(z_{nk}) (x_n-\mu_k)(x_n-\mu_k)^T}{\sum_{n=1}^{N} \gamma(z_{nk})} \tag{9.19} \]
最大化公式 (9.14)需考虑限制条件，混合系数之和等于1. 使用拉格朗日乘法，最大化下面的量 $$ lnp(X|\pi, \mu, \Sigma) + \lambda (\sum_k^K \pi_k-1) \tag{9.20}$$
求导得到
\[0=\sum_{n=1}^N \frac{\mathcal N(x_n|\mu_k, \Sigma_k)}{\sum_{j=1}^K\pi_j \mathcal N (x_n|\mu_j,\Sigma_j)} + \lambda \tag{9.21} \]
上式两边同时乘以\(\pi_k\)，然后再对k 求和得到
\[0=\sum_{k=1}^K \sum_{n=1}^N \frac{\pi_k \mathcal N(x_n|\mu_k, \Sigma_k)}{\sum_{j=1}^K\pi_j \mathcal N (x_n|\mu_j,\Sigma_j)} + \lambda \sum_{k=1}^K \pi_k \]
\[ 0=\sum_{n=1}^N \frac{\sum_{k=1}^K \pi_k \mathcal N(x_n|\mu_k, \Sigma_k)}{\sum_{j=1}^K\pi_j \mathcal N (x_n|\mu_j,\Sigma_j)} + \lambda\]
从而推导得到 $$\lambda = -N$$
再代入公式（9.21）并同时乘以\(\pi_k\)得到
\[\pi_k = \frac{1}{N}\sum_{n=1}^N \frac{\pi_k \mathcal N(x_n|\mu_k, \Sigma_k)}{\sum_{j=1}^K\pi_j \mathcal N (x_n|\mu_j,\Sigma_j)}=\frac{1}{N}\sum_{n=1}^N \gamma(z_{nk})=\frac{N_k}{N} \tag{9.22} \]
附录
[注0] 为了便于大家与原文进行对比，这里的公式索引，仍采用和书本一样的索引号。
[注1]
\[\begin{aligned} p(z)p(x|z) & = \prod_k^{K} \pi_{k}^{z_k} \mathcal N(x|\mu_k, \Sigma_k)^{z_k} \\ & = \prod_k^{K} \{\pi_{k}\mathcal N(x|\mu_k, \Sigma_k)\}^{z_k} \end{aligned}\]
由于z取值是"1-of-K"，因此\(z_k\)只1有个是1其余都等于零，即
\[\begin{aligned} p(z|_{z_k=1})p(x|z|_{z_k=1}) & = \pi_{k}\mathcal N(x|\mu_k, \Sigma_k) \end{aligned}\]
对向量\(z\)的所有可能求和就可以得到
\[p(x)=\sum_z p(z_{z_k=1})p(x|z_{z_k=1})=\sum_k^{K}\pi_k \mathcal N(x|\mu_k, \Sigma_k) \]
[注2]
原文中是这么描述的：\(\gamma(z_k)\) can also be viewed as the responsibility that component \(k\) takes for 'explaining' the observation \(x\). 直白点讲，能观察到数据点\(x\)，高斯混合分布中第k个分布贡献占比是\(\gamma(z_k)\)，所有的分量贡献比之和等于1.
\(z=(1,0,0),z=(0,1,0),z=(0,0,1)\)，依次对应着颜色red,green,blue。变量z的概率分布，以及x的条件分布依次为
\[p(z)|_{z_1=1}=0.1, p(z)|_{z_2=1}=0.3, p(z)|_{z_3=1}=0.6 \]
\[p(x|z_1=1) =\mathcal N(\mu_1, \Sigma_1), z_1=\binom{-2}{-1}, \Sigma_1 = \begin{pmatrix} 1 & -1 \\ -1 & 2 \\ \end{pmatrix} \]
\[p(x|z_2=1) =\mathcal N(\mu_2, \Sigma_2), z_2=\binom{1}{1}, \Sigma_2 = \begin{pmatrix} 1 & 1 \\ 1 & 2 \\ \end{pmatrix} \]
\[p(x|z_3=1) =\mathcal N(\mu_3, \Sigma_3), z_1=\binom{4.5}{2.5}, \Sigma_3 = \begin{pmatrix} 1 & -1 \\ -1 & 2 \\ \end{pmatrix} \]
使用如下代码画出的类似书中图9.5图形如下：
- 图左：随机从概率分布\(p(z)\)取样得到\(z_*\)，再由条件概率分布\(p(x|z_*)\)取样得到二维点\(x_*\)。二维点\(x_*\)可以直接在直角坐标系用点画出，但是\(z_*\)不好直接在绘制在坐标系中，于是映射为一种颜色，然后涂在二维坐标点\(x_*\)上，即一一对应又完整了展示了联合变量\((x,z)\)的取样。
- 图中：忽略变量\(z\)的取值。（实际自然界现象可能很多都是复杂的多元分布，能直观观察到可能就是单个变量的取样，而背后隐藏着无法观察或者未知的，但与当前观察一一对应的某个变量值。）
- 图右：观察到数据点后，按照后验概率\(\gamma(z_k|x)\)公式求出，各个分量贡献占比，然后给数据点着色。可以看出位置点越靠近左右两侧的高斯分布的，其颜色越是接近红色或者蓝色。说明对应分量的影响越大。
import matplotlib.pyplot as plt
import numpy as np
from numpy.linalg import det, inv
global colors
colors = ['red', 'green', 'blue'] # 颜色值：红、绿、蓝
global prob_colors
prob_colors = np.array([0.1, 0.3, 0.6]) # 颜色变量Z取值的概率
# 条件概率 prob(x|z=red)
global mu_red
mu_red = np.array([-2, -1])
global sigma_red
sigma_red = np.array([[1,-1],[-1,2]])
# 条件概率 prob(x|z=green)
global mu_green
mu_green = np.array([1, 1])
global sigma_green
sigma_green = np.array([[1,1],[1,2]])
# 条件概率 prob(x|z=green)
global mu_blue
mu_blue = np.array([4.5, 2.5])
global sigma_blue
sigma_blue = np.array([[1,-1],[-1,2]])
def prob_value(x, mu, sigma):
diff_term = (x - mu).reshape(2,1)
D = len(x)
den = (2 * np.pi) ** (D / 2) * np.sqrt(det(sigma))
num = np.einsum("ij,jk->ik", np.einsum("ij,jk->ik", diff_term.T, inv(sigma)), diff_term)
num = np.exp(-num / 2).squeeze()
return num / den
def gamma_z(x):
x_red = prob_value(x, mu_red, sigma_red)
x_green = prob_value(x, mu_green, sigma_green)
x_blue = prob_value(x, mu_blue, sigma_blue)
den = prob_colors[0] * x_red + prob_colors[1] * x_green + prob_colors[2] * x_blue
return (prob_colors[0]*x_red/den, prob_colors[1]*x_green/den, prob_colors[2]*x_blue/den)
num = 1000 # 总样本点数
x_red = []
x_green = []
x_blue = []
for i in range(num):
rands = np.random.choice(colors, 1, p=prob_colors)
if rands[0] == 'red':
prob_norm = np.random.multivariate_normal(mu_red, sigma_red, size=1)
x_red.append(prob_norm[0])
if rands[0] == 'green':
prob_norm = np.random.multivariate_normal(mu_green, sigma_green, size=1)
x_green.append(prob_norm[0])
if rands[0] == 'blue':
prob_norm = np.random.multivariate_normal(mu_blue, sigma_blue, size=1)
x_blue.append(prob_norm[0])
x_red=np.array(x_red)
x_green=np.array(x_green)
x_blue=np.array(x_blue)
x_all = np.concatenate((x_red, x_green, x_blue), axis=0)
gamma_zk = []
for i in range(num):
xi = x_all[i]
gamma_zk.append(gamma_z(xi))
plt.figure(figsize=(9,3))
plt.subplot(1,3,1)
plt.plot(x_red[:,0], x_red[:,1], '.', alpha=1, color='r')
plt.plot(x_green[:,0], x_green[:,1], '.', alpha=1, color='g')
plt.plot(x_blue[:,0], x_blue[:,1], '.', alpha=1, color='b')
plt.subplot(1,3,2)
x_all = np.concatenate((x_red, x_green, x_blue), axis=0)
plt.plot(x_all[:,0], x_all[:,1], '.', alpha=1, color='pink')
tt = np.array([[1,0,0], [0,1,1]])
plt.subplot(1,3,3)
for ind in range(num):
x = x_all[ind]
plt.plot(x[0], x[1], '.', alpha=1, c=gamma_zk[ind])
plt.show()[注3]：为什么单一高斯分布不会出现奇异性问题？
[注4]：如何推导这个（9.14）关于高斯分量均值\(\mu_k\)的导数
首先，设$$f(X)=X^TAX, X \in R^{n\times1}, A \in R^{n \times n}$$
\[X=\begin{pmatrix} x_1 \\ x_2 \\ \vdots \\ x_n \end{pmatrix}, A=\left ( \begin{matrix} a_{11} & a_{12} & \cdots & a_{1n} \\ a_{21} & a_{22} & \cdots & a_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ a_{n1} & a_{n2} & \cdots & a_{nn} \end{matrix} \right) \]
\[f(X)=(x_1, x_2, \cdots, x_n) \left ( \begin{matrix} a_{11} & a_{12} & \cdots & a_{1n} \\ a_{21} & a_{22} & \cdots & a_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ a_{n1} & a_{n2} & \cdots & a_{nn} \end{matrix} \right) \begin{pmatrix} x_1 \\ x_2 \\ \vdots \\ x_n \end{pmatrix} \]
\[f(X)=(x_1, x_2, \cdots, x_n) \begin{pmatrix} a_{11}x_1 + a_{12}x_2 + \cdots + a_{1n}x_n \\ a_{21}x_1 + a_{22}x_2 + \cdots + a_{2n}x_n \\ \vdots \\ a_{n1}x_1 + a_{n2}x_2 + \cdots + a_{nn}x_n \\ \end{pmatrix} \]
因此\(f(X)\)等于下面矩阵的各个元素相加
\[ \left ( \begin{matrix} a_{11}x_1x_1 & a_{12}x_1x_2 & \cdots & a_{1n}x_1x_n \\ a_{21}x_2x_1 & a_{22}x_2x_2 & \cdots & a_{2n}x_2x_n \\ \vdots & \vdots & \ddots & \vdots \\ a_{n1}x_nx_1 & a_{n2}x_nx_2 & \cdots & a_{nn}x_nx_n \end{matrix} \right) \]
此时假设矩阵\(A\)是对称矩阵, \(A_j\)表示矩阵A的第j列, 现在求\(f(X)\)关于\(X\)中第k个分量的导数，而\(f(X)\)中关于包含与\(x_k\)的，对应上述矩阵的第k行和第k列。再考虑A的对称性，可以得到
\[a_{kk}x_kx_k + 2\sum_{j \neq k}^{n}a_{kj}x_kx_j \]
上式关于\(x_k\)求导得到
\[\frac{\partial}{\partial x_k} = 2\sum_{j=1}^{n} a_{j,k}x_k = 2X^TA_j \]
从而求得
\[\frac{\partial f}{\partial X^T} = (\frac{\partial f}{\partial x_1}, \frac{\partial f}{\partial x_2}, \cdots, \frac{\partial f}{\partial x_n}) = (2X^TA_1, 2X^TA_2, \cdots, 2X^TA_n) = 2X^T(A_1, A_2, \cdots, A_n) = 2X^TA \]
因此,$$f(X)=X^TAX$$关于向量\(X\)的导数为
\[\frac{\partial f}{\partial X^T} = 2X^TA \\ 或\frac{\partial f}{\partial X} = 2AX \]
多元高斯分布的概率密度函数为
\[\mathcal N(X| \mu, \Sigma) = \frac{1}{\sqrt{(2\pi)^n |\Sigma|} } exp\{ -\frac{1}{2}(X-\mu)^T \Sigma^{-1} (X-\mu)\} \]
因此关于\(X\)的导数为
\[\frac{\partial{\mathcal N(X| \mu, \Sigma)}}{\partial{X}} = -\mathcal N(X| \mu, \Sigma)\Sigma^{-1} (X-\mu) \]
\[\]
