隐变量自回归模型是什么 潜变量和隐变量

目录

• 潜（隐)变量模型
• K-means
• GMM模型

• GMM模型参数估计的EM算法总结
• GMM模型和K-means的联系

• EM算法

• 使用EM算法通用步骤重新考虑GMM参数估计
• EM算法通用解释

• python代码实现

潜（隐)变量模型

观测变量：直接观测的数据
潜变量：

• 无法直接被观测到，需要通过模型和观测变量进行推断
• 利用潜变量来判断观测变量的模型，GMM HMM都是潜变量模型
• 潜变量模型将不完数据（只有观测数据）的边缘分布转化为容易处理的完全数据（观测数据+潜变量）的联合分布。

如：聚类问题，潜变量是类别（未知），观测变量是数据点，给定观测变量，如果推断哪些是同一类？K-means

K-means

给定一个含有N个数据点的集合$\left\{\mathbf{x}_{1}, \mathbf{x}_{2}, \ldots, \mathbf{x}_{N}\right\}$,$\mathbf{x}_{\mathbf{n}} \in R^{D}$，聚类的目标是将此N个数据点聚类到K个类别中，且假设K给定。

K-means思路：

1. 引入K个D维均值向量$\mathbf{\mu}_k,k=1,2,...,K$，$\mathbf{\mu}_k$是第k个类别的聚类中心。
2. 计算数据点$\mathbf{x}_{\mathbf{n}}$和所有类中心$\mathbf{\mu}_k$的距离，类中心距离此数据点最近的类别，即为当前数据点的类别。
3. 根据新的聚类结果，使用当前聚集到各个类别的数据的均值来更新当前类别的类中心。
4. 返回第2步，直到满足一定的停止准则。

引入潜变量

• 对于每一个数据点$\mathbf{x}_{n}$引入一个指示因子$r_{nk} \in \{0,1\}$，如果$\mathbf{x}_{n}$属于第k类，则$r_{nk}=1$，否则$r_{nk}=0，$$r_{nk}$即为潜变量
• 定义目标函数$J=\sum_{n=1}^{N} \sum_{k=1}^{K} r_{n k}\left\|\mathbf{x}_{n}-\boldsymbol{\mu}_{k}\right\|^{2}$
• 优化目标：寻找合适的$r_{nk} $和$\boldsymbol{\mu}_{k}$使目标函数最小。

模型优化：两阶段迭代优化（简单EM）

• 选择初始化的$\boldsymbol{\mu}_{k}$值，并保持$\boldsymbol{\mu}_{k}$固定，关于$r_{nk} $最小化$J$(E步)
  $r_{n k}=\left\{\begin{array}{ll} 1 & \text { if } k=\arg \min _{j}\left\|\mathbf{x}_{n}-\boldsymbol{\mu}_{j}\right\|^{2} \\ 0 & \text { otherwise } \end{array}\right.$
• 保持$r_{nk}$固定，关于$\boldsymbol{\mu}_{k}$最小化$J$(M步)
  $\frac{\partial J}{\partial \boldsymbol{\mu}_{k}}=2 \sum_{n=1}^{N} r_{n k}\left(\mathbf{x}_{n}-\boldsymbol{\mu}_{k}\right)=0 \rightarrow \boldsymbol{\mu}_{k}=\frac{\sum_{n} r_{n k} \mathbf{x}_{n}}{\sum_{n} r_{n k}}$

K-means应用：图像分割和压缩

GMM模型

高斯分布

• D维高斯变量的高斯分布：
  $\mathcal{N}(\mathbf{x} \mid \boldsymbol{\mu}, \mathbf{\Sigma})=\frac{1}{(2 \pi)^{D / 2}} \frac{1}{|\mathbf{\Sigma}|^{1 / 2}} \exp \left\{-\frac{1}{2}(\mathbf{x}-\boldsymbol{\mu})^{\mathrm{T}} \boldsymbol{\Sigma}^{-1}(\mathbf{x}-\boldsymbol{\mu})\right\} \tag{1}$
  其中，$\boldsymbol{\mu} \in R^{D}$，高斯分布的均值向量，$\mathbf{\Sigma} \in R^{D \times D}$，高斯分布的协方差矩阵。

为什么使用高斯分布？1. 高斯分布在自然界的数据中广泛存在，2.中心极限定理：在适当条件下，大量相互独立随机变量的均值经适当标准化后依分布收敛于正太分布。

最大似然估计

假设随机变量 $X$ 服从分布 $p(X \mid \boldsymbol{\theta})$, 即 $X \sim p(X \mid \boldsymbol{\theta})$, 其中, $\boldsymbol{\theta}$ 为待估计的参数，如果可以获得N个互相独立的 $X$ 的采样点 $\left\{\mathbf{x}_{1}, \mathbf{x}_{2}, \ldots, \mathbf{x}_{N}\right\}$, 则似然函数的定义为
$p\left(\mathbf{x}_{1}, \mathbf{x}_{2}, \ldots, \mathbf{x}_{N}\right)=\prod_{n=1}^{N} p\left(\mathbf{x}_{n} \mid \boldsymbol{\theta}\right)$
在实际使用中，一般采用对数似然函数：
$\ln p\left(\mathbf{x}_{1}, \mathbf{x}_{2}, \ldots, \mathbf{x}_{N} \mid \boldsymbol{\theta}\right)=\sum_{n=1}^{N} \ln p\left(\mathbf{x}_{n} \mid \boldsymbol{\theta}\right)$
参数 $\boldsymbol{\theta}$ 的最大似然估计为：
$\boldsymbol{\theta}=\arg \max _{\boldsymbol{\theta}} \ln p\left(\mathbf{x}_{1}, \mathbf{x}_{2}, \ldots, \mathbf{x}_{N} \mid \boldsymbol{\theta}\right)$
高斯模型的最大似然估计

$\boldsymbol{\mu}_{M L}=\frac{1}{N} \sum_{n} \mathbf{x}_{n}, \boldsymbol{\Sigma}_{M L}=\frac{1}{N} \sum_{n}\left(\mathbf{x}_{n}-\boldsymbol{\mu}_{M L}\right)\left(\mathbf{x}_{n}-\boldsymbol{\mu}_{M L}\right)^{T}$

高斯混合分布
$p(\mathbf{x})=\sum_{k=1}^{K} \pi_{k} \mathcal{N} {\left(\mathbf{x} \mid \boldsymbol{\mu}_{k}, \mathbf{\Sigma}_{k}\right)} \tag{2}$

$\mathcal{N}(\mathbf{x} \mid \boldsymbol{\mu}, \mathbf{\Sigma})=\frac{1}{(2 \pi)^{D / 2}} \frac{1}{|\boldsymbol{\Sigma}|^{1 / 2}} \exp \left\{-\frac{1}{2}(\mathbf{x}-\boldsymbol{\mu})^{\mathrm{T}} \boldsymbol{\Sigma}^{-1}(\mathbf{x}-\boldsymbol{\mu})\right\}$

$0 \leqslant \pi_{k} \leqslant 1 \quad \sum_{k=1}^{K} \pi_{k}=1$
$\pi _{k}$，$\boldsymbol{\mu}_k$，$\boldsymbol{\Sigma}_k$为待估计参数。

• $\pi _{k}$的解释(直观理解为第k个高斯所占的比重)
  引入一个K维 one-hot(只有一维为1, 其余维度为0) 向量 $z=\left[z_{1}, \ldots, z_{k}, \ldots, z_{K}\right], z_{k} \in\{0,1\}, \Sigma_{k} z_{k}=1$, 概率 $P\left(z_{k}=1\right)$ 为 向量z的第k维为的先验概率 $p\left(z_{k}=1\right)=\pi_{k}$，向量z的分布可以表示为 $p(\mathbf{z})=\prod_{k=1}^{K} \pi_{k}^{z_{k}}$, 等价于 $\mathrm{p}(\mathrm{z})=\pi_{\mathrm{k}}$, where $\mathrm{z}_{\mathrm{k}}=1$
• 条件分布$\quad p\left(\mathbf{x} \mid z_{k}=1\right)=\mathcal{N}\left(\mathbf{x} \mid \boldsymbol{\mu}_{k}, \boldsymbol{\Sigma}_{k}\right) \quad p(\mathbf{x} \mid \mathbf{z})=\prod_{k=1}^{K} \mathcal{N}\left(\mathbf{x} \mid \boldsymbol{\mu}_{k}, \boldsymbol{\Sigma}_{k}\right)^{z_{k}} \tag{4}$
• 联合分布

$p(\boldsymbol{x}, \boldsymbol{z})=\prod_{k=1}^{K} \pi_{k}^{z_{k}} N\left(x \mid u_{k}, \Sigma_{k}\right)^{{z}_{k}}$

• 边缘分布

• 使用贝叶斯公式对潜变量求和：$p(\mathbf{x})=\sum_{\mathbf{z}} p(\mathbf{x}, \mathbf{z})=\sum_{\mathbf{z}} p(\mathbf{z}) p(\mathbf{x} \mid \mathbf{z})=\sum_{k=1}^{K} \pi_{k} \mathcal{N}\left(\mathbf{x} \mid \boldsymbol{\mu}_{k}, \mathbf{\Sigma}_{k}\right)$
• 对每一个观测$\mathbf{x}_n$，都有一个潜变量$\mathbf{z}_n$和其对应，上述公式将变量$\mathbf{x}$和潜变量$\mathbf{z}$联系起来，并且引入了联合分布$p(\boldsymbol{x}, \boldsymbol{z})$，完成了将观测数据的边缘分布转换成观测和潜变量的联合分布

• 后验分布
  $\begin{aligned} \gamma\left(z_{k}\right) \equiv p\left(z_{k}=1 \mid \mathbf{x}\right) &=\frac{p\left(z_{k}=1\right) p\left(\mathbf{x} \mid z_{k}=1\right)}{\sum_{j=1}^{K} p\left(z_{j}=1\right) p\left(\mathbf{x} \mid z_{j}=1\right)} \\ &=\frac{\pi_{k} \mathcal{N}\left(\mathbf{x} \mid \boldsymbol{\mu}_{k}, \boldsymbol{\Sigma}_{k}\right)}{\sum_{j=1}^{K} \pi_{j} \mathcal{N}\left(\mathbf{x} \mid \boldsymbol{\mu}_{j}, \boldsymbol{\Sigma}_{j}\right)} \end{aligned} \tag{3}$

$\gamma\left(z_{k}\right)$为得到观测$x$后，$z_{k}=1$的后验概率，理解为第k个高斯成分对于生成观测$x$的贡献值。

GMM对数似然函数
$\ln p(\mathbf{X} \mid \boldsymbol{\pi}, \boldsymbol{\mu}, \boldsymbol{\Sigma})=\sum_{n=1}^{N} \ln \left\{\sum_{k=1}^{K} \pi_{k} \mathcal{N}\left(\mathbf{x}_{n} \mid \boldsymbol{\mu}_{k}, \boldsymbol{\Sigma}_{k}\right)\right\}$
其中，$\mathbf{X}=\left[\begin{array}{c}\mathbf{x}_{1}^{T} \\ \vdots \\ \mathbf{x}_{N}^{T}\end{array}\right]$，给出潜变量矩阵定义$\mathbf{Z}=\left[\begin{array}{c}\mathbf{z}_{1}^{T} \\ \vdots \\ \mathbf{z}_{N}^{T}\end{array}\right]$。$\mathcal{N}\left(\mathbf{x}_{n} \mid \boldsymbol{\mu}_{k}, \boldsymbol{\Sigma}_{k}\right)$由公式(1)计算。

GMM模型参数估计的EM算法（最大似然准则）

计算似然函数$\ln p(\mathbf{X} \mid \boldsymbol{\pi}, \boldsymbol{\mu}, \boldsymbol{\Sigma})$分别对参数$\boldsymbol{\pi}, \boldsymbol{\mu}, \boldsymbol{\Sigma}$求导：

1. 对$\boldsymbol{\mu}$求导：
   $0=\sum_{n=1}^{N} \underbrace{\frac{\pi_{k} \mathcal{N}\left(\mathbf{x}_{n} \mid \boldsymbol{\mu}_{k}, \boldsymbol{\Sigma}_{k}\right)}{\sum_{j} \pi_{j} \mathcal{N}\left(\mathbf{x}_{n} \mid \boldsymbol{\mu}_{j}, \mathbf{\Sigma}_{j}\right)}}_{\gamma\left(z_{n k}\right)} \boldsymbol{\Sigma}_{k}^{-1}\left(\mathbf{x}_{n}-\boldsymbol{\mu}_{k}\right)$
   $\boldsymbol{\mu}_{k}=\frac{1}{N_{k}} \sum_{n=1}^{N} \gamma\left(z_{n k}\right) \mathbf{x}_{n} \quad N_{k}=\sum_{n=1}^{N} \gamma\left(z_{n k}\right)$
2. 对$\boldsymbol{\Sigma}$求导：
   $0=\sum_{n=1}^{N} \underbrace{\frac{\pi_{k} \mathcal{N}\left(\mathbf{x}_{n} \mid \boldsymbol{\mu}_{k}, \boldsymbol{\Sigma}_{k}\right)}{\sum_{j} \pi_{j} \mathcal{N}\left(\mathbf{x}_{n} \mid \boldsymbol{\mu}_{j}, \mathbf{\Sigma}_{j}\right)}}_{\gamma\left(z_{n k}\right)} \left(1-|\boldsymbol{\Sigma_k}|\right)$
   $\boldsymbol{\Sigma}_{k}=\frac{1}{N_{k}} \sum_{n=1}^{N} \gamma\left(z_{n k}\right)\left(\mathbf{x}_{n}-\boldsymbol{\mu}_{k}\right)\left(\mathbf{x}_{n}-\boldsymbol{\mu}_{k}\right)^{\mathrm{T}}$
3. 对$\boldsymbol{\pi}$求导：
   $\begin{array}{l} \ln p(\mathbf{X} \mid \boldsymbol{\pi}, \boldsymbol{\mu}, \mathbf{\Sigma})+\lambda\left(\sum_{k=1}^{K} \pi_{k}-1\right) \\ 0=\sum_{n=1}^{N} \frac{\mathcal{N}\left(\mathbf{x}_{n} \mid \boldsymbol{\mu}_{k}, \boldsymbol{\Sigma}_{k}\right)}{\sum_{j} \pi_{j} \mathcal{N}\left(\mathbf{x}_{n} \mid \boldsymbol{\mu}_{j}, \boldsymbol{\Sigma}_{j}\right)}+\lambda \\ \pi_{k}=\frac{N_{k}}{N} \end{array}$
   用拉格朗日乘子法，求解$\lambda$:
   $0=\sum_{n=1}^{N} \frac{\mathcal{N}\left(\mathbf{x}_{n} \mid \boldsymbol{\mu}_{k}, \boldsymbol{\Sigma}_{k}\right)}{\sum_{j} \pi_{j} \mathcal{N}\left(\mathbf{x}_{n} \mid \boldsymbol{\mu}_{j}, \boldsymbol{\Sigma}_{j}\right)}+\lambda$
   等号两边同时乘以 $\pi_{k}$, 并对k求和
   $\begin{array}{c} 0=\sum_{n=1}^{N} \sum_{k} \frac{\pi_{k} \mathcal{N}\left(\mathbf{x}_{n} \mid \boldsymbol{\mu}_{k}, \mathbf{\Sigma}_{k}\right)}{\sum_{j} \pi_{j} \mathcal{N}\left(\mathbf{x}_{n} \mid \boldsymbol{\mu}_{j}, \boldsymbol{\Sigma}_{j}\right)}+\lambda \sum_{k} \pi_{k} \\ 0=\sum_{n=1}^{N} \frac{\sum_{k} \pi_{k} \mathcal{N}\left(\mathbf{x}_{\mathrm{n}} \mid \boldsymbol{\mu}_{k}, \boldsymbol{\Sigma}_{k}\right)}{\sum_{j} \pi_{j} \mathcal{N}\left(\mathbf{x}_{n} \mid \boldsymbol{\mu}_{j}, \mathbf{\Sigma}_{j}\right)}+\lambda \\ \lambda=-N \end{array}$

GMM模型参数估计的EM算法总结

上述参数估计方法并不是一个严格的解析解，因为公式中有后验概率$\gamma\left(z_{n k}\right)$，依赖于每个高斯的代估计参数。上述推导过程给出了一个迭代的估计参数的过程，并能保证似然逐步增加。

给定一个GMM模型, 优化目标是寻找使似然函数最大的各个高斯成分的均值向量、协方差矩阵和混合系数

1. 初始化 初始化参数 $\boldsymbol{\mu}_{k}, \boldsymbol{\Sigma}_{k}, \pi_{k}$
2. E步 使用当前参数计算后验概率
   $\gamma\left(z_{n k}\right)=\frac{\pi_{k} \mathcal{N}\left(\mathbf{x}_{n} \mid \boldsymbol{\mu}_{k}, \boldsymbol{\Sigma}_{k}\right)}{\sum_{j} \pi_{j} \mathcal{N}\left(\mathbf{x}_{n} \mid \boldsymbol{\mu}_{j}, \boldsymbol{\Sigma}_{j}\right)}$
3. M步 使用后验重新估计参数
   $\begin{array}{c} \boldsymbol{\mu}_{k}^{n e w}=\frac{1}{N_{\mathrm{k}}} \sum_{n=1}^{N} \gamma\left(z_{n k}\right) \mathbf{x}_{n} \\ \Sigma_{k}^{n e w}=\frac{1}{N_{k}} \sum_{n=1}^{N} \gamma\left(z_{n k}\right)\left(\mathbf{x}_{n}-\boldsymbol{\mu}_{k}^{ {new }}\right)\left(\mathbf{x}_{\mathrm{n}}-\boldsymbol{\mu}_{k}^{ {new }}\right)^{T} \\ \pi_{k}^{ new }=\frac{N_{k}}{N}, \quad N_{k}=\sum_{n=1}^{N} \gamma\left(z_{n k}\right) \end{array}$
4. 重新计算似然函数, 重复2-4，直至满足收敘条件

GMM模型和K-means的联系

K-means可以看作GMM模型的一个特殊情况，假设公式$(2)$中，每个单高斯的分布都具有相同的协方差矩阵，并且有$\boldsymbol{\Sigma}=\epsilon \mathbf{I}$，$\mathbf{I}$是单位矩阵，高斯分布可以简化为：
$p\left(\mathbf{x} \mid \boldsymbol{\mu}_{k}, \mathbf{\Sigma}_{k}\right)=\frac{1}{(2 \pi \epsilon)^{\frac{1}{2}}} \exp \left\{-\frac{1}{2 \epsilon}\left\|\mathbf{x}-\boldsymbol{\mu}_{k}\right\|^{2}\right\}$
公式$(3)$中的后验概率变为：
$\gamma\left(z_{n k}\right)=\frac{\pi_{k} \exp \left\{-\frac{1}{2 \epsilon}\left\|\mathbf{x}_{\mathrm{n}}-\boldsymbol{\mu}_{k}\right\|^{2}\right\}}{\sum_{j} \pi_{j} \exp \left\{-\frac{1}{2 \epsilon}\left\|\mathbf{x}_{\mathrm{n}}-\boldsymbol{\mu}_{\mathrm{j}}\right\|^{2}\right\}}$
当$\epsilon \rightarrow 0$，$-\frac{1}{2 \epsilon}\left\|\mathbf{x}_{\mathrm{n}}-\boldsymbol{\mu}_{\mathrm{j}}\right\|^{2} \rightarrow-\infty, \quad \exp \left\{-\frac{1}{2 \epsilon}\left\|\mathbf{x}_{\mathrm{n}}-\boldsymbol{\mu}_{\mathrm{j}}\right\|^{2}\right\} \rightarrow 0$，对于分母中，假设第m项$\left\|\mathbf{x}_{\mathrm{n}}-\boldsymbol{\mu}_{\mathrm{m}}\right\|^{2}$最小，那么分母上j=m这一项将在$\epsilon \rightarrow 0$的时候以最慢的速度趋于0，因此只有分子上k=m时，$\gamma\left(z_{n k}\right) \rightarrow 1, k \neq m$，$\gamma\left(z_{n k}\right) \rightarrow 0$，显然$\gamma\left(z_{n k}\right) \rightarrow r_{nk}$。

K-means是一种硬对齐方式，某个数据点只能对应在某个类别上，GMM是一种软对齐方式，使用后验概率来表示某个数据点由某个类别产生的概率。

EM算法

上述GMM参数估计过程中已经使用到GMM算法，现在具体看一下EM算法。首先通过一个例子来体会。

问题：假设随机抽取100个男生和女生的身高数据，假设男生和女生的身高分布分别服从高斯分布$\mathcal{N}\left(\mathrm{x} \mid \mu_{M}, \Sigma_{M}\right)$和$\mathcal{N}\left(\mathrm{x} \mid \mu_{W}, \Sigma_{W}\right)$，请用最大似然法估计男生和女生身高分布的均值和方差。

这里分为两种情况：

情况1：已经知道每个数据对应的性别

男	女	女	男	女	女
178	175	170	175	168	169

$\mu_{M}=\frac{1}{\mathrm{~N}_{\mathrm{M}}} \sum_{m=1}^{N_{M}} x_{m}, \Sigma_{M}=\frac{1}{\mathrm{~N}_{\mathrm{M}}} \sum_{m=1}^{N_{M}}\left(x_{m}-\mu_{M}\right)$

其中$N_M$为男生总人数，$x_m$为每个男生的身高数据。同理，可求出女生的数据。

情况2：果冻大意了，只统计了身高数据

-	-	-	-	-	-
178	175	170	175	168	169

这种情况，我们引入一个新的变量$z$，$p(z_i = 1)$表示第i个数据为男生身高的概率，$p(z_i = 0)$表示第i个数据为女生的身高的概率。

$\mu_{M}=\frac{1}{\sum_{1}^{N} p\left(z_{i}=1\right)} \sum_{m=1}^{N_{M}} p\left(z_{i}=1\right) x_{m}$

如何获得$p(z_i)$？可以采用迭代估计法：

1. $\quad$先给出一组随机的参数取值 $\left(\mu_{M}, \Sigma_{M}, \mu_{W}, \Sigma_{W}\right)$
2. $\quad$更新 $p\left(z_{i}=1\right)=\frac{\mathcal{N}\left(\mathrm{x}_{i} \mid \mu_{M}, \Sigma_{M}\right)}{\mathcal{N}\left(\mathrm{x}_{i} \mid \mu_{M}, \Sigma_{M}\right)+\mathcal{N}\left(\mathrm{x}_{i} \mid \mu_{W}, \Sigma_{W}\right)}$
3. 使用新的p更新 $\left(\mu_{M}, \Sigma_{M}, \mu_{W}, \Sigma_{W}\right)$, 重复2-3

含隐变量模型的最大似然估计->EM算法，隐变量与模型参数相互影响，分两步一静一动交替迭代。通用步骤：E步(求期望)、M步(最大化)、重复E、M

给定完全数据 $\{\mathrm{X}, \mathrm{Z}\}$ 的联合概率分布 $\mathrm{p}(\mathrm{X}, \mathrm{Z} \mid \theta)$, 待学习参数 $\theta$, 优化的目标是寻找$\theta$来最大化似然函数p $(\mathrm{X} \mid \theta)$

1. 初始化 初始化参数 $\theta^{old}$
2. E步 计算潜变量的后验概率p(Z|X, $\theta^{{old}}$)
3. M步 使用后验重新估计参数

$Q\left(\theta, \theta^{{old }}\right)=\sum_{Z} \mathrm{p}\left(\mathrm{Z} \mid \mathrm{X}, \theta^{{old }}\right) \ln \mathrm{p}(\mathrm{X}, \mathrm{Z} \mid \theta)=\mathrm{E}_{\mathrm{Z} \sim \mathrm{p}\left(\mathrm{Z} \mid \mathrm{X}, \theta^{ {old }}\right)}[\ln \mathrm{p}(\mathrm{X}, \mathrm{Z} \mid \theta)]$
$\theta^{\text {new }}=\underset{\theta}{\arg \max } {Q}\left(\boldsymbol{\theta}, \boldsymbol{\theta}^{\text {old }}\right)$
4. 重算似然 重新计算似然函数，重复2-4, 更新参数 $\theta^{\text {old }} \leftarrow \theta^{\text {new }}$直至满足收敘条件

将不完全数据（只有观测数据）的边缘分布转换成容易处理的完全数据（观测数据+潜变量）的联合分布

深入理解EM算法

• EM算法的目标是寻找潜变量模型的最大似然解。后边我们假定待估计的参数用$\theta$
• 对数似然函数用完全数据的联合概率表示为：
  $\ln p(\mathbf{X} \mid \boldsymbol{\theta})=\ln \left\{\sum_{\mathbf{Z}} p(\mathbf{X}, \mathbf{Z} \mid \boldsymbol{\theta})\right\}$
  使用EM算法，一般认为完全数据的联合概率分布的似然$\ln p(\mathbf{X}, \mathbf{Z} \mid \boldsymbol{\theta})$容易计算。实际上，完全数据{X,Z}无法获取，但是潜变量Z的后验概率分布$p(\mathbf{Z} \mid \mathbf{X}, \boldsymbol{\theta})$可以进行估计。
• 计算完全数据的似然$\ln p(\mathbf{X}, \mathbf{Z} \mid \boldsymbol{\theta})$在$\mathbf{Z} \sim p\left(\mathbf{Z} \mid \mathbf{X}, \boldsymbol{\theta}^{\text {old }}\right)$时，关于变量Z的期望：

$Q\left(\theta, \theta^{{old }}\right)=\sum_{Z} \mathrm{p}\left(\mathrm{Z} \mid \mathrm{X}, \theta^{{old }}\right) \ln \mathrm{p}(\mathrm{X}, \mathrm{Z} \mid \theta)=\mathrm{E}_{\mathrm{Z} \sim \mathrm{p}\left(\mathrm{Z} \mid \mathrm{X}, \theta^{ {old }}\right)}[\ln \mathrm{p}(\mathrm{X}, \mathrm{Z} \mid \theta)]$

• 寻找使Q函数最大的新参数：
  $\theta^{\text {new }}=\underset{\theta}{\arg \max } {Q}\left(\boldsymbol{\theta}, \boldsymbol{\theta}^{\text {old }}\right)$

当尝试使用EM算法来解决自己的问题时，需要明确Q函数。

使用EM算法通用步骤重新考虑GMM参数估计

在之前的推导中，我们引入了潜变量Z，但并没用用到完全数据的联合概率分布，而是直接对不完全数据X的对数似然进行了求解。根据EM算法的通用步骤，首先考虑完全数据的似然函数：
$\begin{array}{c} p(\mathbf{X}, \mathbf{Z} \mid \boldsymbol{\mu}, \mathbf{\Sigma}, \boldsymbol{\pi})=\prod_{n=1}^{N} \prod_{k=1}^{K} \pi_{k}^{z_{n k}} \mathcal{N}\left(\mathbf{x}_{n} \mid \boldsymbol{\mu}_{k}, \boldsymbol{\Sigma}_{k}\right)^{z_{n k}} \\ \ln p(\mathbf{X}, \mathbf{Z} \mid \boldsymbol{\mu}, \boldsymbol{\Sigma}, \boldsymbol{\pi})=\sum_{n=1}^{N} \sum_{k=1}^{K} z_{n k}\left\{\ln \pi_{k}+\ln \mathcal{N}\left(\mathbf{x}_{n} \mid \boldsymbol{\mu}_{k}, \boldsymbol{\Sigma}_{k}\right)\right\} \end{array}$
根据公式(4)，计算Z的后验概率: $\color{red}{不理解如何推导可以跳过}$
$p(\mathbf{Z} \mid \mathbf{X}, \boldsymbol{\mu}, \boldsymbol{\Sigma}, \boldsymbol{\pi}) \propto \prod_{n=1}^{N} \prod_{k=1}^{K}\left[\pi_{k} \mathcal{N}\left(\mathbf{x}_{n} \mid \boldsymbol{\mu}_{k}, \boldsymbol{\Sigma}_{k}\right)\right]^{z_{n k}}$
完全数据的对数似然关于潜变量的期望值：
$\mathbb{E}_{\mathbf{Z}}[\ln p(\mathbf{X}, \mathbf{Z} \mid \boldsymbol{\mu}, \boldsymbol{\Sigma}, \boldsymbol{\pi})]=\sum_{n=1}^{N} \sum_{k=1}^{K} \gamma\left(z_{n k}\right)\left\{\ln \pi_{k}+\ln \mathcal{N}\left(\mathbf{x}_{n} \mid \boldsymbol{\mu}_{k}, \boldsymbol{\Sigma}_{k}\right)\right\}$

EM算法通用解释

EM算法的一般假设是直接优化观测数据的似然$p(X|\theta)$十分复杂，但是优化完全数据的似然$p(X,Z|\theta)$比较容易。引入一个关于变量Z的任意分布$q(Z)$:

为什么引用任意一个Z的分布？我们只知道Z服从某一个分布，但是并不知道具体服从什么分布，但是后面的推导会证明，无论Z的真实分布是什么都不影响推导。

$\begin{array}{l} \ln p(\mathbf{X} \mid \boldsymbol{\theta})=\sum_{Z} q(\mathbf{Z}) \ln p(\mathbf{X} \mid \boldsymbol{\theta})=\sum_{\mathbf{Z}} q(\mathbf{Z}) \ln \left\{\frac{p(\mathbf{X}, \mathbf{Z} \mid \boldsymbol{\theta})}{q(\mathbf{Z})} \frac{q(\mathbf{Z})}{p(\mathbf{Z} \mid \mathbf{X}, \theta)}\right\}\\ =\sum_{\mathbf{Z}} q(\mathbf{Z}) \ln \frac{p(\mathbf{X}, \mathbf{Z} \mid \boldsymbol{\theta})}{q(\mathbf{Z})}-\sum_{\mathbf{Z}} q(\mathbf{Z}) \ln \frac{p(\mathbf{Z} \mid \mathbf{X}, \boldsymbol{\theta})}{q(\mathbf{Z})} \end{array}$
令：
$\mathcal{L}(q, \boldsymbol{\theta})=\sum_{\mathbf{Z}} q(\mathbf{Z}) \ln \frac{p(\mathbf{X}, \mathbf{Z} \mid \boldsymbol{\theta})}{q(\mathbf{Z})} \\ \mathrm{KL}(q \| p)=-\sum_{\mathbf{Z}} q(\mathbf{Z}) \ln \frac{p(\mathbf{Z} \mid \mathbf{X}, \boldsymbol{\theta})}{q(\mathbf{Z})} \\$

KL散度是衡量两个概率分布之间差异的一个度量。

则：

$\ln p(\mathbf{X} \mid \boldsymbol{\theta})=\mathcal{L}(q, \boldsymbol{\theta})+\operatorname{KL}(q \| p)$

其中：$\mathcal{L}(q, \boldsymbol{\theta})$是$q(\mathbf{Z})$和$\boldsymbol{\theta}$的泛函，$\mathrm{KL}(q \| p) \geq 0$,当且仅当$q=p$时等号成立。

• $\mathrm{KL}(q \| p) \geq 0$，所以$\ln p(\mathbf{X} \mid \boldsymbol{\theta}) \geq \mathcal{L}(q, \boldsymbol{\theta})$，只有当后验分布$p(\mathbf{Z} \mid \mathbf{X}, \boldsymbol{\theta})$和$q(\mathbf{Z})$相等时，等号成立。
• $\mathcal{L}(q, \boldsymbol{\theta})$可以看作是$\ln p(\mathbf{X} \mid \boldsymbol{\theta})$的下界，如果我们无法直接提升$\ln p(\mathbf{X} \mid \boldsymbol{\theta})$的准确值，我们可以提升其下界。

E步：寻找使 $\mathcal{L}(q, \boldsymbol{\theta})$ 最大的 $q(\mathbf{Z})$, 当 $q(\mathbf{Z})=$$p(\mathbf{Z} \mid \mathbf{X}, \boldsymbol{\theta})$ 时， $\mathrm{KL}(q \| p)=0$, 此时 $\mathcal{L}(q, \boldsymbol{\theta})$

M步: 固定 $q(\mathbf{Z})$, 寻找使 $\mathcal{L}(q, \boldsymbol{\theta})$ 增加 的新参数 $\boldsymbol{\theta}^{\text {new }}$,因为参数更新， $q(\mathbf{Z})$ 和 $p\left(\mathbf{Z} \mid \mathbf{X}, \boldsymbol{\theta}^{\text {new }}\right)$ 不再相等，此时$\mathrm{KL}(q \| p)>0$, 因此导致$\ln p\left(\mathbf{X} \mid \boldsymbol{\theta}^{\text {new }}\right)$

$\mathcal{L}(\mathbf{q}, \boldsymbol{\theta})$ 与 $\boldsymbol{Q}\left(\boldsymbol{\theta}, \boldsymbol{\theta}^{\text {old }}\right)$

$q(\mathbf{Z})=$$p(\mathbf{Z} \mid \mathbf{X}, \boldsymbol{\theta}^{old})$ ，$\boldsymbol{Q}\left(\boldsymbol{\theta}, \boldsymbol{\theta}^{\text {old }}\right)$ 等价于$\mathcal{L}(\mathbf{q}, \boldsymbol{\theta})$，推导如下：

$\mathcal{L}(q, \boldsymbol{\theta})=\sum_{\mathbf{Z}} \boldsymbol{p}\left(\mathbf{Z} \mid \mathbf{X}, \boldsymbol{\theta}^{\boldsymbol{o l d}}\right) \ln \boldsymbol{p}(\mathbf{X}, \mathbf{Z} \mid \boldsymbol{\theta})-\sum_{\mathbf{Z}} \boldsymbol{p}\left(\mathbf{Z} \mid \mathbf{X}, \boldsymbol{\theta}^{\boldsymbol{o l d}}\right) \ln \boldsymbol{p}\left(\mathbf{Z} \mid \mathbf{X}, \boldsymbol{\theta}^{\text {old }}\right) \\ =\boldsymbol{Q}\left(\boldsymbol{\theta}, \boldsymbol{\theta}^{\text {old }}\right)+\mathrm{const}$
因此，EM算法中，E步所计算的$\boldsymbol{Q}\left(\boldsymbol{\theta}, \boldsymbol{\theta}^{\text {old }}\right)$ 实际上等价于计算$\mathcal{L}(\mathbf{q}, \boldsymbol{\theta})$。

python代码实现

使用GMM-EM实现0-9孤立词识别，介绍GMM对数似然函数、E步、M步的实现，详细参考：

GMM的对数似然函数：
公式：
$\ln p(\mathbf{X} \mid \boldsymbol{\pi}, \boldsymbol{\mu}, \boldsymbol{\Sigma})=\sum_{n=1}^{N} \ln \left\{\sum_{k=1}^{K} \pi_{k} \mathcal{N}\left(\mathbf{x}_{n} \mid \boldsymbol{\mu}_{k}, \boldsymbol{\Sigma}_{k}\right)\right\}$
$\mathcal{N}(\mathbf{x} \mid \boldsymbol{\mu}, \mathbf{\Sigma})=\frac{1}{(2 \pi)^{D / 2}} \frac{1}{|\mathbf{\Sigma}|^{1 / 2}} \exp \left\{-\frac{1}{2}(\mathbf{x}-\boldsymbol{\mu})^{\mathrm{T}} \boldsymbol{\Sigma}^{-1}(\mathbf{x}-\boldsymbol{\mu})\right\} \tag{1}$

python代码：

def calc_log_likelihood(self , X):
        """Calculate log likelihood of GMM

            param: X: A matrix including data samples, num_samples * D
            return: log likelihood of current model 
        """

        log_llh = 0.0
        N = X.shape[0]
        log_llh = np.sum([np.log(np.sum([self.pi[k] * self.gaussian(X[n], self.mu[k], self.sigma[k]) 
                    for k in range(self.K)])) for n in range(N)])
        return log_llh

E步：

公式：

$\gamma\left(z_{n k}\right)=\frac{\pi_{k} \mathcal{N}\left(\mathbf{x}_{n} \mid \boldsymbol{\mu}_{k}, \boldsymbol{\Sigma}_{k}\right)}{\sum_{j} \pi_{j} \mathcal{N}\left(\mathbf{x}_{n} \mid \boldsymbol{\mu}_{j}, \boldsymbol{\Sigma}_{j}\right)}$

N = X.shape[0]
gamma = [np.zeros(self.K) for i in range(N)] # 后验概率

# E-step
for n in range(N):
    post_prob = [self.pi[k] * self.gaussian(X[n], self.mu[k], self.sigma[k]) for k in range(self.K)]
    post_prob = np.array(post_prob)
    post_prob_sum = np.sum(post_prob)
    gamma[n] = post_prob / post_prob_sum

M步:

公式：

$\begin{array}{c} \boldsymbol{\mu}_{k}^{n e w}=\frac{1}{N_{\mathrm{k}}} \sum_{n=1}^{N} \gamma\left(z_{n k}\right) \mathbf{x}_{n} \\ \Sigma_{k}^{n e w}=\frac{1}{N_{k}} \sum_{n=1}^{N} \gamma\left(z_{n k}\right)\left(\mathbf{x}_{n}-\boldsymbol{\mu}_{k}^{ {new }}\right)\left(\mathbf{x}_{\mathrm{n}}-\boldsymbol{\mu}_{k}^{ {new }}\right)^{T} \\ \pi_{k}^{ new }=\frac{N_{k}}{N}, \quad N_{k}=\sum_{n=1}^{N} \gamma\left(z_{n k}\right) \end{array}$

# M-step
for k in range(self.K):
    Nk = np.sum([gamma[n][k] for n in range(N)])
    if Nk == 0:
        continue
    self.pi[k] = Nk / N
    self.mu[k] = 1.0 / Nk * np.sum([gamma[n][k] * X[n] for n in range(N)], axis=0)
    diffs = X - self.mu[k]
    self.sigma[k] = 1.0 / Nk * np.sum([gamma[n][k] * diffs[n].reshape(self.dim, 1) * diffs[n] 
                    for n in range(N)], axis=0)
                                            
    log_llh = self.calc_log_likelihood(X)