Gaussian Mixture Model(GMM)：高斯混合模型

Gaussian Mixture Model.

高斯混合模型(Gaussian Mixture Model,GMM)是一类概率生成模型，表示为 $K$ 个高斯分布的加权平均。其隐变量(latent variable)是服从多项分布的离散型随机变量，表示属于每个高斯分布的概率，记作 $z$ ，表示如下：

P (z = k) = p^{k}, k = 1, 2, . . ., K . K \sum k = 1 p^{k} = 1

$P(z=k) = p^k, \quad k=1,2,...,K. \quad \sum_{k=1}^{K}p^k = 1$

记观测变量(observable variable)为 $x$ ，则高斯混合模型可表示为：

P (x) = \sum z P (x, z) = \sum z P (z) P (x | z) = K \sum k = 1 P (z = k) P (x | z = k) = K \sum k = 1 p^{k} N (x; μ^{k}, σ^{k})

$P(x) = \sum_{z}^{} P(x,z) = \sum_{z}^{} P(z) P(x|z) \\ = \sum_{k=1}^{K} P(z=k) P(x|z=k) \\ = \sum_{k=1}^{K} p^k \mathcal{N}(x; \mu^k,\sigma^k)$

假设观测数据为 $\{x_1,x_2,...,x_n\}$ ，高斯混合模型的模型参数为 $\theta = \{ p^1,p^2,...,p^K, \mu^1,\mu^2,...,\mu^K,\sigma^1,\sigma^2,...,\sigma^K\}$ 。

尝试直接用极大似然估计方法求解该问题。求解如下：

^θ = arg max θ l o g p (x; θ) = arg max θ l o g n \prod i = 1 P (x_{i}; θ) = arg max θ n \sum i = 1 l o g P (x_{i}; θ) arg max θ n \sum i = 1 l o g K \sum k = 1 P^{k} N (x_{i}; μ^{k}, σ^{k})

$\hat{\theta} = \mathop{\arg \max}_{\theta} logp(x;\theta) = \mathop{\arg \max}_{\theta} log \prod_{i=1}^{n} P(x_i;\theta) = \mathop{\arg \max}_{\theta} \sum_{i=1}^{n} log P(x_i;\theta) \\ \mathop{\arg \max}_{\theta} \sum_{i=1}^{n} log \sum_{k=1}^{K} P^k \mathcal{N}(x_i; \mu^k,\sigma^k)$

上式在对数函数中存在求和项，直接求解是intractable的。采用期望最大算法求解该问题。

E-step

EM算法的E-step计算如下：

E-step： P (z | x; θ^{(t)}) \to E_{P (z | x; θ^{(t)})} [l o g P (x, z; θ)]

$\text{E-step：} P(z|x ; \theta^{(t)}) → \Bbb{E}_{P(z|x ; \theta^{(t)})}[logP(x,z ; \theta)]$

期望计算如下：

E_{P (z | x; θ^{(t)})} [l o g P (x, z; θ)] = \int_{z} P (z | x; θ^{(t)}) l o g P (x, z; θ) d z = \sum z_{1} . . . \sum z_{n} P (z | x; θ^{(t)}) l o g P (x, z; θ) = \sum z_{1} . . . \sum z_{n} n \prod i = 1 P (z_{i} | x_{i}; θ^{(t)}) l o g n \prod j = 1 P (x_{j}, z_{j}; θ) = \sum z_{1} . . . \sum z_{n} n \prod i = 1 P (z_{i} | x_{i}; θ^{(t)}) n \sum j = 1 l o g P (x_{j}, z_{j}; θ) = \sum z_{1} . . . \sum z_{n} n \prod i = 1 P (z_{i} | x_{i}; θ^{(t)}) l o g P (x_{1}, z_{1}; θ) + . . . + \sum z_{1} . . . \sum z_{n} n \prod i = 1 P (z_{i} | x_{i}; θ^{(t)}) l o g P (x_{n}, z_{n}; θ) = \sum z_{1} P (z_{1} | x_{1}; θ^{(t)}) l o g P (x_{1}, z_{1}; θ) n \prod j = 2 \sum z_{j} P (z_{j} | x_{j}; θ^{(t)}) + . . . + \sum z_{n} P (z_{n} | x_{n}; θ^{(t)}) l o g P (x_{n}, z_{n}; θ) n - 1 \prod j = 1 \sum z_{j} P (z_{j} | x_{j}; θ^{(t)}) (由于 \sum z_{j} P (z_{j} | x_{j}; θ^{(t)}) = 1) = n \sum i = 1 \sum z_{i} P (z_{i} | x_{i}; θ^{(t)}) l o g P (x_{i}, z_{i}; θ) = n \sum i = 1 K \sum k = 1 P (z_{i} = k | x_{i}; θ^{(t)}) l o g P (x_{i}, z_{i} = k; θ)

$\Bbb{E}_{P(z|x ; \theta^{(t)})}[logP(x,z ; \theta)] = \int_{z}^{} P(z|x ; \theta^{(t)}) logP(x,z ; \theta) dz \\ = \sum_{z_1}^{}...\sum_{z_n}^{} P(z|x ; \theta^{(t)}) logP(x,z ; \theta) \\ = \sum_{z_1}^{}...\sum_{z_n}^{} \prod_{i=1}^{n} P(z_i|x_i ; \theta^{(t)}) log \prod_{j=1}^{n} P(x_j,z_j ; \theta) \\ = \sum_{z_1}^{}...\sum_{z_n}^{} \prod_{i=1}^{n} P(z_i|x_i ; \theta^{(t)}) \sum_{j=1}^{n} log P(x_j,z_j ; \theta) \\ = \sum_{z_1}^{}...\sum_{z_n}^{} \prod_{i=1}^{n} P(z_i|x_i ; \theta^{(t)}) log P(x_1,z_1 ; \theta) + ... \\ + \sum_{z_1}^{}...\sum_{z_n}^{} \prod_{i=1}^{n} P(z_i|x_i ; \theta^{(t)}) log P(x_n,z_n ; \theta) \\ = \sum_{z_1}^{} P(z_1|x_1 ; \theta^{(t)}) log P(x_1,z_1 ; \theta) \prod_{j=2}^{n} \sum_{z_j}^{} P(z_j|x_j ; \theta^{(t)}) + ... \\ + \sum_{z_n}^{} P(z_n|x_n ; \theta^{(t)}) log P(x_n,z_n ; \theta) \prod_{j=1}^{n-1} \sum_{z_j}^{} P(z_j|x_j ; \theta^{(t)}) \\ (\text{由于}\sum_{z_j}^{} P(z_j|x_j ; \theta^{(t)})=1) \\ = \sum_{i=1}^{n} \sum_{z_i}^{} P(z_i|x_i ; \theta^{(t)}) log P(x_i,z_i ; \theta) \\ = \sum_{i=1}^{n} \sum_{k=1}^{K} P(z_i=k|x_i ; \theta^{(t)}) log P(x_i,z_i=k ; \theta)$

联合概率 $P(x,z)$ 计算如下：

P (x, z) = P (z) P (x | z) = P (z) N (x; μ_{z}, σ_{z})

$P(x,z) = P(z)P(x|z) = P(z)\mathcal{N}(x; \mu_z,\sigma_z)$

条件概率 $P(z|x)$ 计算如下：

P (z | x) = \frac{P (x, z)}{P (x)} = \frac{P (z) N (x; μ_{z}, σ_{z})}{\sum_{k = 1}^{K} p^{k} N (x; μ^{k}, σ^{k})}

$P(z|x) = \frac{P(x,z)}{P(x)} = \frac{P(z)\mathcal{N}(x; \mu_z,\sigma_z)}{\sum_{k=1}^{K} p^k\mathcal{N}(x; \mu^k,\sigma^k)}$

记在当前模型参数 $\theta^{(t)}$ 下第 $i$ 个观测数据来自第 $k$ 个分高斯模型的概率为 $\gamma_{ik}$ ，称为分模型 $k$ 对观测数据 $x_i$ 的响应度。即：

γ_{i k} = P (z_{i} = k | x_{i}; θ^{(t)}) = \frac{p^{k} N (x_{i}; μ^{k}, σ^{k})}{\sum_{k = 1}^{K} p^{k} N (x_{i}; μ^{k}, σ^{k})}

$\gamma_{ik} = P(z_i=k|x_i ; \theta^{(t)}) = \frac{p^k\mathcal{N}(x_i; \mu^k,\sigma^k)}{\sum_{k=1}^{K} p^k\mathcal{N}(x_i; \mu^k,\sigma^k)}$

则原期望计算公式可以表示为：

n \sum i = 1 K \sum k = 1 P (z_{i} = k | x_{i}; θ^{(t)}) l o g P (x_{i}, z_{i} = k; θ) = n \sum i = 1 K \sum k = 1 γ_{i k} l o g P (z_{i} = k) P (x_{i} | z_{i} = k) = n \sum i = 1 K \sum k = 1 γ_{i k} l o g p^{k} N (x_{i}; μ^{k}, σ^{k}) = n \sum i = 1 K \sum k = 1 γ_{i k} [l o g p^{k} + l o g N (x_{i}; μ^{k}, σ^{k})]

$\sum_{i=1}^{n} \sum_{k=1}^{K} P(z_i=k|x_i ; \theta^{(t)}) log P(x_i,z_i=k ; \theta) \\ = \sum_{i=1}^{n} \sum_{k=1}^{K} \gamma_{ik} log P(z_i=k)P(x_i|z_i=k) = \sum_{i=1}^{n} \sum_{k=1}^{K} \gamma_{ik} log p^k\mathcal{N}(x_i; \mu^k,\sigma^k) \\ = \sum_{i=1}^{n} \sum_{k=1}^{K} \gamma_{ik} [log p^k+log\mathcal{N}(x_i; \mu^k,\sigma^k)]$

M-step

EM算法的M-step计算如下：

M-step： θ^{(t + 1)} = arg max θ E_{P (z | x; θ^{(t)})} [l o g P (x, z; θ)]

$\text{M-step：} \theta^{(t+1)} = \mathop{\arg \max}_{\theta} \Bbb{E}_{P(z|x ; \theta^{(t)})}[logP(x,z ; \theta)]$

计算 ${p^k}^{(t+1)}$

{p^{k}}^{(t + 1)} = arg max p^{k} n \sum i = 1 K \sum k = 1 γ_{i k} [l o g p^{k} + l o g N (x_{i}; μ^{k}, σ^{k})] s.t. K \sum k = 1 p^{k} = 1

${p^k}^{(t+1)} = \mathop{\arg \max}_{p^k} \sum_{i=1}^{n} \sum_{k=1}^{K} \gamma_{ik} [log p^k+log\mathcal{N}(x_i; \mu^k,\sigma^k)] \\ \text{s.t. } \sum_{k=1}^{K} p^k=1$

采用拉格朗日乘子法解决上述约束最优化问题。建立拉格朗日函数：

L (p^{k}, λ) = n \sum i = 1 K \sum k = 1 γ_{i k} [l o g p^{k} + l o g N (x_{i}; μ^{k}, σ^{k})] + λ (K \sum k = 1 p^{k} - 1)

$\mathop{L}(p^k,\lambda) = \sum_{i=1}^{n} \sum_{k=1}^{K} \gamma_{ik} [log p^k+log\mathcal{N}(x_i; \mu^k,\sigma^k)] + \lambda (\sum_{k=1}^{K} p^k-1)$

令拉格朗日函数关于参数 $p^k$ 的导数为零：

$\frac{\partial \mathop{L}(p^k,\lambda)}{\partial p^k} = \sum_{i=1}^{n} \gamma_{ik} \cdot \frac{1}{p^k} + \lambda = 0$

解得 $p^k = -\frac{\sum_{i=1}^{n} \gamma_{ik}}{\lambda}$ 。又由 $\sum_{k=1}^{K} p^k=1$ ，可得 $\lambda = \sum_{i=1}^{n}\sum_{k=1}^{K} \gamma_{ik} = n$ ，因此参数 $p^k$ 的估计值为：

${p^k}^{(t+1)} = \frac{\sum_{i=1}^{n} \gamma_{ik}}{n}$

计算 ${\mu^k}^{(t+1)}$

${\mu^k}^{(t+1)} = \mathop{\arg \max}_{\mu^k} \sum_{i=1}^{n} \sum_{k=1}^{K} \gamma_{ik} [log p^k+log\mathcal{N}(x_i; \mu^k,\sigma^k)] \\ = \mathop{\arg \max}_{\mu^k} \sum_{i=1}^{n} \sum_{k=1}^{K} \gamma_{ik} log\mathcal{N}(x_i; \mu^k,\sigma^k) \\ = \mathop{\arg \max}_{\mu^k} \sum_{i=1}^{n} \sum_{k=1}^{K} \gamma_{ik} log \frac{1}{\sqrt{2\pi}\sigma^k} exp \{ -\frac{(x_i-\mu^k)^2}{ 2(\sigma^k)^2} \} \\ = \mathop{\arg \max}_{\mu^k} \sum_{i=1}^{n} \sum_{k=1}^{K} \gamma_{ik} [ log\frac{1}{\sqrt{2\pi}} - log\sigma^k -\frac{(x_i-\mu^k)^2}{ 2(\sigma^k)^2} ]$

令上述表达式关于参数 $\mu^k$ 的导数为零：

$\frac{\partial}{\partial \mu^k} \sum_{i=1}^{n} \sum_{k=1}^{K} - \gamma_{ik} \frac{(x_i-\mu^k)^2}{ 2(\sigma^k)^2} = \sum_{i=1}^{n} \gamma_{ik} \frac{2(x_i-\mu^k)}{ 2(\sigma^k)^2} = 0$

解得：

${\mu^k}^{(t+1)} = \frac{\sum_{i=1}^{n} \gamma_{ik} x_i}{\sum_{i=1}^{n} \gamma_{ik}}$

计算 ${\sigma^k}^{(t+1)}$

${\sigma^k}^{(t+1)} = \mathop{\arg \max}_{\sigma^k} \sum_{i=1}^{n} \sum_{k=1}^{K} \gamma_{ik} [ log\frac{1}{\sqrt{2\pi}} - log\sigma^k -\frac{(x_i-\mu^k)^2}{ 2(\sigma^k)^2} ]$

令上述表达式关于参数 $\sigma^k$ 的导数为零：

$\frac{\partial}{\partial \sigma^k} \sum_{i=1}^{n} \sum_{k=1}^{K} \gamma_{ik} [ - log\sigma^k -\frac{(x_i-\mu^k)^2}{ 2(\sigma^k)^2} ] = \sum_{i=1}^{n} \gamma_{ik} [-\frac{1}{\sigma^k} + \cdot \frac{(x_i-\mu^k)^2}{ (\sigma^k)^3}] = 0$

解得：

$({\sigma^k}^{(t+1)})^2 = \frac{\sum_{i=1}^{n} \gamma_{ik}(x_i-{\mu^k}^{(t+1)})^2}{\sum_{i=1}^{n} \gamma_{ik}}$

算法总结

数据样本集为 $\{x_1,x_2,...,x_n\}$ ，建立高斯混合模型：

$P(x) = \sum_{k=1}^{K} p^k \mathcal{N}(x; \mu^k,\sigma^k)$

待求模型参数为 $\theta = \{ p^1,p^2,...,p^K, \mu^1,\mu^2,...,\mu^K,\sigma^1,\sigma^2,...,\sigma^K\}$ ，随机取参数的初始值开始迭代。

根据当前模型参数，计算分模型 $k$ 对观测数据 $x_i$ 的响应度：

$\gamma_{ik} = P(z_i=k|x_i ; \theta^{(t)}) = \frac{p^k\mathcal{N}(x_i; \mu^k,\sigma^k)}{\sum_{k=1}^{K} p^k\mathcal{N}(x_i; \mu^k,\sigma^k)}$

计算新一轮迭代的模型参数：

${p^k}^{(t+1)} = \frac{\sum_{i=1}^{n} \gamma_{ik}}{n}$

${\mu^k}^{(t+1)} = \frac{\sum_{i=1}^{n} \gamma_{ik} x_i}{\sum_{i=1}^{n} \gamma_{ik}}$

$({\sigma^k}^{(t+1)})^2 = \frac{\sum_{i=1}^{n} \gamma_{ik}(x_i-{\mu^k}^{(t+1)})^2}{\sum_{i=1}^{n} \gamma_{ik}}$

重复迭代直至收敛。

E-step

M-step

计算pk(t+1)pk(t+1){p^k}^{(t+1)}

计算μk(t+1){\mu^k}^{(t+1)}

计算σk(t+1){\sigma^k}^{(t+1)}

算法总结

分享到微信朋友圈

计算 ${p^k}^{(t+1)}$

计算 ${\mu^k}^{(t+1)}$

计算 ${\sigma^k}^{(t+1)}$