数理统计：极大似然估计与EM算法

       极大似然估计（MLE）是统计学习，机器学习中常用的参数估计方法。当然有时候MLE会碰到复杂情况，比如比如数据不完整的时候，或者参数没有闭合解且情况特别复杂时，这时不直接求MLE，而采用确保收敛到MLE的办法，也就是期望极大化准则（EM准则），俗称EM算法。

一 极大似然估计（MLE）

       已知随机变量$X$的密度表达式，其中包含来自$\Theta$的位置参数$\theta$,密度函数可写为$f(x;\theta),\theta \in \Theta$是对应函数族$\mathscr P = \{ f(x;\theta) \mid \theta \in \Theta\}$中的一个函数。 令随机样本$\underline x=(x_1,...,x_n)^T$是来自密度函数$f(x;\theta)$的随机样本$\underline X=(X_1,...,X_n)^T$的观测值。我们希望用极大似然法估计$\theta$的函数$g(\theta)$.

       极大似然估计的基本思想是:在随机采样实验中，小概率事件基本不发生，而已经发生的应该是大概率事件。假定随机样本之间相互独立且分布相同;则$\underline x$发生的概率是$L(\theta;\underline x)= \mathop{\prod}_{i=1}^{n}f(x_i;\theta)$，称为似然函数。

       极大似然估计就是要最大化已经发生的事件的概率$L(\theta;\underline x)$;另外，由于自然对数函数是单调递增的，所以最大化对数似然函数$l(\theta;\underline x)=ln(L(\theta;\underline x))$等价于最大化似然函数$L(\theta;\underline x)$。

1.1 极大似然估计原理

       若样本来自于概率分布$f(x;\theta^*),\theta^* \in \Theta$,则有：

$$P_{\theta}[L(\theta^*;\underline x) \ge L(\theta;\underline x)] \to1;n \to \infty \qquad (1)$$

证明：

$$\begin{array}{left} P_{\theta}[L(\theta^*;\underline X) \ge L(\theta;\underline X)] &= P_{\theta}\bigg[\frac{L(\theta;\underline X) }{ L(\theta^*;\underline X)} \le 1\bigg] \\ & = P_{\theta}\bigg[\frac{1}{n}\sum_{i=1}^{n} ln\bigg(\frac{f(X_i;\theta) }{ f(X_i;\theta^*)}\bigg) \le 0\bigg] \\ \end{array}$$
由强大数定理和Jessen不等式有：
$$\frac{1}{n}\sum_{i=1}^{n} ln\bigg(\frac{f(X_i;\theta) }{ f(X_i;\theta^*)}\bigg) \overset{a.s.}{\rightarrow} E_{\theta}\Bigg(ln\bigg(\frac{f(X_i;\theta) }{ f(X_i;\theta^*)}\bigg) \Bigg) \le lnE_{\theta}\bigg(\frac{f(X_i;\theta) }{ f(X_i;\theta^*)}\bigg)=0$$
得证。

1.2 例子：贝努利分布的MLE

       已知随机样本随机样本$\underline X=(X_1,...,X_n)^T$独立同分布与$Bernoulli(\theta)$，则对数似然函数为：

$$l(\theta;\underline x)=(\sum_{i=1}^nx_i)ln(\theta)+(n-\sum_{i=1}^nx_i)ln(1-\theta)$$
一阶条件为：
$$\frac{\partial l(\theta;\underline x)}{\partial \theta}=\frac{\sum_{i=1}^nx_i}{\theta}+\frac{n-\sum_{i=1}^nx_i}{1-\theta}=0$$
求解得到极大似然估计量$\hat\theta_{ML}=\frac{\sum_{i=1}^nx_i}{n}=\bar x$.

二 EM算法

       有些时候，可能存在一些隐藏随机变量无法直接观测（隐变量），这时候观测不到完整的数据。在求解不完整数据比较复杂，而求解完整数据比较简单的时候，可以采用EM准则。EM准则不直接求出MLE，而是确保可以收敛到MLE。EM准则经典模型有高斯混合模型，另一个经典的模型是受限玻尔兹曼机（RBM）

2.1 EM算法描述

       设可观测的不完整变量为$\underline V=(V_1,...,V_n)^T$，隐藏的变量为$\underline H=(H_1,...,H_m)^T$，完整数据为$\underline X=(\underline V,\underline H)$；$\underline X;\underline V$的密度函数关系如下：

$$f_{\underline v}(\underline v;\theta)=\int f_{\underline X}(\underline v,\underline h;\theta)d\underline h$$

       EM 准则的目标是基于不完整的数据的似然函数$L(\theta;\underline V)(=f_{\underline v}(\underline v;\theta))$求出$\theta$的MLE:$\hat\theta=\hat\theta(\underline V)$.完整的数据的似然函数是$L(\theta;\underline X)=f_{\underline X}(\underline X;\theta)$.
       求解$\hat\theta$分两步，从初始化$\theta^{(0)}$开始：

1. E步：计算期望条件
   $$Q(\theta,\theta^{(0)})=Q(\theta \mid \underline V,\theta^{(0)})\equiv E[lnL(\theta;\underline X)\mid \underline V,\theta^{(0)}] \tag 2$$
2. M 步，极大化$Q(\theta)$.得到更新的$\theta^{(1)}$；
3. 回到E步，$\theta^{(0)}=\theta^{(1)}$，继续迭代直到收敛。

        进一步解释：

E 步：$Q(\theta,\theta^{(0)})$的积分表示如下：

$$E[lnL(\theta;\underline X)\mid \underline V,\theta^{(0)}]\equiv\int f(h\mid v,\theta^{(0)})lnL(\theta;\underline X)\big|_{v}dh$$
M步有： $$Q(\theta^{(j+1)},\theta^{(j)})= \underset{\theta \in \Theta}{\max }Q(\theta,\theta^{(j)})\ge Q(\theta,\theta^{(j)}) \tag 3$$

2.2 EM算法性质

2.2.1 Jenson不等式

离散形式
$\Omega$是有限集合$\{x_1,x_2,...x_n\}$,$f(x),g(x)$是$\Omega$上的两个概率测度，$\phi:R\to R$是凹函数，则有：

$$\phi \bigg(\sum_{i=1}^{n}f(x_i)g(x_i)\bigg)\ge \sum_{i=1}^{n} \bigg(\phi\big(f(x_i)\big)g(x_i)\bigg) \tag 4$$
等式成立的充要条件是：$f(x_i)\equiv c,\forall i\in{1...n}$

---

推论1：在给定概率分布$f(x)$的条件下，任意概率分布$g(x)$，使得

$$\sum_{i=1}^{n} \Bigg(\phi\bigg(\frac{f(x_i)}{g(x_i)}\bigg)g(x_i)\Bigg)\tag 5$$ 取得最大值的充要条件是：
$$f(x_i)=g(x_i)\;,\forall i \in (1...n)$$
证明：根据Jenson不等式，有
$$\sum_{i=1}^{n} \Bigg(\phi\bigg(\frac{f(x_i)}{g(x_i)}\bigg)g(x_i)\Bigg)\le \phi \Bigg(\sum_{i=1}^{n}\bigg(\frac{f(x_i)}{g(x_i)}g(x_i)\bigg)\Bigg)=\phi(1)$$
可知(5)式的极大值是$\phi(1)$，等式成立的充要条件是：$\frac{f(x_i)}{g(x_i)}=c$
并且有$\sum_i(f(x_i)-cg(x_i))=0 \Rightarrow c=1$
所以：$f(x_i)=g(x_i)\;,\forall i \in (1...n)$
证毕。

---

连续形式
$\Omega \subseteq R^n$,$f(x),g(x)$是$\Omega$上的两个概率测度，$\phi:R\to R$是凹函数，则有：

$$\phi \bigg(\int_{x\in \Omega}f(x)g(x)dx\bigg)\ge \int_{x\in \Omega} \phi\big(f(x)\big)g(x)dx \tag {4.1}$$
等式成立的充要条件是：$f(x)\equiv c,\forall{x\in \Omega}$

---

定理 2：在给定概率密度$f(x)$的条件下，任意概率密度$g(x)$，$\phi:R\to R$是凹函数，使得

$$J(g(x))\equiv\int_{x\in \Omega}\phi\bigg(\frac{f(x)}{g(x)}\bigg)g(x)dx\tag {5.1}$$ 取得最大值的充要条件是：
$$f(x)=g(x)\;,\forall{x\in \Omega}$$

2.2.2 可观察变量的似然函数

$$\begin{array}{left} ln L(\theta;\underline V) &=ln\int f_{\underline X}(\underline v,\underline h;\theta)d\underline h\\ &=ln\int f_{\underline H}(\underline h\mid v;\theta)\frac{f_{\underline X}(\underline v \underline h;\theta)}{ f_{\underline H}(\underline h\mid v;\theta)}d\underline h\\ &\ge \int f_{\underline H}(\underline h\mid v;\theta)ln\bigg(\frac{f_{\underline X}(\underline v \underline h;\theta)}{ f_{\underline H}(\underline h\mid v;\theta)}\bigg)d\underline h\\ \end{array} \tag 6$$
       EM 算法不直最大化对数似然函数$ln L(\theta;\underline V)$，而是不断最大化它的下界：
$$l(\theta)=\int f_{\underline H}(\underline h\mid v;\theta)ln\bigg(\frac{f_{\underline X}(\underline v \underline h;\theta)}{f_{\underline H}(\underline h\mid v;\theta)}\bigg)d\underline h \tag7$$

2.2.3 E步的作用
       假设若干步迭代后，得到$\theta^{(j)}$计算

$$Q(\theta,\theta^{(j)})=E[lnL(\theta;\underline X)\mid \underline V,\theta^{(j)}]\equiv\int f(h\mid v,\theta^{(j)})lnL(\theta;\underline X)\big|_{v}dh$$
       实际上有先要计算$f(h\mid v,\theta^{(j)})$,然后跟下式：
$$\begin{align} L(\theta^{(j)};\underline X)\big|_{v}&=f_{\underline X}(\underline v, \underline h;\theta^{(j)})\big|_{v}\\ &=f(h\mid v,\theta^{(j)})*f( v,\theta^{(j)})\big|_{v}\\[2ex] \end{align}$$
得到：
$$\begin{align} \frac{f_{\underline X}(\underline v, \underline h;\theta^{(j)})}{f(h\mid v,\theta^{(j)})}\bigg|_{v}&=f( v,\theta^{(j)})\bigg|_{v}\equiv c\tag8\\ \end{align}$$
       可以看出，E步根据$\theta^{(j)}$更新$Q(\theta,\theta^{(j)})$,实际上同时更新隐变量在给定隐变量的密度函数：$f(h\mid v,\theta^{(j)})\mid_{\underline v}$
根据jesson不等式和(8)式，实际是在给定$f_{\underline X}(\underline v, \underline h;\theta^{(j)})\big|_{v}$的条件下，极大化$l(\theta^{(j)})= L(\theta^{(j)};\underline V)$
       因此有：
$$\begin{align} L(\theta^{(j)};\underline V)&=l(\theta^{(j)})\\ &=\int f_{\underline H}(\underline h\mid v;\theta^{(j)})ln\bigg(\frac{f_{\underline X}(\underline v, \underline h;\theta^{(j)})}{f_{\underline H}(\underline h\mid v;\theta^{(j)})}\bigg)d\underline h\\ &\ge\int f_{\underline H}(\underline h\mid v;\theta^{(j-1)})ln\bigg(\frac{f_{\underline X}(\underline v, \underline h;\theta^{(j)})}{f_{\underline H}(\underline h\mid v;\theta^{(j-1)})}\bigg)d\underline h\\ &=Q(\theta^{(j)},\theta^{(j-1)})-\int f_{\underline H}(\underline h\mid v;\theta^{(j-1)})ln\big({f_{\underline H}(\underline h\mid v;\theta^{(j-1)})}\big)d\underline h\\\tag 9 \end{align}$$

2.2.4 M步的作用

       M步是极大化$Q(\theta,\theta^{(j-1)})$得到$\theta^{(j)}$，即：

$$\theta^{(j)}=\arg \underset{\theta\in \Theta}\max{Q(\theta,\theta^{(j-1)})}$$
       因此有：
$$\;Q(\theta^{(j)},\theta^{(j-1)})\ge Q(\theta,\theta^{(j-1)}) \;,\forall \theta\in \Theta\\[1em] \Rightarrow Q(\theta^{(j)},\theta^{(j-1)})\ge Q(\theta^{(j-1)},\theta^{(j-1)})$$
       结合(9)式有：
$$\begin{align} L(\theta^{(j)};\underline V) &\ge Q(\theta^{(j-1)},\theta^{(j-1)})-\int f_{\underline H}(\underline h\mid v;\theta^{(j-1)})ln\big({f_{\underline H}(\underline h\mid v;\theta^{(j-1)})}\big)d\underline h\\ &=\int f_{\underline H}(\underline h\mid v;\theta^{(j-1)})ln\bigg(\frac{f_{\underline X}(\underline v, \underline h;\theta^{(j-1)})}{f_{\underline H}(\underline h\mid v;\theta^{(j-1)})}\bigg)d\underline h\\ &=l(\theta^{(j-1)})\\ &=L(\theta^{(j-1)};\underline V)\\ \tag {10} \end{align}$$

       由此可以看出，EM每完成一次迭代，都有：$L(\theta^{(j)};\underline V)\ge L(\theta^{(j-1)};\underline V)$，EM算法通过迭不断增大$L(\theta^{(j)};\underline V)$的下界$l(\theta)$来估计参数$\hat\theta$;并不能保证$L(\hat\theta;\underline V)$达到极大值。