EM算法的推出

考虑观测数据$Y=\{y_1, y_2, \dots, y_m\}$,其中不可观测数据为$Z=\{z_1, z_2, \dots, z_k\}$,需要估计的参数为$\theta=\{\theta_1, \theta_2, \dots, \theta_t\}$。$Z$可以是离散或连续型随机变量,以下过程中假设$Z$为离散型($Z$为连续型,则全概率公式由求和改为积分)。则观测数据的对数似然函数为:

为了使$L(\theta)$最大,很容易想到对$\theta$偏导求极值。但有一个困难,即$Z$不可观测,导致的直接后果是对数套求和,计算难度增加。直觉是,否能找到一个方法,将求和放在对数外面?一个常用的技巧是转化为不等式,Jensen’s inequality描述了积分的函数与函数的积分的关系。由于$\log(x)$是凹函数,尝试考察某一次特定的$\theta^{(n)}$取值后,$L(\theta)$与$L(\theta^{(n)})$的差:

从$\eqref{eq:2}$可以得到$L(\theta)$的一个下限:

当$\theta=\theta^{(n)}$时,等号成立。有意思的是,由于$\theta^{(n)}$已知,$\eqref{eq:3}$的右半部分极大化要比$\eqref{eq:1}$简单很多,尽管右边的极值是一个下限,此时的$\theta$的估计值为:

EM算法的核心是通过迭代$L(\theta)$的下限极大化,进而逼近极值。EM算法每次迭代,都会靠近$L(\theta)$,但是最终可能得到局部最优。因此,EM算法对初始值敏感,可以多选几个初始值,对比结果。

双硬币模型

参考资料2中使用EM算法解决双硬币模型:

每次从两个硬币中选择一个,抛$N$次并记录正反面,同样的“选硬币-抛$N$次”进行$m$次,估计两个硬币正面出现的概率。观测数据$Y=\{y_1, y_2, \dots, y_m\}$为每次正面朝上的次数,每次抽样的不可观测数据$Z=\{z_1, z_2\}$为选择的硬币,估计参数为$\theta=\{\theta_1, \theta_2\}$。

对于第$n$次迭代,观察$\P(z_1|y_j, \theta^{(n)})$为:

易得$\P(z_2|y_j, \theta^{(n)}) = 1-\mu_{1j} = \mu_{2j}$

根据$\eqref{eq:4}$得,

分别对$\theta_1$和$\theta_2$求偏导取极限,得第$n+1$次$\theta$估计为:

混合高斯模型

混合高斯分布模型:

每次取样,依概率$\alpha_k$选择第$k$个高斯分布,根据相应概率分布生成观测数据$y_i$,由此组成观测数据$Y=\{y_1, y_2, \dots, y_m\}$。每次取样不可观测数据$Z=\{z_1, z_2, \dots, z_K\}$为抽取的高斯分布。估计参数为$\{\theta_1, \theta_2, \dots, \theta_K, \alpha_1, \alpha_2, \dots, \alpha_K\}$。

对于第$n$次迭代,观察$\P(z_1|y_j, \theta^{(n)})$为:

易得$\sum\limits_{k=1}^{K} \lambda_{kj} = 1$

根据$\eqref{eq:4}$得,

对各个参数求偏导取极限,第$n+1$次参数估计为:

混合高斯模型在生物信息学中有广泛的应用,例如参考资料4用于家族拷贝数变异(copy number variations,CNV)检测。

参考资料

  • 李航《统计学习方法》

  • Do CB, Batzoglou S: What is the expectation maximization algorithm? Nat Biotechnol. 2008;26(8):897-9.

  • The Expectation Maximization Algorithm A short tutorial

  • Chu JH, Rogers A, Ionita-Laza I, Darvishi K, Mills RE, Lee C, Raby BA: Copy number variation genotyping using family information. BMC Bioinformatics. 2013;14:157.

更新记录

2018年2月4日

Comments