EM算法求解三硬币模型参数推导

EM算法

EM算法是一种迭代算法，用于含有隐变量的概率模型的参数估计。该算法通过交替进行两个步骤：E步（Expectation Step，期望步骤）和M步（Maximization Step，最大化步骤）。在E步中，利用当前参数估计值，计算隐变量的后验概率；在M步中，利用E步中计算得到的隐变量后验概率，最大化完整数据的对数似然函数，得到新的参数估计值。不断迭代这两个步骤，直到收敛为止。

EM算法一般流程

EM算法通过迭代求 $L(\theta )=\log P(Y\mid \theta )$ 的极大似然估计。每次迭代包含两步， $E$步, 求期望，$M$ 步, 求极大化。

输入: 观测变量数据 $Y$, 隐变量数据 $Z$, 联合分布 $P(Y,Z\mid \theta )$, 条件分布 $P(Z\mid Y,\theta )$;
目标输出: 模型参数 $\theta$ 。

以下为迭代步骤：
(1) 选择参数的初值 ${\theta }^{(i)}$, 开始迭代;
(2) $\mathrm{E}$ 步: 记 ${\theta }^{(i)}$ 为第 $i$ 次迭代参数 $\theta$ 的估计值, 在第 $i+1$ 次迭代的 $\mathrm{E}$ 步, 计算
$$\begin{array}{rl}Q\left(\theta ,{\theta }^{(i)}\right)& =E_Z\left[\log P(Y,Z\mid \theta )\mid Y,{\theta }^{(i)}\right]\\ & =\sum _Z\log P(Y,Z\mid \theta )P\left(Z\mid Y,{\theta }^{(i)}\right)\end{array}$$

这里, $P\left(Z\mid Y,{\theta }^{(i)}\right)$ 是在给定观测数据 $Y$ 和当前的参数估计 ${\theta }^{(i)}$ 下隐变量数据 $Z$ 的条件概率分布;
(3) $\mathrm{M}$ 步: 求使 $Q\left(\theta ,{\theta }^{(i)}\right)$ 极大化的 $\theta$, 确定第 $i+1$ 次迭代的参数的估计值 ${\theta }^{(i+1)}$
$${\theta }^{(i+1)}=\arg \underset{\theta }{\max }Q\left(\theta ,{\theta }^{(i)}\right)$$
(4) 用得出的${\theta }^{(i+1)}$带入，重复第 (2) 步和第 (3) 步, 直到收敛。

三硬币模型

首先介绍一个使用 EM 算法的例子（引自李航教授《统计学习方法（第2版 175页）》），书上并没有写这个例子的参数是如何得出的，所以我在此进行了更为基础和详细的补充。

例 9.1 (三硬币模型) 假设有 3 枚硬币, 分别记作 A, B, C。这些硬币正面出现的概率分别是 $\pi ,p$ 和 $q$ 。进行如下掷硬币试验: 先掷硬币 $\mathrm{A}$, 根据其结果选出硬币 $\mathrm{B}$或硬币 $C$, 正面选硬币 $B$, 反面选硬币 $C$; 然后掷选出的硬币, 掷硬币的结果, 出现正面记作 1 , 出现反面记作 0 ; 独立地重复 $n$ 次试验 (这里, $n=10$ ), 观测结果如下:
$$1,1,0,1,0,0,1,0,1,1$$

假设只能观测到掷硬币的结果，不能观测掷硬币的过程，也就是不知道一次实验中A硬币的正反面情况的。问如何根据一批实验数据，估计三硬币正面出现的概率, 即三硬币模型的参数。

用数学公式表达出来这个模型

$$\begin{array}{r}P(y\mid \theta )=\sum _{z}P(y,z\mid \theta )\\ =\sum _{z}P(z\mid \theta )P(y\mid z,\theta )\end{array}$$

{ p ( y = 1 ∣ θ ) = π p + ( 1 − π ) q ,  观测到正面  p ( y = 0 ∣ θ ) = π ( 1 − p ) + ( 1 − π ) ( 1 − q ) ,  观测到反面  于是 p ( y ∣ θ ) = π p y ( 1 − p ) 1 − y + ( 1 − π ) y q y ( 1 − q ) 1 − y , y ∈ { 0 , 1 } \begin{aligned} & \left\{\begin{array}{l} p(y=1 \mid \theta)=\pi p+(1-\pi) q, \text { 观测到正面 } \\ p(y=0 \mid \theta)=\pi(1-p)+(1-\pi)(1-q), \text { 观测到反面 } \end{array}\right. \\ &于是 p(y \mid \theta)=\pi p^y(1-p)^{1-y}+(1-\pi)^y q^y(1-q)^{1-y}, y \in\{0,1\} \end{aligned} ​{p(y=1∣θ)=πp+(1−π)q, 观测到正面 p(y=0∣θ)=π(1−p)+(1−π)(1−q), 观测到反面 ​于是p(y∣θ)=πpy(1−p)1−y+(1−π)yqy(1−q)1−y,y∈{0,1}​

上面的模型中$y$表示一次实验最终观测结果是1或者0，随机变量$z$是一个隐变量，表示未观测到的硬币$A$的结果为正面或是反面；$\theta =(\pi ,p,q)$ 为这个模型的参数。

我们写出观测数据$Y$的似然函数
$$P(Y\mid \theta )=\sum _{Z}P(Z\mid \theta )P(Y\mid Z,\theta )$$

$$P(Y\mid \theta )=\prod _{j=1}^n\left[\pi p^{y_j}(1-p{)}^{1-y_j}+(1-\pi )q^{y_j}(1-q{)}^{1-y_j}\right]$$

我们的目标就是求模型参数$\theta =(\pi ,p,q)$的极大似然估计
$$\hat{\theta }=\underset{\theta }{\mathrm{argmax}}[\log P(Y\mid \theta )]$$

针对这个问题，我们使用$EM$算法对上面的三硬币模型进行求解

回顾条件期望知识点

条件期望是指在给定某些条件下，对随机变量的期望进行计算。具体地，设$X$和$Y$是两个随机变量，$Y$的取值范围为$y_1,y_2,\dots ,y_n$，则在给定$Y$的条件下，$X$的条件期望$E(X\mid Y)$定义为：

$$E(X\mid Y)=\sum _{i=1}^n{X}_{i}P(X_i\mid Y)$$

其中，$X_i$表示$X$在$Y=y_i$时的取值，$P(X_i\mid Y)$表示在$Y=y_i$的条件下，$X$取$X_i$的概率。这个式子的意义是，对于每个$y_i$，计算在$Y=y_i$的条件下，$X$的期望值，然后对所有的$i$进行求和，得到$X$在给定$Y$的条件下的期望值。

EM算法

这里我们重点对E步和M步的公式进行拆解，经过M步以后拿到的参数$\theta$就可以进行下一轮迭代

E步

我们引入上面流程中的 $Q$ 函数（这里先不证明怎么来的，先拿来用，详细推导Q函数的过程可以看文末参考文章）

$$\begin{array}{rl}Q\left(\theta ,{\theta }^{(i)}\right)& =E_Z\left[\log P(Y,Z\mid \theta )\mid Y,{\theta }^{(i)}\right]\\ & =\sum _Z\log P(Y,Z\mid \theta )P\left(Z\mid Y,{\theta }^{(i)}\right)\end{array}$$

先看后半部分$P\left(Z\mid Y,{\theta }^{(i)}\right)$按照不同的$z$拆解
对其中$z=1$的情况，利用贝叶斯公式进行推导（$A$为正面，记隐变量$z=1$，表示$y$观测结果来自$B$硬币）
$$\begin{array}{rl}P(z=1\mid y,{\theta }^{(i)})=\frac{P(z=1,y,{\theta }^{(i)})}{P(y,{\theta }^{(i)})}& =\frac{P(z=1)P(y,{\theta }^{(i)}\mid z=1)}{{\sum }_{z}P(z)P(y,{\theta }^{(i)}\mid z)}\\ & =\frac{P(z=1)p(y,{\theta }^{(i)}\mid z=1)}{P(z=1)p(y,{\theta }^{(i)}\mid z=1)+P(z=0)p(y,{\theta }^{(i)}\mid z=0)}\\ & =\frac{\pi p^y(1-p{)}^{1-y}}{\pi p^y(1-p{)}^{1-y}+(1-\pi )q^y(1-q{)}^{1-y}}\end{array}$$
同理可得$z=0$时
$$\begin{array}{rl}P(z=0\mid y,{\theta }^{(i)})=\frac{P(z=0,y,{\theta }^{(i)})}{P(y,{\theta }^{(i)})}& =\frac{P(z=0)P(y,{\theta }^{(i)}\mid z=0)}{{\sum }_{z}P(z)P(y,{\theta }^{(i)}\mid z)}\\ & =\frac{P(z=0)p(y,{\theta }^{(i)}\mid z=0)}{P(z=1)p(y,{\theta }^{(i)}\mid z=1)+P(z=0)p(y,{\theta }^{(i)}\mid z=0)}\\ & =\frac{(1-\pi )q^y(1-q{)}^{1-y}}{\pi p^y(1-p{)}^{1-y}+(1-\pi )q^y(1-q{)}^{1-y}}\end{array}$$
我们设${\mu }^{(i)}=p(z=1\mid y,{\theta }^{(i)})$表示y观测结果来自B硬币的概率（A正面），则硬币来自C的概率（A结果反面）为$1-{\mu }^{(i)}$
此时可以对$Q\left(\theta ,{\theta }^{(i)}\right)$继续拆解
$$\begin{array}{rl}Q\left(\theta ,{\theta }^{(i)}\right)& =E_Z\left[\log P(Y,Z\mid \theta )\mid Y,{\theta }^{(i)}\right]\\ & =\sum _Z\log P(Y,Z\mid \theta )P\left(Z\mid Y,{\theta }^{(i)}\right)\\ & =\log P(Y,Z=1\mid \theta )P\left(Z=1\mid Y,{\theta }^{(i)}\right)+\log P(Y,Z=0\mid \theta )P\left(Z=0\mid Y,{\theta }^{(i)}\right)\\ & =\log P(Z=1\mid \theta )P(Y\mid Z=1,\theta )\cdot {\mu }^{(i)}+\log P(Z=0\mid \theta )P(Y\mid Z=0,\theta )\cdot \left(1-{\mu }^{(i)}\right)\\ & ={\mu }^{(i)}\cdot \log \left[\pi p^y(1-p{)}^{1-y}\right]+\left(1-{\mu }^{(i)}\right)\cdot \log \left[(1-\pi )q^y(1-q{)}^{1-y}\right]\\ & ={\mu }^{(i)}\left[\log \pi +\log p^y(1-p{)}^{1-y}\right]+\left(1-{\mu }^{(i)}\right)\cdot \left[\log (1-\pi )+\log q^y(1-q{)}^{1-y}\right]\end{array}$$

M步

求使 $Q\left(\theta ,{\theta }^{(i)}\right)$ 极大化的 $\theta$，确定第 $i$ 次迭代的参数的估计值 ${\theta }^{(i+1)}$。设共有 $n$ 组观测值 $y_j,1\le j\le n$，则 $Q$ 函数可以表述为 $Q^{\prime }={\sum }_j^{n}Q\left(\theta ,{\theta }^{(i)}\right)$。求使 $Q^{\prime }$ 最大化的参数 $\theta =(\pi ,p,q)$，我们令$Q^{\prime }$的偏导为 $0$，分别求出参数 $\theta =(\pi ,p,q)$。

以下展示求解 $\pi$ 的过程：
$$\begin{array}{rl}\frac{\partial Q^{\prime }}{\partial \pi }=\sum _j^n\frac{\partial Q}{\partial \pi }& =\sum _j^n\frac{\partial }{\partial \pi }\left[{\mu }_j^{(i)}\left[\log \pi +\log p^{y_j}(1-p{)}^{1-y_j}\right]+\left(1-{\mu }_j^{(i)}\right)\cdot \left[\log (1-\pi )+\log q^{y_j}(1-q{)}^{1-y_j}\right]\right]\\ & =\sum _j^n\frac{\partial }{\partial \pi }\left[{\mu }_j^{(i)}\log \pi +(1-{\mu }_j^{(i)})\log (1-\pi )\right]\\ & =\sum _j^n\frac{{\mu }_j^{(i)}}{\pi }+\frac{{\mu }_j^{(i)}-1}{1-\pi }=0\\ & =\sum _j^n(1-\pi ){\mu }_j^{(i)}+\pi \left({\mu }_j^{(i)}-1\right)=0\\ & \sum _j^n{\mu }_j^{(i)}=\sum _j^n\pi \end{array}$$
即 $$\sum _j^n{\mu }_j^{(i)}=n\pi$$ 于是$$\pi =\frac{1}{n}\sum _j^n{\mu }_j^{(i)}$$
同样的步骤，我们分别令${\sum }_j^n\frac{\partial Q}{\partial p}=0$和${\sum }_j^n\frac{\partial Q}{\partial q}=0$，即可求出下次迭代的全部参数 $\theta =(\pi ,p,q)$ 可以进行下一轮迭代优化。

$$\begin{array}{rl}& \text{ 令 }\sum _j^n\frac{\partial Q}{\partial p}=0\\ & =\sum _j^n\frac{\partial }{\partial p}\left[{\mu }_j^{(i)}\left[\log \pi +\log p^{y_j}(1-p{)}^{1-y_j}\right]+\left(1-{\mu }_j^{(i)}\right)\cdot \left[\log (1-\pi )+\log q^{y_j}(1-q{)}^{1-y_j}\right]\right]\\ & =\sum _j^n\frac{\partial }{\partial p}\left[{\mu }_j^{(i)}\cdot \log p^{y_j}(1-p{)}^{(1-y_j)}\right]\\ & =\sum _j^n\frac{\partial }{\partial p}[y_j\cdot {\mu }_j^{(i)}\log p+{\mu }_j^{(i)}(1-y_j)\log (1-p)]\\ & =\sum _j^n\left[\frac{y_j\cdot {\mu }_j^{(i)}}{p}+\frac{{\mu }_j^{(i)}\cdot (1-y_j)\cdot (-1)}{1-p}\right]=0\\ & =\sum _j^n{y}_j\cdot {\mu }_j^{(i)}-{\mu }_j^{(i)}\cdot y_j\cdot p+p\cdot y_j\cdot {\mu }_j^{(i)}-p\cdot {\mu }_j^{(i)}=0\end{array}$$
即
$$\begin{gathered} \sum _j^n{y}_j\cdot {\mu }_j^{(i)}=\sum _j^{n}p\cdot {\mu }_j^{(i)} \\ p=\frac{{\sum }_j^n{y}_j\cdot {\mu }_j^{(i)}}{{\sum }_j^n{\mu }_j^{(i)}} \end{gathered}$$

接下来
$$\begin{array}{rl}& \text{ 令 }\sum _j^n\frac{\partial Q}{q}=0\\ & =\sum _j^n\frac{\partial }{\partial q}\left[{\mu }_j^{(i)}\left[\log \pi +\log p^{y_j}(1-p{)}^{1-y_j}\right]+\left(1-{\mu }_j^{(i)}\right)\cdot \left[\log (1-\pi )+\log q^{y_j}(1-q{)}^{1-y_j}\right]\right]\\ & =\sum _j^n\frac{\partial }{\partial q}\left[\left(1-{\mu }_j^{(i)}\right)\cdot \log q^{y_j}+\left(1-{\mu }_j^{(i)}\right)\cdot \left(1-y_j\right)\cdot \log (1-q)\right]\\ & =\sum _j^n\frac{\left(1-{\mu }_j^{(i)}\right)\cdot y_j}{q}+\frac{\left(1-{\mu }_j^{(i)}\right)\cdot \left(1-y_j\right)\cdot (-1)}{1-q}\\ & =\sum _j^n\left(y_j\cdot (1-q)-{\mu }_j^{(i)}\cdot y_j(1-q)+q\cdot \left(y_j-1\right)-q\cdot {\mu }_j^{(i)}\cdot \left(y_j-1\right)\right)\\ & =\sum _j^n\left(y_j-y_j\cdot q-{\mu }_j^{(i)}\cdot y_j+{\mu }_j^{(i)}\cdot y_j\cdot q+y_{j}q-q-{\mu }_j^{(i)}\cdot y_j\cdot q+{\mu }_j^{(i)}\cdot q\right)\\ & =\sum _j^n\left(y_j-{\mu }_j^{(i)}\cdot y_j-q+{\mu }_j^{(i)}\cdot q\right)\\ & \end{array}$$
即
$$\begin{array}{r}\sum _j^n{y}_j-{\mu }_j^{(i)}\cdot y_j=\sum _j^{n}q\left(1-{\mu }_j^{(i)}\right)\\ q=\frac{{\sum }_j^n\left(1-{\mu }_j^{(i)}\right)\cdot y_j}{{\sum }_j^n\left(1-{\mu }_j^{(i)}\right)}\end{array}$$

至此，我们已经求得参数 $\theta =(\pi ,p,q)$分别为
$$\begin{gathered} \pi =\frac{1}{n}\sum _j^n{\mu }_j^{(i)} \\ p=\frac{{\sum }_j^n{y}_j\cdot {\mu }_j^{(i)}}{{\sum }_j^n{\mu }_j^{(i)}} \\ q=\frac{{\sum }_j^n\left(1-{\mu }_j^{(i)}\right)\cdot y_j}{{\sum }_j^n\left(1-{\mu }_j^{(i)}\right)} \end{gathered}$$

有了上面参数的迭代关系及样本
$$1,1,0,1,0,0,1,0,1,1$$
我们设初始参数
$${\pi }^{(0)}=0.5,p^{(0)}=0.5,q^{(0)}=0.5$$
由${\mu }_j^{(i)}$的表达式得到无论$y_j=0$还是$y_j=1$，均有${\mu }_j^{(1)}=0.5$
于是下一轮
$$\begin{array}{rl}{\pi }^{(1)}& =\frac{1}{10}\left({\mu }_1^{(0)}+{\mu }_2^{(0)}+{\mu }_3^{(0)}\cdots +{\mu }_{10}^{(0)}\right)\\ & =0.5\\ p^{(1)}& =\frac{y_1\cdot {\mu }_1^{(0)}+y_2\cdot {\mu }_2^{(0)}+y_3\cdot {\mu }_3^{(0)}+\cdots +y_{10}\cdot {\mu }_{10}^{(0)}}{\left({\mu }_1^{(0)}+{\mu }_2^{(0)}+{\mu }_3^{(0)}\cdots +{\mu }_{10}^{(0)}\right)}\\ & =\frac{6\times 0.5}{10\times 0.5}\\ & =0.6\\ q^{(1)}& =\frac{\left(1-{\mu }_1^{(0)}\right)\cdot y_1+\left(1-{\mu }_2^{(0)}\right)\cdot y_2\cdots +\left(1-{\mu }_{10}^{(0)}\right)\cdot y_{10}}{\left(1-{\mu }_1^{(0)}\right)+\left(1-{\mu }_2^{(0)}\right)+\cdots +\left(1-{\mu }_{10}^{(0)}\right)}\\ & =\frac{6\times 0.5}{10-0.5\times 10}=0.6\end{array}$$
一直迭代下去，我们可以得到最终参数 $\theta =(\pi ,p,q)$的一个最终值，即为三硬币模型的参数。

使用Python代码进行实验模拟，检测是否有效

以下为一段代码，用于生成实验数据，并用EM算法预估参数，用得到的预估参数进行多次模拟实验，检验模型预测能力

import random
import numpy as np
import time
def coin_experiment(pi, p, q, n):
    """
    模拟掷硬币实验
    pi: 硬币 A 正面出现的概率
    p: 硬币 B 正面出现的概率
    q: 硬币 C 正面出现的概率
    n: 实验次数
    """
    results = []
    results_A = []
    results_B = []
    results_C = []
    for i in range(n):
        # 先掷硬币 A
        if random.random() < pi:
            # 选硬币 B
            coin = 'B'
            p_head = p
        else:
            # 选硬币 C
            coin = 'C'
            p_head = q
        
        # 接着掷选出的硬币
        if random.random() < p_head:
            results.append(1)
        else:
            results.append(0)
        
        # 记录每个硬币的正反面
        if coin == 'B':
            if random.random() < p:
                results_B.append(1)
            else:
                results_B.append(0)
            results_A.append(1)
        else:
            if random.random() < q:
                results_C.append(1)
            else:
                results_C.append(0)
            results_A.append(0)
    
    # 计算 A、B、C 硬币的正面概率
    p_A = sum(results_A) / len((results_A))
    p_B = sum(results_B) / len(results_B)
    p_C = sum(results_C) / len(results_C)
    
    return results, p_A, p_B, p_C


pi = 0.2
p = 0.3
q = 0.8
n = 100000
s=time.time()
print(f'开始模拟，模拟参数为pi={pi},p={p},q={q}……')
Y,pi,p,q=coin_experiment(pi, p, q, n)
Y= np.array(Y)
print(f'模拟结束，共模拟{n}次，耗时{time.time()-s}……')
print(f"模拟结果，pi={pi},p={p},q={q},整体实验Y=1概率={sum(Y)/len(Y)}")

#EM参数反推，任意设定初始参数
pi_0 = 0.9
p_0 = 0.8
q_0 = 0.9
epsiodes=1000 #迭代次数
count=1
while count<=epsiodes:
    mu = pi_0 * p_0**Y * (1-p_0)**(1-Y) / (pi_0 * p_0**Y * (1-p_0)**(1-Y) + (1-pi_0) * q_0**Y * (1-q_0)**(1-Y))
    pi=(1/n)*sum(mu)
    p=sum(Y*mu)/sum(mu)
    q=sum((1-mu)*Y)/sum(1-mu)
    if count%100==0:
        print(f"第{count}次迭代，估算参数分别为：pi={pi},p={p},q={q}")
    pi_0 = pi
    p_0 = p
    q_0 = q
    count+=1

#用拿到的估计参数重新模拟下
s=time.time()
print(f'二次模拟检验，模拟参数为pi={pi},p={p},q={q}……')
Y,pi,p,q=coin_experiment(pi, p, q, n)
Y= np.array(Y)
print(f'模拟结束，共模拟{n}次，耗时{time.time()-s}……')
print(f"模拟结果，pi={pi},p={p},q={q},整体实验Y=1概率={sum(Y)/len(Y)}")

观察实验，按照我们设定的初始pi = 0.2，p = 0.3，q = 0.8，进行100000次实验

模拟过程中pi,p,q和我们设定的值非常接近，实验确实是按这个参数进行，整体实验结果y=1的概率为0.696，最终EM算法迭代1000次求出来的参数pi = 0.91，p = 0.68，q = 0.83，虽然和我们模型实际参数相差很大，但是我们用这个参数进行了新的100000次实验，最终整体实验结果y=1的概率为0.6983，和我们实际模型的结果很接近，实现了比较精准的模拟，模型是有效的。

参考资料

[1].EM算法Q函数推导过程详解