最大似然估计与EM算法

最大似然估计(1922 by Filch)

例1 估计两点分布 $b(1,p)$ 的参数 $p$. 其分布列为 $P(X = x) = p^x(1-p)^{1-x},x=0,1.$

解 抽取n个样本$x_1,\cdots,x_n$, 这批观测值发生的概率为

$$P(X_1=x_1,\cdots,X_n=x_n;p) \overset{\text{独立性}}{=} \prod_{i=1}^n p^{x_i}(1-p)^{1-x_i} = p^{\sum_{i=1}^n x_i}(1-p)^{n-\sum_{i=1}^n x_i} \quad(eq1)$$

Note1: 注意这里是n个两点分布$b(1,p)$, 而不是1个二项分布$b(n,p)$. 因为我们观测的是n次抽取的结果, 即n个$b(1,p)$, 而不是整体n个产品中的合格数, 即1个$b(n,p)$. 两者计算过程中相差一个二项系数$C_n^k, k=\sum_{i=1}^n x_i$.

Note2: 对于事件$X_i=x_i$, $p^{x_i}(1-p)^{1-x_i}$为其发生的概率, 因为

$$P(X_i = x_i) = p^{x_i}(1-p)^{1-x_i} = \begin{cases} p^1(1-p)^0=p, & x_i = 1;\\ p^0(1-p)^1=1-p, & x_i = 0. \end{cases} \notag$$

由于$p$是未知的, 根据最大似然原理, 我们应该选择$p$使得$eq1$式中的概率尽可能大. 将$eq1$式看作位置参数$p$的函数, 用$L(p)$表示, 称作似然函数

$$L(p) = p^{\sum_{i=1}^n x_i}(1-p)^{n-\sum_{i=1}^n x_i} \quad(eq2)$$

求似然函数$eq2$式的最大值并不难, 它是一个幂函数的形式, 对其关于$p$求导令其为0, 得到下面方程, 称其为似然方程:

令$\sum_{i=1}^n(x_i)=k,$ 于是$L(p) = p^k(1-p)^{n-k}$

$$\begin{aligned} \frac{\partial L(p)}{\partial p} &= kp^{k-1}(1-p)^{n-k}-(n-k)p^k(1-p)^{n-k-1} = 0\\ \Leftrightarrow k(1-p) &= (n-k)p \\ \Leftrightarrow p &= k/n = \frac{1}{n}\sum_{i=1}^n x_i \qquad(\text{恰好是样本均值}) \end{aligned}$$

除了上面的直接求导还有另外一种常用的方式, 对$eq2$式两边取对数再求导, 这种方式计算导数要更简单

$$\begin{aligned} \ln L(p) &= k\ln(p) + (n-k)\ln(1-p), \qquad k=\sum_{i=1}^n x_i\\ \text{两边对 }p\text{ 求导}\rightarrow \frac{\partial \ln L(p)}{\partial p} &= \frac{k}{p} - \frac{n-k}{1-p} = 0 \\ \text{求解得到: } p &= k/n = \frac{1}{n}\sum_{i=1}^n x_i \end{aligned}$$

定义 设总体的概率函数为$p(x;\theta), \theta\in\Theta$, 其中$\theta$是一个未知参数或几个未知参数组成的参数向量, $\Theta$是参数空间, $x_1,\cdots,x_n$是来自该总体的独立的样本, 将样本的联合概率函数看成$\theta$的函数, 用$L(\theta;x_1,\cdots,x_n)$表示, 简记为$L(\theta)$,

$$L(\theta) = L(\theta;x_1,\cdots,x_n) = p(x_1;\theta)p(x_2;\theta)\cdots p(x_n;\theta),$$

$L(\theta)$称为样本的似然函数. 如果某统计量$\hat{\theta}=\hat{\theta}(x_1,\cdots,x_n)$满足

$$L(\hat{\theta}) = \max_{\theta\in\Theta} L(\theta),$$

则称$\hat{\theta}$是$\theta$的最大似然估计, 简记为MLE (Maximum Likelihood Estimate).

由于$\ln x$ 是 $x$ 的单调增函数, 因此, 使对数似然函数 $\ln L(\theta)$ 达到最大与使 $L(\theta)$ 达到最大是等价的. 当 $L(\theta)$ 可微时, 求导是求最大似然估计的常用方法, 而对数似然函数的求导更加简单, 所以大家通常更习惯于由 $\ln L(\theta)$ 出发寻找 $\theta$ 的最大似然估计.

Note: 从最大似然估计的定义可以看出, 若 $L(\theta)$ 与联合概率函数相差一个与 $\theta$ 无关的比例因子, 不会影响最大似然估计, 因此, 可以在$L(\theta)$ 中剔去与 $\theta$ 无关的因子

EM算法(1977 by Dempster)

MLE是一种非常有效的参数估计方法, 但当分布中有多余参数或者数据为截尾或者缺失时, 求其MLE是比较困难的. 于是Dempster等人于1977年提出了EM算法, 其出发点是要把求MLE的过程分两步走, 第一步求期望, 以便把多余的部分去掉, 第二步求极大值.

例2 设有多项分布 $M(1,\frac{1}{2}-\frac{\theta}{4},\frac{1-\theta}{4},\frac{1+\theta}{4},\frac{\theta}{4}),\theta\in(0,1)$, 现进行了197次试验, 四种结果出现的次数分别为75, 18, 70, 34. 试估计参数 $\theta$.

解 以 $y_1,y_2,y_3,y_4$表示四种结果发生的次数, 多项分布的似然函数为:

$$\begin{split} L(\theta;y) &\propto (\frac{1}{2}-\frac{\theta}{4})^{y_1}(\frac{1-\theta}{4})^{y_2}(\frac{1+\theta}{4})^{y_3}(\frac{\theta}{4})^{y_4}\\ &\propto (2-\theta)^{y_1}(1-\theta)^{y_2}(1+\theta)^{y_3}\theta^{y_4}. \end{split}$$

其对数似然函数为一个关于 $\theta$ 的三次多项式

$$\ln L\propto(y_1+y_2+y_3+y_4)\theta^3 +(-y_2-3y_3-2y_4)\theta^2+(-y_1+2y_2+2y_3-y_4)\theta+2y_4 \quad(eq8)$$

通过它求 $\theta$ 的MLE就是要求 $eq8$ 式在 $(0,1)$ 上的零点.

​我们通过引入2个变量 $z_1,z_2$ 后, 使得求解变得比较容易. 现假设第一种结果可以分成两部分, 其发生的概率分别为 $\frac{1-\theta}{4}$ 和 $\frac{1}{4}$, 令 $z_1$ 和 $y_1-z_1$ 分别表示落入这两部分的次数. 再假设第三种结果也分成两部分, 其发生的概率分别为 $\frac{\theta}{4}$ 和 $\frac{1}{4}$, 令 $z_2$ 和 $y_3-z_2$ 分别表示落入这两部分的次数. 由于 $z_1,z_2$ 是人为引入的, 是不可观测的, 我们称它们为 潜变量(latent variable). 称数据 $(y,z)$ 为完全数据(complete data), 而观测到的数据 $y$ 称为不完全数据.

​ 引入潜变量 $z_1,z_2$ 后, 此时, 完全数据的似然函数为

$$\begin{split} L(\theta;y,z) &\propto (\frac{1}{4})^{y_1-z_1}(\frac{1-\theta}{4})^{z_1+y_2}(\frac{1}{4})^{y_3-z_2}(\frac{\theta}{4})^{z_2+y_4}\\ &\propto \theta^{z_2+y_4}(1-\theta)^{z_1+y_2}, \end{split}$$

其对数似然函数为

$$\ln L(\theta;y,z) = (z_2+y_4)\ln \theta + (z_1+y_2)\ln(1-\theta).$$

如果 $(y,z)$ 均已知, 则由上式很容易求得 $\theta$ 的MLE, 但是很遗憾我们仅知道 $y$, 而不知道 $z$ 的值. 虽然我们不能知道 $z$ 的值, 但是当 $y$ 和 $\theta$ 已知时, 我们可以知道 $z$ 的分布

$$z_1\sim b\left(y_1,\frac{1-\theta}{2-\theta}\right),\qquad z_2\sim b\left(y_3,\frac{\theta}{1+\theta}\right)$$

Note: $P(Z_1|Y_1) \overset{Z_1\subset Y_1}{=} \frac{P(Z_1)}{P(Y_1)} = \frac{(1-\theta)/4}{(2-\theta)/4} = \frac{1-\theta}{2-\theta},\quad P(Z_2|Y_3) \overset{Z_2\subset Y_3}{=} \frac{P(Z_2)}{P(Y_3)} = \frac{\theta/4}{(1+\theta)/4} = \frac{\theta}{1+\theta}$

于是, Dempster等人建议分两步迭代求解 $\theta$

E步: 在已有观测数据 $y$ 及第 $i$ 步 $\theta$ 的估计值 $\theta^{(i)}$ 的条件下, 求 $\ln L(\theta;y,z)$ 关于 $z$ 的期望

$$E_z \left[\ln L(\theta;y,z)\right] ,\qquad z\sim F(z|y,\theta^{(i)}) \quad (eq11)$$

Note1: 这一步的目的是把其中关于 $z$ 的部分积分掉.

Note2: 求 $z$ 的期望需要知道 $z$ 的分布, 而 $z$ 的分布是与 $y$ 和 $\theta$ 相关的, 其中 $y$ 是可观测的已知的, 而 $\theta$ 是未知的, 所以这里采用第 $i$ 步的估计值 $\theta^{(i)}$ 做为 $\theta$ 的值.

M步: 求 $eq11$ 式关于 $\theta$ 的最大值 $\theta^{(i+1)}$, 即找 $\theta^{(i+1)}$ 使得

$$\theta^{(i+1)} = \underset{\theta}{\arg\max} \left\{E_z \left[\ln L(\theta;y,z)\right]\right\}. \quad (eq12)$$

这样就完成了由 $\theta^{(i)}$ 到 $\theta^{(i+1)}$ 的一次迭代. 重复 $eq11$ 式和 $eq12$ 式, 直至收敛, 即可得到 $\theta$ 的 MLE 的一个数值解.

​对于本例, 其 E 步为

$$\begin{split} E_z\ln L &= E\Big[(z_2+y_4)\ln \theta + (z_1+y_2)\ln(1-\theta)\Big]\\ &= (E\big[z_2\big|y_3,\theta^{(i)}\big]+y_4)\ln\theta + (E\big[z_1\big|y_1,\theta^{(i)}\big]+y_2)\ln(1-\theta)\\ &= \left(\frac{\theta^{(i)}}{1+\theta^{(i)}}y_3+y_4\right)\ln\theta + \left(\frac{1-\theta^{(i)}}{2-\theta^{(i)}}y_1+y_2\right)\ln(1-\theta) \end{split}$$

M 步为对 $E_z\ln L$ 关于 $\theta$ 求导, 并令其等于 $0$ , 得到

$$\frac{\frac{\theta^{(i)}}{1+\theta^{(i)}}y_3+y_4}{\theta} - \frac{\frac{1-\theta^{(i)}}{2-\theta^{(i)}}y_1+y_2}{1-\theta} = 0$$

解得如下迭代公式

$$\theta^{(i+1)} = \frac{\frac{\theta^{(i)}}{1+\theta^{(i)}}y_3+y_4}{\frac{\theta^{(i)}}{1+\theta^{(i)}}y_3+y_4+\frac{1-\theta^{(i)}}{2-\theta^{(i)}}y_1+y_2} \quad (eq16)$$

实际上 $eq16$ 式等号右边就是 $eq8$ 式的一个不动点迭代函数. $eq16$ 式可以整理为下面的形式

$$\theta = \frac{(y_3+y_4)\theta^2-(2y_3+y_4)\theta-2y_4}{(y_1+y_2+y_3+y_4)\theta^2-(y_2+2y_3+y_4)\theta-(y_1+2y_2+2y_4)}\notag$$

最后使用题设的数据75, 18, 70, 34计算 $\theta$ 的MLE的近似值

$i$	0	1	2	3	4	5	6	7	8	9	10
$\theta^{i}$	0.5	0.571 492	0.594 816	0.602 681	0.605 357	0.6067 271	0.606 584	0.606 697	0.606 728	0.606 746	0.606 744