EM算法解决含隐变量的掷硬币问题

数据挖掘十大经典算法–EM算法（详解+代码）

适用场景：用于含有隐变量的概率模型参数的极大似然估计，或者极大后验概率估计

先对隐变量和极大似然估计做出解释

考虑这样两个问题：

1、抛一枚硬币，假设正面向上的概率是$\theta$，先对其进行了10次试验，结果为8次正面，2次反面。

2、有三枚硬币A，B，C，其正面向上的概率分别为$\pi$，$p$，$q$，现进行10次实验，先掷硬币A，若A为正面，则掷出B，否者掷出C，记录最终的掷硬币结果，现已知实验的观测结果为6次正面，4次反面。

背景知识

隐变量

可以看出，问题2比问题1多出一个中间过程，以投掷A硬币的结果取确定投掷B，C，并对最终结果产生影响。但是在实验过程中只知道最终结果，不知道这一中间结果。隐变量就是这一会对最终结果产生影响的中间结果。

极大似然估计

利用已知样本的结果信息，反推最具有可能导致这些样本结果出现的模型的参数值。步骤为：

1. 写出似然函数
2. 对似然函数取对数并整理
3. 求导数，并令其等于0
4. 解似然方程

极大似然估计的使用在总体为离散型和连续性时有一定的差别

1、离散型

若总体$x$为离散型，其概率分布列为
$$P\left\{\begin{array}{c}X=x\end{array}\right\}=p(x_i,\theta )$$
其中$\theta$为未知参数，设$(X_1,X_2\cdots ,X_n)$为取自总体的样本容量为$n$的一个样本，则$(X_1,X_2\cdots ,X_n)$的联合分布律为${\prod }_{i=1}^{n}p(x_i,\theta )$。

联合分布律：多维随机变量的分布律

例：$F(x,y)=P\left\{\begin{array}{c}(X\le x)\bigcap (Y\le y)\end{array}\right\}$ 分布律的表达形式为 $p(X\le x,Y\le y)$

设$(X_1,X_2\cdots ,X_n)$的一组观测值为$(x_1,x_2\cdots ,x_n)$，则样本$(X_1,X_2\cdots ,X_n)$取观测值$(x_1,x_2\cdots ,x_n)$的概率为：
$$L(\theta )=L(x_1,x_2\cdots x_n;\theta )=\prod _{i=1}^{n}p(x_i,\theta )$$
这一概率随$\theta$的取值未改变，它是$\theta$的函数，称$L(\theta )$为样本的似然函数

2、连续性

若总体X为连续型，其概率密度函数为$f(x;\theta )$，其中$\theta$为未知参数。，设$(X_1,X_2\cdots ,X_n)$为取自总体的样本容量为$n$的一个样本，则$(X_1,X_2\cdots ,X_n)$的联合联合概率密度为${\prod }_{i=1}^{n}f(x_i,\theta )$。

联合概率密度：对于二元随机变量$(X,Y)$的分布函数$F(x,y)$，如果存在非负函数$f(x,y)$，使得对于任意$x$，$y$，都有
$$F(x,y)={\int }_{-\infty }^x{\int }_{-\infty }^{y}f(u,v)dudv$$
称$(X,Y)$为二元连续型随机变量，$f(x,y)$为二元随机变量$(X,Y)$的联合概率密度函数

设$(X_1,X_2\cdots ,X_n)$的一组观测值为$(x_1,x_2\cdots ,x_n)$，则随机点$(X_1,X_2\cdots ,X_n)$落在点$(x_1,x_2\cdots ,x_n)$的邻边内的概率近似地为${\prod }_{i=1}^{n}f(x_i,\theta )d{x}_i$
$$L(\theta )=L(x_1,x_2\cdots x_n;\theta )=\prod _{i=1}^{n}f(x_i,\theta )$$
$L(\theta )$为样本的似然函数

极大似然估计就是，固定样本观测值$(x_1,x_2\cdots ,x_n)$，挑选$\theta$使得：
$$L(x_1,x_2\cdots x_n;\theta )=\max L(x_1,x_2\cdots x_n;\theta )$$
这样的$\hat{\theta }$与样本值有关，$\hat{\theta }(x_1,x_2\cdots ,x_n)$称为参数$\theta$的极大似然估计，其对应的统计量$\hat{\theta }(X_1,X_2\cdots ,X_n)$称为$\theta$的极大似然估计量。

求解$\hat{\theta }$：利用$\ln L(\theta )$是$L(\theta )$的增函数，故$\ln L(\theta )$和$L(\theta )$在同一点处达到最大值，于是对于似然函数$L(\theta )$取对数，利用微分学知识转换为求解对数似然方程
$$\frac{\partial \ln L(\theta )}{\partial {\theta }_j}=0,j=1,2\cdots k$$
求解$\theta$的值

极大似然估计求解问题1

例：抛一枚硬币，假设正面向上的概率是$\theta$，先对其进行了10次试验，结果为8次正面，2次反面。求$\theta$的极大似然估计。

解：先写出似然函数：
$$L(\theta )={\theta }^8(1-\theta {)}^2$$
似然函数对$\theta$求导数，并令其等于0：
$$\frac{\partial \ln L(\theta )}{\partial \theta }={\theta }^7(1-\theta )(8-10\theta )=0$$
解得：
$$\hat{\theta }=0.8$$
在该问题中，解$\hat{\theta }$的过程可以简化为如下：
$$\hat{\theta }=\arg \underset{\theta }{\max }L(\theta )$$

argmax：设$y=f(t)$在对应一个$t$值时会唯一确定一个$y$值

则：

$y=\max f(t)$，此时$y$的取值为$f(t)$的最大的output

$y=\arg \max f(t)$，此时$y$的取值为$f(t)$的最大的input

例：

设函数$f(t)，t\in (0,1,2)$，此时$f(t=0)=10$，$f(t=1)=20$，$f(t=2)=7$

则$y=\max f(t)=20$，$y=\arg \max f(t)=1$

凹函数与凸函数

中国大陆数学界某些机构关于函数凹凸性定义和国外的定义是相反的。这里使用向下方向为凸，来定义凸函数。

凸函数的主要性质有：

1.若$f$为定义在凸集$S$上的凸函数，则对任意实数$\beta \ge 0$，函数$\beta f$也是定义在$S$上的凸函数;

2.若$f_1$和$f_2$为定义在凸集$S$上的两个凸函数，则其和$f=f_1+f_2$仍为定义在$S$上的凸函数;

3.若$f_i(i=1，2\cdots m)$为定义在凸集$S$上的凸函数，则对任意实数${\beta }_i\ge 0$，函数${\beta }_i{f}_i$也是定义在$S$上的凸函数;

4.若$f$为定义在凸集$S$上的凸函数，则对每一实数$c$，水平集 S c = { x ∣ x ∈ S ， f ( x ) ≤ c } Sc=\{x|x\in S，f(x)\le c\} Sc={x∣x∈S，f(x)≤c}是凸集.

5.对于凸函数来说，局部最小值就是全局最小值。

对数似然函数

因为似然函数$L(\theta )$不是凹函数，求解其极大值时较为困难，所以一般使用与其具有相同单调性的log-likelihood

即原似然方程：
$$\hat{\theta }=\arg \underset{\theta }{\max }L(\theta )$$
转换为对数似然方程：
$$\hat{\theta }=\arg \underset{\theta }{\max }\log L(\theta )$$

EM算法

EM算法（Expectation-Maximum）分为两步，即

1. 期望（E）
2. 最大化（M）

输入：观测变量数据$Y$，隐变量数据$Z$，联合分布$P(Y,Z|\theta )$，条件分布$P(Z|Y,\theta )$

输出：模型参数$\theta$

• 选择参数处置${\theta }^{(0)}$开始迭代

• E步骤：记${\theta }^{(i)}$为第$i$次迭代参数$\theta$的估计值，在第$i+1$次迭代的E步，计算：
  $$\begin{gathered} Q(\theta ,{\theta }^{(i)})=E_z[\log P(Y,Z|\theta )|Y,{\theta }^{(i)}] \\ =\sum _z\log P(Y,Z|\theta )P(z|Y,{\theta }^{(i)}) \end{gathered}$$

这里$P(z|Y,{\theta }^{(i)})$是在给定观测数据$Y$和当前的参数估计${\theta }^{(i)}$下隐变量数据$z$的条件概率分布。

• M步骤：求使得$Q(\theta ,{\theta }^{(i)})$极大化的$\theta$，确定第$i+1$次迭代的参数的估计值${\theta }^{(i+1)}$
  $${\theta }^{(i+1)}=\arg \underset{\theta }{\max }Q(\theta ,{\theta }^{(i)})$$

• 重复上两步直至收敛

其中$Q$函数为EM算法的核心

$Q$函数：完全数据的对数似然函数$\log P(Y,Z|\theta )$关于在给定观测数据$Y$和当前参数${\theta }^{(i)}$下对未观测数据$Z$的条件概率分布$P(Y,Z|{\theta }^{(i)})$的期望称为$Q$函数，即：
$$Q(\theta ,{\theta }^{(i)})=E_z[\log P(Y,Z|\theta )|Y,{\theta }^{(i)}]$$
EM算法的几点说明：

1. 参数初值可以任意选择，但需注意EM算法对初值是敏感的。

2. E步骤求$Q(\theta ,{\theta }^{(i)})$，$Q$函数式中$Z$是未观测数据，$Y$是观测数据。注意，$Q(\theta ,{\theta }^{(i)})$的第一个变元表示要极大化的参数，第二个变元表示参数的当前估计值。每次迭代实际上在求Q函数及其极大。

3. M步骤求$Q(\theta ,{\theta }^{(i)})$的极大化，得到${\theta }^{(i+1)}$，完成一次${\theta }^{(i)}\to {\theta }^{(i+1)}$。

4. 各处停止迭代的条件，一般是对较小的正数${\xi }_1$，${\xi }_2$，若满足
   $$||{\theta }^{(i+1)}-{\theta }^{(i)}||<{\xi }_1或||Q({\theta }^{(i+1)},{\theta }^{(i)})-Q({\theta }^{(i)},{\theta }^{(i)})||<{\xi }_2$$
   则停止迭代

---

例：有三枚硬币A，B，C，其正面向上的概率分别为$\pi$，$p$，$q$，现进行10次实验，先掷硬币A，若A为正面，则掷出B，否者掷出C，记录最终的掷硬币结果，现已知实验的观测结果为1,1,0,1,0,0,1,0,1,1。估算三枚硬币正面出现的概率。

解：设$y$未观测变量（结果），$z$未隐变量（A的正反面），$\theta$为模型参数$（\pi ,p,q）$
$$\begin{gathered} L(\theta )=\sum _{z}P(y,z|\theta )=\sum _{z}P(z|\theta )P(y|z,\theta ) \\ =\pi p^y(1-p{)}^{(1-y)}+(1-\pi )q^y(1-q{)}^{(1-y)} \end{gathered}$$
将观测数据表示为$Y=(Y_1,Y_2\cdots Y_n{)}^T$，为观测到的数据表示为$Z=(Z_1,Z_2\cdots Z_n{)}^T$，则观测数据的似然函数为：
$$L(\theta )=P(Y,\theta )=\sum _{z}P(z|\theta )P(Y|Z,\theta )$$
即：
$$P(Y,\theta )=\prod _{j=1}^n[\pi p^{y_i}(1-p{)}^{1-y_i}+(1-\pi )q^{y_i}(1-q{)}^{1-y_i}]$$
求模型参数$\theta =(\pi ,p,q)$的极大似然估计即：
$$\hat{\theta }=\arg \underset{\theta }{\max }\log P(Y|\theta )$$
迭代求解上述过程：

选取模型参数初值${\theta }^{(0)}=({\pi }^{(0)},p^{(0)},q^{(0)})$，第$i$次迭代时参数为${\theta }^{(i)}=({\pi }^{(i)},p^{(i)},q^{(i)})$，则第$i+1$次迭代：

E步骤：计算在模型参数${\pi }^{(i)},p^{(i)},q^{(i)}$下观测数据$y_j$来自掷硬币B的概率

设：
$$\begin{gathered} A={\pi }^{(i)}(p^{(i)}{)}_j^y(1-p^{(i)}{)}^{(1-y_j)} \\ B=(1-{\pi }^{(i)})(q^{(i)}{)}_j^y(1-q^{(i)}{)}^{(1-y_j)} \end{gathered}$$
则E步骤：
$${\mu }_j^{(i+1)}=\frac{A}{A+B}$$
M步骤：计算模型参数的新估计值：
$${\pi }^{(i+1)}=\frac{1}{n}\sum _{j=1}^n{\mu }_j^{(i+1)}$$

$$p^{(i+1)}=\frac{{\sum }_{j=1}^n{\mu }_j^{(i+1)}{y}_j}{{\sum }_{j=1}^n{\mu }_j^{(i+1)}}$$

$$q^{(i+1)}=\frac{{\sum }_{j=1}^n(1-{\mu }_j^{(i+1)})y_j}{{\sum }_{j=1}^n(1-{\mu }_j^{(i+1)})}$$

迭代到满足设置条件

代码

EM算法的思路即为上面的内容，现使用python实现问题二的求解：

"""
@File:EM.py
@Date:2021/5/9 16:39
@Author:博0_oer~
"""


def EM(observation, theta):
    new_theta = []  # 模型参数结果
    theta_pi = theta[0]  # A的几率
    theta_p = theta[1]  # B的几率
    theta_q = theta[2]  # C的几率

    for obs in observation:
        mu_j_list = []
        for y in obs:
            A = theta_pi * theta_p ** y * (1 - theta_p) ** (1 - y)
            B = (1 - theta_pi) * theta_q ** y * (1 - theta_q) ** (1 - y)
            L_theta = A + B
            mu_j = A / L_theta
            mu_j_list.append(mu_j)
        theta_pi = 1 / len(obs) * sum(mu_j_list)

        p_mol = 0
        q_mol = 0
        p_den = 0
        q_den = 0
        for i in range(len(obs)):
            p_mol += int(obs[i]) * mu_j_list[i]
            p_den += mu_j_list[i]
            q_mol += int(obs[i]) * (1 - mu_j_list[i])
            q_den += (1 - mu_j_list[i])

        theta_p = p_mol / p_den
        theta_q = q_mol / q_den
    new_theta.append(round(theta_pi, 4))
    new_theta.append(round(theta_p, 4))
    new_theta.append(round(theta_q, 4))

    return new_theta


def iteration(observation, theta, iter_time):
    theta_new = theta
    for i in range(iter_time):
        theta_new = EM(observation, theta_new)
        print("The result of the {}-th iteration is: pi={}, p={}, q={}"
              .format(i + 1, theta_new[0], theta_new[1], theta_new[2]))
    return theta_new


if __name__ == '__main__':
    observations = [[1, 1, 0, 1, 0, 0, 1, 0, 1, 1]]
    thetas = [0.4, 0.6, 0.7]
    iter_times = 5
    print(iteration(observations, thetas, iter_times))