标签归档：统计

统计算法的理论与实例之一 最大似然估计与EM算法

发表评论

版权声明：原创作品，欢迎转载，但转载请以超链接形式注明文章来源(planckscale.info)、作者信息和本声明，否则将追究法律责任。

去年和朋友做了一些统计方面的东西，一些笔记发在这儿。统计我是外行，贻笑大方之处请指正。

\section{记号约定}
本文中用大写字母X,Y,Z...表示随机变量，小写字母x,y,z...表示随机变量的观测值，粗体表示向量/矩阵。\mathcal{X},\mathcal{Y},\mathcal{Z}...表示向量值随机变量或随机变量X,Y,Z...的观测样本。

\section{估计理论}
给随机变量一个待定参数的概率模型，通过一系列样本观测，萃取出概率模型最可能参数的理论。

\subsection{一些数学}
考虑随机变量X,Y的概率模型p(X,Y),我们经常需要考虑一系列观测样本\mathcal{X}=\{x_1,...x_N\},\mathcal{Y}=\{y_1,...y_N\}的分布

(1)   \begin{eqnarray*} p(\mathcal{X},\mathcal{Y})=\prod_i p(x_i,y_i) \end{eqnarray*}

观测样本之间的独立性意味着边缘分布p(\mathcal{X})（或p(\mathcal{Y})），条件分布p(\mathcal{X}\vert\mathcal{Y})（或p(\mathcal{Y}\vert\mathcal{X})）也满足类似的因子分解性质。严格说来，假设Y取值于\Delta,则\mathcal{Y}取值于\otimes\Delta^N,我们有

(2)   \begin{eqnarray*} p(\mathcal{X})&=&\sum_{\mathcal{Y}}p(\mathcal{X},\mathcal{Y})\\ &=&\sum_{\mathcal{Y}\in \otimes\Delta^N}\prod_i p(x_i,y_i)\\ &=&\prod_i\sum_{y_i\in\Delta} p(x_i,y_i)\\ &=&\prod_i p(x_i)\\ p(\mathcal{X}\vert\mathcal{Y})&=&\frac{p(\mathcal{X},\mathcal{Y})}{p(\mathcal{Y})}\\ &=&\prod_i \frac{p(x_i,y_i)}{p(y_i)}\\ &=&\prod_i p(x_i\vert y_i) \end{eqnarray*}

\subsection{最大似然估计}
假设随机变量X的概率模型为p(X\vert\theta),其中\theta为模型参数。我们有一系列独立观测样本\mathcal{X}=\{x_1,...x_N\},欲定出最可能的模型参数\theta_{\star}.

最大似然估计即假设观测结果是最大可能的结果，从而定出分布参数的估计方法。定义对数似然函数

(3)   \begin{eqnarray*} \log L(\theta;\mathcal{X})&=&\log p(\mathcal{X}\vert\theta)\\ &=&\log\prod_i p(x_i\vert\theta)\\ &=&\sum_{i}\log p(x_i\vert\theta) \end{eqnarray*}

则模型的最大似然估计为

(4)   \begin{eqnarray*} \theta_{\star}=\max_{\theta} \log L(\theta) \end{eqnarray*}

注意-\log p(x_i,\theta)是第i次观测所获取到的信息量，此目标函数一个含义是\textbf{调整参数使得观测到的信息量取极小}。

\subsection{概率分布的描述能力}
估计理论使得我们能够以估计分布的形式得到样本空间中的几何对象。概率分布的描述能力是广泛的，这里举几个例子。

\textbf{简单椭球区域}。可以简单的由Gauss分布p(\boldsymbol{X})=\mathcal{N}(\boldsymbol{\mu},\boldsymbol{\Sigma})刻画，\boldsymbol{\mu}给出了位置而\boldsymbol{\Sigma}则给出了椭球在各方向上的轴长。

\textbf{多个椭球区域}。由多个Gauss分布的叠加，即\textbf{混合Gauss分布}描述。概率密度分布\mathcal{M}(\{ \boldsymbol{\mu}_a,\boldsymbol{\Sigma}_a,\boldsymbol{\pi}_a\})为

(5)   \begin{eqnarray*} p(\boldsymbol{X})=\mathcal{M}(\{ \boldsymbol{\mu}_a,\boldsymbol{\Sigma}_a,\boldsymbol{\pi}_a\}) \sim \sum_a \pi_a \mathcal{N}(\boldsymbol{\mu}_a,\boldsymbol{\Sigma}_a) \end{eqnarray*}

分布描述了样本空间里多个简单椭球区域，每个区域可以代表一个数据类。通常在应用中，我们把观测不到的隐含变量，类别标签Z加入到统计量集合中，于是概率密度函数成为

(6)   \begin{eqnarray*} p(\boldsymbol{X},Z)=\sum_a\delta_{Za}\pi_a\mathcal{N}(\boldsymbol{\mu}_a,\boldsymbol{\Sigma}_a) \end{eqnarray*}

\textbf{函数图形}。函数f(x)在样本空间里的图形可以由概率分布

(7)   \begin{eqnarray*} p(X,Y)=\mathcal{N}(f(X),\sigma^2) \end{eqnarray*}

刻画。给定任意X,Y都满足一个以f(X)为中心的Gauss分布，因而这里刻画了一个对Y=f(X)这一背后本质的数量关系所做的测量，而该测量是带有Gauss噪声的。

\subsection{最小二乘法}
考虑函数族f_{\boldsymbol{\alpha}}(x),我们通过一系列测量数据\mathcal{X}\equiv\{x_{1},...x_N\},\mathcal{Y}\equiv\{y_{1},...y_N\}去估计背后真实的函数关系f_*(x).
假设观测噪声为Gauss分布，于是观测的概率分布为

(8)   \begin{eqnarray*} p(X,Y)=\mathcal{N}(f_*(X),\sigma^2) \end{eqnarray*}

构造对数似然函数

(9)   \begin{eqnarray*} \log L(\boldsymbol{\theta};\mathcal{X},\mathcal{Y})&=&\log p(\mathcal{X},\mathcal{Y}\vert\boldsymbol{\theta})\\ &=&\log\prod_i p(x_i,y_i\vert\boldsymbol{\theta})\\ &=&\sum_{i}\log \mathcal{N}(f_{\boldsymbol{\alpha}}(x_i),\sigma)\\ &=&\sum_{i} (-\log\sigma-\frac{1}{2}\log\pi -\frac{(y_i-f(x_i,\boldsymbol{\alpha}))^2}{2\sigma^2}) \end{eqnarray*}

于是，\max_{\theta} \log L(\theta;\mathcal{X},\mathcal{Y})给出

(10)   \begin{eqnarray*} &\max_{\boldsymbol{\alpha},\sigma}&\{-N\log\sigma-\frac{\sum_i(y_i-f_{\boldsymbol{\alpha}}(x_i))^2}{2\sigma^2}\}\\ \Longleftrightarrow&\min_{\boldsymbol{\alpha},\sigma}&\{N\log\sigma+\frac{\sum_i(y_i-f_{\boldsymbol{\alpha}}(x_i))^2}{2\sigma^2}\}\\ \Longleftrightarrow&\min_{\boldsymbol{\alpha}}&\{\frac{N}{2}\log\frac{\sum_i(y_i-f_{\boldsymbol{\alpha}}(x_i))^2}{N}+\frac{N}{2}\}\\ \Longleftrightarrow&\min_{\boldsymbol{\alpha}}&\{\sum_i(y_i-f_{\boldsymbol{\alpha}}(x_i))^2\} \end{eqnarray*}

其中\sigma=\sqrt{\sum_i(y_i-f_{*}(x_i))^2/N}.

于是我们得到了最小二乘法。

\subsection{EM算法}
对于一个联合分布，实际应用中并非每一个随机变量都可以被观测到。比如上述混合Gauss分布中，通常类别标签Z是无法观测到的。这种不完备观测情况下直接使用最大似然估计会出现计算量过大，实际不可解的问题。于是对于这类估计问题，我们使用\textbf{EM算法}来求解。

考虑混合Gauss分布

(11)   \begin{eqnarray*} p(X,Z\vert\boldsymbol{\theta})=\sum_a\delta_{Za}\pi_a\mathcal{N}(\mu_a,\sigma^2_a) \end{eqnarray*}

设\mathcal{X}\equiv\{x_{1},...x_N\}为观测样本,\mathcal{Z}\equiv\{z_1,...z_N\}为观测样本对应的隐变量。对数似然函数为

(12)   \begin{eqnarray*} \log L(\boldsymbol{\theta};\mathcal{X})&=&\log\sum_{\mathcal{Z}} p(\mathcal{X},\mathcal{Z}\vert\boldsymbol{\theta})\\ &=&\log\sum_{\mathcal{Z}}\prod_i p(x_i,\mathcal{Z}\vert\boldsymbol{\theta}) \end{eqnarray*}

其中这一对数似然函数在计算上有诸多不便，首先加和号阻止了对数运算作用在p上，再者在z变量取值返回较广时计算代价也会很大。因此传统的最大似然估计在此处不适用。

在EM算法中，我们在每一个迭代步都构造一个新的函数Q^{(t)},用\max Q^{(t)}来近似\log L(\boldsymbol{\theta};\mathcal{X})在当前解\boldsymbol{\theta}^{(t)}附近邻域上的变化规律,迭代收敛到最终解

(13)   \begin{eqnarray*} Q^{(t)}(\boldsymbol{\theta})\equiv Q(\boldsymbol{\theta}\vert\boldsymbol{\theta}^{(t)})&\equiv& E_{\mathcal{Z}\vert \mathcal{X},\boldsymbol{\theta}^{(t)}}[\log L(\boldsymbol{\theta};\mathcal{X},\mathcal{Z})]\\ \boldsymbol{\theta}^{(t+1)}&=&\mathrm{argmax}_{\boldsymbol{\theta}} Q(\boldsymbol{\theta}\vert\boldsymbol{\theta}^{(t)}) \end{eqnarray*}

其中

(14)   \begin{eqnarray*} p(\mathcal{Z}\vert \mathcal{X},\boldsymbol{\theta}^{(t)})= \frac{p(\mathcal{X},\mathcal{Z}\vert \boldsymbol{\theta}^{(t)})}{\sum_{\mathcal{Z}} p(\mathcal{X},\mathcal{Z}\vert \boldsymbol{\theta}^{(t)})} \end{eqnarray*}

该算法的正确性证明留到下一小节处理。

\subsection{EM算法的理论基础}
在目标函数较为复杂的情况下，极大/极小值的寻找通常采用迭代算法来处理。计算目标函数精确的形状是困难的，因而我们需要\textbf{在每一迭代步通过有限的计算，探知目标函数在当前迭代解之邻域上的变化规律，用一个简单的函数来刻画这一局域变化趋势，然后求解这一局部的极小/极大化问题，最后进入下一轮迭代}。这样，一个复杂目标函数的极小/极大问题就被分解成了一系列容易求解的简单问题，通过迭代不断靠近极优解。

具体说来，不同的算法会为目标函数的局域变化建立不同的模型，比如通过计算导数，用平面来建模该处的目标函数的变化；再如计算到二阶导数，用抛物面来近似目标函数，等等。

EM算法也是这样一种算法。\textbf{它的特殊之处在于，由于目标函数，即概然函数的特殊结构，我们可以找到比一次/二次近似更好的局域变化模型}，从而获得比常规算法更好的性能。这一局域变化模型就是Q^{(t)}(\boldsymbol{\theta}).

首先我们来寻找Q^{(t)}(\boldsymbol{\theta})和\log p(\mathcal{X}\vert \boldsymbol{\theta})之间的关系

(15)   \begin{eqnarray*} \log p(\mathcal{X}\vert \boldsymbol{\theta})&=&\log p(\mathcal{X},\mathcal{Z}\vert\boldsymbol{\theta}) - \log p(\mathcal{Z}\vert \mathcal{X},\boldsymbol{\theta})\\ \log p(\mathcal{X}\vert \boldsymbol{\theta})&=&\sum_{\mathcal{Z}} p(\mathcal{Z}\vert \mathcal{X},\boldsymbol{\theta}^{(t)})\log p(\mathcal{X},\mathcal{Z}\vert\boldsymbol{\theta}) - \sum_{\mathcal{Z}} p(\mathcal{Z}\vert \mathcal{X},\boldsymbol{\theta}^{(t)}) \log p(\mathcal{Z}\vert \mathcal{X},\boldsymbol{\theta})\\ \log p(\mathcal{X}\vert \boldsymbol{\theta})&=&Q(\boldsymbol{\theta} \vert \boldsymbol{\theta}^{(t)})+ H(\boldsymbol{\theta} \vert \boldsymbol{\theta}^{(t)}) \end{eqnarray*}

易证H(\boldsymbol{\theta} \vert \boldsymbol{\theta}^{(t)})\geq 0,可见

(16)   \begin{eqnarray*} \log p(\mathcal{X}\vert \boldsymbol{\theta})&\geq&Q(\boldsymbol{\theta} \vert \boldsymbol{\theta}^{(t)}) \end{eqnarray*}

即\textbf{Q(\boldsymbol{\theta} \vert \boldsymbol{\theta}^{(t)})构成对数似然函数的下界}。将15减去其在\boldsymbol{\theta}^{(t)}处的值，又有

(17)   \begin{eqnarray*} \Delta \log p(\mathcal{X}\vert \boldsymbol{\theta}) - \Delta Q(\boldsymbol{\theta} \vert \boldsymbol{\theta}^{(t)}) &=& H(\boldsymbol{\theta}\vert \boldsymbol{\theta}^{(t)}) - H(\boldsymbol{\theta}^{(t)}\vert \boldsymbol{\theta}^{(t)})\\ &=&D(p(\mathcal{Z}\vert \mathcal{X},\boldsymbol{\theta}^{(t)}) \mid p(\mathcal{Z}\vert \mathcal{X},\boldsymbol{\theta}))\\ &\geq& 0 \end{eqnarray*}

即\textbf{每一迭代步中必有\Delta \log L \geq \Delta Q}.

可见Q^{(t)}(\boldsymbol{\theta})可以反映\log p(\mathcal{X}\vert \boldsymbol{\theta})在\boldsymbol{\theta}^{(t)}处的变化规律。更重要的是，Q^{(t)}(\boldsymbol{\theta})构成了对数概然函数的一个下界，这样以来我们就可以放心的求解\max Q^{(t)}(\boldsymbol{\theta})这一问题来得到一个更好的迭代解，而不用像LM算法中那样担心二阶模型只在一个局域的范围内构成目标函数的有效近似，需要用置信域方法避免大步长的优化。

算法的收敛性由这一点来保证。由于D(p(\mathcal{Z}\vert \mathcal{X},\boldsymbol{\theta}^{(t)}) \mid p(\mathcal{Z}\vert \mathcal{X},\boldsymbol{\theta}))在\boldsymbol{\theta}=\boldsymbol{\theta}^{(t)}附近的一阶展开项为零，于是17的一阶展开给出

(18)   \begin{eqnarray*} \frac{\mathrm{d}}{\mathrm{d} \boldsymbol{\theta}} \log p(\mathcal{X}\vert \boldsymbol{\theta})\vert_{\boldsymbol{\theta}=\boldsymbol{\theta}^{(t)}} &=& \frac{\mathrm{d}}{\mathrm{d}\boldsymbol{\theta}} Q(\boldsymbol{\theta} \vert \boldsymbol{\theta}^{(t)})\vert_{\boldsymbol{\theta}=\boldsymbol{\theta}^{(t)}} \end{eqnarray*}

算法迭代到Q不再有明显变化(\boldsymbol{\theta}^{(t+1)}\approx\boldsymbol{\theta}^{(t)})时，上式右端为零，式18意味着\boldsymbol{\theta}=\boldsymbol{\theta}^{(t)}同样也是对数似然函数的极值点。于是Q的收敛意味着对数似然函数的收敛，算法必然收敛于一个极优解。