概率模型（一）：极大似然估计与EM算法

极大似然估计

对于随机事件，在确定了概率模型，但是不确定模型参数的时候，可以用最大似然估计法来估计模型的参数。

比如，对于离散型随机事件，已知事件的概率分布函数为$P(x;\theta )$（分号左边是随机变量，分号右边是模型的参数，这个符号的意思是随机变量x的概率分布函数是以$\theta$为参数的函数），我们随机采样得到$n$个样本值$x_1,x_2,...,x_n$，这n个样本值满足独立同分布，则这些样本值的联合概率为：

$$P(x_1,x_2,...,x_n;\theta )=P(x_1;\theta )\ast P(x_2;\theta )\ast ...P(x_n;\theta )=\prod _{i=1}^{n}P(x_i;\theta )$$

但是对于该事件的概率分布函数$P(x;\theta )$的参数$\theta$的取值我们并不知道，如何估计一个合理的$\theta$值呢？

最大似然估计法是可以这样取$\theta$的值:

$$\theta =argmax\prod _{i=1}^{n}P(x_i;\theta )$$

背后的思想就是在众多的$\theta$取值中，取使得样本出现概率最大的$\theta$值。

我们定义似然函数为：

$$L(\theta )=\prod _{i=1}^{n}P(x_i;\theta )$$

有时为了计算方便，会加上对数函数，得到对数似然函数：

$$lnL(\theta )=ln\prod _{i=1}^{n}P(x_i;\theta )=\sum _{i=1}^{n}lnP(x_i;\theta )$$

对数似然函数取平均之后得到对数平均似然函数：

$$\hat{\ell }=\frac{1}{n}\sum _{i=1}^{n}lnP(x_i;\theta )$$

下面通过一个例子来说明最大似然估计方法的用途：例如有一大箱苹果，有部分苹果坏了，但是坏苹果的比例未知，我们想知道有多少苹果是坏的，但是不方便把所有的苹果都倒出来数一遍。现在我们从箱子中任意拿出一个苹果来，记录苹果是否坏了，然后放回去，把苹果摇匀后再拿出来一个，记录苹果是否坏了。如果反复，假设我们记录了100次，其中有10次拿到的苹果是坏的，那么箱子中坏的苹果的比例是多少呢？是10%吗？为什么呢？

我们假设箱子中坏苹果的比例是p，每次从箱子中取出一个苹果，记录是否为坏苹果，然后再返回去摇匀，说明每次取苹果是独立同分布事件，100次中取出10次的坏苹果的似然函数是：

$$P=p^{10}(1-p{)}^{90}$$

为了是似然函数取值最大，p该去什么值呢？

$$\frac{\partial P}{\partial p}=10p^9(1-p{)}^{90}-90p^{10}(1-p{)}^{89}=0$$

解出方程的解：

$$p=0.1$$

所以10%是对坏苹果比例的极大似然估计值。

类似的，对于连续型随机变量，已知概率密度函数为$f(x;\theta )$，对于随机采样得到的n个样本$x_1,x_2,...,x_n$，满足独立同分布的条件，但是对于连续性变量求取某个特定值的概率没有意义，我们需要求落在某个值附近区域的概率，随意n次采样落到$x_1,x_2,...,x_n$附近区域的联合概率为：

$$f(x_1;\theta )\Delta xf(x_2;\theta )\Delta x...f(x_n;\theta )\Delta x=\Delta x^n\prod _{i=1}^{n}f(x_i;\theta )$$

因为$\Delta x^n$是个固定的值，所以我们定义连续型随机事件的似然函数为：

$$L(\theta )=\prod _{i=1}^{n}f(x_i;\theta )$$

类似的，也可以定义它的对数平均似然函数：

$$\hat{\ell }=\frac{1}{n}\sum _{i=1}^{n}lnf(x_i;\theta )$$

下面我们来通过一个例子看一下连续性随机变量的极大似然估计：例如我们想要知道全国所有男性中，身高在1.7米到1.8之间的有多少人。我们没有办法统计出所有男性的身高，怎么办呢？

我们有理由相信全国男性的身高符合正太分布，也就是身高的概率密度函数为：

$$f(x)=\frac{1}{\sigma \sqrt{2\pi }}{e}^{-\frac{(x-\mu {)}^2}{2{\sigma }^2}}$$

假如我们采样得到了n份身高样本$[x_1,x_2,...,x_n]$,如果估计出参数$\mu ,{\sigma }^2$呢？

采用极大似然估计，似然函数为：

$$L(\mu ,{\sigma }^2)=\prod _{i=1}^n\frac{1}{\sigma \sqrt{2\pi }}{e}^{-\frac{(x_i-\mu {)}^2}{2{\sigma }^2}}=(2\pi \sigma {)}^{-\frac{n}{2}}{e}^{-\frac{1}{2{\sigma }^2}{\sum }_{i=1}^n(x_i-\mu {)}^2}$$

对数似然函数为：

$$lnL(\mu ,{\sigma }^2)=-\frac{n}{2}ln(2\pi )-\frac{n}{2}ln({\sigma }^2)-\frac{1}{2{\sigma }^2}\sum _{i=1}^n(x_i-\mu {)}^2$$

求对数似然函数极值，得到一下方程组：

$$\frac{\partial lnL(\mu ,{\sigma }^2)}{\partial \mu }=\frac{1}{{\sigma }^2}\sum _{i=1}^n(x_i-\mu )=0$$
$$\frac{\partial lnL(\mu ,{\sigma }^2)}{\partial \sigma }=-\frac{n}{2{\sigma }^2}+\frac{1}{2{\sigma }^4}\sum _{i=1}^n(x_i-\mu {)}^2=0$$

分别得出解为：

$$\mu =\frac{1}{n}\sum _{i=1}^n{x}_i$$

$${\sigma }^2=\frac{1}{n}\sum _{i=1}^n(x_i-\mu {)}^2$$

估计出正太分布的参数后，在知道人数总量的情况下，就可以计算身高落在1.7到1.8米之间的人数了,假设男性人口总量为M：

$$N=M\ast P(1.7<h<1.8)=M\ast {\int }_{1.7}^{1.8}f(x)dx$$

EM算法

我们来看一个比较经典的例子：有A,B两块硬币，抛投落地后正面朝上的概率分别记为$P_A,P_B$（当然，背面朝上的概率分别为$1-P_A,1-P_B$），为了估计两个概率值，我们连续做5次实验，每次实验随机取一个硬币（等概率），随机抛投5次，记录每次抛投的结果：

实验编号	硬币	实验结果
1	A	正正反正反
2	B	反反正正反
3	A	正反反反反
4	B	正反反正正
5	A	反正正反反

采用极大似然估计方法，很容易计算出$P_A，P_B$的估计值，A硬币的的似然函数为：

$$L(P_A)=P_A^6(1-P_A{)}^9$$

求极值得以下方程：

$$\frac{\partial L}{\partial P_A}=6P_A^5(1-P_A{)}^9-9P_A^6(1-P_A{)}^8=0$$

得出解：$P_A=\frac{6}{15}=0.4$

同理，得出：$P_B=\frac{5}{10}=0.5$

下面加大一下问题的难度，假如每次实验我们没办法知道所选的硬币是A还是B，还有办法估计A和B的正面朝上的概率$P_A，P_B$吗？

实验编号	硬币	实验结果
1	？	正正反正反
2	？	反反正正反
3	？	正反反反反
4	？	正反反正正
5	？	反正正反反

EM算法就排上用场了，可以按照下列步骤估计$P_A,P_B$：

• 设置$P_A$和$P_B$的估计初始值，比如$P_A=0.2,P_B=0.7$
• 设置每次实验选择A，B硬币是等概率的，记为$P(z=A)=P(z=B)=0.5$，设随机变量$x_i$表示每次实验中，硬币的朝向，目前为止我们可以得到下列方程：

$$P(x_i|z=A)=P_A=0.2$$
$$P(x_i|z=B)=P_B=0.7$$
$$P(x_1,x_2,...,x_5|z=A)=\prod _{i=1}^{5}P(x_i|z=A)$$
$$P(x_1,x_2,...,x_5|z=B)=\prod _{i=1}^{5}P(x_i|z=B)$$
$$P(x_1,x_2,...,x_5)=P(x_1,x_2,...,x_5|z=A)P(z=A)+P(x_1,x_2,...,x_5|z=B)P(z=B)$$

• 进入E步，我们计算一下再当前初始值情况下第一次实验选择了A，B硬币的后验概率,计算方式如下：

由贝叶斯公式：

$$P(A|B)=\frac{P(B|A)P(A)}{P(B)}$$

得到：

$$P(z=A|x_1,x_2,...,x_5)=\frac{P(x_1,x_2,...,x_5|z=A)P(z=A)}{P(x_1,x_2,...,x_5)}$$

$$=\frac{P(x_1,x_2,...,x_5|z=A)P(z=A)}{P(x_1,x_2,...,x_5|z=A)P(z=A)+P(x_1,x_2,...,x_5|z=B)P(z=B)}$$

$$P(z=B|x_1,x_2,...,x_5)=\frac{P(x_1,x_2,...,x_5|z=A)P(z=B)}{P(x_1,x_2,...,x_5)}$$

$$=\frac{P(x_1,x_2,...,x_5|z=A)P(z=A)}{P(x_1,x_2,...,x_5|z=A)P(z=A)+P(x_1,x_2,...,x_5|z=B)P(z=B)}$$

对所有的实验，采用相同的计算法方式，计算结果如下表：

实验编号	硬币	实验结果	硬币为A概率	硬币为B概率
1	？	正正反正反	0.14	0.86
2	？	反反正正反	0.61	0.39
3	？	正反反反反	0.94	0.06
4	？	正反反正正	0.14	0.86
5	？	反正正反反	0.61	0.39

• 进入M步，估计$P_A,P_B$的值，最大化下面函数$Q(\theta ,{\theta }^{(i)})$：

$$Q(\theta ,{\theta }^{(i)})=\sum _{i=1}^{5}P(z^{(i)}|x_1^{(i)},...,x_5^{(i)},{\theta }^{(i)})logP(x_1^{(i)},...,x_5^{(i)},z^{(i)}|\theta )$$

$$=\sum _{i=1}^{5}P(z^{(i)}|x_1^{(i)},...,x_5^{(i)},{\theta }^{(i)})logP(x_1^{(i)},...,x_5^{(i)}|z^{(i)},\theta )P(z^{(i)}|\theta )$$

$$=0.14log{P}_A^3(1-P_A{)}^2{P}_z+0.86log{P}_B^3(1-P_B{)}^2(1-P_z)$$
$$+0.61log{P}_A^2(1-P_A{)}^3{P}_z+0.39log{P}_B^2(1-P_B{)}^3(1-P_z)$$
$$+0.94log{P}_A(1-P_A{)}^4{P}_z+0.06log{P}_B(1-P_B{)}^4(1-P_z)$$
$$+0.14log{P}_A^3(1-P_A{)}^2{P}_z+0.86log{P}_B^3(1-P_B{)}^2(1-P_z)$$
$$+0.61log{P}_A^2(1-P_A{)}^3{P}_z+0.39log{P}_B^2(1-P_B{)}^3(1-P_z)$$

由取极值条件得方程组：

$$\frac{\partial Q}{\partial P_z}=0$$
$$\frac{\partial Q}{\partial P_A}=0$$
$$\frac{\partial Q}{\partial P_B}=0$$

求得$P_z=0.49,P_A=0.35,P_B=0.53$

• 重复以上的E和M步，直到$P_z,P_A,P_B$收敛或$Q$函数收敛为止。

EM算法推导

EM算法希望通过迭代计算，最大化可观测变量的似然函数：

$$L(\theta )=logP(Y|\theta )=log\sum _{Z}P(Y,Z|\theta )$$
$$=log\sum _{Z}P(Y|Z,\theta )P(Z|\theta )$$

但是直接最大化此函数比较困难，函数中存在隐含变量Z.

假设我们初始的估计值为${\theta }^{(0)}$,第$i$次的估计值为${\theta }^{(i)}$

$$L({\theta }^{(i)})=logP(Y|{\theta }^{(i)})$$

$$L(\theta )-L({\theta }^{(i)})=log\sum _{Z}P(Y,Z|\theta )-logP(Y|{\theta }^{(i)})$$

$$=log\sum _{Z}P(Y|Z,\theta )P(Z|\theta )-logP(Y|{\theta }^{(i)})$$

$$=log\sum _{Z}P(Z|Y,{\theta }^{(i)})\frac{P(Y|Z,\theta )P(Z|\theta )}{P(Z|Y,{\theta }^{(i)})}-logP(Y|{\theta }^{(i)})$$

$$\ge \sum _{Z}P(Z|Y,{\theta }^{(i)})log\frac{P(Y|Z,\theta )P(Z|\theta )}{P(Z|Y,{\theta }^{(i)})}-logP(Y|{\theta }^{(i)})$$

$$=\sum _{Z}P(Z|Y,{\theta }^{(i)})log\frac{P(Y|Z,\theta )P(Z|\theta )}{P(Z|Y,{\theta }^{(i)})}-\sum _{Z}P(Z|Y,{\theta }^{(i)})logP(Y|{\theta }^{(i)})$$

$$=\sum _{Z}P(Z|Y,{\theta }^{(i)})log\frac{P(Y|Z,\theta )P(Z|\theta )}{P(Z|Y,{\theta }^{(i)})P(Y|{\theta }^{(i)})}$$

综上得出：

$$L(\theta )\ge L({\theta }^{(i)})+\sum _{Z}P(Z|Y,{\theta }^{(i)})log\frac{P(Y|Z,\theta )P(Z|\theta )}{P(Z|Y,{\theta }^{(i)})P(Y|{\theta }^{(i)})}$$

上式中右侧部分记为$B(\theta ,{\theta }^{(i)})$，则有：

$$L(\theta )\ge B(\theta ,{\theta }^{(i)})$$

$${\theta }^{(i+1)}=argma{x}_{\theta }B(\theta ,{\theta }^{(i)})$$
$$=argma{x}_{\theta }(logP(Y|{\theta }^{(i)})+\sum _{Z}P(Z|Y,{\theta }^{(i)})log\frac{P(Y|Z,\theta )P(Z|\theta )}{P(Z|Y,{\theta }^{(i)})P(Y|{\theta }^{(i)})})$$

$$=argma{x}_{\theta }\sum _{Z}P(Z|Y,{\theta }^{(i)})log\frac{P(Y|Z,\theta )P(Z|\theta )}{P(Z|Y,{\theta }^{(i)})}$$

$$=argma{x}_{\theta }\sum _{Z}P(Z|Y,{\theta }^{(i)})logP(Y|Z,\theta )P(Z|\theta )-\sum _{Z}P(Z|Y,{\theta }^{(i)})logP(Z|Y,{\theta }^{(i)})$$

$$=argma{x}_{\theta }\sum _{Z}P(Z|Y,{\theta }^{(i)})logP(Y,Z|\theta )-\sum _{Z}P(Z|Y,{\theta }^{(i)})logP(Z|Y,{\theta }^{(i)})$$

省略常数项得：

$${\theta }^{(i+1)}=argma{x}_{\theta }\sum _{Z}P(Z|Y,{\theta }^{(i)})logP(Y,Z|\theta )$$

这里就是上面例子中提到的$Q$函数的定义

总结来说，EM算法的过程为：

1. 模型包含可观察变量$Y$和隐藏变量$Z$,模型参数$\theta$初始值为${\theta }^{(0)}$，第$i$次迭代的值记为${\theta }^{(i)}$
2. E步:计算函数$Q(\theta ,{\theta }^{(i)})={\sum }_{Z}logP(Y,Z|\theta )P(Z|{\theta }^{(i)})$
3. M步:求参数${\theta }^{(i+1)}=argma{x}_{\theta }Q(\theta ,{\theta }^{(i)})$
4. 重复2，3两步，直到Q函数或参数收敛。

另外还需要注意的是，EM算法对初始化的参数$\theta$是敏感的，不同的初始化值可能会造成最后的收敛值的不同