EM算法学习笔记
引入
设袋子里有红白两种球,比例为\(p:1-p\)。现在不知道比例值是多少,因此可以随机抽样:\(x_1,x_2...x_n\),每一个\(x_i\)都是随机变量的随机抽样结果,并进行如下估计:\(p=\frac{|\{x_i=红\}|}{n}\)
如果有两个袋子,各自都有红色的球和白色的球,抽取的时候先随机挑选一个袋子,设选第一个袋子的概率为\(w\),第一个袋子中\(红球:白球=p:1-p\),第二个为\(红球:白球=q:1-q\)。
如果此时进行抽取,抽取到红球的概率为:\(P\{x=红\}=w\times p+(1-w)\times q\),抽到白球的概率同理。
如果想要知道从第一个袋子中抽到红球的概率:\(P\{第一个|x_i=红\}\),可以把概率乘以权重放到矩阵中:
$$\begin{pmatrix}
wp&w(1-p)\
(1-w)q&(1-w)(1-q)\
\end{pmatrix}$$
取矩阵的第一列分析,不难得出\(P\{第一个|x_i=红\}=\frac{wp}{wp+(1-w)q}\)。
此时模型的\(w,p,q\)三个参数均未知,可以使用极大似然法对其进行估计。令1代表红色,0代表白色,每个\(x_i\)的密度函数为\(wp^{x_i}(1-p)^{1-x_i}+(1-w)q^{x_i}(1-q)^{1-x_i}\),如果\(x_i=1\)(红色),带入上式可以发现结果就是\(wp+(1-w)q\),也就是矩阵第一列的和。
根据极大似然估计,需要把所有样本的密度函数做乘积,似然函数就是:
\(\Pi_{i=1}^n [wp^{x_i}(1-p)^{1-x_i}+(1-w)q^{x_i}(1-q)^{1-x_i}]\),需要对其求极大值。如果按照极大似然的思路,会发现这个问题很复杂(中括号里的部分不易处理,不信你试试)。不过这个函数可以转换为如下的形式:
\(\Pi_{i=1}^n [(wp+(1-w)q)^{x_i}(w(1-p)+(1-w)(1-q))^{1-x_i}]\),这直观上也很容易理解,因为\(x_i\)只有0和1两种形式(带入\(x_i=1 \ or\ 0\)试试)。对其处理的话,可以令\(P=wp+(1-w)q\),\(Q=w(1-p)+(1-w)(1-q)\),函数继续变形为\(\Pi_{i=1}^nP^{x_i}Q^{1-x_i}=\Pi_{i=1}^nP^{x_i}(1-P)^{1-x_i}\),对其取对数后得\((\Sigma_{i=1}^nx_i)lnP+(n-\Sigma_{i=1}^nx_i)ln(1-P)\),对P求导并令导数为0,可得\(P=\frac{\Sigma_{i=1}^nx_i}{n}=\frac{k}{n}\),\(k\)为样本中1(红球)的个数,且有\(w(1-p)+(1-w)(1-q)=\frac{n-k}{n}\)(本质上等于上面的方程)。
总结一下,刚才设计的模型(先随机选袋子,再随机选球)从极大似然估计的角度,\(wp+(1-w)q=\frac{k}{n}\)就是最佳估计,但是一个方程里含有三个未知数,所以只要三个参数满足这个方程就是一个合理的估计,即给出两个参数就可以把另一个参数解出来。
从另一个角度看问题
现在给出三个参数\(p^{(0)}, q^{(0)}, w^{(0)}\),对于\(x_i\),要根据它的颜色判断它来自于哪个袋子(利用上面的矩阵来计算):
红球来自第一个袋子的概率:\(\frac{wp}{wp+(1-w)q}=a_i\),白球来自第一个袋子的概率:\(\frac{w(1-p)}{w(1-p)+(1-w)(1-q)}=b_i\)...同理设红球和白球来自第二个袋子的概率为\(c_i,d_i\),那么对于每个球,这四个概率都可以求出来,列成一张表:
\(x_1\) | \(x_2\) | .. | ||
---|---|---|---|---|
第一个-红 | \(\frac{wp}{wp+(1-w)q}=a_1\) | |||
第一个-白 | \(\frac{w(1-p)}{wp+(1-w)(1-q)}=b_1\) | |||
第二个-红 | \(\frac{(1-w)p}{wp+(1-w)q}=c_1\) | |||
第二个-白 | \(\frac{(1-w)(1-p)}{wp+(1-w)(1-q)}=d_1\) |
每一列的总和都是一样的(都是1),因为这里考虑的是后验概率,即每个球已知是红色还是白色,所以对于每一列只是把和已知颜色对应的行加起来(而非把所有的行加起来)。
因此来自于第一个袋子的概率\(w\)可以给出一个重新的估计:\(w^{(1)}=\frac{\Sigma_{x_i=1}a_i+\Sigma_{x_i=1}b_i}{n}\),即表中前两行相加再除以n。这个过程就是:选定\(w^{(0)}\)为初值,根据初值计算出一个个概率值,再根据上表的概率值计算出新的\(w^{(1)}\)。
同理\(p^{(1)}=\frac{\Sigma_{x_i=1}a_i}{\Sigma_{x_i=1}a_i+\Sigma_{x_i=0}b_i}\),\(q^{(1)}=\frac{\Sigma_{x_i=1}c_i}{\Sigma_{x_i=1}c_i+\Sigma_{x_i=0}d_i}\)。如此可以继续迭代下去得到新的参数值,能够证明这个算法是收敛的,至少可以收敛到局部极大值(不一定是全局)。
整个过程就叫做EM算法。
推广
对于高斯分布\(N(\mu,\sigma^2)\),设有\(K\)个高斯分布,从中随机抽取到的概率分别为\(w_1,w_2...w_k\),且满足\(\Sigma_{i=1}^k=1\)。
\(X\)的密度函数\(w_1f_1(x)+....+w_kf_k(x)\),\(f(x)=\frac{1}{\sqrt{2\pi}\sigma}e^{-\frac{(x-\mu)^2}{2\sigma^2}}\)
现在给出\(x_1,...x_n\),如何估计出所有的\(\mu_j,\sigma_j^2,w_i\)?
按照极大似然的过程,\(max\ \Pi_{i=1}^n[\Sigma_{j=1}^Kw_jf_j(x_i)]\)等价于\(max\ \Sigma_{i=1}^nlog(\Sigma_{j=1}^Kw_jf_j(x_i))\)。但是这样并没有对问题进行化简,所以根据前面的结果,给出一组参数\(\mu_j^{(0)},\sigma_j^{2(0)},w_i^{(0)}\),考虑建一个表:
高斯分布 | \(x_1\) | \(x_2\) | .. | |
---|---|---|---|---|
1 | \(\frac{w_1f_1(x_1)}{w_1f_1(x_1)+...+w_kf_k(x_1)}=w_{11}\) | \(\frac{w_1f_1(x_2)}{w_1f_1(x_2)+...}=w_{21}\) | .. | |
2 | \(\frac{w_2f_2(x_1)}{w_1f_1(x_1)+...}=w_{12}\) | .. | .. | |
.. | .. | .. | .. | |
k | .. | .. | .. |
这样就又可以进行估计了。同样,每一列的总和是1,则对于\(w_1\)有一个很自然的估计:\(w_1^{(1)}=\frac{\Sigma_{i=1}^nw_{i1}}{n}\),其余的w值同理。现在看均值和方差:
\(\mu_1^{(1)}=\frac{\Sigma_{i=1}^nw_{i1}x_i}{\Sigma_{i=1}^nw_{i1}}\),注意这里分布是对每一行的\(w_{i1}\)求和,并不是1。\(\sigma_1^{(1)}=\frac{\Sigma_{i=1}^nw_{i1}(x_i-\mu_1^{(1)})^2}{\Sigma_{i=1}^nw_{i1}}\)
迭代若干次以后就会得到比较合适的高斯分布了。
对于二维高斯分布,其参数也都可以用EM算法进行估计。
需要注意的是,这里的\(K\)是超参数,需要人为给定。